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$:
\[
\begin{array}{rlrl}
\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},
& & \\
\end{array}
\]
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).

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 (https://github.com/project-fluxo/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. https://arxiv.org/abs/2202.00576.

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 (https://github.com/project-fluxo/fluxo).

References
[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.

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.