Utilizing a mapping $(\xi, \eta) \in \mathbb{R}^2 \rightarrow (x, y, z) \in \mathbb{R}^3$ and specifically tailored tensor-product Legendre–Gauss–Lobatto basis functions [1,2], discontinuous Galerkin (DG) simulations can be performed on a curved surface. In particular, we are interested in running DG simulations on a spherical shell.

To achieve this, we create a two-dimensional cubed-sphere mesh for tessellating the sphere’s surface. This mesh proves advantageous as it avoids singularities at the poles present in latitude-longitude grids, while still facilitating a highly regular tessellation of the simulation domain through the use of quadrilaterals.

Within this mesh, we address the linear advection equations incorporating position-dependent advection velocity:
\[ \frac{\partial \rho}{\partial t} + \nabla \cdot \left( \vec{v} (x, y, z) \rho \right) = 0. \]
Additionally, a solid-body rotation velocity tangential to the spherical surface is imposed.

The presented video showcases the advection of a Gaussian pulse across the globe using Trixi.jl (https://github.com/trixi-framework/Trixi.jl) with the two-dimensional p4est solver. The simulation also incorporates a custom implementation of the three-dimensional linear advection equations.

References

[1] Song, C. & Wolf, J. P. (1999). The scaled boundary finite element method—alias consistent infinitesimal finite element cell method—for diffusion. International Journal for Numerical Methods in Engineering, 45(10), 1403-1431.
[2] Giraldo, F. X. (2001). A spectral element shallow water model on spherical geodesic grids. International Journal for Numerical Methods in Fluids, 35(8), 869-901.

When solving PDEs with the finite volume method, one must choose a numerical flux function. We use a convex combination of the central flux $F^C$ and local Lax Friedrich flux $F^{LLF}$ to solve the Burgers equation.

\begin{align*}
&u_t + \left(\frac{u^2}{2}\right)_x =0\\
&\frac{du_j}{dt} =\frac{1}{\Delta x} \left[ \alpha_j \left(F_{j-1/2}^{LLF} -F_{j+1/2}^{LLF}\right) + (1-\alpha_j)\left(F_{j-1/2}^{C} – F_{j+1/2}^{C}\right) \right]
\end{align*}

The central flux is highly accurate but can lead to oscillating and thus unstable solutions. The local Lax Friedrich flux is stable but also dissipative. Therefore, we want to choose alpha so that the solution is stable and at the same time as accurate as possible. We train a Reinforcement Learning Agent to choose alpha. One agent performs locally in one grid cell. The action of agent j is $\alpha_j$ and the state is given by

\begin{equation*}
r_j=\frac{u_j – u_{j-1}}{u_{j+1} – u_{j}}.
\end{equation*}

The policy is a Neural Network with one hidden layer of size 10, a relu activation function in the first layer, and a hard sigmoid function in the output layer. The deep deterministic policy gradient (DDPG) approach is used to train the agent.

The video below shows the chosen amount of local Lax Friedrich flux and the resulting solution of the Burgers equation. The initial solution is a sine wave with two discontinuities. Randomly changed variants of this were used for training.

The agent performs well in unseen situations, e.g. negative sine wave, as shown below.

This FV simulation uses linear reconstruction to reach second order. The minmod slope limiter is used. The mesh contains triangles and quads and is a high-resolution version (with 327680 cells) of this one:

It is created and organized using t8code.

The simulation uses an initial condition with a pressure blast wave in the center and periodic boundaries.

This simulation uses weakly compressible smoothed particle hydrodynamics for the fluid and a total Lagrangian SPH formulation for the elastic solid.

DGSEM simulations of the two-dimensional Euler equations considering the astrophysical jet with Mach number $\textrm{Ma} \approx 2000$ [1].
The computational domain, $\Omega = [-0.5,0.5]^2$, is filled with a mono-atomic gas ($\gamma = 5/3$) at rest with

\begin{equation}
\rho(x,y) = 0.5, \qquad
p(x,y) = 0.4127, \qquad
v_1(x,y) = 0, \qquad
v_2(x,y) = 0,
\end{equation}

and on the left boundary there is a hypersonic inflow with

\begin{equation}
\rho(-0.5,y_B) = 5, ~
p(-0.5,y_B) = 0.4127, ~
v_1(-0.5,y_B) = 800, ~
v_2(-0.5,y_B) = 0,
\end{equation}

for $y_B \in [-0.05, 0.05]$, which corresponds to a Mach number of $\textrm{Ma}=2156.91$ with respect to the speed of sound of the jet gas, and $\textrm{Ma}=682.08$ with respect to the speed of sound of the ambient gas.

A spatial resolution of $256\times 256$ quadrilateral elements, a polynomial degree of $N=3$, and four different CFL numbers are used.

The simulations are stabilized by using two different subcell limiting approaches, the monolithic convex limiting (MCL) and the FCT/IDP limiting.

The resulting density contours at $t=10^{-3}$:

The dependence of the spatial discretization on the time-step size for FCT/IDP methods causes the number of vortical structures to be highly dependent on the CFL number and the total number of time steps to be not inversely proportional to the CFL number.

In fact, the amount of dissipation is reduced for small CFL numbers, which leads to lower minimum densities and higher maximum pressures in FCT/IDP, as can be seen in this Table.

References
[1] Rueda-Ramírez, A., Bolm, B., Kuzmin, D. & Gassner, G. (2023). Monolithic Convex Limiting for Legendre-Gauss-Lobatto Discontinuous Galerkin Spectral Element Methods. arXiv preprint arXiv:2303.00374

The discontinuous Galerkin spectral element method (DGSEM) is a nodal discontinuous Galerkin collocation scheme that uses a set of Legendre-Gauss or Legendre-Gauss-Lobatto nodes to represent the solution with Lagrange interpolating polynomials and to compute integrals with a quadrature rule.

It was recently shown [1] that the DGSEM semi-discretization of a conservation law using Legendre-Gauss nodes can be written in finite-volume form as
\begin{equation*}
J_{j} \dot{\mathbf{u}}^{DG}_{j} =
\frac{1}{\omega_j}
\left(
\hat{\mathbf{f}}^{DG}_{(j-1,j)}
– \hat{\mathbf{f}}^{DG}_{(j,j+1)}
\right),
\end{equation*}
where $J_j$ is the mapping Jacobian, $\omega_j$ are the Legendre-Gauss quadrature weights, and $\hat{\mathbf{f}}^{DG}_{(\cdot,\cdot)}$ are high-order DGSEM fluxes between adjacent nodes.

The existence of a FV form for the Legendre–Gauss DGSEM enables the use of the subcell limiting strategies described in [2]. In these methods, a hybrid DGSEM/FV method is used, where the semi-discretization at each node reads
\begin{equation*}
J_{j} \dot{\mathbf{u}}_{j} =
\frac{1}{\omega_j}
\left(
\hat{\mathbf{f}}_{(j-1,j)}
– \hat{\mathbf{f}}_{(j,j+1)}
\right),
\end{equation*}
and the local fluxes are obtained as a convex combination of high-order DGSEM and low-order FV fluxes,
\begin{equation*}
\hat{\mathbf{f}}_{(i,j)} = (1-\alpha_{(i,j)}) \hat{\mathbf{f}}^{DG}_{(i,j)}
+
\alpha_{(i,j)}
\hat{\mathbf{f}}^{DG}_{(i,j)},
~~~~
\alpha_{(i,j)} \in [0,1].
\end{equation*}

To illustrate the shock-capturing capacity of the hybrid DGSEM/FV method, we simulate a Sedov blast problem describing the evolution of a blast wave expanding from an initial concentration of density and pressure, and adjust the blending coefficient $\alpha_{(i,j)}$ to avoid non-physical density oscillations in the vicinity of shocks.

The video shows the distribution of density and the blending coefficient obtained for entropy-stable variants of the Legendre-Gauss and Legendre-Gauss-Lobatto DGSEM methods.

The simulation runs stably without any spurious oscillations for both variants of the DGSEM method.
The figures below show the mean blending coefficient and the number of time steps for the simulation of the blast wave.
Although the Legendre–Gauss DGSEM requires more limiting than the Legendre-Gauss-Lobatto DGSEM to avoid spurious density oscillations, it completes the simulation in fewer time steps.
The Legendre-Gauss subcell distribution allows for longer time-step sizes than the Legendre-Gauss-Lobatto subcell distribution for the same CFL number.

Mean blending coefficient

Number of time steps

Figure 1: Evolution of the mean blending coefficient and number of time steps taken for the simulation of the blast wave with CFL=0.9.

References
[1] Mateo-Gabín, A., Rueda-Ramírez, A. M., Valero, E., & Rubio, G. (2022). Entropy-stable flux-differencing formulation with Gauss nodes for the DGSEM. arXiv preprint arXiv:2211.05066.
[2] Rueda-Ramírez, A. M., Pazner, W., & Gassner, G. J. (2022). Subcell limiting strategies for discontinuous Galerkin spectral element methods. Computers & Fluids, 247, 105627. arXiv preprint arXiv:2202.00576.

SPH simulation of two elastic spheres with different densities falling into a tank of water:

We solve the shallow water equations with DGSEM on a mesh with 103 elements and polynomial degree 12 in combination with subcell shock capturing and limiting, positivity preservation and handling of wet/dry regions.