Constraint Preserving Discontinuous Galerkin Method for Ideal Compressible MHD on 2-D Cartesian Grids

Abstract

We propose a constraint preserving discontinuous Galerkin method for ideal compressible MHD in two dimensions and using Cartesian grids, which automatically maintains the global divergence-free property. The approximation of the magnetic field is achieved using Raviart–Thomas polynomials and the DG scheme is based on evolving certain moments of these polynomials which automatically guarantees divergence-free property. We also develop HLL-type multi-dimensional Riemann solvers to estimate the electric field at vertices which are consistent with the 1-D Riemann solvers. When limiters are used, the divergence-free property may be lost and it is recovered by a divergence-free reconstruction step. We show the performance of the method on a range of test cases up to fourth order of accuracy.

Appendices

A Limiting and Divergence-Free Reconstruction

When the solution on the faces $b_x$, $b_y$ is limited as explained in Sect. 7, we lose the divergence-free property of the magnetic field. To recover this property, we have to perform a divergence-free reconstruction step. We explain this reconstruction process for second, third and fourth order accuracy. The fifth order version is given in [30] together with more details on the reconstruction idea. The resulting polynomial has the structure of the BDM polynomial on rectangles, see [12], Equation (3.29). Here we explain how the RT polynomial can be modified to recover divergence-free property. In two dimensions, the BDM polynomial has $(k+1)(k+2)+2$ degrees of freedom while its divergence has ${\frac{1}{2}}k(k+1)$ coefficients. Up to third order accuracy, the reconstruction can be performed using only the solution on the faces $(b_x, b_y)$, but at fourth order and higher, we need additional information which is supplied in the form of the curl of the magnetic field. This additional information is available to us via the cell moments $\alpha , \beta$.

For example, at fourth order, the BDM polynomial has $(3+1)(3+2)+2 =22$ degrees of freedom, while the face solution $(b_x,b_y)$, which are polynomials of degree 3, provide $(4+4+4+4)-1=15$ degrees of freedom, where one piece of information is redundant since the face solution satisfies $\int _{\partial K} \varvec{B}\cdot \varvec{n}\text{ d }s = 0$ on each cell K. The divergence-free condition on the BDM polynomial yields ${\frac{1}{2}}(3)(3+1)=6$ conditions. So we have a total of $15+6 = 21$ equations but 22 coefficients to be determined. Hence we need to supply one additional piece of information to completely determine the BDM polynomial.

We have the following inclusions $\mathbb {P}_k^2 \subset \text {BDM}(k) \subset \text {RT}(k)$ and the RT polynomial has many more basis functions than the BDM polynomial. In the reconstruction process, we set some of the coefficients $\{a,b\}$ in the RT polynomial to zero but this does not affect the accuracy since only coefficients $a_{ij}, b_{ij}$ with $i+j > k$ are set to zero, and we retain the $\mathbb {P}_k$ part of the solution.

1.1 A.1 Degree $k=1$

The divergence of the vector field $\varvec{B}\in \text {RT}(1)$ is given by

$$\begin{aligned} \begin{aligned} \nabla \cdot \varvec{B}&= \frac{a_{10}}{\Delta x} + \frac{b_{01}}{\Delta y} + \left( \frac{2 a_{20}}{\Delta x} + \frac{b_{11}}{\Delta y} \right) \phi _1(\xi ) + \left( \frac{a_{11}}{\Delta x} + \frac{2 b_{02}}{\Delta y} \right) \phi _1(\eta ) \\&\quad + 2 \left( \frac{a_{21}}{\Delta x} + \frac{b_{12}}{\Delta y} \right) \phi _1(\xi ) \phi _1(\eta ) \end{aligned} \end{aligned}$$

The coefficients $a_{ij}, b_{ij}$ are related to the face solution and cell moments according to Table 1. The constant term is already zero. The linear terms can be made zero by setting

$$\begin{aligned} \alpha _{00}= & {} {\frac{1}{2}}(a_0^- + a_0^+) + \frac{1}{12} (b_1^+ - b_1^-) \frac{\Delta x}{\Delta y}\\ \beta _{00}= & {} {\frac{1}{2}}(b_0^- + b_0^+) + \frac{1}{12} (a_1^+ - a_1^-) \frac{\Delta y}{\Delta x} \end{aligned}$$

This will however destroy the conservation property since $\alpha _{00}$, $\beta _{00}$ are cell averages of $B_x$, $B_y$ respectively. The bilinear term can be made zero by individually setting $a_{21} = b_{12} = 0$ which yields

$$\begin{aligned} \alpha _{01} = {\frac{1}{2}}(a_1^- + a_1^+), \qquad \beta _{10} = {\frac{1}{2}}(b_1^- + b_1^+) \end{aligned}$$

By this process we would have modified all the cell moments and the resulting reconstruction coincides with that of Balsara.

1.2 A.2 Degree $k=2$

The divergence of the vector field $\varvec{B}\in \text {RT}(2)$ is given by

$$\begin{aligned} \nabla \cdot \varvec{B}&= \left[ \frac{1}{\Delta x}\left( a_{10}+\frac{a_{30}}{10}\right) + \frac{1}{\Delta y}\left( b_{01} + \frac{b_{03}}{10}\right) \right] + \left[ \frac{2}{\Delta x}a_{20} +\frac{1}{\Delta y}\left( b_{11}+\frac{b_{13}}{10}\right) \right] \phi _1(\xi ) \\&+ \left[ \frac{1}{\Delta x}\left( a_{11}+\frac{a_{31}}{10}\right) +\frac{2}{\Delta y}b_{02}\right] \phi _1(\eta ) + \left[ \frac{3}{\Delta x}a_{30}+\frac{1}{\Delta y}\left( b_{21}+\frac{b_{23}}{10}\right) \right] \phi _2(\xi ) \\&+ 2\left[ \frac{a_{21}}{\Delta x}+\frac{b_{12}}{\Delta y}\right] \phi _1(\xi ) \phi _1(\eta ) + \left[ \frac{1}{\Delta x}\left( a_{12}+\frac{a_{32}}{10}\right) +\frac{3}{\Delta y}b_{03}\right] \phi _2(\eta ) \\&+ \left[ \frac{2}{\Delta x} a_{22} +\frac{3}{\Delta y}b_{13}\right] \phi _1(\xi )\phi _2(\eta )+\left[ \frac{3}{\Delta x}a_{31}+\frac{2}{\Delta y}b_{22} \right] \phi _2(\xi )\phi _1(\eta )\\&+ 3\left[ \frac{a_{32}}{\Delta x} + \frac{b_{23}}{\Delta y}\right] \phi _2(\xi )\phi _2(\eta ) \end{aligned}$$

The constant term is already zero. The linear terms can be made zero by setting

$$\begin{aligned} \alpha _{00} = {\frac{1}{2}}(a_0^- + a_0^+) + \frac{1}{12} (b_1^+ - b_1^-) \frac{\Delta x}{\Delta y}, \qquad \beta _{00} = {\frac{1}{2}}(b_0^- + b_0^+) + \frac{1}{12} (a_1^+ - a_1^-) \frac{\Delta y}{\Delta x} \end{aligned}$$

The quadratic terms are zero by choosing

$$\begin{aligned} \alpha _{10} = a_0^+ - a_0^- + \frac{1}{30}(b_2^+ - b_2^-) \frac{\Delta x}{\Delta y}, \qquad \beta _{01} = b_0^+ - b_0^- + \frac{1}{30}(a_2^+ - a_2^-) \frac{\Delta y}{\Delta x} \end{aligned}$$

and also setting $a_{21} = b_{12} = 0$ which yields

$$\begin{aligned} \alpha _{01} = {\tfrac{1}{2}}(a_1^- + a_1^+), \qquad \beta _{10} = {\tfrac{1}{2}}(b_1^- + b_1^+) \end{aligned}$$

In the cubic terms, we set each coefficient to zero, $a_{22} = a_{31} = a_{32} = b_{22} = b_{13} = b_{23} = 0$, which yields

$$\begin{aligned} \begin{array}{ll} \alpha _{02} = {\tfrac{1}{2}}(a_2^- + a_2^+) &{} \qquad \qquad \beta _{20} = {\tfrac{1}{2}}(b_2^- + b_2^+) \\ \alpha _{11} = a_1^+ - a_1^- &{} \qquad \qquad \beta _{11} = b_1^+ - b_1^- \\ \end{array} \end{aligned}$$

Finally, in the biquadratic term, we set $a_{32} = b_{23} = 0$ to obtain

$$\begin{aligned} \alpha _{12} = a_2^+ - a_2^-, \qquad \beta _{21} = b_2^+ - b_2^- \end{aligned}$$

1.3 A.3 Degree $k=3$

The divergence of the vector field $\varvec{B}\in \text {RT}(3)$ is given by

$$\begin{aligned} \nabla \cdot \varvec{B}=&\left[ \frac{1}{\Delta x}\left( a_{10}+\frac{a_{30}}{10}\right) + \frac{1}{\Delta y}\left( b_{01} + \frac{b_{03}}{10}\right) \right] \\ +&\left[ \frac{1}{\Delta x}\left( 2a_{20}+\frac{6}{35}a_{40}\right) +\frac{1}{\Delta y}\left( b_{11}+\frac{b_{13}}{10}\right) \right] \phi _1(\xi ) \\ +&\left[ \frac{1}{\Delta x}\left( a_{11}+\frac{a_{31}}{10}\right) +\frac{1}{\Delta y}\left( 2b_{02}+\frac{6}{35}b_{04}\right) \right] \phi _1(\eta )\\ +&\left[ \frac{3}{\Delta x}a_{30}+\frac{1}{\Delta y}\left( b_{21}+\frac{b_{23}}{10}\right) \right] \phi _2(\xi )+ \left[ \frac{1}{\Delta x}\left( a_{12}+\frac{a_{32}}{10}\right) +\frac{3}{\Delta y}b_{03}\right] \phi _2(\eta ) \\ +&\left[ \frac{1}{\Delta x}\left( 2a_{21}+\frac{6}{35}a_{41}\right) +\frac{1}{\Delta y}\left( 2b_{12}+\frac{6}{35}b_{14}\right) \right] \phi _1(\xi ) \phi _1(\eta ) \\ +&\left[ \frac{4}{\Delta x}a_{40} + \frac{1}{\Delta y}\left( b_{31}+\frac{b_{33}}{10}\right) \right] \phi _3(\xi )+\left[ \frac{1}{\Delta x}\left( 2 a_{22}+\frac{6}{35}a_{42}\right) +\frac{3}{\Delta y}b_{13}\right] \phi _1(\xi )\phi _2(\eta )\\ +&\left[ \frac{3}{\Delta x}a_{31}+\frac{1}{\Delta y}\left( 2b_{22}+\frac{6}{35}b_{24}\right) \right] \phi _2(\xi )\phi _1(\eta ) + \left[ \frac{1}{\Delta x}\left( a_{13}+\frac{a_{33}}{10}\right) +\frac{4}{\Delta y} b_{04}\right] \phi _3(\eta )\\ +\,&3\left[ \frac{a_{32}}{\Delta x} + \frac{b_{23}}{\Delta y}\right] \phi _2(\xi )\phi _2(\eta ) + \left[ \frac{1}{\Delta x}\left( 2 a_{23}+\frac{6}{35}a_{43}\right) +\frac{4}{\Delta y} b_{14}\right] \phi _1(\xi ) \phi _3(\eta ) \\ +&\left[ \frac{4}{\Delta x}a_{41} +\frac{1}{\Delta y}\left( 2b_{32} +\frac{6}{35}b_{34} \right) \right] \phi _3(\xi )\phi _1(\eta ) \\ +&\left[ \frac{3}{\Delta x} a_{33} +\frac{4}{\Delta y}b_{24}\right] \phi _2(\xi ) \phi _3(\eta ) +\left[ \frac{4}{\Delta x}a_{42} + \frac{3}{\Delta y}b_{33}\right] \phi _3(\xi ) \phi _2(\eta ) \\ +&\left[ \frac{4}{\Delta x} a_{43}+\frac{4}{\Delta y}b_{34} \right] \phi _3(\xi ) \phi _3(\eta ) \end{aligned}$$

The constant term is already zero. The linear terms can be made zero by setting

$$\begin{aligned} \alpha _{00} = {\frac{1}{2}}(a_0^- + a_0^+) + \frac{1}{12} (b_1^+ - b_1^-) \frac{\Delta x}{\Delta y}, \qquad \beta _{00} = {\frac{1}{2}}(b_0^- + b_0^+) + \frac{1}{12} (a_1^+ - a_1^-) \frac{\Delta y}{\Delta x} \end{aligned}$$

The quadratic terms which are coefficients of $\phi _2(\xi )$, $\phi _2(\eta )$ become zero by choosing

$$\begin{aligned} \alpha _{10} = a_0^+ - a_0^- + \frac{1}{30}(b_2^+ - b_2^-) \frac{\Delta x}{\Delta y}, \qquad \beta _{01} = b_0^+ - b_0^- + \frac{1}{30}(a_2^+ - a_2^-) \frac{\Delta y}{\Delta x} \end{aligned}$$

The coefficient of $\phi _1(\xi )\phi _1(\eta )$ gives only one equation but there are two unknowns; adding an extra equation $\omega = \beta _{10}-\alpha _{01}$, we can solve for the two coefficients

$$\begin{aligned} \alpha _{01} = \frac{1}{\left( 1+\frac{\Delta y}{\Delta x}\right) }\left( r_2-\omega +r_1\frac{\Delta y}{\Delta x}\right) , \qquad \beta _{10} = \omega +\alpha _{01} \end{aligned}$$

where $r_1={\frac{1}{2}}(a_1^{-}+a_1^{+})$ and $r_2={\frac{1}{2}}(b_1^{-}+b_1^{+})$. The cubic terms are zeros by choosing

$$\begin{aligned}&\alpha _{20} = -{\frac{1}{2}}(b_1^{+}-b_1^{-})\frac{\Delta x}{\Delta y}+\frac{3}{140}(b_3^{+}-b_3^{-})\frac{\Delta x}{\Delta y},\\&\beta _{02}=-\frac{1}{2}(a_1^{+}-a_1^{-})\frac{\Delta y}{\Delta x}+\frac{3}{140}(a_3^{+}-a_3^{-})\frac{\Delta y}{\Delta x} \end{aligned}$$

and setting $2 a_{22}+\frac{6}{35}a_{42} = 2b_{22}+\frac{6}{35}b_{24} = a_{31} = b_{13} = 0$, which yields

$$\begin{aligned} \begin{array}{ll} \alpha _{02} = {\tfrac{1}{2}}(a_2^- + a_2^+) &{} \qquad \qquad \beta _{20} = {\tfrac{1}{2}}(b_2^- + b_2^+) \\ \alpha _{11} = a_1^+ - a_1^- &{} \qquad \qquad \beta _{11} = b_1^+ - b_1^- \\ \end{array} \end{aligned}$$

In the higher order term greater then three, we set each coefficient to zero, $a_{32}=b_{23} = 2a_{23}+\frac{6}{35}a_{43} = b_{14}=a_{41} = 2b_{32}+\frac{6}{35}b_{34} = a_{33}=b_{24}=a_{42}=b_{33}=a_{43}=b_{34}=0$, which yields

$$\begin{aligned} \begin{array}{ll} \alpha _{12} = a_2^+ - a_2^- &{} \qquad \qquad \beta _{21} = b_2^+ - b_2^-\\ \alpha _{03} = {\tfrac{1}{2}}(a_3^+ +a_3^-) &{} \qquad \qquad \beta _{12} = 6(r_2-\beta _{10})\\ \alpha _{21} = 6(r_1-\alpha _{01}) &{}\qquad \qquad \beta _{30} = {\tfrac{1}{2}}(b_3^+ +b_3^-)\\ \alpha _{13} = \alpha _3^+ -\alpha _3^- &{}\qquad \qquad \beta _{22} = 0\\ \alpha _{22} = 0 &{}\qquad \qquad \beta _{31} = b_3^+-b_3^-\\ \alpha _{23} = 0 &{} \qquad \qquad \beta _{32} = 0 \end{array} \end{aligned}$$

We see that at fourth order, we need an extra information which we took in the form of the quantity $\omega$ in order to complete the divergence-free reconstruction. Note $\omega$ gives us information about the curl of the magnetic field.

B Setting the Initial Condition

Let $\psi _h \in \mathbb {Q}_{k+1,k+1}$ be a continuous interpolation of the magnetic potential $\psi$ which can be achieved using $(k+2)\times (k+2)$ GLL nodes. Then we can set the magnetic field as $B_x = \frac{\partial \psi _h}{\partial y}$, $B_y = -\frac{\partial \psi _h}{\partial x}$ which will be exactly divergence-free. But in our work, we want to set the initial condition in terms of the polynomials $b_x,b_y$ and the moments $\alpha ,\beta$. We can perform an $L^2$ projection of $\nabla \times (\psi _h e_z)$ to initialize $b_x,b_y$ which will be exact, and the moments can be computed using the same GLL nodes for quadrature as are used to define $\psi _h$. Let $\xi _i, i=1,2,\ldots ,k+2$ denote the GLL nodes and let $\ell _i(\xi ), i=1,2,\ldots ,k+2$ be the Lagrange polynomials. Define the barycentric weights

$$\begin{aligned} w_j = \frac{1}{\prod _{i=1, i\ne j}^{k+2} (\xi _j - \xi _i)} \end{aligned}$$

Then the derivatives of Lagrange polynomials at the GLL nodes are given by

$$\begin{aligned} D_{ij} = \ell _j'(\xi _i) = \frac{w_j}{w_i} \frac{1}{\xi _i - \xi _j}, \quad i \ne j, \qquad D_{ii} = \ell _i'(\xi _i) = -\sum _{j=1, j \ne i}^{k+2} D_{ij} \end{aligned}$$

The derivatives of the potential at the GLL nodes are given by

$$\begin{aligned} \frac{\partial \psi _h}{\partial x}(\xi _i,\xi _j) = \frac{1}{\Delta x} [D \cdot \psi (:,j)]_i, \qquad \frac{\partial \psi _h}{\partial y}(\xi _i,\xi _j) = \frac{1}{\Delta y} [D \cdot \psi (i,:)]_j \end{aligned}$$

The cell moments are initialized as $\alpha = \alpha ^h(\partial _y \psi _h)$ and $\beta = \beta ^h(-\partial _x\psi _h)$ where the superscript h denotes that we compute the integrals using $(k+2)^2$-point GLL quadrature which is exact for the integrands involved in the cell moments.