An entropy stable nodal discontinuous Galerkin method for the resistive MHD equations: Continuous analysis and GLM divergence cleaning
This work presents an extension of discretely entropy stable discontinuous Galerkin (DG) methods to the resistive magnetohydrodynamics (MHD) equations. Although similar to the compressible Navier-Stokes 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 divergence-free constraint on the magnetic field components must be incorporated as a non-conservative term in a form either proposed by Powell or Janhunen. Consequently, this non-conservative 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 semi-definite 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 divergence-free 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 built-in 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 high-order DG method for the resistive MHD equations.
keywords:magnetohydrodynamics, resistivity, entropy stability, discontinuous Galerkin spectral element method, hyperbolic divergence cleaning, summation-by-parts
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 time-dependent 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 hyperbolic-parabolic 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. , 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 divergence-free constraint of the magnetic field [2, 25]
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 non-conservative term that is proportional to the divergence of the magnetic field, as first introduced by Godunov . In the context of numerical schemes, several such non-conservative terms have been proposed over the years by Brackbill and Barnes , Powell  and Janhunen . On the continuous level, adding a non-conservative 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 non-conservative terms .
Such differences in stability and accuracy are largely due to the discrete satisfaction of the divergence-free 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 divergence-free in the magnetic field. Therefore, many numerical techniques have been devised to control errors introduced into the divergence-free constraint by a numerical approximation. These include adding the aforementioned non-conservative term, the projection approach of Brackbill and Barnes , the method of constrained transport introduced by Evans and Hawley  and the generalized Lagrange multiplier (GLM) hyperbolic divergence cleaning technique originally proposed for the ideal MHD equations by Dedner et al. . A thorough review of the behavior and efficacy of these techniques, except hyperbolic divergence cleaning, is provided by Tóth . 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.  modified the additional GLM divergence cleaning system in such a way that the resulting ideal GLM-MHD 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 high-order discontinuous Galerkin (DG) method for the modified GLM-MHD system augmented with resistive terms to obtain a robust, high-order accurate numerical discretization with inbuilt divergence cleaning. Recently, the construction and analysis of discretely entropy stable DG methods for non-linear conservation laws have made major progress [6, 24, 35]. The key to discrete entropy stability is to mimic the integration-by-parts property in the DG operators with the so-called summation-by-parts (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 Navier-Stokes equations an entropy stable DG method has drastically increased robustness in comparison to the standard DG and even in comparison to DG methods with de-aliasing, e.g. [6, 24]. Solutions of the resistive MHD equations can be even more complicated than solutions to the compressible Navier-Stokes equations. As such, we aim for the construction of a high-order method that is able to capture the wide range of scales in a robust and stable way.
However, the construction of entropy stable high-order 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 non-conservative 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 non-conservative term is necessary. We note that there is recent work on entropy stable DG methods applied to the ideal MHD equations by Rossmanith , Liu et al.  and Gallego-Valencia . In particular, we extend the ideas presented in Liu et al.  and provide clarification and motivation to the structure of the presented discretization of the non-conservative term and add detailed proofs of the 2D extension. As already mentioned, another important difference is the numerical control of the divergence-free constraint. To do so, we incorporate the modified GLM divergence cleaning approach in our full high-order DG entropy analysis. Furthermore, we build on the recent work of Gassner et al.  who showed that using a Bassi-Rebay (BR1) type scheme  is entropy stable for the viscous terms of the compressible Navier-Stokes 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 semi-definite.
The remainder of this paper is organized as follows: Our first main result, in Sec. 2, is the continuous entropy analysis of the three-dimensional resistive GLM-MHD equations, which demonstrates that the model indeed satisfies the entropy inequality and that the resistive terms can be recast into a symmetric and positive semi-definite 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 high-order 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 GLM-MHD equations. In general, we consider systems of conservation laws in a domain defined as
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
and the spatial gradient of a state as
The gradient of a spatial vector is a second order tensor, written in matrix form as
and the divergence of a flux written as a block vector is defined as
2.1 Resistive GLM-MHD 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.  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
with as well as the advective and viscous fluxes:
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 
and the heat flux is defined as
The introduced constants describe the viscosity from the Navier-Stokes equations, resistivity of the plasma, thermal conductivity and the universal gas constant, respectively. In particular, the constants and are first-order transport coefficients that describe the kinematic viscosity and the diffusivity of the magnetic field . We close the system with the ideal gas assumption, which relates the total energy and pressure
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 non-conservative 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 . In this work we restrict to the alternative Janhunen source term , for two reasons:
It is sufficient to restore the satisfaction of the second law of thermodynamics even when , e.g. , 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
We note that because on the continuous level the modification of to the resistive GLM-MHD system is consistent.
2.2 Thermodynamic properties of the system
To discuss the thermodynamic properties of the resistive GLM-MHD 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 non-conservative term proportional to (1.1). The ideal GLM-MHD 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 well-developed 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 GLM-MHD equations a suitable entropy function is the entropy density divided by the constant for convenience as
where is the physical entropy . From the entropy function we define the entropy variables to be
and we introduce , which is proportional to the inverse temperature.
It is known, that for smooth solutions when we contract the ideal GLM-MHD 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 
where we define the entropy fluxes to be
Additionally, as it will be necessary in later derivations and the proof of discrete entropy stability, we compute the entropy flux potential to be
where we introduce notation for the non-conservative state vector (2.11) contracted into entropy space
Furthermore, e.g. in case of shock discontinuities, the solution satisfies the following entropy inequality
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 GLM-MHD 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 GLM-MHD equations in (2.7) can be expressed by gradients of the entropy variables as
with a block matrix that is symmetric and positive semi-definite matrix, i.e,
First, we consider the viscous fluxes of the resistive GLM-MHD system in (2.7) separately
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:
We collect all these block matrices into the matrix
which clearly yields
For clarification, we present the first matrix
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
To show that the matrix is positive semi-definite is more involved. We first note that it is possible to split the matrix (2.26) into the viscous terms associated with the Navier-Stokes equations and the resistive terms of the magnetic fields that arise in the resistive GLM-MHD equations. We exploit this fact and rewrite the total diffusion matrix into two pieces
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 Navier-Stokes part, , is known to be positive semi-definite 
Thus, all that remains is to demonstrate that the additional resistive dissipation matrix, , is positive semi-definite. To do so we examine the eigenvalues of the system. We use the computer algebra system Maxima  to find an explicit expression of the eigenvalues to be
Under the physical assumptions that and we see that the eigenvalues (2.32) of the matrix are all non-negative. Hence, the block matrix is symmetric and positive semi-definite. ∎
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 GLM-MHD equations).
We start by contracting the resistive GLM-MHD system (2.6) with the entropy variables and integrate over the domain:
Here, denotes the inner product, e.g.,
From the definition of the entropy variables we have
Next, for clarity, we separate the advective flux into Euler, ideal MHD and GLM parts
The Euler terms generate the divergence of the entropy flux, e.g., 
and the GLM terms vanish directly with the modifications introduced in 
The damping source term for the GLM divergence cleaning is zero in all but its ninth component, so we see
Thus, we have
Next, to demonstrate the entropy stability for the resistive GLM-MHD equations, we address the viscous flux components. So, we integrate by parts for the remaining inner product including the viscous fluxes and obtain
where denotes the outward pointing normal vector at the domain boundaries. Assuming these to be periodic and using Lemma 1, we arrive at
Since is symmetric and positive semi-definite and , all inner products on the right hand side are positive and thus
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 high-order DG approximation. As such, the entropy flux potential is split into Euler, ideal MHD and GLM components
In summary, we have demonstrated that the resistive GLM-MHD 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 three-dimensional 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 GLM-MHD equations.
3 Split form discontinous Galerkin approximation
The goal of this paper is to construct an entropy stable high-order DG method for the resistive GLM-MHD 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 Legendre-Gauss-Lobatto nodes, the discrete derivative matrix and the discrete mass matrix satisfy the summation-by-parts (SBP) property for any polynomial order . This is a key property as it allows us to use results from the work of Fisher et al.  and Fisher and Carpenter . We follow the notation introduced in  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 high-order SBP method, i.e. to the higher order DG method.
In the following we introduce the so-called 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 two-dimensional Cartesian mesh with uniform squared elements. So, the first step in the discretization is to subdivide the computational domain into non-overlapping square elements and map each of these elements to the reference element . The bilinear mappings are defined as
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
for equal element side lengths , .
Next, on every element, each conserved variable is approximated by a polynomial of degree in each spatial direction on