Emil M. Constantinescu

Shrirang Abhyankar, Jed Brown, Emil M. Constantinescu et al.

High-quality ordinary differential equation (ODE) solver libraries have a long history, going back to the 1970s. Over the past several years we have implemented, on top of the PETSc linear and nonlinear solver package, a new general-purpose, extensive, extensible library for solving ODEs and differential algebraic equations (DAEs). Package includes support for both forward and adjoint sensitivities that can be easily utilized by the TAO optimization package, which is also part of PETSc. The ODE/DAE integrator library strives to be highly scalable but also to deliver high efficiency for modest-sized problems. The library includes explicit solvers, implicit solvers, and a collection of implicit-explicit solvers, all with a common user interface and runtime selection of solver types, adaptive error control, and monitoring of solution progress. The library also offers enormous flexibility in selection of nonlinear and linear solvers, including the entire suite of PETSc iterative solvers, as well as several parallel direct solvers.

Daniel S. Abdi, Francis X. Giraldo, Emil M. Constantinescu et al.

We present the acceleration of an IMplicit-EXplicit (IMEX) non-hydrostatic atmospheric model on manycore processors such as GPUs and Intel's MIC architecture. IMEX time integration methods sidestep the constraint imposed by the Courant-Friedrichs-Lewy condition on explicit methods through corrective implicit solves within each time step. In this work, we implement and evaluate the performance of IMEX on manycore processors relative to explicit methods. Using 3D-IMEX at Courant number C=15 , we obtained a speedup of about 4X relative to an explicit time stepping method run with the maximum allowable C=1. In addition, we demonstrate a much larger speedup of 100X at C=150 using 1D-IMEX due to the unconditional stability of the method in the vertical direction. Several improvements on the IMEX procedure were necessary in order to outperform our results with explicit methods: a) reducing the number of degrees of freedom of the IMEX formulation by forming the Schur complement; b) formulating a horizontally-explicit vertically-implicit (HEVI) 1D-IMEX scheme that has a lower workload and potentially better scalability than 3D-IMEX; c) using high-order polynomial preconditioners to reduce the condition number of the resulting system; d) using a direct solver for the 1D-IMEX method by performing and storing LU factorizations once to obtain a constant cost for any Courant number. Without all of these improvements, explicit time integration methods turned out to be difficult to beat. We discuss in detail the IMEX infrastructure required for formulating and implementing efficient methods on manycore processors. Finally, we validate our results with standard benchmark problems in NWP and evaluate the performance and scalability of the IMEX method using up to 4192 GPUs and 16 Knights Landing processors.

Debojyoti Ghosh, Emil M. Constantinescu

This paper presents a characteristic-based flux partitioning for the semi-implicit time integration of atmospheric flows. Nonhydrostatic models require the solution of the compressible Euler equations. The acoustic time-scale is significantly faster than the advective scale, yet it is typically not relevant to atmospheric and weather phenomena. The acoustic and advective components of the hyperbolic flux are separated in the characteristic space. High-order, conservative additive Runge-Kutta methods are applied to the partitioned equations so that the acoustic component is integrated in time implicitly with an unconditionally stable method, while the advective component is integrated explicitly. The time step of the overall algorithm is thus determined by the advective scale. Benchmark flow problems are used to demonstrate the accuracy, stability, and convergence of the proposed algorithm. The computational cost of the partitioned semi-implicit approach is compared with that of explicit time integration.

Simone Marras, Michal A. Kopera, Emil M. Constantinescu et al.

The high-order numerical solution of the non-linear shallow water equations (and of hyperbolic systems in general) is susceptible to unphysical Gibbs oscillations that form in the proximity of strong gradients. The solution to this problem is still an active field of research as no general cure has been found yet. In this paper, we tackle this issue by presenting a dynamically adaptive viscosity based on a residual-based sub-grid scale model that has the following properties: $(i)$ it removes the spurious oscillations in the proximity of strong wave fronts while preserving the overall accuracy and sharpness of the solution. This is possible because of the residual-based definition of the dynamic diffusion coefficient. $(ii)$ For coarse grids, it prevents energy from building up at small wave-numbers. $(iii)$ The model has no tunable parameter. Our interest in the shallow water equations is tied to the simulation of coastal inundation, where a careful handling of the transition between dry and wet surfaces is particularly challenging for high-order Galerkin approximations. In this paper, we extend to a unified continuous/discontinuous Galerkin (CG/DG) framework a very simple, yet effective wetting and drying algorithm originally designed for DG [Xing, Zhang, Shu (2010)]. We show its effectiveness for problems in one and two dimensions on domains of increasing characteristic lengths varying from centimeters to kilometers. Finally, to overcome the time-step restriction incurred by the high-order Galerkin approximation, we advance the equations forward in time via a three stage, second order explicit-first-stage, singly diagonally implicit Runge-Kutta (ESDIRK) time integration scheme. Via ESDIRK, we are able to preserve numerical stability for an advective CFL number 10 times larger than its explicit counterpart.

Shinhoo Kang, Emil M. Constantinescu

We propose a new approach to learning the subgrid-scale model when simulating partial differential equations (PDEs) solved by the method of lines and their representation in chaotic ordinary differential equations, based on neural ordinary differential equations (NODEs). Solving systems with fine temporal and spatial grid scales is an ongoing computational challenge, and closure models are generally difficult to tune. Machine learning approaches have increased the accuracy and efficiency of computational fluid dynamics solvers. In this approach neural networks are used to learn the coarse- to fine-grid map, which can be viewed as subgrid-scale parameterization. We propose a strategy that uses the NODE and partial knowledge to learn the source dynamics at a continuous level. Our method inherits the advantages of NODEs and can be used to parameterize subgrid scales, approximate coupling operators, and improve the efficiency of low-order solvers. Numerical results with the two-scale Lorenz 96 ODE, the convection-diffusion PDE, and the viscous Burgers' PDE are used to illustrate this approach.

Shinhoo Kang, Emil M. Constantinescu

Computational advances have fundamentally transformed the landscape of numerical simulations, enabling unprecedented levels of complexity and precision in modeling physical phenomena. While these high-fidelity simulations offer invaluable insights for scientific discovery and problem solving, they impose substantial computational requirements. Consequently, low-fidelity models augmented with subgrid-scale parameterizations are employed to achieve computational feasibility. We introduce an end-to-end differentiable framework for solving the compressible Navier--Stokes equations. This integrated approach combines a differentiable discontinuous Galerkin (DG) solver with a neural network source term. Through the implementation of neural ordinary differential equations (NODEs) for network parameter optimization, our methodology ensures continuous interaction with the governing equations throughout the training process. We refer to this approach as NODE-DG. This hybrid approach combines the accuracy of numerical methods with the efficiency of machine learning, offering the following key advantages: (1) improved accuracy of low-order DG approximations by capturing subgrid-scale dynamics; (2) robustness against nonuniform or missing temporal data; (3) elimination of operator-splitting errors; (3) total mass conservation; and (4) a continuous-in-time operator that enables variable time step predictions, which accelerate projected high-order DG simulations. We demonstrate the performance of the proposed framework through two examples: two-dimensional Kelvin--Helmholtz instability and three-dimensional Taylor--Green vortex examples.