We build a multi-element variant of the smoothness increasing accuracy conserving (SIAC) shock capturing technique proposed for single element spectral methods by Wissink et al. (J Sci Comput 77:579–596, 2018). In particular, the baseline scheme of our method is the nodal discontinuous Galerkin spectral element method (DGSEM) for approximating the solution of systems of conservation laws. It is well known that high-order methods generate spurious oscillations near discontinuities which can develop in the solution for nonlinear problems, even when the initial data is smooth. We propose a novel multi-element SIAC filtering technique applied to the DGSEM as a shock capturing method. We design the SIAC filtering such that the numerical scheme remains high-order accurate and that the shock capturing is applied adaptively throughout the domain. The shock capturing method is derived for general systems of conservation laws. We apply the novel SIAC filter to the two-dimensional Euler and ideal magnetohydrodynamics equations to several standard test problems with a variety of boundary conditions.
The first paper of this series presents a discretely entropy stable discontinuous Galerkin (DG) method for the resistive magnetohydrodynamics (MHD) equations on three-dimensional curvilinear unstructured hexahedral meshes. Compared to other fluid dynamics systems such as the shallow water equations or the compressible Navier-Stokes equations, the resistive MHD equations need special considerations because of the divergence-free constraint on the magnetic field. For instance, it is well known that for the symmetrization of the ideal MHD system as well as the continuous entropy analysis a non-conservative term proportional to the divergence of the magnetic field, typically referred to as the Powell term, must be included. As a consequence, the mimicry of the continuous entropy analysis in the discrete sense demands a suitable DG approximation of the non-conservative terms in addition to the ideal MHD terms.
This paper focuses on the resistive MHD equations: Our first contribution is a proof that the resistive terms are symmetric and positive-definite when formulated in entropy space as gradients of the entropy variables, which enables us to show that the entropy inequality holds for the resistive MHD equations. This continuous analysis is the key for our DG discretization and guides the path for the construction of an approximation that discretely mimics the entropy inequality, typically termed entropy stability. Our second contribution is a detailed derivation and analysis of the discretization on three-dimensional curvilinear meshes. The discrete analysis relies on the summation-by-parts property, which is satisfied by the DG spectral element method (DGSEM) with Legendre-Gauss-Lobatto (LGL) nodes. Although the divergence-free constraint is included in the non-conservative terms, the resulting method has no particular treatment of the magnetic field divergence errors, which might pollute the solution quality. Our final contribution is the extension of the standard resistive MHD equations and our DG approximation with a divergence cleaning mechanism that is based on a generalized Lagrange multiplier (GLM).
As a conclusion to the first part of this series, we provide detailed numerical validations of our DGSEM method that underline our theoretical derivations. In addition, we show a numerical example where the entropy stable DGSEM demonstrates increased robustness compared to the standard DGSEM.
This article serves as a summary outlining the mathematical entropy analysis of the ideal magnetohydrodynamic (MHD) equations. We select the ideal MHD equations as they are particularly useful for mathematically modeling a wide variety of magnetized fluids. In order to be self-contained we first motivate the physical properties of a magnetic fluid and how it should behave under the laws of thermodynamics. Next, we introduce a mathematical model built from hyperbolic partial differential equations (PDEs) that translate physical laws into mathematical equations. After an overview of the continuous analysis, we thoroughly describe the derivation of a numerical approximation of the ideal MHD system that remains consistent to the continuous thermodynamic principles. The derivation of the method and the theorems contained within serve as the bulk of the review article. We demonstrate that the derived numerical approximation retains the correct entropic properties of the continuous model and show its applicability to a variety of standard numerical test cases for MHD schemes. We close with our conclusions and a brief discussion on future work in the area of entropy consistent numerical methods and the modeling of plasmas.
Entropy stable schemes can be constructed with a specific choice of the numerical flux function. First, an entropy conserving flux is constructed. Secondly, an entropy stable dissipation term is added to this flux to guarantee dissipation of the discrete entropy. Present works in the field of entropy stable numerical schemes are concerned with thorough derivations of entropy conservative fluxes for ideal MHD. However, as we show in this work, if the dissipation operator is not constructed in a very specific way, it cannot lead to a generally stable numerical scheme. The two main findings presented in this paper are that the entropy conserving flux of Ismail & Roe can easily break down for certain initial conditions commonly found in astrophysical simulations, and that special care must be taken in the derivation of a discrete dissipation matrix for an entropy stable numerical scheme to be robust. We present a convenient novel averaging procedure to evaluate the entropy Jacobians of the ideal MHD and the compressible Euler equations that yields a discretization with favorable robustness properties.
We describe a high-order numerical magnetohydrodynamics (MHD) solver built upon a novel non-linear entropy stable numerical flux function that supports eight travelling wave solutions. By construction the solver conserves mass, momentum, and energy and is entropy stable. The method is designed to treat the divergence-free constraint on the magnetic field in a similar fashion to a hyperbolic divergence cleaning technique. The solver described herein is especially well-suited for flows involving strong discontinuities. Furthermore, we present a new formulation to guarantee positivity of the pressure. We present the underlying theory and implementation of the new solver into the multi-physics, multi-scale adaptive mesh refinement (AMR) simulation code FLASH (http://flash.uchicago.edu). The accuracy, robustness and computational efficiency is demonstrated with a number of tests, including comparisons to available MHD implementations in FLASH.
The paper presents two contributions in the context of the numerical simulation of magnetized fluid dynamics. First, we show how to extend the ideal magnetohydrodynamics (MHD) equations with an inbuilt magnetic field divergence cleaning mechanism in such a way that the resulting model is consistent with the second law of thermodynamics. As a byproduct of these derivations, we show that not all of the commonly used divergence cleaning extensions of the ideal MHD equations are thermodynamically consistent. Secondly, we present a numerical scheme obtained by constructing a specific finite volume discretization that is consistent with the discrete thermodynamic entropy. It includes a mechanism to control the discrete divergence error of the magnetic field by construction and is Galilean invariant. We implement the new high-order MHD solver in the adaptive mesh refinement code FLASH where we compare the divergence cleaning efficiency to the constrained transport solver available in FLASH (unsplit staggered mesh scheme).
Non-conforming numerical approximations offer increased flexibility for applications that require high resolution in a localized area of the computational domain or near complex geometries. Two key properties for non-conforming methods to be applicable to real world applications are conservation and energy stability. The summation-by-parts (SBP) property, which certain finite-difference and discontinuous Galerkin methods have, finds success for the numerical approximation of hyperbolic conservation laws, because the proofs of energy stability and conservation can discretely mimic the continuous analysis of partial differential equations. In addition, SBP methods can be developed with high-order accuracy, which is useful for simulations that contain multiple spatial and temporal scales. However, existing non-conforming SBP schemes result in a reduction of the overall degree of the scheme, which leads to a reduction in the order of the solution error. This loss of degree is due to the particular interface coupling through a simultaneous-approximation-term (SAT). We present in this work a novel class of SBP-SAT operators that maintain conservation, energy stability, and have no loss of the degree of the scheme for non-conforming approximations. The new degree preserving discretizations require an ansatz that the norm matrix of the SBP operator is of a degree ≥ 2p, in contrast to, for example, existing finite difference SBP operators, where the norm matrix is 2p − 1 accurate. We demonstrate the fundamental properties of the new scheme with rigorous mathematical analysis as well as numerical verification.
This work examines the development of an entropy conservative (for smooth solutions) or entropy stable (for discontinuous solutions) space-time discontinuous Galerkin (DG) method for systems of nonlinear hyperbolic conservation laws. The resulting numerical scheme is fully discrete and provides a bound on the mathematical entropy at any time according to its initial condition and boundary conditions. The crux of the method is that discrete derivative approximations in space and time are summation-by-parts (SBP) operators. This allows the discrete method to mimic results from the continuous entropy analysis and ensures that the complete numerical scheme obeys the second law of thermodynamics. Importantly, the novel method described herein does not assume any exactness of quadrature in the variational forms that naturally arise in the context of DG methods. Typically, the development of entropy stable schemes is done on the semidiscrete level ignoring the temporal dependence. In this work, we demonstrate that creating an entropy stable DG method in time is similar to the spatial discrete entropy analysis, but there are important (and subtle) differences. Therefore, we highlight the temporal entropy analysis throughout this work. For the compressible Euler equations, the preservation of kinetic energy is of interest besides entropy stability. The construction of kinetic energy preserving (KEP) schemes is, again, typically done on the semidiscrete level similar to the construction of entropy stable schemes. We present a generalization of the KEP condition from Jameson to the space-time framework and provide the temporal components for both entropy stability and kinetic energy preservation. The properties of the space-time DG method derived herein are validated through numerical tests for the compressible Euler equations. Additionally, we provide, in appendices, how to construct the temporal entropy stable components for the shallow water or ideal magnetohydrodynamic (MHD) equations.
This work presents an entropy stable discontinuous Galerkin (DG) spectral element approximation for systems of non-linear conservation laws with general geometric (h) and polynomial order (p) non-conforming rectangular meshes. The crux of the proofs presented is that the nodal DG method is constructed with the collocated Legendre–Gauss–Lobatto nodes. This choice ensures that the derivative/mass matrix pair is a summation-by-parts (SBP) operator such that entropy stability proofs from the continuous analysis are discretely mimicked. Special attention is given to the coupling between non-conforming elements as we demonstrate that the standard mortar approach for DG methods does not guarantee entropy stability for non-linear problems, which can lead to instabilities. As such, we describe a precise procedure and modify the mortar method to guarantee entropy stability for general non-linear hyperbolic systems on h / p non-conforming meshes. We verify the high-order accuracy and the entropy conservation/stability of fully non-conforming approximation with numerical examples.
We show how to modify the original Bassi and Rebay scheme (BR1) [F. Bassi and S. Rebay, A High Order Accurate Discontinuous Finite Element Method for the Numerical Solution of the Compressible Navier-Stokes Equations, Journal of Computational Physics, 131:267–279, 1997] to get a provably stable discontinuous Galerkin collocation spectral element method (DGSEM) with Gauss-Lobatto (GL) nodes for the compressible Navier-Stokes equations (NSE) on three dimensional curvilinear meshes.
Specifically, we show that the BR1 scheme can be provably stable if the metric identities are discretely satisfied, a two-point average for the metric terms is used for the contravariant fluxes in the volume, an entropy conserving split form is used for the advective volume integrals, the auxiliary gradients for the viscous terms are computed from gradients of entropy variables, and the BR1 scheme is used for the interface fluxes.
Our analysis shows that even with three dimensional curvilinear grids, the BR1 fluxes do not add artificial dissipation at the interior element faces. Thus, the BR1 interface fluxes preserve the stability of the discretization of the advection terms and we get either energy stability or entropy-stability for the linear or nonlinear compressible NSE, respectively.
In this work, we design an arbitrary high order accurate nodal discontinuous Galerkin spectral element type method for the one dimensional shallow water equations. The novel method uses a skew-symmetric formulation of the continuous problem. We prove that this discretisation exactly preserves the local mass and momentum. Furthermore, we show that combined with a special numerical interface flux function, the method exactly preserves the entropy, which is also the total energy for the shallow water equations. Finally, we prove that the surface fluxes, the skew-symmetric volume integrals, and the source term are well balanced. Numerical tests are performed to demonstrate the theoretical findings.
Fisher and Carpenter (High-order entropy stable finite difference schemes for non-linear conservation laws: Finite domains, Journal of Computational Physics, 252:518–557, 2013) found a remarkable equivalence of general diagonal norm high-order summation-by- parts operators to a subcell based high-order finite volume formulation. This equivalence enables the construction of provably entropy stable schemes by a specific choice of the sub-cell finite volume flux. We show that besides the construction of entropy stable high order schemes, a careful choice of subcell finite volume fluxes generates split formulations of quadratic or cubic terms. Thus, by changing the subcell finite volume flux to a specific choice, we are able to generate, in a systematic way, all common split forms of the compressible Euler advection terms, such as the Ducros splitting and the Kennedy and Gruber splitting. Although these split forms are not entropy stable, we present a systematic way to prove which of those split forms are at least kinetic energy preserving. With this, we show we construct a unified high-order split form DG framework. We investigate with three dimensional numerical simulations of the inviscid Taylor-Green vortex and show that the new split forms enhance the robustness of high order simulations in comparison to the standard scheme when solving turbulent vortex dominated flows. In fact, we show that for certain test cases, the novel split form discontinuous Galerkin schemes are more robust than the discontinuous Galerkin scheme with over-integration.
We design a novel provably stable discontinuous Galerkin spectral element (DGSEM) approximation to solve systems of conservation laws on moving domains. To incorporate the motion of the domain, we use an arbitrary Lagrangian-Eulerian formulation to map the governing equations to a fixed reference domain. The approximation is made stable by a discretization of a skew-symmetric formulation of the problem. We prove that the discrete approximation is stable, conservative and, for constant coefficient problems, maintains the free- stream preservation property. We also provide details on how to add the new skew-symmetric ALE approximation to an existing discontinuous Galerkin spectral element code. Lastly, we provide numerical support of the theoretical results.
We compare and contrast information provided by the energy analysis of Kreiss and the entropy theory of Tadmor for systems of nonlinear hyperbolic conservation laws. The two-dimensional nonlinear shallow water equations are used to highlight the similarities and differences since the total energy of the system is a mathematical entropy function. We demonstrate that the classical energy method is consistent with the entropy analysis, but significantly more fundamental as it guides proper boundary treatments. In particular, the energy analysis provides information on what type of and how many boundary conditions are required, which is lacking in the entropy analysis. For the shallow water system we determine the number and the type of boundary conditions needed for subcritical and supercritical flows on a general domain. As eigenvalues are augmented in the nonlinear analysis, we find that a flow may be classified as subcritical, but the treatment of the boundary resembles that of a supercritical flow. Because of this, we show that the nonlinear energy analysis leads to a different number of boundary conditions compared with the linear energy analysis. We also demonstrate that the entropy estimate leads to erroneous boundary treatments by over specifying and/or under specifying boundary data causing the loss of existence and/or energy bound, respectively. Our analysis reveals that the nonlinear energy analysis is the only one that provides an estimate for open boundaries. Both the entropy and linear energy analysis fail.
It is known that HLL-type schemes are more dissipative than schemes based on characteristic decompositions. However, HLL-type methods offer greater flexibility to large systems of hyperbolic conservation laws because the eigenstructure of the flux Jacobian is not needed. We demonstrate in the present work that several HLL-type Riemann solvers are provably entropy stable. Further, we provide convex combinations of standard dissipation terms to create hybrid HLL-type methods that have less dissipation while retaining entropy stability. The decrease in dissipation is demonstrated for the ideal MHD equations with a numerical example.
We design an arbitrary high-order accurate nodal discontinuous Galerkin spectral element approximation for the non-linear two dimensional shallow water equations with non- constant, possibly discontinuous, bathymetry on unstructured, possibly curved, quadrilateral meshes. The scheme is derived from an equivalent flux differencing formulation of the split form of the equations. We prove that this discretisation exactly preserves the local mass and momentum. Furthermore, combined with a special numerical interface flux function, the method exactly preserves the mathematical entropy, which is the total energy for the shallow water equations. By adding a specific form of interface dissipation to the baseline entropy conserving scheme we create a provably entropy stable scheme. That is, the numerical scheme discretely satisfies the second law of thermodynamics. Finally, with a particular discretisation of the bathymetry source term we prove that the numerical approximation is well-balanced. We provide numerical examples that verify the theoretical findings and furthermore provide an application of the scheme for a partial break of a curved dam test problem.
We extend the entropy stable high order nodal discontinuous Galerkin spectral element approximation for the non-linear two dimensional shallow water equations presented by Wintermeyer et al. [N. Wintermeyer, A. R. Winters, G. J. Gassner, and D. A. Kopriva. An entropy stable nodal discontinuous Galerkin method for the two dimensional shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry. Journal of Computational Physics, 340:200-242, 2017] with a shock capturing technique and a positivity preservation capability to handle dry areas. The scheme preserves the entropy inequality, is well-balanced and works on unstructured, possibly curved, quadrilateral meshes. For the shock capturing, we introduce an artificial viscosity to the equations and prove that the numerical scheme remains entropy stable. We add a positivity preserving limiter to guarantee non-negative water heights as long as the mean water height is non-negative. We prove that non-negative mean water heights are guaranteed under a certain additional time step restriction for the entropy stable numerical interface flux. We implement the method on GPU architectures using the abstract language OCCA, a unified approach to multi-threading languages. We show that the entropy stable scheme is well suited to GPUs as the necessary extra calculations do not negatively impact the runtime up to reasonably high polynomial degrees (around N = 7). We provide numerical examples that challenge the shock capturing and positivity properties of our scheme to verify our theoretical findings.
We outline and extend results for an explicit local time stepping (LTS) strategy designed to operate with the discontinuous Galerkin spectral element method (DGSEM). The LTS procedure is derived from Adams-Bashforth multirate time integration methods. The new results of the LTS method focus on parallelization and reformulation of the LTS integrator to maintain conservation. Discussion is focused on a moving mesh implementation, but the procedures remain applicable to static meshes. In numerical tests, we demonstrate the strong scaling of a parallel, LTS implementation and compare the scaling properties to a parallel, global time stepping (GTS) Runge-Kutta implementation. We also present time-step refinement studies to show that the redesigned, conservative LTS approximations are spectrally accurate in space and have design temporal accuracy.
We describe a unique averaging procedure to design an entropy stable dissipation operator for the ideal magnetohydrodynamic (MHD) and compressible Euler equations. Often in the derivation of an entropy conservative numerical flux function much care is taken in the design and averaging of the entropy conservative numerical flux. We demonstrate in this work that if the discrete dissipation operator is not carefully chosen as well it can have deleterious effects on the numerical approximation. This is particularly true for very strong shocks or high Mach number flows present, for example, in astrophysical simulations. We present the underlying technique of how to construct a unique averaging technique for the discrete dissipation operator. We also demonstrate numerically the increased robustness of the approximation.
In this work, we compare and contrast two provably entropy stable and high-order accurate nodal discontinuous Galerkin spectral element methods applied to the one dimensional shallow water equations for problems with non-constant bottom topography. Of particular importance for numerical approximations of the shallow water equations is the well-balanced property. The well-balanced property is an attribute that a numerical approximation can preserve a steady-state solution of constant water height in the presence of a bottom topography. Numerical tests are performed to explore similarities and differences in the two high-order schemes.
In this work, we design an entropy stable, finite volume approximation for the ideal magnetohydrodynamics (MHD) equations. The method is novel as we design an affordable analytical expression of the numerical interface flux function that discretely preserves the entropy of the system. To guarantee the discrete conservation of entropy requires the addition of a particular source term to the ideal MHD system. Exact entropy conserving schemes cannot dissipate energy at shocks, thus to compute accurate solutions to problems that may develop shocks, we determine a dissipation term to guarantee entropy stability for the numerical scheme. Numerical tests are performed to demonstrate the theoretical findings of entropy conservation and robustness.
In this work, we design an entropy stable, finite volume approximation for the shallow water magnetohydrodynamics (SWMHD) equations. The method is novel as we design an affordable analytical expression of the numerical interface flux function that exactly preserves the entropy, which is also the total energy for the SWMHD equations. To guarantee the discrete conservation of entropy requires a special treatment of a consistent source term for the SWMHD equations. With the goal of solving problems that may develop shocks, we determine a dissipation term to guarantee entropy stability for the numerical scheme. Numerical tests are performed to demonstrate the theoretical findings of entropy conservation and robustness.
We derive a spectrally accurate moving mesh method for mixed material interface problems modeled by Maxwell's or the classical wave equation. We use a discontinuous Galerkin spectral element approximation with an arbitrary Lagrangian-Eulerian mapping and derive the exact upwind numerical fluxes to model the physics of wave reflection and transmission at jumps in material properties. Spectral accuracy is obtained by placing moving material interfaces at element boundaries and solving the appropriate Riemann problem. We present numerical examples showing the performance of the method for plane wave reflection and transmission at dielectric and acoustic interfaces.
We derive and evaluate an explicit local time stepping (LTS) integration for the discontinuous Galerkin spectral element method (DGSEM) on moving meshes. The LTS procedure is derived from Adams-Bashforth multirate time integration methods. We also present speedup and memory estimates, which show that the explicit LTS integration scales well with problem size. Time-step refinement studies with static and moving meshes show that the approximations are spectrally accurate in space and have design temporal accuracy. The numerical tests validate theoretical estimates that the LTS procedure can reduce computational cost by as much as an order of magnitude for time accurate problems.
This work focuses on the accuracy and stability of high-order nodal discontinuous Galerkin (DG) methods for under-resolved turbulence computations. In particular we consider the inviscid Taylor-Green vortex (TGV) flow to analyse the implicit large eddy simulation (iLES) capabilities of DG methods at very high Reynolds numbers. The governing equations are discretised in two ways in order to suppress aliasing errors introduced into the discrete variational forms due to the under-integration of non-linear terms. The first, more straightforward way relies on consistent/over-integration, where quadrature accuracy is improved by using a larger number of integration points, consistent with the degree of the non-linearities. The second strategy, originally applied in the high-order finite difference community, relies on a split (or skew-symmetric) form of the governing equations. Different split forms are available depending on how the variables in the non-linear terms are grouped. The desired split form is then built by averaging conservative and non-conservative forms of the governing equations, although conservativity of the DG scheme is fully preserved. A preliminary analysis based on Burgers’ turbulence in one spatial dimension is conducted and shows the potential of split forms in keeping the energy of higher-order polynomial modes close to the expected levels. This indicates that the favourable dealiasing properties observed from split-form approaches in more classical schemes seem to hold for DG. The remainder of the study considers a comprehensive set of (under-resolved) computations of the inviscid TGV flow and compares the accuracy and robustness of consistent/over-integration and split form discretisations based on the local Lax-Friedrichs and Roe-type Riemann solvers. Recent works showed that relevant split forms can stabilize higher-order inviscid TGV test cases otherwise unstable even with consistent integration. Here we show that stable high-order cases achievable with both strategies have comparable accuracy, further supporting the good dealiasing properties of split form DG. The higher-order cases achieved only with split form schemes also displayed all the main features expected from consistent/over-integration. Among test cases with the same number of degrees of freedom, best solution quality is obtained with Roe-type fluxes at moderately high orders (around sixth order). Solutions obtained with very high polynomial orders displayed spurious features attributed to a sharper dissipation in wavenumber space. Accuracy differences between the two dealiasing strategies considered were, however, observed for the low-order cases, which also yielded reduced solution quality compared to high-order results.