We present scalable iterative solvers and preconditioning strategies for Hybridizable Discontinuous Galerkin (HDG) discretizations of partial differential equations (PDEs) on graphics processing units (GPUs). The key idea lies in the systematic reformulation of the HDG system to expose dense face-level block structure, enabling the entire linear algebra pipeline to be executed using strided batched dense kernels rather than conventional sparse matrix formats. By eliminating the local element degrees of freedom in parallel and assembling the global condensed system on the device, the approach avoids sparse data structures, increases arithmetic intensity, and achieves sustained high throughput on modern GPU architectures. We develop and compare a family of additive-type preconditioners such as block-Jacobi, overlapping additive Schwarz, and polynomial precondioning, which are designed for GPU efficiency and portability. All preconditioners are implemented in batched form using vendor-optimized dense linear algebra libraries, ensuring high concurrency while maintaining architectural portability across NVIDIA and AMD platforms. Comprehensive performance studies are conducted for a broad class of elliptic and compressible flow problems, including the Poisson equation, Burgers equation, linear and nonlinear elasticity, Euler equations, Navier-Stokes equations, and Reynolds-Averaged Navier-Stokes equations, on structured and unstructured meshes with varying element types and polynomial orders. Detailed profiling and benchmarking demonstrate that the dominant kernels operate near architectural limits and that additive and polynomial-enhanced strategies provide stable convergence across diverse physical regimes. The results demonstrate a portable and high-performance pathway for HDG solvers on GPU architectures.

Advances in machine-learned interatomic potentials have enabled the prediction of complex material properties with accuracy approaching that of ab initio methods. However, it is unclear how the finite capacity of such models affects their ability to achieve consistent accuracy across diverse thermodynamic conditions without introducing trade-offs. In this paper, we present two computationally efficient interatomic potentials capable of accurately simulating the behavior of hafnium and hafnium dioxide across a very wide variety of thermodynamic conditions. Our approach combines Latin Hypercube and Monte Carlo Sampling for generating diverse data sets, with an extended formulation of the recently-developed environment-adaptive proper orthogonal descriptors. Molecular dynamics simulations show that the resulting potentials accurately reproduce density functional theory results and experimental data for pressure- and temperature-induced phase transitions as well as other properties associated with the materials’ polymorphs and liquid phases. We further showcase the versatility of the environment-adaptive formulation by using our potential to compute the shock Hugoniot of hafnium up to temperatures and pressures of 1 MK and 1 TPa, respectively; good agreement with available experimental data is observed.

We present novel model reduction methods for rapid solution of parametrized nonlinear partial differential equations (PDEs) in real-time or many-query contexts. Our approach combines reduced basis (RB) space for global approximation of the parametric solution manifold, Galerkin projection of the underlying PDEs onto the RB space for dimensionality reduction, and high-order empirical interpolation for efficient treatment of the nonlinear terms. We propose a class of high-order empirical interpolation methods to derive basis functions and interpolation points by using high-order partial derivatives of the nonlinear terms. We develop error indicator to estimate the interpolation errors and determine parameter points via greedy sampling. Furthermore, we introduce two hyperreduction schemes to construct reduced-order models: one that applies the hyperreduction technique before Newton’s method and another after. The latter scheme significantly reduces hyperreduction errors while maintaining computational efficiency. Numerical results are presented to demonstrate the accuracy and efficiency of our approach.

Direct numerical simulation (DNS) of compressible transitional and turbulent flows requires numerical methods that combine high-order accuracy, robustness, and computational efficiency to resolve a broad range of spatial and temporal scales. This paper presents a massively parallel hybridizable discontinuous Galerkin (HDG) solver for DNS of the compressible Navier-Stokes equations on GPU-accelerated high-performance computing systems. The proposed solver combines high-order HDG discretization with robust shock capturing, diagonally implicit Runge-Kutta (DIRK) time integration, and an efficient Newton-GMRES solution strategy accelerated by additive Schwarz preconditioning and reduced-basis approximation. A distributed implementation of these methods based on GPU-aware MPI, Kokkos, and CUDA/HIP libraries enables scalable execution on heterogeneous computing platforms. The solver is demonstrated on three canonical benchmark problems covering a wide range of Mach-number flow regimes: subsonic transitional flow over the Eppler 387 airfoil, the supersonic Taylor-Green vortex, and hypersonic boundary-layer transition. Numerical results are compared with available experimental measurements and published DNS data, showing good agreement across distinct flow regimes. The results demonstrate the ability of the proposed solver to resolve laminar-turbulent transition, strong compressibility effects, shock-associated flow structures, and fully three-dimensional turbulent dynamics.

This work presents an adaptive high-order method for the simulation of parametrized high-speed flows. To obviate the proliferation of degrees of freedom inherent in the discontinuous Galerkin discretization chosen here, we use a hybridized method that decreases the size of the global system along with r-adaptive meshing. For efficiency and robustness, we use a strategy of artificial viscosity continuation, wherein the amount of added dissipation is ramped down until physicality and smoothness constraints are violated. We demonstrate this procedure on problems with greater physical complexity than considered in prior works by incorporating coupled aerothermal interactions between fluid and solid domains. In addition, we have previously found that this choice of mesh adaptation method aids in the creation of surrogate models for flows with shocks. We evaluate whether these advantages persist for viscous flows with multiphysics couplings, demonstrating the use of this method for the forward propagation of uncertain boundary conditions for high-speed flows.

We present a model reduction approach for the real-time solution of time-dependent nonlinear partial differential equations (PDEs) with parametric dependencies. A major challenge in constructing efficient and accurate reduced-order models for nonlinear PDEs is the efficient treatment of nonlinear terms. We address this by unifying the implementation of hyperreduction methods to deal with nonlinear terms. Furthermore, we introduce a first-order empirical interpolation method (EIM) to provide an efficient approximation of the nonlinear terms in time-dependent PDEs. We demonstrate the effectiveness of our approach on the Allen-Cahn equation, which models phase separation, and the Buckley-Leverett equation, which describes two-phase fluid flow in porous media. Numerical results highlight the accuracy, efficiency, and stability of the proposed method compared with both the Galerkin-Newton approach and hyper-reduced models using the standard EIM.

We implemented a scalable hybridizable discontinuous Galerkin (HDG) code to model gas discharges. Protection against lightning-induced ignition in aircraft fuel tanks is a critical design and certification consideration, especially for next-generation aircraft utilizing hydrogen fuels, which have ignition energies an order of magnitude lower than those of conventional kerosene. Fuel tanks, with exposed fasteners and complex geometries, present multiple pathways for energy deposition during a lightning strike, including sparking, arcing, and edge glow. To quantify ignition risk and assess these pathways, we require robust numerical models capable of handling various gas compositions, states, and geometries. In this work, we develop and validate a high-fidelity numerical framework for modeling gaseous breakdown and low-temperature plasmas (LTPs) using Exasim software. Exasim is a scalable DG/HDG finite element solver that supports high-order approximations on unstructured grids and allows for modular, multiphysics coupling. HDG methods, with their reduced global degrees of freedom and optimal convergence properties, are particularly suited for discharge modeling, which depends on accurate electric field computations for species transport and reaction rates. Our framework employs a fluid approximation for the transport of charged species and the Poisson equation for a self-consistent electric field description. The framework is designed to be scalable and flexible, enabling prototyping on local systems and deployment on high-performance computing clusters with graphics processors. We present three simulation cases to demonstrate the framework's capabilities, ranging from basic validation to increasingly complex configurations. This work provides a validated, extensible and scalable tool for studying the physics of gas discharges and the energy deposition pathways critical to lightning-induced ignition risk assessment in advanced aircraft designs.

We present a generative reduced basis (RB) approach for the rapid and reliable solution parametrized linear partial differential equations. Central to this approach is the construction generative RB spaces that provide rapidly convergent approximations of the solution manifold. We propose a generative snapshot method to generate significantly larger sets of snapshots from a small initial set of solution snapshots. This method leverages multivariate nonlinear transformations to enrich the RB spaces, thereby enabling a more accurate approximation of the solution manifold than commonly used dimensionality reduction techniques such proper orthogonal decomposition and greedy sampling. We employ the generative RB spaces to construct reduced order models and compute aposteriori error estimates. The error estimates allow us to efficiently explore the parameter space and select parameter points that improve the efficiency and accuracy of the reduced order model. Through numerical experiments, demonstrate that the generative RB method not only improves the accuracy of the reduced order model but also provides tight error estimates.

We introduce the proper orthogonal descriptors for efficient and accurate interatomic potentials of multi-element chemical systems. The potential energy surface of a multi-element system is represented as a many-body expansion of parametrized potentials which are functions of atom positions, atom types, and parameters. The proper orthogonal decomposition is employed to decompose the parametrized potentials as a linear combination of orthogonal basis functions. The orthogonal basis functions are used to construct proper orthogonal descriptors based on the elements of atoms, thus leading to multi-element descriptors. We compose the multi-element proper orthogonal descriptors to develop linear and quadratic interatomic potentials. We devise an algorithm to efficiently compute the total energy and forces of the interatomic potentials constructed from the proper orthogonal descriptors. The potentials are demonstrated for indium phosphide and titanium dioxide in comparison with the spectral neighbor analysis potential (SNAP) and atomic cluster expansion (ACE) potentials.

The development of interatomic potentials that can accurately capture a wide range of phys ical phenomena and diverse environments is of significant interest, but it presents a formidable challenge. This challenge arises from the numerous structural forms, multiple phases, complex intramolecular and intermolecular interactions, and varying external conditions. In this paper, we present a method to construct environment-adaptive interatomic potentials by adapting to the local atomic environment of each atom within a system. The collection of atomic environments of interest is partitioned into several clusters of atomic environments. Each cluster represents a distinctive local environment and is used to define a corresponding local potential. We introduce a many-body many-potential expansion to smoothly blend these local potentials to ensure global continuity of the potential energy surface. This is achieved by computing the probability functions that determine the likelihood of an atom belonging to each cluster. We apply the environment-adaptive machine learning potentials to predict observable properties for Ta element and InP compound, and compare them with density functional theory calculations.

The leading edge region plays a crucial role in the receptivity of hypersonic boundary layers to freestream disturbances due to the amplification of flow instabilities downstream. However, the vast majority of large eddy simulations (LES) and direct numerical simulations (DNS) exclude this very important region from the computational domain. Accurately resolving the shock wave caused by the leading edge can require a mesh resolution that is an order of magnitude finer than that which is needed to resolve the turbulent structures downstream. This results in a dramatic reduction of the maximum allowable timestep that can be used by the solver, ultimately making the simulation take too long to run. In this paper, we explore the nuances of leading edge region mesh generation as well as synthetic turbulence generation. Our DNS results for a Mach 6 hypersonic flat plate with an impinging shock wave are compared against experimental data and previous DNS studies.

In recent years, high-order discontinuous Galerkin (DG) methods have emerged as an attractive approach for numerical simulations of compressible flows. This paper presents an overview of the recent development of DG methods for compressible flows with particular focus on hypersonic flows. First, we survey state-of-the-art DG methods for computational fluid dynamics. Next, we discuss both matrix-based and matrix-free iterative methods for the solution of discrete systems stemming from the spatial DG discretizations of the compressible Navier–Stokes equations. We then describe various shock capturing methods to deal with strong shock waves in hypersonic flows. We discuss adaptivity techniques to refine high-order meshes, and synthetic boundary conditions to simulate free-stream disturbances in hypersonic boundary layers. We present a few examples to demonstrate the ability of high-order DG methods to provide accurate solutions of hypersonic laminar flows. Furthermore, we present direct numerical simulations of hypersonic transitional flow past a flared cone at Reynolds number 10.8×106, and hypersonic transitional shock wave boundary layer interaction flow over a flat plate at Reynolds number 3.97×106. These simulations run entirely on hundreds of graphics processing units (GPUs) and demonstrate the ability of DG methods to directly resolve hypersonic transitional flows, even at high Reynolds numbers, without relying on transition or turbulence models. We end the paper by offering our perspectives on error estimation, turbulence modeling, and real gas effects in hypersonic flows.

We present an optimal transport approach for mesh adaptivity and shock capturing of compressible flows. Shock capturing is based on a viscosity regularization of the governing equations by introducing an artificial viscosity field as solution of the modified Helmholtz equation. Mesh adaptation is based on the optimal transport theory by formulating a mesh mapping as solution of Monge-Ampère equation. The marriage of optimal transport and viscosity regularization for compressible flows leads to a coupled system of the compressible Euler/Navier-Stokes equations, the Helmholtz equation, and the Monge-Ampère equation. We propose an iterative procedure to solve the coupled system in a sequential fashion using homotopy continuation to minimize the amount of artificial viscosity while enforcing positivity-preserving and smoothness constraints on the numerical solution. We explore various mesh monitor functions for computing r-adaptive meshes in order to reduce the amount of artificial dissipation and improve the accuracy of the numerical solution. The hybridizable discontinuous Galerkin method is used for the spatial discretization of the governing equations to obtain high-order accurate solutions. Extensive numerical results are presented to demonstrate the optimal transport approach on transonic, supersonic, hypersonic flows in two dimensions. The approach is found to yield accurate, sharp yet smooth solutions within a few mesh adaptation iterations.

We present model reduction techniques for parametrized nonlinear partial differential equations (PDEs). The main ingredients of our approach are reduced basis (RB) spaces to provide rapidly convergent approximations to the parametric manifold; Galerkin projection of the underlying PDEs onto the RB space to provide reduction in the number of degrees of freedom; and empirical interpolation schemes to provide rapid evaluation of the nonlinear terms associated with the Galerkin projection. We devise a first-order empirical interpolation method to construct an inexpensive and stable interpolation of the nonlinear terms. We consider two different hyper-reduction strategies: hyper-reduction followed by linearization, and linearization followed by hyper-reduction. We extend empirical interpolation to nonintrusive model reduction and apply it to compressible flows in both supersonic and hypersonic regimes. Numerical results are presented to illustrate the accuracy, efficiency, and stability of the reduced-order models.

We introduce two hybridizable discontinuous Galerkin (HDG) methods for numerically solving the Monge-Ampere equation. The first HDG method is devised to solve the nonlinear elliptic Monge-Ampere equation by using Newton's method. The second HDG method is devised to solve a sequence of the Poisson equation until convergence to a fixed-point solution of the Monge-Ampere equation is reached. Numerical examples are presented to demonstrate the convergence and accuracy of the HDG methods. Furthermore, the HDG methods are applied to r-adaptive mesh generation by redistributing a given scalar density function via the optimal transport theory. This r-adaptivity methodology leads to the Monge-Ampere equation with a nonlinear Neumann boundary condition arising from the optimal transport of the density function to conform the resulting high-order mesh to the boundary. Hence, we extend the HDG methods to treat the nonlinear Neumann boundary condition. Numerical experiments are presented to illustrate the generation of r-adaptive high-order meshes on planar and curved domains.

The accurate prediction of aerodynamic drag on satellites orbiting in the upper atmosphere is critical to the operational success of modern space technologies, such as satellite-based communication or navigation systems, which have become increasingly popular in the last few years due to the deployment of constellations of satellites in low-Earth orbit. As a result, physics-based models of the ionosphere and thermosphere have emerged as a necessary tool for the prediction of atmospheric outputs under highly variable space weather conditions. This paper proposes a high-fidelity approach for physics-based space weather modeling based on the solution of the Navier-Stokes equations using a high-order discontinuous Galerkin method, combined with a matrix-free strategy suitable for high-performance computing on GPU architectures. The approach consists of a thermospheric model that describes a chemically frozen neutral atmosphere in non-hydrostatic equilibrium driven by the external excitation of the Sun. A novel set of variables is considered to treat the low densities present in the upper atmosphere and to accommodate the wide range of scales present in the problem. At the same time, and unlike most existing approaches, radial and angular directions are treated in a non-segregated approach. The study presents a set of numerical examples that demonstrate the accuracy of the approximation and validate the current approach against observational data along a satellite orbit, including estimates of established empirical and physics-based models of the ionosphere-thermosphere system. Finally, a 1D radial derivation of the physics-based model is presented and utilized for conducting a parametric study of the main thermal quantities under various solar conditions.

We present the proper orthogonal descriptors for efficient and accuracy representation of the potential energy surface. The potential energy surface is represented as a many-body expansion of parametrized potentials in which the potentials are functions of atom positions and parameters. The proper orthogonal decomposition is employed to decompose the parametrized potentials into a set of proper orthogonal descriptors (PODs). Because of the rapid convergence of the proper orthogonal decomposition, relevant snapshots can be sampled exhaustively to represent the atomic neighborhood environment accurately with a small number of descriptors. The proper orthogonal descriptors are used to develop interatomic potentials by using a linear expansion of the descriptors and determining the expansion coefficients from a weighted least-squares regression against a density functional theory (DFT) training set. We present a comprehensive evaluation of the POD potentials on previously published DFT data sets comprising Li, Mo, Cu, Ni, Si, Ge, and Ta elements. The data sets represent a diverse pool of metals, transition metals, and semiconductors. The accuracy of the POD potentials are comparable to that of state-of-the-art machine learning potentials such as the spectral neighbor analysis potential (SNAP) and the atomic cluster expansion (ACE).

This contribution motivates the need to investigate new methods to address the lightning protection of novel transport aircraft, that look radically different to the conventional wing-tube-configuration, and for which historical information is not available. Physics-based models can not only fill this gap but also present an opportunity to include the lightning protection constraints early on in the design phase. In addition, they can lead to innovations in the field. Along these lines, a physics-based numerical tool is presented to automatically estimate the first and second leader attachment points for any aircraft geometry. The model is also used to investigate a charge control strategy to reduce the risk of a strike. Making this vision a reality will require revisiting some of the fundamentals for leader inception and propagation, as well as the differences between positive and negative discharges.

This paper describes the application of an implicit matrix-free discontinuous Galerkin method for hypersonic flow in chemical non-equilibrium. An implicit high-order method is coupled with an open-source library providing the physical modeling terms relevant to hypersonic flow. The main focus of study is the development of robust artificial viscosity methods to provide stability when sharp features appear in the solution. These techniques have been extensively studied and refined for ideal gas flows, and this work adapts two different approaches for non-equilibrium reacting flows. Initial validation results of the coupling between the high-order method and hypersonics modeling codes are presented and results are presented on a 1D shock tube and 2D cylinder problem.