In this paper, we present a discontinuous Galerkin solver based on previous work by Markert et al. (2021) for magneto-hydrodynamics in form of a new fluid solver module integrated into the established and well-known multi-physics simulation code FLASH. Our goal is to enable future research on the capabilities and potential advantages of discontinuous Galerkin methods for complex multi-physics simulations in astrophysical settings. We give specific details and adjustments of our implementation within the FLASH framework and present extensive validations and test cases, specifically its interaction with several other physics modules such as (self-)gravity and radiative transfer. We conclude that the new DG solver module in FLASH is ready for use in astrophysics simulations and thus ready for assessments and investigations.
Author Archives: Christof Czernik
New paper submitted: Efficient implementation of modern entropy stable and kinetic energy preserving discontinuous Galerkin methods for conservation laws
Many modern discontinuous Galerkin (DG) methods for conservation laws make use of summation by parts operators and flux differencing to achieve kinetic energy preservation or entropy stability. While these techniques increase the robustness of DG methods significantly, they are also computationally more demanding than standard weak form nodal DG methods. We present several implementation techniques to improve the efficiency of flux differencing DG methods that use tensor product quadrilateral or hexahedral elements, in 2D or 3D respectively. Focus is mostly given to CPUs and DG methods for the compressible Euler equations, although these techniques are generally also useful for GPU computing and other physical systems including the compressible Navier-Stokes and magnetohydrodynamics equations. We present results using two open source codes, Trixi.jl written in Julia and FLUXO written in Fortran, to demonstrate that our proposed implementation techniques are applicable to different code bases and programming languages.
New paper published: On the theoretical foundation of overset grid methods for hyperbolic problems: Well-posedness and conservation
We use the energy method to study the well-posedness of initial-boundary value problems approximated by overset mesh methods in one and two space dimensions for linear constant-coefficient hyperbolic systems. We show that in one space dimension, for both scalar equations and systems of equations, the problem where one domain partially oversets another is well-posed when characteristic coupling conditions are used. If a system cannot be diagonalized, as is usually the case in multiple space dimensions, then the energy method does not give proper bounds in terms of initial and boundary data. For those problems, we propose a novel penalty approach. We show, by using a global energy that accounts for the energy in the overlap region of the domains, that under well-defined conditions on the coupling matrices the penalized overset domain problems are energy bounded, conservative, well-posed and have solutions equivalent to the original single domain problem.
Published in Journal of Computational Physics (ScienceDirect)
Snapshot: Simulation of a stiff multi-species reacting flow
Simulation of an asymmetric two-dimensional stiff multi-species reacting flow for the compressible Euler equations. The computation domain [0,6]x[0,2] of this detonation problem is split into three zones. Zone A [0,0.5]x[0,2] as well as zone B (0.5,6]x[1.2,2] is filled with burnt gas, namely O2,OH and H2O. The last zone C (0.5,6]x[0,1.2) contains unburnt gas, namely H2 and O2.
The reacting model consists of two reactions H2+O2->2OH, 2OH+H2->2H2O and five species H2,O2,OH,H2O,N2 whereas reaction 1 has a smaller ignition temperature as well as a much faster chemical rate than reaction 2. Species N2 acts as a dilute catalyst here.
This example has been calculated using Trixi.jl with a P4est mesh using 3072x1024=3.145.728 elements, Dirichlet boundary conditions in x-direction as well as solid wall boundaries in y-direction. The final time is T = 0.1. The convection part is solved with FV whereas the chemical network is solved with KROME. (Note: Not all of these features are currently available in the main code of Trixi.jl.)
Pressure distribution at final time T=0.1:
Temperature distribution at final time T=0.1:
Density distribution at final time T=0.1:
OH-fraction distribution at final time T=0.1:
Snapshot: Simulation of a colliding flow setup with DGSEM
Simulation of a colliding flow setup with DGSEM (polynomial degree 3,
compatible sub-cell LGL type FV for shock capturing) of the compressible
Euler equations, where at x = +/- 64 inflow boundary conditions with
Mach number Ma = 70 is set (with periodic boundary conditions in y
direction).
Density distribution, where the color bar is plotted in logarithmic scale, at
final simulation time t=25:
Simulation with AMR (level 3 up to level 9), showing the AMR indicator at
final time t=25:
Plot of the DGSEM/FV blending factor for shock capturing at final time
t=25:
The combination of the AMR indicator and the shock capturing blending
factor triggers the final mesh refinement, as can be seen here at the final time t=25:
New paper published: An entropy stable nodal discontinuous Galerkin method for the resistive MHD equations. Part II: Subcell finite volume shock capturing
The second paper of this series presents two robust entropy stable shock-capturing methods for discontinuous Galerkin spectral element (DGSEM) discretizations of the compressible magneto-hydrodynamics (MHD) equations. Specifically, we use the resistive GLM-MHD equations, which include a divergence cleaning mechanism that is based on a generalized Lagrange multiplier (GLM). For the continuous entropy analysis to hold, and due to the divergence-free constraint on the magnetic field, the GLM-MHD system requires the use of non-conservative terms, which need special treatment.
Hennemann et al. (2020) [25] recently presented an entropy stable shock-capturing strategy for DGSEM discretizations of the Euler equations that blends the DGSEM scheme with a subcell first-order finite volume (FV) method. Our first contribution is the extension of the method of Hennemann et al. to systems with non-conservative terms, such as the GLM-MHD equations. In our approach, the advective and non-conservative terms of the equations are discretized with a hybrid FV/DGSEM scheme, whereas the visco-resistive terms are discretized only with the high-order DGSEM method. We prove that the extended method is semi-discretely entropy stable on three-dimensional unstructured curvilinear meshes. Our second contribution is the derivation and analysis of a second entropy stable shock-capturing method that provides enhanced resolution by using a subcell reconstruction procedure that is carefully built to ensure entropy stability.
We provide a numerical verification of the properties of the hybrid FV/DGSEM schemes on curvilinear meshes and show their robustness and accuracy with common benchmark cases, such as the Orszag-Tang vortex and the GEM (Geospace Environmental Modeling) reconnection challenge. Finally, we simulate a space physics application: the interaction of Jupiter’s magnetic field with the plasma torus generated by the moon Io.
Published in Journal of Computational Physics (ScienceDirect)
Snapshot: Evolution of a magnetized torus differentially rotating around a static point mass
Evolution of a magnetized torus differentially rotating around a static point mass. [1] The torus gradually forms into an accretion disk with expanding filaments into the surrounding space due to magneto-rotational instabilities similar to the rich structure of the corona of our sun. The simulation was done with a new Discontinuous-Galerkin based MHD solver [2] in FLASH [3]. Fig. 1 and Fig. 2 show the log. scale density cross section in the x-y and x-z plane. Fig 3 highlights the rich structure of the magnetic field (white stream lines) overlayed on top of the log. scale magnetic pressure in the x-z plane.
[1] Machida, Mami, Mitsuru R. Hayashi, and R. Matsumoto.
“Global simulations of differentially rotating magnetized disks:
Formation of low-beta filaments and structured coronae.”
The Astrophysical Journal Letters 532.1 (2000): L67.
[2] Markert et al. “A Discontinuous Galerkin Method with Sub‑Cell Adaptive Shock Capturing for
FLASH” (in preparation)
[3] http://flash.uchicago.edu/
Snapshot: Gingerbread man simulation with unstructured quadrilateral elements
Here we are simulating a gingerbread man, whereby the gingerbread man itself has 861 unstructured quadrilateral elements. For the run parameters, we used polydeg = 10 with the compressible Euler equations and the manufactured (convergence test) solution. The visualization is done on nvisnodes=12 uniform points in each direction, on each element.
We thank Prof. Dr. Andrew Winters (https://liu.se/en/employee/andwi94) for providing this video/simulation.
New paper submitted: On the Theoretical Foundation of Overset Grid Methods for Hyperbolic Problems: Well-Posedness and Conservation
We use the energy method to study the well-posedness of initial-boundary value problems approximated by overset mesh methods in one and two space dimensions for linear constant-coefficient hyperbolic systems. We show that in one space dimension, for both scalar equations and systems of equations, the problem where one domain partially oversets another is well-posed when characteristic coupling conditions are used. If a system cannot be diagonalized, as is ususally the case in multiple space dimensions, then the energy method does not give proper bounds in terms of initial and boundary data. For those problems, we propose a novel penalty approach. We show, by using a global energy that accounts for the energy in the overlap region of the domains, that under well-defined conditions on the coupling matrices the penalized overset domain problems are energy bounded, conservative, well-posed and have solutions equivalent to the original single domain problem.
Preprint available at: https://arxiv.org/abs/2105.04664
New paper published: Stability of Discontinuous Galerkin Spectral Element Schemes for Wave Propagation when the Coefficient Matrices have Jumps
Abstract:
We use the behavior of the L2 norm of the solutions of linear hyperbolic equations with discontinuous coefficient matrices as a surrogate to infer stability of discontinuous Galerkin spectral element methods (DGSEM). Although the L2 norm is not bounded by the initial data for homogeneous and dissipative boundary conditions for such systems, the L2 norm is easier to work with than a norm that discounts growth due to the discontinuities. We show that the DGSEM with an upwind numerical flux that satisfies the Rankine-Hugoniot (or conservation) condition has the same energy bound as the partial differential equation does in the L2 norm, plus an added dissipation that depends on how much the approximate solution fails to satisfy the Rankine-Hugoniot jump.
Accepted in Journal of Scientific Computing.
arXiv link: https://arxiv.org/abs/2011.11746