Search

We study a temporal step size control of explicit Runge-Kutta (RK) methods for compressible computational fluid dynamics (CFD), including the Navier-Stokes equations and hyperbolic systems of conservation laws such as the Euler equations. We demonstrate that error-based approaches are convenient in a wide range of applications and compare them to more classical step size control based on a Courant-Friedrichs-Lewy (CFL) number. Our numerical examples show that the error-based step size control is easy to use, robust, and efficient, e.g., for (initial) transient periods, complex geometries, nonlinear shock capturing approaches, and schemes that use nonlinear entropy projections. We demonstrate these properties for problems ranging from well-understood academic test cases to industrially relevant large-scale computations with two disjoint code bases, the open source Julia packages Trixi.jl with OrdinaryDiffEq.jl and the C/Fortran code SSDC based on PETSc.

In this paper we will review a recent emerging paradigm shift in the construction and analysis of high order Discontinuous Galerkin methods applied to approximate solutions of hyperbolic or mixed hyperbolic-parabolic partial differential equations (PDEs) in computational physics. There is a long history using DG methods to approximate the solution of partial differential equations in computational physics with successful applications in linear wave propagation, like those governed by Maxwell’s equations, incompressible and compressible fluid and plasma dynamics governed by the Navier-Stokes and the Magnetohydrodynamics equations, or as a solver for ordinary differential equations (ODEs), e.g., in structural mechanics. The DG method amalgamates ideas from several existing methods such as the Finite Element Galerkin method (FEM) and the Finite Volume method (FVM) and is specifically applied to problems with advection dominated properties, such as fast moving fluids or wave propagation. In the numerics community, DG methods are infamous for being computationally complex and, due to their high order nature, as having issues with robustness, i.e., these methods are sometimes prone to crashing easily. In this article we will focus on efficient nodal versions of the DG scheme and present recent ideas to restore its robustness, its connections to and influence by other sectors of the numerical community, such as the finite difference community, and further discuss this young, but rapidly developing research topic by highlighting the main contributions and a closing discussion about possible next lines of research.

Discontinuous Galerkin (DG) methods have a long history in computational physics and engineering to approximate solutions of partial differential equations due to their high-order accuracy and geometric flexibility. However, DG is not perfect and there remain some issues. Concerning robustness, DG has undergone an extensive transformation over the past seven years into its modern form that provides statements on solution boundedness for linear and nonlinear problems. This chapter takes a constructive approach to introduce a modern incarnation of the DG spectral element method for the compressible Navier-Stokes equations in a three-dimensional curvilinear context. The groundwork of the numerical scheme comes from classic principles of spectral methods including polynomial approximations and Gauss-type quadratures. We identify aliasing as one underlying cause of the robustness issues for classical DG spectral methods. Removing said aliasing errors requires a particular differentiation matrix and careful discretization of the advective flux terms in the governing equations.

The second paper of this series presents two robust entropy stable shock-capturing methods for discontinuous Galerkin spectral element (DGSEM) discretizations of the compressible magneto-hydrodynamics (MHD) equations. Specifically, we use the resistive GLM-MHD equations, which include a divergence cleaning mechanism that is based on a generalized Lagrange multiplier (GLM). For the continuous entropy analysis to hold, and due to the divergence-free constraint on the magnetic field, the GLM-MHD system requires the use of non-conservative terms, which need special treatment. Hennemann et al. (2020) [25] recently presented an entropy stable shock-capturing strategy for DGSEM discretizations of the Euler equations that blends the DGSEM scheme with a subcell first-order finite volume (FV) method. Our first contribution is the extension of the method of Hennemann et al. to systems with non-conservative terms, such as the GLM-MHD equations. In our approach, the advective and non-conservative terms of the equations are discretized with a hybrid FV/DGSEM scheme, whereas the visco-resistive terms are discretized only with the high-order DGSEM method. We prove that the extended method is semi-discretely entropy stable on three-dimensional unstructured curvilinear meshes. Our second contribution is the derivation and analysis of a second entropy stable shock-capturing method that provides enhanced resolution by using a subcell reconstruction procedure that is carefully built to ensure entropy stability. We provide a numerical verification of the properties of the hybrid FV/DGSEM schemes on curvilinear meshes and show their robustness and accuracy with common benchmark cases, such as the Orszag-Tang vortex and the GEM (Geospace Environmental Modeling) reconnection challenge. Finally, we simulate a space physics application: the interaction of Jupiter's magnetic field with the plasma torus generated by the moon Io.

One of the challenges when simulating astrophysical flows with self-gravity is to compute the gravitational forces. In contrast to the hyperbolic hydrodynamic equations, the gravity field is described by an elliptic Poisson equation. We present a purely hyperbolic approach by reformulating the elliptic problem into a hyperbolic diffusion problem, which is solved in pseudotime using the same explicit high-order discontinuous Galerkin method we use for the flow solution. The flow and the gravity solvers operate on a joint hierarchical Cartesian mesh and are two-way coupled via the source terms. A key benefit of our approach is that it allows the reuse of existing explicit hyperbolic solvers without modifications, while retaining their advanced features such as non-conforming and solution-adaptive grids. By updating the gravitational field in each Runge-Kutta stage of the hydrodynamics solver, high-order convergence is achieved even in coupled multi-physics simulations. After verifying the expected order of convergence for single-physics and multi-physics setups, we validate our approach by a simulation of the Jeans gravitational instability. Furthermore, we demonstrate the full capabilities of our numerical framework by computing a self-gravitating Sedov blast with shock capturing in the flow solver and adaptive mesh refinement for the entire coupled system., Comment: 30 pages, 7 figures, second revision

In this work we analyze the entropic properties of the Euler equations when the system is closed with the assumption of a polytropic gas. In this case, the pressure solely depends upon the density of the fluid and the energy equation is not necessary anymore as the mass conservation and momentum conservation then form a closed system. Further, the total energy acts as a convex mathematical entropy function for the polytropic Euler equations. The polytropic equation of state gives the pressure as a scaled power law of the density in terms of the adiabatic index $$\gamma $$ γ . As such, there are important limiting cases contained within the polytropic model like the isothermal Euler equations ($$\gamma {=}1$$ γ = 1 ) and the shallow water equations ($$\gamma {=}2$$ γ = 2 ). We first mimic the continuous entropy analysis on the discrete level in a finite volume context to get special numerical flux functions. Next, these numerical fluxes are incorporated into a particular discontinuous Galerkin (DG) spectral element framework where derivatives are approximated with summation-by-parts operators. This guarantees a high-order accurate DG numerical approximation to the polytropic Euler equations that is also consistent to its auxiliary total energy behavior. Numerical examples are provided to verify the theoretical derivations, i.e., the entropic properties of the high order DG scheme.

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.

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., arXiv admin note: substantial text overlap with arXiv:1509.06902; text overlap with arXiv:1007.2606 by other authors

In this paper, we extend the entropy conservative/stable algorithms presented by Del Rey Fernandez et al. (2019) for the compressible Euler and Navier–Stokes equations on nonconforming p-refined/coarsened curvilinear grids to h/p refinement/coarsening. The main difficulty in developing nonconforming algorithms is the construction of appropriate coupling procedures across nonconforming interfaces. Here, we utilize a computationally simple and efficient approach based upon using decoupled interpolation operators. The resulting scheme is entropy conservative/stable and element-wise conservative. Numerical simulations of the isentropic vortex and viscous shock propagation confirm the entropy conservation/stability and accuracy properties of the method (achieving $$\sim p+1$$ convergence), which are comparable to those of the original conforming scheme (Carpenter et al. in SIAM J Sci Comput 36(5):B835–B867, 2014; Parsani et al. in SIAM J Sci Comput 38(5):A3129–A3162, 2016). Simulations of the Taylor–Green vortex at $$\hbox {Re}=1600$$ and turbulent flow past a sphere at $$\hbox {Re}_{\infty }=2000$$ show the robustness and stability properties of the overall spatial discretization for unstructured grids. Finally, to demonstrate the entropy conservation property of a fully-discrete explicit entropy stable algorithm with h/p refinement/coarsening, we present the time evolution of the entropy function obtained by simulating the propagation of the isentropic vortex using a relaxation Runge–Kutta scheme.

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.

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.

An open-source code that implements the entropy stable discontinuous Galerkin scheme with Legendere–Gauss–Lobatto collocation (DGSEM) on curved unstructured hexahedral grids for compressible Navier–Stokes equations (NSE) is available at github.com/project-fluxo.

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)., Comment: 54 pages

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.

When we published this article, the biographies of Gregor J Gassner and Stefanie Walch have been mixed up during the type-setting process. Sadly, this mistake remained unnoticed. The correct biographies are: Gregor Gassner is a professor in numerical analysis/scientific computing at the Mathematical Institute at the University of Cologne. Gregor obtained a diploma degree in aerospace engineering (2004) and a diploma degree in mathematics (2005) from the University of Stuttgart. In 2009 he received his doctoral degree at the Institute of Aerodynamics and Gasdynamics, Stuttgart. He got an offer for a Junior professorship at TU Braunschweig in computational fluid dynamics, an offer for scientific computing at the TU Hamburg-Harburg and an offer from University of Cologne, where he joined in 2013. In 2016 he received an ERC starting grant on the topic of adaptive, robust and highly efficient numerical methods for the simulation of non-linear conservation laws. Stefanie Walch is a full professor for theoretical astrophysics at the University of Cologne since 2013. After the diploma in physics in 2004, Stefanie received her doctoral degree in theoretical physics from the Ludwig-Maximilians-Universitat Munchen in 2008. Afterwards she has spent 3 years at Cardiff University in Wales, and approximately 2 years at the Max-Planck-Institute for Astrophysics in Garching as a postdoctoral researcher before being appointed at the University of Cologne. She specialises in the physics of the interstellar medium including star formation and stellar feedback in the form of ionizing radiation, stellar winds, and supernovae. The evolution of gas in galaxies is studied using massively parallel, three-dimensional magneto-hydrodynamics simulations with additional (self-)gravity, radiative transfer, and chemical networks to follow the formation and dissociation of molecules and to bridge the gap to observations in different tracers. Synthetic observations ranging from X-ray emission to CO and dust polarization maps are now used to put modern theories of star formation and feedback to the test. In 2015, she received an ERC Starting Grant on the topic of “The radiative interstellar medium”. The original article has been corrected. The publisher apologises for this mistake.

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., 6 pages

Fisher and Carpenter (\textit{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 subcell 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 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.

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., in Journal of Scientific Computing, 2015

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.

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.