We show that even though the Discontinuous Galerkin Spectral Element Method is stable for hyperbolic boundary-value problems, and the overset domain problem is well-posed in an appropriate norm, the energy of the approximation of the latter is bounded by data only for fixed polynomial order, mesh, and time. In the absence of dissipation, coupling of the overlapping domains is destabilizing by allowing positive eigenvalues in the system to be integrated in time. This coupling can be stabilized in one space dimension by using the upwind numerical flux. To help provide additional dissipation, we introduce a novel penalty method that applies dissipation at arbitrary points within the overlap region and depends only on the difference between the solutions. We present numerical experiments in one space dimension to illustrate the implementation of the well-posed penalty formulation, and show spectral convergence of the approximations when sufficient dissipation is applied.

In this work, we develop a new hydrostatic reconstruction procedure to construct well-balanced schemes for one and multilayer shallow water flows, including wetting and drying. Initially, we derive the method for a path-conservative finite volume scheme and combine it with entropy conservative fluxes and suitable numerical dissipation to preserve an entropy inequality in the semi-discrete case. We then combine the novel hydrostatic reconstruction with a collocated nodal split-form discontinuous Galerkin spectral element method, extending the method to high-order and curvilinear meshes. The high-order method incorporates an additional positivity-limiter and is blended with a compatible subcell finite volume method to maintain well-balancedness at wet/dry fronts. We prove entropy stability, well-balancedness, and positivity-preservation for both methods. Numerical results for the high-order method validate the theoretical findings and demonstrate the robustness of the scheme.

High-order numerical methods for conservation laws are highly sought after due to their potential efficiency. However, it is challenging to ensure their robustness, particularly for under-resolved flows. Baseline high-order methods often incorporate stabilization techniques that must be applied judiciously—sufficient to ensure simulation stability but restrained enough to prevent excessive dissipation and loss of resolution. Recent studies have demonstrated that combining upwind summation-by-parts (USBP) operators with flux vector splitting can increase the robustness of finite difference (FD) schemes without introducing excessive artificial dissipation. This work investigates whether the same approach can be applied to nodal discontinuous Galerkin (DG) methods. To this end, we demonstrate the existence of USBP operators on arbitrary grid points and provide a straightforward procedure for their construction. Our discussion encompasses a broad class of USBP operators, not limited to equidistant grid points, and enables the development of novel USBP operators on Legendre–Gauss–Lobatto (LGL) points that are well-suited for nodal DG methods. We then examine the robustness properties of the resulting DG-USBP methods for challenging examples of the compressible Euler equations, such as the Kelvin–Helmholtz instability. Similar to high-order FD-USBP schemes, we find that combining flux vector splitting techniques with DG-USBP operators does not lead to excessive artificial dissipation. Furthermore, we find that combining lower-order DG-USBP operators on three LGL points with flux vector splitting indeed increases the robustness of nodal DG methods. However, we also observe that higher-order USBP operators offer less improvement in robustness for DG methods compared to FD schemes. We provide evidence that this can be attributed to USBP methods adding dissipation only to unresolved modes, as FD schemes typically have more unresolved modes than nodal DG methods.

We use the framework of upwind summation-by-parts (SBP) operators developed by Mattsson (2017) and study different flux vector splittings in this context. To do so, we introduce discontinuous-Galerkin-like interface terms for multi-block upwind SBP methods applied to nonlinear conservation laws. We investigate the behavior of the upwind SBP methods for flux vector splittings of varying complexity on Cartesian as well as unstructured curvilinear multi-block meshes. Moreover, we analyze the local linear/energy stability of these methods following Gassner, Svärd, and Hindenlang (2022). Finally, we investigate the robustness of upwind SBP methods for challenging examples of shock-free flows of the compressible Euler equations such as a Kelvin-Helmholtz instability and the inviscid Taylor-Green vortex.

We present an entropy stable nodal discontinuous Galerkin spectral element method (DGSEM) for the two-layer shallow water equations on two dimensional curvilinear meshes. We mimic the continuous entropy analysis on the semi-discrete level with the DGSEM constructed on Legendre–Gauss–Lobatto (LGL) nodes. The use of LGL nodes endows the collocated nodal DGSEM with the summation-by-parts property that is key in the discrete analysis. The approximation exploits an equivalent flux differencing formulation for the volume contributions, which generate an entropy conservative split-form of the governing equations. A specific combination of a numerical surface flux and discretization of the nonconservative terms is then applied to obtain a high-order path-conservative scheme that is entropy conservative. Furthermore, we find that this combination yields an analogous discretization for the pressure and nonconservative terms such that the numerical method is well-balanced for discontinuous bathymetry on curvilinear domains. Dissipation is added at the interfaces to create an entropy stable approximation that satisfies the second law of thermodynamics in the discrete case, while maintaining the well-balanced property. We conclude with verification of the theoretical findings through numerical tests and demonstrate results about convergence, entropy stability and well-balancedness of the scheme.

We extend the construction of so-called encapsulated global summation-by-parts operators to the general case of a mesh which is not boundary conforming. Owing to this development, energy stable discretizations of nonlinear and variable coefficient initial boundary value problems can be formulated in simple and straightforward ways using high-order accurate operators of generalized summation-by-parts type. Encapsulated features on a single computational block or element may include polynomial bases, tensor products as well as curvilinear coordinate transformations. Moreover, through the use of inner product preserving interpolation or projection, the global summation-by-parts property is extended to arbitrary multi-block or multi-element meshes with non-conforming nodal interfaces.

HOHQMesh generates unstructured all-quadrilateral and hexahedral meshes with high order boundaryinformation for use with spectral element solvers. Model input by the user requires only anoptional outer boundary curve plus any number of inner boundary curves that are built aschains of simple geometric entities (lines and circles), user defined equations, and cubic splines.Inner boundary curves can be designated as interface boundaries to force element edges alongthem. Quadrilateral meshes are generated automatically with the mesh sizes guided by abackground grid and the model, without additional input by the user. Hexahedral meshesare generated by extrusions of a quadrilateral mesh, including sweeping along a curve, andcan follow bottom topography. The mesh files that HOHQMesh generates include high orderpolynomial interpolation points of arbitrary order.

Many modern discontinuous Galerkin (DG) methods for conservation laws make use of summation by parts operators and flux differencing to achieve kinetic energy preservation or entropy stability. While these techniques increase the robustness of DG methods significantly, they are also computationally more demanding than standard weak form nodal DG methods. We present several implementation techniques to improve the efficiency of flux differencing DG methods that use tensor product quadrilateral or hexahedral elements, in 2D or 3D respectively. Focus is mostly given to CPUs and DG methods for the compressible Euler equations, although these techniques are generally also useful for other physical systems including the compressible Navier-Stokes and magnetohydrodynamics equations. We present results using two open source codes, Trixi.jl written in Julia and FLUXO written in Fortran, to demonstrate that our proposed implementation techniques are applicable to different code bases and programming languages.

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.

We derive boundary conditions and estimates based on the energy and entropy analysis of systems of the nonlinear shallow water equations in two spatial dimensions. It is shown that the energy method provides more details, but is fully consistent with the entropy analysis. The details brought forward by the nonlinear energy analysis allow us to pinpoint where the difference between the linear and nonlinear analysis originate. We find that the result from the linear analysis does not necessarily hold in the nonlinear case. The nonlinear analysis leads in general to a different minimal number of boundary conditions compared with the linear analysis. In particular, and contrary to the linear case, the magnitude of the flow does not influence the number of required boundary conditions.