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