An entropy stable nodal discontinuous Galerkin method for the resistive MHD equations: Continuous analysis and GLM divergence cleaning
Abstract
This work presents an extension of discretely entropy stable discontinuous Galerkin (DG) methods to the resistive magnetohydrodynamics (MHD) equations. Although similar to the compressible NavierStokes equations at first sight, there are some important differences concerning the resistive MHD equations that need special focus. The continuous entropy analysis of the ideal MHD equations, which are the advective parts of the resistive MHD equations, shows that the divergencefree constraint on the magnetic field components must be incorporated as a nonconservative term in a form either proposed by Powell or Janhunen. Consequently, this nonconservative term needs to be discretized, such that the approximation is consistent with the entropy. As an extension of the ideal MHD system, we address in this work the continuous analysis of the resistive MHD equations and show that the entropy inequality holds. Thus, our first contribution is the proof that the resistive terms are symmetric and positive semidefinite when formulated in entropy space as gradients of the entropy variables. Moreover, this enables the construction of an entropy stable DG discretization for the resistive MHD equations. However, the resulting method suffers from large errors in the divergencefree constraint, since no particular treatment of divergence errors is included in the standard resistive MHD model. Hence, our second contribution is the extension of the resistive MHD equations with proper divergence cleaning based on a generalized Lagrange multiplier (GLM) strategy. We construct and analyze a DG method that is entropy stable for the resistive MHD equations and has a builtin GLM divergence cleaning mechanism. The theoretical derivations and proofs are then verified by several numerical examples. Specifically, we demonstrate that the entropy stable method with GLM divergence cleaning is more robust than the standard highorder DG method for the resistive MHD equations.
keywords:
magnetohydrodynamics, resistivity, entropy stability, discontinuous Galerkin spectral element method, hyperbolic divergence cleaning, summationbyparts1 Introduction
The resistive magnetohydrodynamic (MHD) equations are of great interest in many areas of plasma, space and astrophysics. This stems from a wide range of applications such as electromagnetic turbulence in conducting fluids, magnetically confined fusion for power generation, modeling the action of dynamos and predicting the interaction of the solar wind with planets or moons. The governing equations are able to describe both dense and sparse plasmas that are timedependent with motions that feature a wide range of temporal and spatial scales, e.g., compressible MHD turbulence. In addition, the resistive MHD equations exhibit a mixed hyperbolicparabolic character depending on the strength of the viscous or resistive effects. Another important property, in a closed physical system is the second law of thermodynamics, i.e., the evolution of the entropy. In the absence of resistivity and viscosity, that is for the ideal MHD model, and if the resulting solutions are smooth, the entropy of the system is an additional conserved quantity, although not explicitly built into the mathematical model. Further, in the presence of shocks, the second law of thermodynamics becomes the entropy inequality, e.g. [26], which guarantees that entropy is always dissipated with the correct sign. It is assumed that the additional resistive terms have a pure entropy dissipative effect as well. But, to the best of our knowledge, no continuous entropy analysis of the resistive MHD equations has been presented in literature and it is unclear if the entropy inequality holds for the full resistive MHD equations.
The discrete analog to the entropy inequality of the continuous partial differential equations (PDE) is entropy stability. Constructing such entropy stable methods for the advective parts of the MHD model has been studied by many authors, e.g., [2, 8, 35, 43, 54]. A complication when discussing entropy stability for MHD models is the necessity of the involution condition, that is, the divergencefree constraint of the magnetic field [2, 25]
(1.1) 
The condition (1.1), like entropy stability, is an additional PDE not explicitly built into the resistive MHD equations. To cancel extraneous terms in the entropy analysis requires the addition of a nonconservative term that is proportional to the divergence of the magnetic field, as first introduced by Godunov [25]. In the context of numerical schemes, several such nonconservative terms have been proposed over the years by Brackbill and Barnes [5], Powell [41] and Janhunen [29]. On the continuous level, adding a nonconservative term scaled by (1.1) is a clever way of adding zero to the model. However, for numerical approximations, there are known stability and accuracy issues that differ between the three types of nonconservative terms [45].
Such differences in stability and accuracy are largely due to the discrete satisfaction of the divergencefree condition. It is well known in the MHD numerics community that even if the initial conditions of a problem satisfy (1.1), it is not guaranteed that the discrete evolution of a plasma will remain divergencefree in the magnetic field. Therefore, many numerical techniques have been devised to control errors introduced into the divergencefree constraint by a numerical approximation. These include adding the aforementioned nonconservative term, the projection approach of Brackbill and Barnes [5], the method of constrained transport introduced by Evans and Hawley [14] and the generalized Lagrange multiplier (GLM) hyperbolic divergence cleaning technique originally proposed for the ideal MHD equations by Dedner et al. [10]. A thorough review of the behavior and efficacy of these techniques, except hyperbolic divergence cleaning, is provided by Tóth [51]. Due to its relative ease of implementation and computational efficiency we are most interested in the method of hyperbolic divergence cleaning. However, the current work is also concerned with constructing entropy stable numerical approximations. Recent work by Derigs et al. [11] modified the additional GLM divergence cleaning system in such a way that the resulting ideal GLMMHD system is consistent with the continuous entropy analysis and provides inbuilt divergence cleaning capabilities.
The goal of this work is to construct an entropy stable highorder discontinuous Galerkin (DG) method for the modified GLMMHD system augmented with resistive terms to obtain a robust, highorder accurate numerical discretization with inbuilt divergence cleaning. Recently, the construction and analysis of discretely entropy stable DG methods for nonlinear conservation laws have made major progress [6, 24, 35]. The key to discrete entropy stability is to mimic the integrationbyparts property in the DG operators with the socalled summationbyparts (SBP) property. This enables the construction of DG methods that are entropy stable without the assumption of exact evaluation of the variational forms. The reason we are particularly interested in entropy stability for our DG method is the expected increase in robustness. We know that for the compressible NavierStokes equations an entropy stable DG method has drastically increased robustness in comparison to the standard DG and even in comparison to DG methods with dealiasing, e.g. [6, 24]. Solutions of the resistive MHD equations can be even more complicated than solutions to the compressible NavierStokes equations. As such, we aim for the construction of a highorder method that is able to capture the wide range of scales in a robust and stable way.
However, the construction of entropy stable highorder DG methods for the resistive MHD equations involves several key differences in comparison to standard conservation laws. One major difference is the addition of a nonconservative PDE term proportional to the discrete divergence of the magnetic field. Typically, entropy stable discontinuous Galerkin spectral element methods (DGSEM) are built for conservation laws. Thus, precise and careful numerical treatment of the nonconservative term is necessary. We note that there is recent work on entropy stable DG methods applied to the ideal MHD equations by Rossmanith [43], Liu et al. [35] and GallegoValencia [20]. In particular, we extend the ideas presented in Liu et al. [35] and provide clarification and motivation to the structure of the presented discretization of the nonconservative term and add detailed proofs of the 2D extension. As already mentioned, another important difference is the numerical control of the divergencefree constraint. To do so, we incorporate the modified GLM divergence cleaning approach in our full highorder DG entropy analysis. Furthermore, we build on the recent work of Gassner et al. [23] who showed that using a BassiRebay (BR1) type scheme [3] is entropy stable for the viscous terms of the compressible NavierStokes equations, if the scheme is formulated with the gradients of the entropy variables instead of the gradients of the conservative or primitive variables. To use this proof, we reformulate the resistive terms by the gradients of the entropy variables and show that the resulting resistive coefficient matrices are symmetric and positive semidefinite.
The remainder of this paper is organized as follows: Our first main result, in Sec. 2, is the continuous entropy analysis of the threedimensional resistive GLMMHD equations, which demonstrates that the model indeed satisfies the entropy inequality and that the resistive terms can be recast into a symmetric and positive semidefinite form. The discontinuous Galerkin spectral element method (DGSEM) is outlined in Sec. 3. The second main result of this work, in Sec. 4, is to demonstrate the entropy stability of the numerical approximation by discretely mimicking the continuous entropy analysis with special attention given to the GLM divergence cleaning and resistive parts. We present numerical results in Sec. 5 to verify the highorder nature and theoretical entropy conservation (without viscous terms) and entropy stability (including viscous terms) of the derived scheme. We furthermore demonstrate the increased robustness of the novel DG method. Conclusion remarks are given in the final section.
2 Continuous entropy analysis
As a starting point for the discussion and numerical approximations in this work we outline the resistive GLMMHD equations. In general, we consider systems of conservation laws in a domain defined as
(2.1) 
where denotes the vector of conserved variables and the multidimensional flux vector. These definitions allow for a compact notation that will simplify the analysis, i.e., we define block vectors with the double arrow as
(2.2) 
and the spatial gradient of a state as
(2.3) 
The gradient of a spatial vector is a second order tensor, written in matrix form as
(2.4) 
and the divergence of a flux written as a block vector is defined as
(2.5) 
2.1 Resistive GLMMHD equations
The equations that govern the evolution of resistive, conducting fluids depend on the solution as well as its gradient [12, 15, 52]. In the advective components we also include the recent modifications proposed by Derigs et al. [11] to incorporate a GLM divergence cleaning framework into the model. Thus, the complete mathematical model considered in this work can be defined in the domain by
(2.6) 
with as well as the advective and viscous fluxes:
(2.7) 
Here, are the mass density, fluid velocities, pressure and total energy, respectively, denotes the magnetic field components and the identity matrix. Furthermore, the viscous stress tensor reads [31]
(2.8) 
and the heat flux is defined as
(2.9) 
The introduced constants describe the viscosity from the NavierStokes equations, resistivity of the plasma, thermal conductivity and the universal gas constant, respectively. In particular, the constants and are firstorder transport coefficients that describe the kinematic viscosity and the diffusivity of the magnetic field [55]. We close the system with the ideal gas assumption, which relates the total energy and pressure
(2.10) 
where denotes the adiabatic coefficient.
The system (2.6) contains the GLM extension indicated by the additional field variable , which controls the divergence error and propagates it out of the physical domain by the constant wave speed . We also introduce two nonconservative terms: The first one is scaled by and, thus, zero in the continuous case. The particular choice of is related to the thermodynamic properties of (2.6) as it ensures the entropy conservation for the advective part of the system outlined in the next section. In the original approach to derive thermodynamically consistent methods for the ideal MHD equations is chosen to be the Powell source term [2, 41], which also recovers the symmetrization of the system [25]. In this work we restrict to the alternative Janhunen source term [29], for two reasons:

It is sufficient to restore the satisfaction of the second law of thermodynamics even when , e.g. [54], which is important in the construction of entropy stable numerical methods.

The equations remain conservative in all the hydrodynamic variables so the satisfaction of the first law of thermodynamics, i.e. conservation of the total energy, is unchanged.
In particular, the components of the Janhunen source term are
(2.11) 
We note that because on the continuous level the modification of to the resistive GLMMHD system is consistent.
2.2 Thermodynamic properties of the system
To discuss the thermodynamic properties of the resistive GLMMHD equations (2.6) we translate the concepts of the first and second law of thermodynamics into a mathematical context. To do so, we first exclusively examine the advective and nonconservative term proportional to (1.1). The ideal GLMMHD equations satisfy the first law of thermodynamics, because the evolution of the total fluid energy is one of the conserved quantities. This is true for any choice of the vector because (1.1) is assumed to hold in the continuous analysis. But, on the discrete level, this is not the case as noted by many authors [11, 41, 29, 51]. However, the mathematical description of the second law of thermodynamics is more subtle, because the entropy is not explicitly built into the system. Thus, we require more formalism and utilize the welldeveloped entropy analysis tools for hyperbolic systems, e.g. [26, 46, 38]. As such, we define a strongly convex entropy function that is then used to define an injective mapping between state space and entropy space. Note that we adopt the mathematical notation that the entropy function is negative as is often done, e.g. [26, 48]. Once the advective terms are accounted for in entropy space, we demonstrate the first main result of this work: the resistive terms do not violate the second law of thermodynamics.
For the ideal and the resistive GLMMHD equations a suitable entropy function is the entropy density divided by the constant for convenience as
(2.13) 
where is the physical entropy [32]. From the entropy function we define the entropy variables to be
(2.14) 
and we introduce , which is proportional to the inverse temperature.
It is known, that for smooth solutions when we contract the ideal GLMMHD equations without viscous fluxes nor any source term on the right hand side by the entropy variables (2.14) we obtain the entropy conservation law [11]
(2.15) 
where we define the entropy fluxes to be
(2.16) 
Additionally, as it will be necessary in later derivations and the proof of discrete entropy stability, we compute the entropy flux potential to be
(2.17) 
where we introduce notation for the nonconservative state vector (2.11) contracted into entropy space
(2.18) 
We note that, if the Powell nonconservative state vector [41] was used, the result (2.18) still holds [11].
Furthermore, e.g. in case of shock discontinuities, the solution satisfies the following entropy inequality
(2.19) 
which is the mathematical description of the second law of thermodynamics for the ideal MHD equations. Next, we account for the resistive terms to demonstrate the entropy inequality for the resistive GLMMHD equations. To do so, we require a suitable representation to discuss how the resistive terms affect (2.19).
Lemma 0 (Entropy representation of viscous fluxes).
The viscous fluxes of the resistive GLMMHD equations in (2.7) can be expressed by gradients of the entropy variables as
(2.20) 
with a block matrix that is symmetric and positive semidefinite matrix, i.e,
(2.21) 
Proof.
First, we consider the viscous fluxes of the resistive GLMMHD system in (2.7) separately
(2.22) 
with
where
Using the vector of entropy variables from (2.14)
we find the following relations:
With some algebraic effort we can determine the matrices to express the viscous fluxes in terms of matrices times the gradients of entropy variables:
(2.23)  
(2.24)  
(2.25) 
We collect all these block matrices into the matrix
(2.26) 
which clearly yields
(2.27) 
For clarification, we present the first matrix
(2.28) 
The other matrices are explicitly stated in A. It is straightforward to verify that the matrix is symmetric by inspecting the block matrices listed in (2.28) and (A.1)  (A.8) where the following relationships hold
(2.29) 
To show that the matrix is positive semidefinite is more involved. We first note that it is possible to split the matrix (2.26) into the viscous terms associated with the NavierStokes equations and the resistive terms of the magnetic fields that arise in the resistive GLMMHD equations. We exploit this fact and rewrite the total diffusion matrix into two pieces
(2.30) 
where all terms with are put in and all terms with are in . It is easy to verify that the NS and RMHD block matrices are symmetric, as both satisfy (2.29). A further convenience is that the NavierStokes part, , is known to be positive semidefinite [13]
(2.31) 
Thus, all that remains is to demonstrate that the additional resistive dissipation matrix, , is positive semidefinite. To do so we examine the eigenvalues of the system. We use the computer algebra system Maxima [36] to find an explicit expression of the eigenvalues to be
(2.32) 
Under the physical assumptions that and we see that the eigenvalues (2.32) of the matrix are all nonnegative. Hence, the block matrix is symmetric and positive semidefinite. ∎
With the ability to rewrite the viscous fluxes as a linear combination of the entropy variable gradients, we present our first main result:
Theorem 2 (Entropy inequality for the resistive GLMMHD equations).
Proof.
We start by contracting the resistive GLMMHD system (2.6) with the entropy variables and integrate over the domain:
(2.33) 
Here, denotes the inner product, e.g.,
(2.34) 
From the definition of the entropy variables we have
(2.35) 
Next, for clarity, we separate the advective flux into Euler, ideal MHD and GLM parts
(2.36) 
The Euler terms generate the divergence of the entropy flux, e.g., [26]
(2.37) 
the ideal MHD and nonconservative term cancel, e.g., [2, 35]
(2.38) 
and the GLM terms vanish directly with the modifications introduced in [11]
(2.39) 
The damping source term for the GLM divergence cleaning is zero in all but its ninth component, so we see
(2.40) 
Thus, we have
(2.41) 
Next, to demonstrate the entropy stability for the resistive GLMMHD equations, we address the viscous flux components. So, we integrate by parts for the remaining inner product including the viscous fluxes and obtain
(2.42) 
where denotes the outward pointing normal vector at the domain boundaries. Assuming these to be periodic and using Lemma 1, we arrive at
(2.43) 
Since is symmetric and positive semidefinite and , all inner products on the right hand side are positive and thus
(2.44) 
∎
Remark 0.
Splitting the flux contraction in the continuous entropy analysis is useful to keep track of which terms contribute to the entropy flux and which terms cancel. This will also be the case in the discrete entropy analysis of the highorder DG approximation. As such, the entropy flux potential is split into Euler, ideal MHD and GLM components
(2.45) 
where
(2.46)  
(2.47)  
(2.48) 
In summary, we have demonstrated that the resistive GLMMHD equations satisfy an entropy inequality. To do so, we separated the advective contributions into Euler, ideal MHD and GLM pieces and considered the viscous contributions separately, which served to clarify how each term contributed to the entropy analysis. A major result is that it is possible to rewrite the resistive terms of the threedimensional system in an entropy consistent way to demonstrate that those terms are entropy dissipative. We will use an identical splitting of the advective and diffusive terms in the discrete entropy stability proofs in Sec. 4 to directly mimic the continuous analysis. However, we first must build the necessary components of a discontinuous Galerkin approximation for the resistive GLMMHD equations.
3 Split form discontinous Galerkin approximation
The goal of this paper is to construct an entropy stable highorder DG method for the resistive GLMMHD system (2.6).
In this section, we introduce the building blocks of our nodal DG discretization. Most importantly we highlight that, as long as the nodal DG approximation is built with the LegendreGaussLobatto nodes, the discrete derivative matrix and the discrete mass matrix satisfy the summationbyparts (SBP) property for any polynomial order [22]. This is a key property as it allows us to use results from the work of Fisher et al. [17] and Fisher and Carpenter [16]. We follow the notation introduced in [24] and present a split form DG approximation, where we have two numerical fluxes, one at the surface and one inside the volume for the special split form volume integral. Carpenter and Fisher showed that when using an entropy conservative finite volume flux for the numerical volume flux in a SBP discretization, the property of entropy conservation extends to the highorder SBP method, i.e. to the higher order DG method.
In the following we introduce the socalled nodal discontinuous Galerkin collocation spectral element method (DGSEM). We use tensor product elements and derive the strong form including the SBP property of the volume discretization. Furthermore, we restrict the DGSEM discretization to a twodimensional Cartesian mesh with uniform squared elements. So, the first step in the discretization is to subdivide the computational domain into nonoverlapping square elements and map each of these elements to the reference element . The bilinear mappings are defined as
(3.1) 
for , where are the four corners of the element . Since we consider the uniform Cartesian case, we find the Jacobian of the mappings to be constant
(3.2) 
for equal element side lengths , .
Next, on every element, each conserved variable is approximated by a polynomial of degree in each spatial direction on