Limiting and divergence cleaning for continuous finite element discretizations of the MHD equations

doi:10.1016/j.jcp.2020.109230

Journal of Computational Physics

Volume 407, 15 April 2020, 109230

https://doi.org/10.1016/j.jcp.2020.109230 Get rights and content

Highlights

•
Localized element-based FCT algorithm for continuous finite element discretizations of the MHD system.
•
Synchronized and sequential a priori limiters based on local maximum principles for control variables.
•
Synchronized a posteriori limiting to enforce positivity preservation for the thermodynamic pressure.
•
Built-in penalization of divergence errors without violating maximum principles or conservation laws.

Abstract

This work introduces a new type of constrained algebraic stabilization for continuous piecewise-linear finite element approximations to the equations of ideal magnetohydrodynamics (MHD). At the first step of the proposed flux-corrected transport (FCT) algorithm, the Galerkin element matrices are modified by adding graph viscosity proportional to the fastest characteristic wave speed. At the second step, limited antidiffusive corrections are applied and divergence cleaning is performed for the magnetic field. The limiting procedure developed for this stage is designed to enforce local maximum principles, as well as positivity preservation for the density and thermodynamic pressure. Additionally, it adjusts the magnetic field in a way which penalizes divergence errors without violating conservation laws or positivity constraints. Numerical studies for 2D test problems are performed to demonstrate the ability of the proposed algorithms to accomplish this task in applications to ideal MHD benchmarks.

Introduction

In recent years, significant advances have been made in the analysis and design of physics-compatible finite element methods for nonlinear hyperbolic systems. In particular, the principle of invariant domain preservation was introduced by Guermond and Popov [17] as a useful generalization of scalar discrete maximum principles to systems of conserved quantities [17]. Artificial diffusion operators and localized edge-based flux-corrected transport (FCT) algorithms based on this design criterion were proposed in [18], [19], [20]. Similar algebraic approaches to construction of FCT-like finite element schemes for systems of conservation laws were explored in [8], [13], [25], [28]. The synchronized and sequential limiting techniques developed in these publications make it possible to constrain a set of coupled conserved variables in a manner which guarantees the validity of relevant maximum principles for nonlinear functions of these variables. In the context of ideal MHD, a well-designed limiter must guarantee positivity preservation for the density and pressure (internal energy).

Custom-made positivity-preserving limiters for the ideal MHD system were proposed in [4], [8], [10]. Similarly to limiters for the Euler equations, their design relies on the assumption that the underlying low-order scheme produces physically admissible solutions. As shown by Wu and Shu [42], [43], [44], the validity of this assumption cannot be guaranteed for Riemann problems with discontinuous magnetic fields and numerical approximations violating a discrete form of the divergence-free condition. However, a rigorous proof of positivity preservation could be obtained for the nonconservative Godunov-Powell form [34] of the MHD system discretized using a discontinuous Galerkin (DG) method [43].

Another important recent development is the advent of Multi-dimensional Optimal Order Detection (MOOD) methods based on the idea of a posteriori limiting [30], [41], [45], [46]. Instead of constraining numerical fluxes, element contributions, or gradients to keep the quantities of interest in the admissible range, MOOD approaches evolve the numerical solution using a standard high-order method. If violations of physical or numerical admissibility conditions are detected in any mesh cell, the restriction of the high-order candidate solution to this cell is rejected and an admissible approximation is calculated using a bound-preserving subcell finite volume scheme. Such a posteriori fixes are well suited for finite volume and discontinuous Galerkin methods. The MOOD paradigm can also be used for selective activation of artificial viscosity or adaptive choice of the limiting strategy for continuous finite elements (cf. [2], [26]).

A robust scheme for solving the MHD equations must incorporate a mechanism for avoiding (unbounded growth of) divergence errors in numerical approximations to the magnetic field [9], [39], [40]. This requirement can be satisfied, e.g., by using constrained transport (CT) schemes [6], [14], [21], [35], mixed finite elements [8], [23], penalizing source terms [34], or various divergence cleaning procedures [12]. Many ways of keeping the magnetic field approximately divergence-free are closely related to numerical methods for the incompressible Navier-Stokes equations (projection schemes, grad-div stabilization, artificial compressibility methods). Unfortunately, most of them lead to numerical schemes that are nonconservative or may produce spurious undershoots/overshoots.

In this paper, we present a new localized version of the element-based FCT algorithm developed in [8]. The proposed methodology constrains the numerical solution of the MHD system to stay in the admissible set using a priori limiting to enforce local maximum principles for the density and other quantities of interest. Additional a posteriori limiting is performed if negative values of the thermodynamic pressure are detected. At the a priori limiting stage, changes of the conserved quantities are constrained in a synchronized or sequential manner. The calculation of correction factors for the synchronized FCT limiter is based solely on the density constraints. The sequential FCT algorithm is based on the methodology developed in [13], [22] for continuous and discontinuous Galerkin discretizations of the Euler equations. Imposing local maximum principles on the density, velocity, specific total energy, and magnetic field, it uses different correction factors for different variables. Both versions of the a priori limiter guarantee positivity preservation for the density field. The positivity-preserving pressure fix is based on a simplified version of the synchronized limiter developed in [8]. In contrast to fully synchronized a priori limiters for MHD, the new FCT scheme does not automatically smear the magnetic field in regions of constant density and vice versa. Moreover, it features embedded divergence cleaning based on localized grad-div penalization. The corresponding penalty terms are built into the antidiffusive element contributions to the magnetic field, leading to a conservative and bound-preserving correction procedure. The numerical results for standard MHD benchmarks illustrate the shock-capturing and divergence cleaning capabilities of the proposed FCT method.

Section snippets

Continuous finite element discretization

We consider the conservative form of the ideal compressible MHD equations $\frac{\partial U}{\partial t} + \nabla \cdot F (U) = 0,$ where $U = [\begin{matrix} ρ \\ ρ u \\ ρ E \\ B \end{matrix}], F (U) = [\begin{matrix} ρ u \\ ρ u \otimes u + p_{tot} I - B \otimes B \\ ρ E u + p_{tot} u - B (u \cdot B) \\ u \otimes B - B \otimes u \end{matrix}] .$ The conserved quantities are the density ρ, momentum $ρ u$ , total energy ρE, and the magnetic field B. The velocity $u = \frac{ρ u}{ρ}$ and total pressure $p_{tot} = p + \frac{1}{2} | B |^{2}$ are derived quantities. The thermal pressure p is related to the density and internal energy by an equation of state. The EOS for ideal gases reads $p = (γ - 1) (ρ E - \frac{ρ | u |^{2}}{2} - \frac{1}{2} | B |^{2}),$ where γ is the

Low-order artificial viscosity

Since $U_{h}$ and $F_{h}$ are continuous on the internal faces of the finite element mesh, upwinding techniques based on MHD Riemann solvers [3], [4], [5] cannot be used to manipulate the numerical fluxes directly. Following [8], [18], [26], we will stabilize the lumped-mass Galerkin system (18) by adding some artificial viscosity. The resulting low-order approximation will be of the form $M_{L} \frac{d U}{d t} = - C \cdot F + D U + G,$ where $D = {d_{i j}}_{i, j = 1}^{N_{h}}$ is a discrete diffusion operator with entries $d_{i j} = \sum_{e \in E_{i} \cap E_{j}} d_{i j}^{e} .$

Let $i \in {1, \dots, N_{h}}$ be

FCT correction and divergence cleaning

At the second stage of our flux-corrected transport (FCT) algorithm, we will correct the low-order approximation $U_{i}^{L} \in G$ by adding limited element contributions which should remove excessive artificial diffusion, reduce phase errors due to mass lumping, and penalize divergence errors in the finite element approximation to the magnetic field. The outcome of this correction will be $U_{i}^{C} = U_{i}^{L} + \frac{1}{m_{i}} \sum_{e \in E_{i}} m_{i}^{e} F_{i}^{e, C} = \frac{1}{\sum_{e \in E_{i}} m_{i}^{e}} \sum_{e \in E_{i}} m_{i}^{e} U_{i}^{e, C},$ where $U_{i}^{e, C} = U_{i}^{L} + F_{i}^{e, C}$ and $F_{i}^{e, C}$ are element contributions satisfying

Limiting procedures for MHD

Let $F_{i}^{e, T} = {[f_{i}^{e, ρ}, f_{i}^{e, ρ u}, f_{i}^{e, ρ E}, f_{i}^{e, B}]}^{T}$ be defined by (42). For $F_{i}^{e, C} = F_{i}^{e, T}$ , the FCT correction step (39) will produce the target approximation $ρ_{i}^{e, T} = ρ_{i}^{L} + f_{i}^{e, ρ},$ ${(ρ u)}_{i}^{e, T} = {(ρ u)}_{i}^{L} + f_{i}^{e, ρ u},$ ${(ρ E)}_{i}^{e, T} = {(ρ E)}_{i}^{L} + f_{i}^{e, ρ E},$ $B_{i}^{e, T} = B_{i}^{L} + f_{i}^{e, B} .$ In general, the FCT limiter for $F_{i}^{e, C}$ should be designed to guarantee that $U_{i}^{e, C} (α^{e}) = {[ρ_{i}^{e}, {(ρ u)}_{i}^{e}, {(ρ E)}_{i}^{e}, B_{i}^{e}]}^{T}$ is physically admissible, i.e., $U_{i}^{e, C} \in G$ . Additionally, scalar quantities of interest may be constrained to satisfy local maximum principles of the form $v_{i}^{\min} \leq v (U_{i}^{e, C} (α^{e})) \leq$

Summary of the FCT algorithm

In summary, the new FCT scheme for calculating constrained approximations $U_{i}^{C} \in G$ to the solution of system (1) involves the following steps:

1.
Advance U in time using the low-order method presented in Section 3.
2.
For each element $e \in {1, \dots, E_{h}}$ , calculate the nodal states ${\tilde{U}}_{i}^{e, C} \in {\tilde{G}}_{i}$ using the synchronized a priori limiter presented in Section 5.1 or the sequential a priori limiter presented in Section 5.2.
3.
If ${\tilde{U}}_{i}^{e, C} \in G$ , then set $U_{i}^{e, C} = {\tilde{U}}_{i}^{e, C}$ . Otherwise, calculate $U_{i}^{e, C} \in G$ using the synchronized a posteriori

Numerical examples

In this section, we perform numerical studies for two standard MHD benchmarks. The methods under investigation are abbreviated as follows:

LO	low-order scheme (Section 3) without any corrections;
FCT-A	localized element-based FCT with synchronized a priori limiting
FCT-A	(Section 5.1) and synchronized a posteriori limiting (Section 5.3);
FCT-B	localized element-based FCT with sequential a priori limiting
FCT-B	(Section 5.2) and synchronized a posteriori limiting (Section 5.3).

The two-dimensional implementation of

Conclusions

Building on recent advances in the field of bound-preserving high-resolution finite element schemes for hyperbolic systems [13], [18], [20], [22], [43], the proposed upgrade of the FCT algorithm developed in [8] equips it with more advanced limiting techniques. The presented new features include a positivity-preserving correction of the low-order predictor and sequential limiting subject to global energy constraints. The use of directional maximum principles for the momentum and magnetic field

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper

Acknowledgements

This research was supported by the German Research Association (DFG) under grant KU 1530/15-2. The authors would like to thank Dr. John N. Shadid and Dr. Sibusiso Mabuza (Sandia National Laboratories) for inspiring discussions of limiting techniques for MHD problems.

References (46)

D.S. Balsara
Multidimensional HLLE Riemann solver: application to Euler and magnetohydrodynamic flows
J. Comput. Phys.
(2010)
D.S. Balsara
Self-adjusting, positivity preserving high order schemes for hydrodynamics and magnetohydrodynamics
J. Comput. Phys.
(2012)
D.S. Balsara et al.
Multidimensional HLLC Riemann solver for unstructured meshes — with application to Euler and MHD flows
J. Comput. Phys.
(2014)
D.S. Balsara et al.
A staggered mesh algorithm using high order Godunov fluxes to ensure solenoidal magnetic fields in magnetohydrodynamic simulations
J. Comput. Phys.
(1999)
M. Basting et al.
An FCT finite element scheme for ideal MHD equations in 1D and 2D
J. Comput. Phys.
(2017)
Y. Cheng et al.
Positivity-preserving DG and central DG methods for ideal MHD equations
J. Comput. Phys.
(2013)
C.J. Cotter et al.
Embedded discontinuous Galerkin transport schemes with localised limiters
J. Comput. Phys.
(2016)
A. Dedner et al.
Hyperbolic divergence cleaning for the MHD equations
J. Comput. Phys.
(2002)
V. Dobrev et al.
Sequential limiting in continuous and discontinuous Galerkin methods for the Euler equations
J. Comput. Phys.
(2018)
F. Frank et al.
Bound-preserving flux limiting schemes for DG discretizations of conservation laws with applications to the Cahn-Hilliard equation
Comput. Methods Appl. Mech. Eng.
(2020)

J.-L. Guermond et al.

Invariant domain preserving discretization-independent schemes and convex limiting for hyperbolic systems

Comput. Methods Appl. Mech. Eng.

(2019)

C. Helzel et al.

An unstaggered constrained transport method for the 3D ideal magnetohydrodynamic equations

J. Comput. Phys.

(2011)

H. Hajduk et al.

New directional vector limiters for discontinuous Galerkin methods

J. Comput. Phys.

(2019)

D. Kuzmin et al.

A flux-corrected transport algorithm for handling the close-packing limit in dense suspensions

J. Comput. Appl. Math.

(2012)

D. Kuzmin et al.

Failsafe flux limiting and constrained data projections for equations of gas dynamics

J. Comput. Phys.

(2010)

C. Lohmann et al.

Synchronized flux limiting for gas dynamics variables

J. Comput. Phys.

(2016)

C. Lohmann et al.

Flux-corrected transport algorithms for continuous Galerkin methods based on high order Bernstein finite elements

J. Comput. Phys.

(2017)

M. Olshanskii et al.

Grad–div stabilization and subgrid pressure models for the incompressible Navier–Stokes equations

Comput. Methods Appl. Mech. Eng.

(2009)

K.G. Powell et al.

A solution-adaptive upwind scheme for ideal magnetohydrodynamics

J. Comput. Phys.

(1999)

V. Selmin

The node-centred finite volume approach: bridge between finite differences and finite elements

Comput. Methods Appl. Mech. Eng.

(1993)

G. Tóth

The $\nabla \cdot B = 0$ constraint in shock-capturing magnetohydrodynamics codes

J. Comput. Phys.

(2000)

G. Tóth et al.

Comparison of some flux corrected transport and total variation diminishing numerical schemes for hydrodynamic and magnetohydrodynamic problems

J. Comput. Phys.

(1996)

F. Vilar

A posteriori correction of high-order discontinuous Galerkin scheme through subcell finite volume formulation and flux reconstruction

J. Comput. Phys.

(2019)

Cited by (20)

Dissipation-based WENO stabilization of high-order finite element methods for scalar conservation laws
2023, Journal of Computational Physics
We present a new perspective on the use of weighted essentially nonoscillatory (WENO) reconstructions in high-order methods for scalar hyperbolic conservation laws. The main focus of this work is on nonlinear stabilization of continuous Galerkin (CG) approximations. The proposed methodology also provides an interesting alternative to WENO-based limiters for discontinuous Galerkin (DG) methods. Unlike Runge–Kutta DG schemes that overwrite finite element solutions with WENO reconstructions, our approach uses a reconstruction-based smoothness sensor to blend the numerical viscosity operators of high- and low-order stabilization terms. The so-defined WENO approximation introduces low-order nonlinear diffusion in the vicinity of shocks, while preserving the high-order accuracy of a linearly stable baseline discretization in regions where the exact solution is sufficiently smooth. The underlying reconstruction procedure performs Hermite interpolation on stencils consisting of a mesh cell and its neighbors. The amount of numerical dissipation depends on the relative differences between partial derivatives of reconstructed candidate polynomials and those of the underlying finite element approximation. All derivatives are taken into account by the employed smoothness sensor. To assess the accuracy of our CG-WENO scheme, we derive error estimates and perform numerical experiments. In particular, we prove that the consistency error of the nonlinear stabilization is of the order $p + 1 / 2$ , where p is the polynomial degree. This estimate is optimal for general meshes. For uniform meshes and smooth exact solutions, the experimentally observed rate of convergence is as high as $p + 1$ .
Subcell limiting strategies for discontinuous Galerkin spectral element methods
2022, Computers and Fluids
Citation Excerpt :
Although the non-conservative form is necessary to show entropy stability [68], we use the conservative form in this work since it is more convenient to cast in the framework of the present collection of subcell limiting strategies. The Orszag–Tang Vortex is an inviscid 2D MHD case that was originally proposed by Orszag and Tang [71], which is widely used to test the robustness of MHD codes [68,72–74]. The simulation starts from a smooth initial condition, which evolves into complex shock patterns with several shock–shock interactions and transitions to MHD turbulence.
We present a general family of subcell limiting strategies to construct robust high-order accurate nodal discontinuous Galerkin (DG) schemes. The main strategy is to construct compatible low order finite volume (FV) type discretizations that allow for convex blending with the high-order variant with the goal of guaranteeing additional properties, such as bounds on physical quantities and/or guaranteed entropy dissipation. For an implementation of this main strategy, four main ingredients are identified that may be combined in a flexible manner: (i) a nodal high-order DG method on Legendre–Gauss–Lobatto nodes, (ii) a compatible robust subcell FV scheme, (iii) a convex combination strategy for the two schemes, which can be element-wise or subcell-wise, and (iv) a strategy to compute the convex blending factors, which can be either based on heuristic troubled-cell indicators, or using ideas from flux-corrected transport methods.
By carefully designing the metric terms of the subcell FV method, the resulting methods can be used on unstructured curvilinear meshes, are locally conservative, can handle strong shocks efficiently while directly guaranteeing physical bounds on quantities such as density, pressure or entropy. We further show that it is possible to choose the four ingredients to recover existing methods such as a provably entropy dissipative subcell shock-capturing approach or a sparse invariant domain preserving approach.
We test the versatility of the presented strategies and mix and match the four ingredients to solve challenging simulation setups, such as the KPP problem (a hyperbolic conservation law with non-convex flux function), turbulent and hypersonic Euler simulations, and MHD problems featuring shocks and turbulence.
An implicit monolithic AFC stabilization method for the CG finite element discretization of the fully-ionized ideal multifluid electromagnetic plasma system
2022, Journal of Computational Physics
Citation Excerpt :
The most common modifications result in a set of purely-hyperbolic equations that are more amenable to fully-explicit or IMEX time discretizations [1,3,30,35,55,83], while others yield a mixed hyperbolic-parabolic or hyperbolic-elliptic structure [22]. An eliminated parabolic divergence cleaning method for CG discretizations has recently been shown to be effective for the solution of ideal MHD equations [41,58]. This approach is leveraged in the current work for the solution of Maxwell's equations using a nodal CG discretization.
This study considers an implicit finite element formulation for an ideal fully-ionized multifluid electromagnetic plasma system. The formulation is based on fully-implicit Runge-Kutta time discretizations and a monolithic discrete algebraic flux corrected (AFC) continuous Galerkin (CG) spatial discretization of the coupled system. The AFC approach adds scalar artificial diffusion to the high-order, semi-discrete Galerkin method and uses mass lumping in the time derivative term. The result is a low-order method that attempts to enforce local-extremum-diminishing properties for the hyperbolic system. An element-based iterative limiter is applied to reduce the amount of artificial diffusion that is used in regions where the solution is smooth and the additional stabilization is not required. Two models are considered for the electromagnetics portion of the system: an electrostatic model, and a full Maxwell system with a parabolic divergence cleaning approach that enforces the required involutions on the electric and magnetic fields. Results are presented that demonstrate the accuracy and robustness of the formulation for smooth and discontinuous solutions to challenging plasma physics problems. This includes a demonstration that the solution of the full multifluid system yields the expected behavior in the ideal shock-MHD limit.
Monolithic parabolic regularization of the MHD equations and entropy principles
2022, Computer Methods in Applied Mechanics and Engineering
We show at the PDE level that the monolithic parabolic regularization of the equations of ideal magnetohydrodynamics (MHD) is compatible with all the generalized entropies, fulfills the minimum entropy principle, and preserves the positivity of density and internal energy. We then numerically investigate this regularization for the MHD equations using continuous finite elements in space and explicit strong stability preserving Runge–Kutta methods in time. The artificial viscosity coefficient of the regularization term is constructed to be proportional to the entropy residual of MHD. It is shown that the method has a high order of accuracy for smooth problems and captures strong shocks and discontinuities accurately for non-smooth problems.
A linearity preserving nodal variation limiting algorithm for continuous Galerkin discretization of ideal MHD equations
2020, Journal of Computational Physics
In this work, a stabilized continuous Galerkin (CG) method for magnetohydrodynamics (MHD) is presented. Ideal, compressible inviscid MHD equations are discretized in space on unstructured meshes using piecewise linear or bilinear finite element bases to get a semi-discrete scheme. Stabilization is then introduced to the semi-discrete method in a strategy that follows the algebraic flux correction paradigm. This involves adding some artificial diffusion to the high order, semi-discrete method and mass lumping in the time derivative term. The result is a low order method that provides local extremum diminishing properties for hyperbolic systems. The difference between the low order method and the high order method is scaled element-wise using a limiter and added to the low order scheme. The limiter is solution dependent and computed via an iterative linearity preserving nodal variation limiting strategy. The stabilization also involves an optional consistent background high order dissipation that reduces phase errors. The resulting stabilized scheme is a semi-discrete method that can be applied to inviscid shock MHD problems and may be even extended to resistive and viscous MHD problems. To satisfy the divergence free constraint of the MHD equations, we add parabolic divergence cleaning to the system. Various time integration methods can be used to discretize the scheme in time. We demonstrate the robustness of the scheme by solving several shock MHD problems.
Monolithic convex limiting for continuous finite element discretizations of hyperbolic conservation laws
2020, Computer Methods in Applied Mechanics and Engineering
Citation Excerpt :
The AFC formalism presented in [16,30] provides algebraic interpretations of local extremum diminishing (LED) schemes [31–33] that use artificial diffusion operators and flux limiters to enforce discrete maximum principles. Localized element-based limiting procedures for scalar conservation laws and hyperbolic systems were proposed in [15,17,34–36]. The way in which the limiter is localized in element-based and edge-based FCT schemes of this kind forms the basis of what Guermond et al. [4,27] named convex limiting in the context of second-order invariant domain preserving schemes for hyperbolic systems.
Using the theoretical framework of algebraic flux correction and invariant domain preserving schemes, we introduce a monolithic approach to convex limiting in continuous finite element schemes for linear advection equations, nonlinear scalar conservation laws, and hyperbolic systems. In contrast to flux-corrected transport (FCT) algorithms that apply limited antidiffusive corrections to bound-preserving low-order solutions, our new limiting strategy exploits the fact that these solutions can be expressed as convex combinations of bar states belonging to a convex invariant set of physically admissible solutions. Each antidiffusive flux is limited in a way which guarantees that the associated bar state remains in the invariant set and preserves appropriate local bounds. There is no free parameter and no need for limit fluxes associated with the consistent mass matrix of time derivative term separately. Moreover, the steady-state limit of the nonlinear discrete problem is well defined and independent of the pseudo-time step. In the case study for the Euler equations, the components of the bar states are constrained sequentially to satisfy local maximum principles for the density, velocity, and specific total energy in addition to positivity preservation for the density and pressure. The results of numerical experiments for standard test problems illustrate the ability of built-in convex limiters to resolve steep fronts in a sharp and nonoscillatory manner.

View all citing articles on Scopus

View full text