Snapshot: Discontinuous Galerkin Solver for the Spherical Shallow Water Equations

A preliminary step in the development of a global atmospheric model is the construction of a suitable horizontal discretization for partial differential equations in spherical geometry. The spherical shallow water equations, which consist of a scalar equation governing mass conservation coupled with a vector equation for the momentum balance under gravitational and Coriolis forces, serve as a simplified model for the horizontal dynamics of the Earth’s atmopheric circulation. Since the spherical shallow water equations exhibit many of the characteristic features and numerical challenges associated with atmospheric fluid flow, they provide a useful testbed for the development and assessment of numerical schemes for weather prediction and climate modelling (see, for example, Williamson et al. [1]).

By extending Trixi.jl to solve hyperbolic partial differential equations on curved manifolds, we are able to simulate shallow-water flows on the surface of the sphere using discontinuous Galerkin methods formulated with respect to the two-dimensional tangent space associated with a cubed-sphere grid [2, 3]. The video below depicts the relative vorticity field for the numerical solution to the spherical shallow water equations in flux form, which we discretize similarly to Bao et al. [4] using a discontinuous Galerkin method employing 5400 curved quadrilateral elements of polynomial degree seven. The initial condition corresponds to a mid-latitude jet with a small perturbation added to initiate a barotropic instability, which was proposed by Galewsky et al. [5] as a test case exhibiting complex nonlinear dynamics characteristic of those present in numerical weather prediction and climate models.


References

[1] D. L. Williamson, J. B. Drake, J. J. Hack, R. Jakob, and P. N. Swarztrauber. A standard test set for numerical approximations to the shallow water equations in spherical geometry. Journal of Computational Physics, 102(1):211-224, 1992.

[2] R. Sadourny. Conservative finite-difference approximations of the primitive equations on quasi-uniform spherical grids. Monthly Weather Review, 100(2):136-144, 1972.

[3] C. Ronchi, R. Iacono, and P. S. Paolucci. The “cubed sphere”: A new method for the solution of partial differential equations in spherical geometry. Journal of Computational Physics, 124(1):93-114, 1996.

[4] L. Bao, R. D. Nair, and H. M. Tufo. A mass and momentum flux-form high-order discontinuous Galerkin shallow water model on the cubed-sphere. Journal of Computational Physics, 271:224-243, 2014.

[5] J. Galewsky, R. K. Scott, and L. M. Polvani. An initial-value problem for testing numerical models of the global shallow-water equations. Tellus A, 56(5):429–440, 2004.

This snapshot was created by Tristan Montoya.

Snapshot: Postdoctoral researcher Dr. Simone Chiocchetti joins our research group

Simone Chiocchetti studied Civil, Environmental, and Mechanical Engineering at the University of Trento, Italy, under the supervision of Prof. Dr.-Ing. Michael Dumbser.

He spent a year in Stuttgart at the IAG (Institute of Aerodynamics and Gasdynamics), funded by the DFG project “DROPIT startup grant” and he is currently working at the University of Cologne in the Numerical Simulation Research group led by Prof. Dr.-Ing. Gregor Gassner, under a Marie Skłodowska-Curie Postdoctoral Fellowship (European Union’s Horizon Europe Research and Innovation Programme, grant agreement No. 101109532).

Simone’s research interests include: high order numerical methods, of the Finite Volume and Discontinuous Galerkin families, in the context of first order hyperbolic partial differential equations; algorithms for generation and optimization of high quality unstructured meshes; treatment of stiff algebraic source terms; modeling of fluid and solid mechanics as well as multiphase flow; high-performance computing; semi-implicit numerical methods.

Simone’s current MSCA project (MoMeNTUM: Modern high order numerical Methods based on No-compromise moving Voronoi Tessellations, a Unified solver for continuum Mechanics) concerns the development of high order numerical methods on Voronoi grids, together with new meshing algorithms that allow to fully take advantage of the versatility provided by unstructured moving polygonal meshes, with special attention to the computational efficiency of the implementation.

Snapshot: Postdoctoral researcher Dr. Aleksey Sikstel joins our research group


Aleksey defended his doctoral thesis on coupling of hyperbolic conservation laws in 2020 at the RWTH Aachen University while working with Prof. Siegfried Müller and Prof. Michael Herty. After that he continued working with Prof. Jan Giesselmann at the TU Darmstadt on a posteriori error estimators for systems of hyperbolic conservation laws before switching to Cologne in 2022.

Alekseys research focus lies on adaptivity for high-order numerical methods for hyperbolic conservation laws, coupling of conservation laws via sharp interfaces, multiphase flows and stochastic-Galerkin PDE. In the Numerical Simulation Research Group Aleksey is working on multi-ion MHD systems and novel coupling and temporal schemes within the Snubic-project.

Snapshot: Postdoctoral researcher Dr. Benedict Geihe joins our research group

Benedict Geihe studied Mathematics and Computer Science at the University of Bonn, Germany. After graduation he joined the research group of Professor Martin Rumpf at the Institute for Numerical Simulation as a research and teaching assistant. In his PhD studies he investigated and implemented finite element based numerical simulation and shape optimization of multi-scale microstructures in elastic composite materials.

Benedict then moved to the German Aerospace Center (DLR) and worked as a researcher at the Numerical Methods department of the Institute of Propulsion Technology. He was part of the development team of the flow solver TRACE, used for academic and industrial simulation and analysis of flows in turbomachines. His research covered resonance phenomena induced by fluid structure interaction and frequency based numerical methods.

Benedict is now part of the Numerical Simulation Research Group at the University of Cologne, where he joined the Trixi.jl development team. He is currently working on project ADAPTEX, funded by the German Federal Ministry of Education and Research within the SCALEXA initiative. The objective is to deploy Trixi.jl in exascale-capable flow simulations on CPUs and GPUs for applications in earth system modeling on dynamic adaptive meshes.

Snapshot: 2D Sloshing simulation

2D Sloshing simulation in a round tank oscillating from left to right. The simulation was done with the CFD code Trixi.jl [1] which was extended with a three-equations model [2]. The model consists of a two-phase flow simplification of the Baer–Nunziato equations of compressible multi-phase flows [3]. The shown simulation is part of combined research efforts of the German Aerospace Center (DLR) and the University of Cologne to create computational models of sloshing dynamics in liquid hydrogen storage systems on mobile platforms such as cars, ships, and planes.


References:

[1] https://trixi-framework.github.io/
[2] Dumbser, Michael. “A simple two-phase method for the simulation of complex free surface flows.” Computer methods in applied mechanics and engineering 200.9-12 (2011): 1204-1219.
[3] M.R. Baer, J.W. Nunziato: “A two-phase mixture theory for the deflagration-todetonation transition (DDT) in reactive granular materials”, J. Multiphase Flow 12 (1986) 861–889.
[4] Project HyTaZer: https://elib.dlr.de/201347/ and https://www.dlr.de/sy/PortalData/17/Resources/dokumente/wissenschaftstag/2022/6_Wissenschaftstag_2022_Hytazer_Freund.pdf

This Snapshot was created by Johannes Markert (http://jmark.de/).

Snapshot: Discontinuous Galerkin simulation on a spherical shell

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.

Snapshot: Blending Finite Volume Fluxes with Reinforcement Learning

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.

Snapshot: Second-order FV method on hybrid mesh using t8code

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.