Workshop: Novel Adaptive Discontinuous Galerkin Approaches for the Simulation of Atmospheric Flows

On November 30th and December 1st, the workshop on “Novel Adaptive Discontinuous Galerkin Approaches for the Simulation of Atmospheric Flows” will take place at the Department of Mathematics and Computer Science with researchers from the German Aerospace Center (DLR) Cologne, Deutscher Wetterdienst (DWD), and University of Cologne (UoC). This workshop is supported by and embedded into the Center for Earth System Observation and Computational Analysis (CESOC).

The idea of the workshop is to discuss novel developments for discontinuous Galerkin methods applied to atmospheric flow simulations, with a focus on implicit time integration, meshing of the sphere, adaptivity, efficient implementations, benchmark test cases for atmospheric flows, entropy-stable split formulations of the nonlinear partial differential equations, and bounds-preserving numerical schemes.

A detailed program with abstracts is available here.

Snapshot: Multi-component magneto-hydrodynamics simulations with a bounds-preserving discontinuous Galerkin method

We solve the multi-component magneto-hydrodynamics (MHD) equations using a hybrid finite volume (FV)/discontinuous Galerkin (DG) method [1], which combines the two discretization approaches at the node level. To test the robustness of the hybrid FV/DG method, we compute a modification of the Orszag-Tang vortex problem [2], where we use two ion species and initialize the flow with a sharp interface at $y=0$:
\rho^1(x,y,t=0) &= \begin{cases}
{25}/{36 \pi} & y \ge 0 \\
10^{-8} & y < 0
\end{cases} ,
&\rho^2(x,y,t=0) &= \begin{cases}
10^{-8} & y < 0\\
{25}/{36 \pi} & y \ge 0
\end{cases} , \\
v_1 (x,y,t=0) &= – \sin (2 \pi y),
&v_2 (x,y,t=0) &= \sin (2 \pi x), \\
B_1 (x,y,t=0) &= -\frac{1}{\sqrt{4 \pi}} \sin (2 \pi y),
&B_2 (x,y,t=0) &= -\frac{1}{\sqrt{4 \pi}} \sin (4 \pi x). \\
p (x,y,t=0) &= \frac{5}{12 \pi},
& & \\
where $\rho^1$ is the density of the first ion species, $\rho^2$ is the density of the second ion species,$\vec{v}=(v_1,v_2)$ is the plasma velocity, $\vec{B}=(B_1,B_2)$ is the magnetic field, and $p$ is the plasma pressure (without magnetic pressure).
The fourth-order DG method delivers high accuracy, but fails to describe shocks and material interfaces. Therefore, we use a modal indicator [3] to detect the elements where the solution changes abruptly. In those elements, we combine the DG method with a robust lower order FV scheme locally (at the node level) to impose TVD-like conditions on the density of the ion species,
\min_{k \in \mathcal{N} (j)} {\rho}^{i,\mathrm{FV}}_{k}
\le \rho^i_{j} \le
\max_{k \in \mathcal{N} (j)} {\rho}^{i,\mathrm{FV}}_{k},
and a local minimum principle on the specific entropy,
s({\mathbf{u}}_{j}) \ge \min_{k \in \mathcal{N} (j)} s(\mathbf{u}^{\mathrm{FV}}_{k}),.
where $\mathcal{N} (j)$ denotes the collection of nodes in the low-order stencil of each node $j$.
The video shows the total density ($\rho=\rho^1+\rho^2$), the density of the first ion species ($\rho^1$), and the so-called blending coefficient ($\alpha$) for the combination of a fourth-order accurate entropy stable DG method with first- and second-order accurate FV methods (Rusanov solver) and $1024^2$ degrees of freedom. A blending coefficient $\alpha=0$ means that only the DG is being used, and a blending coefficient $\alpha=1$ means that only the FV method is being used. The example shows that the hybrid FV/DG scheme is able to capture the shocks of the simulation correctly, deal with near-vacuum conditions and sharp interfaces.

[1] Rueda-Ramírez, A. M.; Pazner, W., & Gassner, G. J. (2022). Subcell limiting strategies for discontinuous Galerkin spectral element methods. Computers & Fluids. arXiv:2202.00576.
[2] Orszag, S.; Tang, C (1979). Small-scale structure of two-dimensional magnetohydrodynamic turbulence. Journal of Fluid Mechanics.
[3] Persson, P. O.; Peraire, J. (2006). Sub-cell shock capturing for discontinuous Galerkin methods. In 44th AIAA Aerospace Sciences Meeting and Exhibit (p. 112).

New paper submitted: Entropy-Stable Gauss Collocation Methods for Ideal Magneto-Hydrodynamics

In this paper, we present an entropy-stable Gauss collocation discontinuous Galerkin (DG) method on 3D curvilinear meshes for the GLM-MHD equations: the single-fluid magneto-hydrodynamics (MHD) equations with a generalized Lagrange multiplier (GLM) divergence cleaning mechanism. For the continuous entropy analysis to hold and to ensure Galilean invariance in the divergence cleaning technique, the GLM-MHD system requires the use of non-conservative terms.
Traditionally, entropy-stable DG discretizations have used a collocated nodal variant of the DG method, also known as the discontinuous Galerkin spectral element method (DGSEM) on Legendre-Gauss-Lobatto (LGL) points. Recently, Chan et al. (“Efficient Entropy Stable Gauss Collocation Methods”. SIAM -2019) presented an entropy-stable DGSEM scheme that uses Legendre-Gauss points (instead of LGL points) for conservation laws. Our main contribution is to extend the discretization technique of Chan et al. to the non-conservative GLM-MHD system.
We provide a numerical verification of the entropy behavior and convergence properties of our novel scheme on 3D curvilinear meshes. Moreover, we test the robustness and accuracy of our scheme with a magneto-hydrodynamic Kelvin-Helmholtz instability problem. The numerical experiments suggest that the entropy-stable DGSEM on Gauss points for the GLM-MHD system is more accurate than the LGL counterpart.


Snapshot: Supersonic turbulence simulations with FLUXO

We show a supersonic turbulence simulation in a periodic box fueled by perpetual kinetic energy injection on the largest three modes. Such setups are subject of active research in astrophysics and are believed to play an important role in the dynamics of galactic clouds and star formation, e.g., [1,2,3]. An initially uniform density distribution (left plot in the video) is gradually driven to a turbulent state till an average velocity of Mach 3.5 is reached (about t=0.19). Locally, speeds can spike above Mach 35 (bottom right plot in the video). The simulation was computed on FLUXO ( with a robust invariant domain preserving Discontinuous Galerkin scheme (DG) [4] solving the compressible, quasi-isothermal ($\gamma =1.001$) Euler equations on a uniform grid with $128^3$ degrees-of-freedom.

Smooth blending with a monotone finite volume scheme stabilizes the high-order DG at strongly shocked regions (filaments of high density concentration) displayed by the blending parameter plot (top right in the video). Blue encodes regions of 100% high order DG while red dots indicate focused blending with FV. Hypersonic turbulence regimes are very challenging for any numerical scheme considering the stark gradients in density, bubbles of near-vacuum and the extremely high speeds involved in such simulations.

[4] Rueda-Ramírez, A. M., Pazner, W., & Gassner, G. J. (2022). Subcell limiting strategies for discontinuous Galerkin spectral element methods. arXiv preprint arXiv:2202.00576.

Snapshot: Hypersonic astrophysical jet (Mach~2000) simulation using a hybrid DG/FV solver

Simulation of a hypersonic astrophysical jet with Mach=2156.91 moving through a medium at rest.

The simulation was performed on a 256 × 256 quadrilateral grid with a nodal fourth-order entropy-stable Discontinuous Galerkin (DG) scheme. Local bounds on the density and specific entropy are imposed by blending the DG scheme locally (at the node level) with a second-order Finite Volume (FV) method. Both the DG and the FV methods use the Harten-Lax-van Leer-Einfeldt (HLLE) Riemann solver.

The simulation results were obtained with the open-source code FLUXO (

The numerical methods used to perform the simulation are described in our preprint:
Rueda-Ramírez, A. M., Pazner, W., Gassner, G.J. (2022). Subcell Limiting Strategies for Discontinuous Galerkin Spectral Element Methods. Submitted.

New paper submitted: Subcell limiting strategies for discontinuous Galerkin spectral element methods

We present a general family of subcell limiting strategies to construct robust high-order accurate nodal discontinuous Galerkin (DG) schemes. The main strategy is to construct compatible low order finite volume (FV) type discretizations that allow for convex blending with the high-order variant with the goal of guaranteeing additional properties, such as bounds on physical quantities and/or guaranteed entropy dissipation. For an implementation of this main strategy, four main ingredients are identified that may be combined in a flexible manner: (i) a nodal high-order DG method on Legendre-Gauss-Lobatto nodes, (ii) a compatible robust subcell FV scheme, (iii) a convex combination strategy for the two schemes, which can be element-wise or subcell-wise, and (iv) a strategy to compute the convex blending factors, which can be either based on heuristic troubled-cell indicators, or using ideas from flux-corrected transport methods.
By carefully designing the metric terms of the subcell FV method, the resulting methods can be used on unstructured curvilinear meshes, are locally conservative, can handle strong shocks efficiently while directly guaranteeing physical bounds on quantities such as density, pressure or entropy. We further show that it is possible to choose the four ingredients to recover existing methods such as a provably entropy dissipative subcell shock-capturing approach or a sparse invariant domain preserving approach.
We test the versatility of the presented strategies and mix and match the four ingredients to solve challenging simulation setups, such as the KPP problem (a hyperbolic conservation law with non-convex flux function), turbulent and hypersonic Euler simulations, and MHD problems featuring shocks and turbulence.


Snapshot: Immersed Boundary Discontinuous Galerkin Simulation of Flow Past a Cylinder at Re=100

Simulation of the viscous flow past a cylinder with a discontinuous Galerkin spectral element method on a hierarchical Cartesian mesh with adaptive mesh refinement (AMR). The polynomial degree is N = 3 and the mesh uses five refinement levels. The cylinder boundary is imposed with an immersed boundary method (IBM) based on volume penalization [1], which adds a source term to the compressible Navier-Stokes equations,
\partial_t \mathbf{u} + \nabla \cdot \vec{\mathbf{f}} (\mathbf{u},\nabla \mathbf{u}) = \underbrace{\beta (\mathbf{u}_s – \mathbf{u})}_{\text{IBM source}},
where $\mathbf{u}$ is the vector of conservative quantities, $\vec{\mathbf{f}}$ is the Navier-Stokes flux, $\beta$ is the penalty parameter, and $\mathbf{u}_s$ is an imposed state. For this simulation, $\mathbf{u}_s$ mimics an isothermal no-slip wall boundary condition.

The simulation was obtained with the open-source code FLUXO (

[1] Kou, J., Joshi, S., Hurtado-de-Mendoza, A., Puri, K., Hirsch, C., & Ferrer, E. (2022). Immersed boundary method for high-order flux reconstruction based on volume penalization. Journal of Computational Physics, 448, 110721.

Snapshot: Numerical simulations of earth’s climate

Earth’s climate is changing due to human activity. Therefore, the ability to predict future climate change is of utmost importance if we want to reduce the human interference with the climate system (climate change mitigation) and our vulnerability to the harmful effects of climate change (climate change adaptation).

It is possible to model the time-dependent local temperatures on the globe using a simple two-dimensional energy balance model [1,2],
$$ C(\vec{x}) \frac{\partial T}{\partial t} + A + BT = \nabla \cdot \left( D(\vec{x}) \nabla T \right) + Q S(\vec{x},t) (1-a(\vec{x})), $$
where $T$ is the local temperature; $C(\vec{x})$ is the local heat capacity for a column of air on land, water or ice; the term $A+BT$ models the outgoing long-wave radiation, where $A$ depends on the amount of green-house gases in the atmosphere and $B$ models several feedback effects; $D$ is a latitude-dependent diffusion coefficient; and the last term includes the heat input in the system, where $Q$ is the solar radiation at earth’s position in the solar system, $S$ is the seasonal and latitudinal distribution of insolation, which depends on orbital parameters, such as the eccentricity of earth’s orbit, and the precession and tilt of earth’s rotation axis, and $a(\vec{x})$ is the local albedo, which measures the reflection of solar radiation and depends on the surface cover.

We developed a julia code that implements the EBM model on the sphere surface and used it to describe the mean temperatures on the globe assuming a constant atmospheric CO$_2$ concentration (of year 1950). The results are presented in the animation.

[1] North, G.R.; Mengel, J.G.; Short, D.A. (1983). Simple Energy Balance Model Resolving the Seasons and the Continents’ Application to the Astronomical Theory of the Ice Ages. Journal of Geophysical Research.
[2] Zhuang, K.; North, G.R.; Stevens, Mark J. (2017). A NetCDF version of the two-dimensional energy balance model based on the full multigrid algorithm. SoftwareX.

Snapshot: Interaction of Jupiter’s moon Io with the Jovian magnetosphere

Io is the most volcanically active body of the solar system. Furthermore, it is embedded in Jupiter’s magnetic field, the largest and most powerful planetary magnetosphere of the solar system. Due to its strong volcanic activity, Io expels ions and neutrals, which are in turn ionized by ultraviolet and electron impact ionization, forming a plasma torus around Jupiter [1,2]. As Io moves inside the plasma torus, elastic collisions of ions and neutrals inside its atmosphere generate a magnetospheric disturbance that propagates away from Io along the background magnetic field lines at the Alfvén wave speed. This phenomenon creates a pair of Alfvén current tubes that are commonly called Alfvén wings, which have been observed by several flybys [1].

The figure shows the momentum (non-dimensional) and magnetic field of the plasma that surrounds Io, obtained with the magneto-hydrodynamic (MHD) plasma model of FLUXO*. The Alfvén wings can be observed as disturbances of both the background magnetic and momentum fields. In yellow is the trajectory of the I31 flyby of the Galileo spacecraft, which visited Io in 2001 [1]. The Galileo spacecraft is also depicted (not to scale). A 3D model of the moon’s surface developed by NASA [3] was superposed on the simulation results.

*FLUXO ( is an MPI parallel high-order Discontinuous Galerkin code, which supports unstructured curvilinear hexahedral grids, and is able to perform Adaptive Mesh Refinement (AMR).

[1] M. Kivelson, K. Khurana, C. Russell, R. Walker, S. Joy, J. Mafi, GALILEO ORBITER AT JUPITER CALIBRATED MAG HIGH RES V1.0, GO-J-MAG-3-RDR-HIGHRES-V1.0, Technical Report, NASA Planetary Data System, 1997.
[2] J. Saur, F. M. Neubauer, J. E. P. Connerney, Plasma interaction of Io with its plasma torus, 2004.

Snapshot: Simulation of a Kelvin-Helmholtz instability using second order Finite Volume schemes and fourth order Discontinuous Galerkin methods

We present in-viscid and viscous simulations of a Kelvin-Helmholtz instability using second a order accurate monotoniced-central finite volume (FV) method and a fourth order accurate discontinuous Galerkin (DG) method. The initial condition is given by [1]:

$$\rho (t=0) = \frac{1}{2}
+ \frac{3}{4} B,
p (t=0) = 1,~~~~~~~~~~

v_1 (t=0) = \frac{1}{2} \left( B-1 \right),
v_2 (t=0) = \frac{1}{10} \sin(2 \pi x),

with $$B=\tanh \left( 15 y + 7.5 \right) – \tanh(15y-7.5).$$

We first present the FV results at end time $t=3.7$, which are computed using a monotoniced-central second order discretization of the Euler equations of gas dynamics on uniform grids.

The next results use a fourth order DG discretization of the Navier-Stokes equations on uniform grids using $Re=320.000$ at end time $t=3.7$. The highest resolution (4096² DOFs) is a direct numerical simulation (DNS) of the problem, where all scales are resolved.

It is remarkable that the numerical dissipation of the second order FV scheme causes the in-viscid simulation with 2048² DOFs to look very similar to the viscous DNS solution at $Re=320.000$.

[1] A.M. Rueda-Ramírez, G.J Gassner (2021). A Subcell Finite Volume Positivity-Preserving Limiter for DGSEM Discretizations of the Euler Equations.