High-order Discontinuous Galerkin Discretizations Research Articles

Domain-specific implementation of high-order Discontinuous Galerkin methods in spherical geometry

In recent years, domain-specific languages (DSLs) have achieved significant success in large-scale efforts to reimplement existing meteorological models in a performance portable manner. The dynamical cores of these models are based on finite difference and finite volume schemes, and existing DSLs are generally limited to supporting only these numerical methods. In the meantime, there have been numerous attempts to use high-order Discontinuous Galerkin (DG) methods for atmospheric dynamics, which are currently largely unsupported in main-stream DSLs. In order to link these developments, we present two domain-specific languages which extend the existing GridTools (GT) ecosystem to high-order DG discretization. The first is a C++-based DSL called G4GT, which, despite being no longer supported, gave us the impetus to implement extensions to the subsequent Python-based production DSL called GT4Py to support the operations needed for DG solvers. As a proof of concept, the shallow water equations in spherical geometry are implemented in both DSLs, thus providing a blueprint for the application of domain-specific languages to the development of global atmospheric models. We believe this is the first GPU-capable DSL implementation of DG in spherical geometry. The results demonstrate that a DSL designed for finite difference/volume methods can be successfully extended to implement a DG solver, while preserving the performance-portability of the DSL.

High-Order Finite Element Second Moment Methods for Linear Transport

We present high-order, finite element–based Second Moment Methods (SMMs) for solving radiation transport problems in two spatial dimensions. We leverage the close connection between the Variable Eddington Factor (VEF) method and SMM to convert existing discretizations of the VEF moment system to discretizations of the SMM moment system. The moment discretizations are coupled to a high-order Discontinuous Galerkin discretization of the SN transport equations. We show that the resulting methods achieve high-order accuracy on high-order (curved) meshes, preserve the thick diffusion limit, and are effective on a challenging multimaterial problem both in outer fixed-point iterations and in inner preconditioned iterative solver iterations for the discrete moment systems. We also present parallel scaling results and provide direct comparisons to the VEF algorithms from which the SMM algorithms were derived.

Analysis on physical-constraint-preserving high-order discontinuous Galerkin method for solving Kapila's five-equation model

This work focuses on the analysis and development of a high-order physical-constraint-preserving discontinuous Galerkin (DG) method for compressible two-fluid flows by solving Kapila's five-equation model with the ideal gas equation of state. The proposed method is proved to possess the following two attractive advantages: first, the high-order DG discretization intrinsically satisfies the mechanical equilibrium criterion around an isolated material interface, and the numerical oscillations due to the use of high-order approximations are well suppressed by using a cell-based characteristic-wise weighted essentially non-oscillatory (WENO) limiter without destroying the mechanical equilibrium criterion; second, it maintains the physical-constraint-preserving property, including the bound-preserving property of volume fractions and positivity-preserving property of partial densities and internal energy. A thorough analysis on the sufficient condition is given to achieve the physical-constraint-preserving property. Several benchmark test cases are examined to validate the accuracy and robustness of the proposed method. Further discussion on the application of the physical-constraint-preserving DG method for solving Kapila's five-equation model with the stiffened gas equation of state is also presented.

Deep reinforcement learning for adaptive mesh refinement

Finite element discretizations of problems in computational physics often rely on adaptive mesh refinement (AMR) to preferentially resolve regions containing important features during simulation. However, these spatial refinement strategies are often heuristic and rely on domain-specific knowledge or trial-and-error. We treat the process of adaptive mesh refinement as a local, sequential decision-making problem under incomplete information, formulating AMR as a partially observable Markov decision process. Using a deep reinforcement learning (RL) approach, we train policy networks for AMR strategy directly from numerical simulation. The training process does not require an exact solution or a high-fidelity ground truth to the partial differential equation (PDE) at hand, nor does it require a pre-computed training dataset. The local nature of our deep RL (DRL) allows the policy network to be trained inexpensively on much smaller problems than those on which they are deployed. The new DRL-AMR method is not specific to any particular PDE, problem dimension, or numerical discretization. The RL policy networks, trained on simple examples, can generalize to more complex problems, and can flexibly incorporate diverse problem physics. To that end, we apply the method to a range of PDEs, using a variety of high-order discontinuous Galerkin and hybridizable discontinuous Galerkin finite element discretizations. We show that the resultant DRL policies are competitive with common AMR heuristics and strike a favorable balance between accuracy and cost such that they often lead to a higher accuracy per problem degree of freedom, and are effective across a wide class of PDEs and problems.

High-Order Mixed Finite Element Variable Eddington Factor Methods

We apply high-order mixed finite element discretization techniques and their associated preconditioned iterative solvers to the Variable Eddington Factor (VEF) equations in two spatial dimensions. The mixed finite element VEF discretizations are coupled to a high-order Discontinuous Galerkin (DG) discretization of the discrete ordinates transport equation to form effective linear transport algorithms that are compatible with high-order (curved) meshes. This combination of VEF and transport discretizations is motivated by the use of high-order mixed finite element methods in hydrodynamics calculations at the Lawrence Livermore National Laboratory (LLNL). Due to the mathematical structure of the VEF equations, the standard Raviart Thomas (RT) mixed finite elements cannot be used to approximate the vector variable in the VEF equations. Instead, we investigate three alternatives based on the use of continuous finite elements for each vector component, a non-conforming RT approach where DG-like techniques are used, and a hybridized RT method. We present numerical results that demonstrate high-order accuracy, compatibility with curved meshes, and robust and efficient convergence in iteratively solving the coupled transport-VEF system and in the preconditioned linear solvers used to invert the discretized VEF equations.

Assessment of an Implicit Discontinuous Galerkin Solver for Incompressible Flow Problems with Variable Density

Multi-component flow problems are typical of many technological and engineering applications. In this work, we propose an implicit high-order discontinuous Galerkin discretization of the variable density incompressible (VDI) flow model for the simulation of multi-component problems. Indeed, the peculiarity of the VDI model is that the density is treated as an advected property, which can be used to possibly track multiple (more than two) components. The interface between fluids is described by a smooth, but sharp, variation in the density field, thus not requiring any geometrical reconstruction. Godunov numerical fluxes, density positivity, mass conservation, and Gibbs-type phenomena at material interfaces are challenges that are considered during the numerical approach development. To avoid Courant-related time step restrictions, high-order single-step multi-stage implicit schemes are applied for the temporal integration. Several test cases with known analytical solutions are used to assess the current approach in terms of space, time, and mass conservation accuracy. As a challenging application, the simulation of a 2D droplet impinging on a thin liquid film is performed and shows the capabilities of the proposed DG approach when dealing with high-density (water–air) multi-component problems.

Comparative study of inner–outer Krylov solvers for linear systems in structured and high-order unstructured CFD problems

Advanced Krylov subspace methods are investigated for the solution of large sparse linear systems arising from stiff adjoint-based aerodynamic shape optimization problems. A special attention is paid to the flexible inner–outer GMRES strategy combined with most relevant preconditioning and deflation techniques. The choice of this specific class of Krylov solvers for challenging problems is based on its outstanding convergence properties. Typically in our implementation the efficiency of the preconditioner is enhanced with a domain decomposition method with overlapping. However, maintaining the performance of the preconditioner may be challenging since scalability and efficiency of a preconditioning technique are properties often antagonistic to each other. In this paper we demonstrate how flexible inner–outer Krylov methods are able to overcome this critical issue. A numerical study is performed considering either a Finite Volume (FV), or a high-order Discontinuous Galerkin (DG) discretization which affect the arithmetic intensity and memory-bandwidth of the algebraic operations. We consider test cases of transonic turbulent flows with RANS modeling over the two-dimensional supercritical OAT15A airfoil and the three-dimensional ONERA M6 wing. Benefits in terms of robustness and convergence compared to standard GMRES solvers are obtained. Strong scalability analysis shows satisfactory results. Based on these representative problems a discussion of the recommended numerical practices is proposed.

Efficient hyperreduction of high-order discontinuous Galerkin methods: Element-wise and point-wise reduced quadrature formulations

We develop and assess projection-based model reduction methods for high-order discontinuous Galerkin (DG) discretizations of parametrized nonlinear partial differential equations with applications in aerodynamics. Our emphasis is on the choice of hyperreduction methods. We analyze computational complexity and use numerical examples to show that typical hyperreduction methods based on element-wise or index-wise sampling of the residual vector cannot effectively reduce high-order DG discretizations, with a large number of degrees of freedom and quadrature points per element and conversely a small number of elements. To overcome this limitation, we devise a hyperreduction method for high-order DG methods that samples individual quadrature points, instead of elements, to provide a finer decomposition of the residual. We compare the formulation, implementation, and computational complexity of the element-wise and point-wise formulations. We then numerically assess the two formulations using Euler flow over a shape-parametrized airfoil and laminar and turbulent Navier-Stokes flow over a three-dimensional aerodynamic body. In each case, we observe that the point-wise hyperreduced model reduces the computational time by several orders of magnitude relative to the high-order DG method and a few orders of magnitude relative to the element-wise hyperreduced model. In addition, the computational performance of the point-wise formulation does not deteriorate with the DG polynomial degree, unlike for the element-wise formulation.

An entropy stable scheme for the non-linear Boltzmann equation

The kinetic equations govern the behavior of gaseous flows, compressibility, turbulence, reactions with internal energy exchange, that form the essence of numerous physical processes. Despite their wide applicability, their six-dimensional (seven including time) nature presents a huge computational challenge. With the advent of the modern high performance computing systems and few recent advances in numerical methods, it is now possible to numerically study the behavior of these systems. However, to understand the instability of rich molecular processes and find a structure in the apparent chaos, a scheme that is efficient in multi-dimensions, exhibits high parallel-efficiency, and foremost produces an entropy solution as per the theory of solutions of hyperbolic systems may prove useful. In this work, first, we construct an entropy stable flux for non-linear inhomogeneous (full) Boltzmann equation. Second, to ensure geometrical flexibility, we couple the scheme with a class of high order discontinuous Galerkin discretization (Jaiswal 2019 [41]) which satisfy summation-by-part (SBP) discretely. Third, utilizing SBP, we prove that the resulting semi-discrete scheme is locally and globally conservative in flat spaces; preserves the entropy-decay (H-theorem) property; efficient and simple; and therefore suitable for dealing with highly complex non-smooth flow problems. Fourth, we show that the fully-discrete kinetic scheme, utilizing an implicit-explicit time-discretization, in the limit of vanishing Knudsen number, becomes an entropy-stable explicit scheme applied to the Euler system. Fifth, we carry out a series of verification tests to illustrate the stability, accuracy, and conservation properties of the proposed method. These tests involve over 85 million degrees of freedom and over 50 billion operations per time-step.

Large-Eddy Simulation of Transonic Buffet Using Matrix-Free Discontinuous Galerkin Method

We present an implicit large-eddy simulation of transonic buffet over the OAT15A supercritical airfoil at Mach number 0.73, angle of attack 3.5 deg, and Reynolds number . The simulation is performed using a matrix-free discontinuous Galerkin (DG) method and a diagonally implicit Runge–Kutta scheme on graphics processor units. We propose a Jacobian-free Newton–Krylov method to solve nonlinear systems arising from the discretization of the Navier–Stokes equations. The method successfully predicts the buffet onset, the buffet frequency, and turbulence statistics owing to the high-order DG discretization and an efficient mesh refinement for the laminar and turbulent boundary layers. A number of physical phenomena present in the experiment are captured in our simulation, including periodical low-frequency oscillations of shock wave in the streamwise direction, strong shear layer detached from the shock wave due to shock-wave/boundary-layer interaction and small-scale structures broken down by the shear-layer instability in the transition region, and shock-induced flow separation. The pressure coefficient, the root mean square of the fluctuating pressure, and the streamwise range of the shock wave oscillation agree well with the experimental data. The results suggest that the proposed method can accurately predict the onset of turbulence and buffet phenomena at high Reynolds numbers without a subgrid scale model or a wall model.

Aeroacoustic Instability Analysis of Helmholtz Resonators Through High-Order Unsteady Reynolds-Averaged Navier–Stokes Simulations

This paper presents a numerical analysis of the aeroacoustic instabilities in a slit-type Helmholtz resonator subjected to a low Mach number grazing flow. This configuration can be susceptible to nonlinear feedback phenomena due to the mutual coupling between the aerodynamic and acoustic fields in the proximity of the slit. The aeroacoustic flowfield is simulated through a direct approach by solving the unsteady compressible Reynolds-averaged Navier–Stokes equations using a high-order discontinuous Galerkin discretization. The response of the Helmholtz resonator to an incident acoustic field is analyzed using a two-port characterization technique. The onset of the aeroacoustic instabilities is investigated by applying an acoustic energy criterion based on the scattering matrix formulation. The Nyquist stability diagram is studied to determine the system stability and identify the whistling occurrence. The two- and three-dimensional numerical results are compared to experimental data obtained from a dedicated measurement campaign. The comparison between the numerical and experimental data demonstrates a very good agreement, highlighting the potential of the applied numerical approach for the investigation of the aeroacoustic instabilities. This allows a deep description of the aeroacoustic feedback mechanisms leading to the whistling onset, proposing a valid engineering tool for the early design stage.

A quadratic programming flux correction method for high-order DG discretizations of SN transport

We present a new flux-fixup approach for arbitrarily high-order discontinuous Galerkin (DG) discretizations of the SN transport equation, and we demonstrate the compatibility of this approach with the Variable Eddington Factor (VEF) method [1,2]. The new fixup approach is sweep-compatible: during a transport sweep (block Gauss-Seidel iteration in which the scattering source is lagged), a local quadratic programming (QP) problem is solved in each spatial element to ensure that the solution satisfies certain physical constraints, including local particle balance. In this paper, we describe two choices of physical constraints, resulting in two variants of the method: QP Zero (QPZ) and QP Maximum Principle (QPMP). In QPZ, the finite element coefficients of the solution are constrained to be nonnegative. In QPMP, they are constrained to adhere to an approximate discrete maximum principle.There are two primary takeaways in this paper. First, when the positive Bernstein basis is used for DG discretization, the QPMP method eliminates negativities, preserves high-order accuracy for smooth problems, and significantly dampens unphysical oscillations in the solution. The latter feature – the dampening of unphysical oscillations – is an improvement upon standard, simpler fixup approaches such as the approach described in [3] (denoted as the “zero and rescale” (ZR) method in this paper). This improvement comes at a moderate computational cost, but it is not prohibitive. Our results show that, even in an unrealistic worst-case scenario where 83% of the spatial elements require a fixup, the computational cost of performing a transport sweep with fixup is only ∼31% greater than performing one without fixup.The second takeaway is that the VEF method can be used to accelerate the convergence of transport sweeps even when a fixup is applied. When optically thick regions are present, transport sweeps converge slowly, regardless of whether a fixup is applied, and acceleration is needed. However, attempting to apply standard diffusion synthetic acceleration (DSA) to fixed-up transport sweeps results in divergence for optically thick problems. Our results show that the same is not true for VEF. When VEF is combined with fixed-up transport sweeps, the result is a scheme that produces a nonnegative solution, converges independently of the mean free path, and, in the case of the QPMP fixup, adheres to an approximate discrete maximum principle.

Hybrid multigrid methods for high-order discontinuous Galerkin discretizations

The present work develops hybrid multigrid methods for high-order discontinuous Galerkin discretizations of elliptic problems, which are, for example, a key ingredient of incompressible flow solvers in the field of computational fluid dynamics. Fast matrix-free operator evaluation on tensor product elements is used to devise a computationally efficient PDE solver. The multigrid hierarchy exploits all possibilities of geometric, polynomial, and algebraic coarsening, targeting engineering applications on complex geometries. Additionally, a transfer from discontinuous to continuous function spaces is performed within the multigrid hierarchy. This does not only further reduce the problem size of the coarse-grid problem, but also leads to a discretization most suitable for state-of-the-art algebraic multigrid methods applied as coarse-grid solver. The relevant design choices regarding the selection of optimal multigrid coarsening strategies among the various possibilities are discussed with the metric of computational costs as the driving force for algorithmic selections. We find that a transfer to a continuous function space at highest polynomial degree (or on the finest mesh), followed by polynomial and geometric coarsening, shows the best overall performance. The success of this particular multigrid strategy is due to a significant reduction in iteration counts as compared to a transfer from discontinuous to continuous function spaces at lowest polynomial degree (or on the coarsest mesh). The coarsening strategy with transfer to a continuous function space on the finest level leads to a multigrid algorithm that is robust with respect to the penalty parameter of the symmetric interior penalty method. Detailed numerical investigations are conducted for a series of examples ranging from academic test cases to more complex, practically relevant geometries. Performance comparisons to state-of-the-art methods from the literature demonstrate the versatility and computational efficiency of the proposed multigrid algorithms.

Implicit shock tracking using an optimization-based high-order discontinuous Galerkin method

A novel framework for resolving discontinuous solutions of conservation laws, e.g., contact lines, shock waves, and interfaces, using implicit tracking and a high-order discontinuous Galerkin (DG) discretization was introduced in [39]. Central to the framework is an optimization problem whose solution is a discontinuity-aligned mesh and the corresponding high-order approximation to the flow that does not require explicit meshing of the unknown discontinuity surface. The method was shown to deliver highly accurate solutions on coarse, high-order discretizations without nonlinear stabilization and recover optimal convergence rates O(hp+1) even for problems with discontinuous solutions. This work extends the implicit tracking framework such that robustness is improved and convergence accelerated. In particular, we introduce an improved formulation of the central optimization problem and an associated sequential quadratic programming (SQP) solver. The new error-based objective function penalizes violation of the DG residual in an enriched test space and is shown to have excellent tracking properties. The SQP solver simultaneously converges the nodal coordinates of the mesh and DG solution to their optimal values and is equipped with a number of features to ensure robust, fast convergence: Levenberg-Marquardt approximation of the Hessian with weighted elliptic regularization, backtracking line search based on the ℓ1 merit function, and rigorous convergence criteria. We use the proposed method to solve a range of inviscid conservation laws of varying difficulty. We show the method is able to deliver accurate solutions on coarse, high-order meshes and the SQP solver is robust and usually able to drive the first-order optimality system to tight tolerances.

Algebraic Multigrid Schemes for High-Order Nodal Discontinuous Galerkin Methods

We present algebraic multigrid (AMG) methods for the efficient solution of the linear system of equations stemming from high-order discontinuous Galerkin (DG) discretizations of second-order elliptic problems. For DG methods, standard multigrid approaches cannot be employed because of redundancy of the degrees of freedom associated to the same grid point. We present new aggregation procedures and test them in extensive two-dimensional numerical experiments to demonstrate that the proposed AMG method is uniformly convergent with respect to all of the discretization parameters, namely the mesh-size and the polynomial approximation degree.

Modern discontinuous Galerkin methods for the simulation of transitional and turbulent flows in biomedical engineering: A comprehensive LES study of the FDA benchmark nozzle model.

This work uses high-order discontinuous Galerkin discretization techniques to simulate transitional and turbulent flows through medical devices. Flows through medical devices are characterized by moderate Reynolds numbers and typically involve different flow regimes such as laminar, transitional, and turbulent flows. Previous studies for the FDA benchmark nozzle model revealed limitations of Reynolds-averaged Navier-Stokes turbulence models when applied to more complex flow scenarios. Recent works based on large-eddy simulation approaches indicate that these limitations can be overcome but also highlight potential limitations due to a high sensitivity with respect to numerical parameters. The methodology presented in this work introduces two novel ingredients compared with previous studies. Firstly, we use high-order discontinuous Galerkin methods for discretization in space. The inherent dissipation and dispersion properties of high-order discontinuous Galerkin discretizations are expected to render this approach well suited for transitional and turbulent flow simulations. Secondly, to mimic blinded CFD studies, we propose to use a precursor simulation approach in order to predict the inflow boundary condition for laminar, transitional, and turbulent flow regimes instead of prescribing analytical velocity profiles at the inflow. We investigate the whole range of Reynolds numbers as suggested by the FDA benchmark nozzle problem and compare the numerical results to experimental data obtained by particle image velocimetry in order to critically assess the predictive capabilities of the solver on the one hand and the suitability of the FDA nozzle problem as a benchmark in computational fluid dynamics on the other hand.

A GPU accelerated discontinuous Galerkin incompressible flow solver

We present a GPU-accelerated version of a high-order discontinuous Galerkin discretization of the unsteady incompressible Navier–Stokes equations. The equations are discretized in time using a semi-implicit scheme with explicit treatment of the nonlinear term and implicit treatment of the split Stokes operators. The pressure system is solved with a conjugate gradient method together with a fully GPU-accelerated multigrid preconditioner which is designed to minimize memory requirements and to increase overall performance. A semi-Lagrangian subcycling advection algorithm is used to shift the computational load per timestep away from the pressure Poisson solve by allowing larger timestep sizes in exchange for an increased number of advection steps. Numerical results confirm we achieve the design order accuracy in time and space. We optimize the performance of the most time-consuming kernels by tuning the fine-grain parallelism, memory utilization, and maximizing bandwidth. To assess overall performance we present an empirically calibrated roofline performance model for a target GPU to explain the achieved efficiency. We demonstrate that, in the most cases, the kernels used in the solver are close to their empirically predicted roofline performance.

Comparison of hybrid and standard discontinuous Galerkin methods in a mesh-optimisation setting

This paper assesses the benefits of hybridization on the accuracy and efficiency of high-order discontinuous Galerkin (DG) discretizations. Two hybridized methods are considered in addition to DG: hybridized DG (HDG) and embedded DG (EDG). These methods offer memory and computational time savings by introducing trace degrees of freedom on faces that become the only globally-coupled unknowns. To mitigate the effects of solution singularities on accuracy, the methods are compared in an adaptive setting on meshes optimised for the accurate prediction of chosen scalar outputs. Compressible flow results for the Euler and Reynolds-averaged Navier-Stokes equations demonstrate that the hybridized methods offer cost savings relative to DG in memory and computational time. In addition, for the cases tested, EDG yields the lowest error levels for a given number of degrees of freedom. These benefits disappear on uniformly-refined meshes, indicating the importance of using order-optimised meshes when comparing the discretizations.

Evaluating high order discontinuous Galerkin discretization of the Boltzmann collision integral in <inline-formula><tex-math id="M1">\begin{document}$ \mathcal{O}(N^2) $\end{document}</tex-math></inline-formula> operations using the discrete fourier transform

We present a numerical algorithm for evaluating the Boltzmann collision operator with \begin{document}$O(N^2)$\end{document} operations based on high order discontinuous Galerkin discretizations in the velocity variable. To formulate the approach, Galerkin projection of the collision operator is written in the form of a bilinear circular convolution. An application of the discrete Fourier transform allows to rewrite the six fold convolution sum as a three fold weighted convolution sum in the frequency space. The new algorithm is implemented and tested in the spatially homogeneous case, and results in a considerable improvement in speed as compared to the direct evaluation. Split and non-split forms of the collision operator are considered, which are forms of the collision operator that have separate and simultaneous evaluations of the gain and loss terms, respectively. Smaller numerical errors are observed in the conserved quantities in simulations using the non-split form.

Discontinuous Galerkin scheme for the spherical shallow water equations with applications to tsunami modeling and prediction

We present a novel high-order discontinuous Galerkin discretization for the spherical shallow water equations, able to handle wetting/drying and non-conforming, curved meshes in a well-balanced manner. This requires a well-balanced discretization, that cannot rely on exact quadrature, due to the curved mesh. Using the strong form of the discontinuous Galerkin discretization, we achieve a splitting of the well-balanced condition into individual problems for the flux and volume terms, which has significant advantages: It allows for the construction of non-conforming, well-balanced flux discretizations, i.e. we can perform non-conforming mesh refinement while preserving the well-balanced property of the scheme. More importantly, this approach enables the development of a new method for handling wet/dry transitions. In contrast to other wetting/drying methods, it is well-balanced and able to handle wetting/drying robustly at any polynomial order, without the introduction of physical model assumptions such as viscosity, artificial porosity or cancellation of gravity.We perform a series of one-dimensional tests and analyze the properties of our scheme. In order to validate our method for the simulation of large-scale tsunami events on the rotating sphere, we perform numerical simulations of the 2011 Tohoku tsunami and compare our results to real-world buoy data. The method is able to predict arrival times and wave amplitudes accurately even over long distances. This indicates that our method accurately captures all physical phenomena relevant to the long-term evolution of tsunami waves.