Structure-Preserving Oscillation-Eliminating Discontinuous Galerkin Schemes for Ideal MHD Equations: Locally Divergence-Free and Positivity-Preserving
Mengqing Liu Kailiang Wu


[figure]labelfont=bf,labelformat=default,labelsep=period,name=Fig.
[table]labelsep=newline, singlelinecheck=false




mymainaddress]Mengqing Liu
liumq@sustech.edu.cn

mymainaddress,mysecondaryaddress,mythirdaddress]Kailiang Wu
mycorrespondingauthor
[mycorrespondingauthor]Corresponding author
wukl@sustech.edu.cn


[mymainaddress]Department of Mathematics and Shenzhen International Center for Mathematics, Southern University of Science and Technology, Shenzhen 518055, China
[mysecondaryaddress]National Center for Applied Mathematics Shenzhen (NCAMS), Shenzhen 518055, China
[mythirdaddress]Guangdong Provincial Key Laboratory of Computational Science and Material Design, Shenzhen 518055, China





   
   This paper develops structure-preserving, oscillation-eliminating discontinuous Galerkin (OEDG) schemes for ideal magnetohydrodynamics (MHD), as a sequel to our recent work <cit.>. The schemes are based on a locally divergence-free (LDF) oscillation-eliminating (OE) procedure to suppress spurious oscillations while maintaining many of the good properties of original DG schemes, such as conservation, local compactness, and optimal convergence rates. The OE procedure is built on a simple damping equation—a linear ordinary differential equation (ODE) whose solution can be exactly formulated. Because this OE procedure does not interfere with DG spatial discretization and RK stage update, it can be easily incorporated to existing DG codes as an independent module. These features make the proposed LDF OEDG schemes highly efficient and easy to implement. 
   In addition, we present a positivity-preserving (PP) analysis of the LDF OEDG schemes on Cartesian meshes via the optimal convex decomposition technique and the geometric quasi-linearization (GQL) approach. Efficient PP LDF OEDG schemes are obtained with the HLL flux under a condition accessible by the simple local scaling PP limiter. 
   Several one- and two-dimensional MHD tests confirm the accuracy, effectiveness, and robustness of the proposed structure-preserving OEDG schemes.



magnetohydrodynamics (MHD)oscillation-eliminating discontinuous Galerkin (OEDG) locally divergence-free (LDF)positivity-preserving (PP) multidimensional







§ INTRODUCTION


The magnetohydrodynamic (MHD) equations are widely utilized for modeling macro-scale behaviors in various domains, including laboratory, space, and astrophysical plasmas. Notable applications of these equations can be found in studies of phenomena such as coronal mass ejections <cit.>, the evolution of solar coronal magnetic fields <cit.>, and changes within Earth's magnetosphere <cit.>.
As a nonlinear hyperbolic system, the ideal compressible MHD equations pose significant challenges in obtaining analytical solutions due to their complexity. Consequently, numerical methods have become the principal approach for deriving solutions in this field. Among these, discontinuous Galerkin (DG) methods are particularly favored for their capacity to achieve high-order accuracy, excellent local compactness, and compatibility with arbitrary mesh configurations (see, e.g., <cit.> and references therein). These attributes have led to their widespread adoption in MHD numerical simulations, as evidenced by various studies (e.g., <cit.>). However, when addressing MHD problems featuring strong discontinuities, such as shocks, conventional high-order DG methods can produce spurious oscillations near these discontinuities, potentially leading to non-physical outcomes and simulation failures.

To suppress spurious oscillations in the DG methods for hyperbolic conservation laws, various techniques have been developed, including the total variation bounded limiter <cit.>, the weighted essentially non-oscillatory limiter <cit.>, and artificial diffusion <cit.>. These strategies effectively control oscillations and enhance the stability of DG schemes. However, they may compromise certain desirable attributes of the original DG schemes (e.g., optimal error estimates or superconvergence) or depend on some problem-specific parameters.
Recently, Lu, Liu, and Shu <cit.> proposed the so-called oscillation-free DG (OFDG) method for scalar conservation laws, which skillfully integrates a damping term into the semi-discrete DG schemes. The OFDG method can mitigate numerical oscillations near strong discontinuities while preserving many advantageous properties of the conventional DG schemes, including conservation, superconvergence, and optimal error estimates.
Subsequently, the OFDG method was successfully extended to hyperbolic systems <cit.> with characteristic decomposition and applied to shallow water equations <cit.>, nonlinear degenerate parabolic equations <cit.>, and multi-component chemically reacting flows <cit.>. However, the characteristic decomposition in the OFDG method may be computationally expensive and require careful adjustments to the eigenvectors due to the lack of scale invariance. Furthermore, as Liu, Lu, and Shu have mentioned <cit.>, the stiffness of the large damping terms near strong discontinuities makes the Courant–Friedrichs–Lewy (CFL) condition more restrictive, necessitating the use of modified exponential Runge–Kutta methods to alleviate this restriction. 
More recently, motivated by <cit.>, a novel oscillation-eliminating DG (OEDG) method was systematically proposed in <cit.> by alternately evolving the conventional semi-discrete DG scheme and a new damping equation. Compared to the OFDG method, the OEDG method decouples the damping terms from the DG formulations, thereby being non-intrusive for integration into existing DG codes and simpler and more efficient for implementation. Moreover, the OEDG method is free from characteristic decomposition for hyperbolic systems, is stable under the normal CFL condition without requiring the exponential Runge–Kutta discretization, and exhibits desirable scale-invariant and evolution-invariant properties for problems spanning various scales and wave speeds.





In ideal MHD, the magnetic field is divergence-free, and both density and internal energy (or pressure) are positive.
Maintaining these conditions—zero divergence of the magnetic field and positivity of density and pressure—poses significant challenges in MHD simulations. Large divergence errors can induce nonphysical features or lead to numerical instabilities <cit.>.
To reduce the divergence errors, various techniques have been developed, including the projection method <cit.>, constrained transport method <cit.>, eight-wave method with the Godunov–Powell source term <cit.>, hyperbolic divergence cleaning method <cit.>, locally divergence-free (LDF) method <cit.>, and globally divergence-free (GDF) method <cit.>.
These approaches have been extensively applied across numerous schemes, thanks to their efficacy in controlling divergence. Positivity of density and pressure can be compromised in challenging MHD simulations involving conditions such as low density or pressure, high Mach numbers, large magnetic energies, or strong discontinuities. To address this, early research primarily focused on developing  one-dimensional robust approximate Riemann solvers; see, e.g., <cit.>. 
Waagan <cit.> introduced a second-order MUSCL-Hancock scheme using relaxation Riemann solvers <cit.> that ensures the positivity of density and pressure; the robustness of this scheme was further validated through extensive testing <cit.>. There has been notable progress in the development of positivity-preserving (PP) limiters for higher-order schemes <cit.>. Balsara <cit.> proposed a self-adjusting PP limiter to maintain the positivity of reconstructed density and pressure in ideal MHD. Cheng et al. <cit.> extended the PP limiter from <cit.> and showed the PP property of one-dimensional (central) DG schemes under some assumptions for ideal MHD. Christlieb et al. <cit.> developed high-order PP finite difference methods for ideal MHD utilizing a parameterized flux limiter. These PP techniques have enhanced the robustness of MHD schemes, as evidenced by numerous numerical tests. However, the intrinsic complexity of MHD equations and the unclear relationship between the PP property and the divergence-free condition of the magnetic field have limited the proposal of numerical schemes that can rigorously prove to preserve positivity for MHD in earlier research.

In recent years, significant advances <cit.> have been made in revealing the theoretical connections between the PP property and the divergence-free condition for MHD systems.  The work <cit.> first rigorously proved that a discrete divergence-free (DDF) condition is intrinsically linked to achieving the PP property for the ideal MHD equations. As discovered in <cit.>, the PP property may be compromised even if there is a minor violation of the DDF condition. Furthermore, it has been found in  <cit.> that even the exact smooth solutions of the conservative MHD system with nonzero magnetic divergence may result in negative pressure. The symmetrizable MHD system (i.e., the MHD system with the Godunov–Powell source term <cit.>) does not exhibit this problem <cit.>, and a provably PP DG scheme for this system was first developed in <cit.>. Subsequently, <cit.> established a framework for provably positive high-order finite volume and DG schemes on general meshes. All these findings are based on an innovative technique known as geometric quasilinearization (GQL), which simplifies nonlinear constraints—such as the positivity of internal energy (or pressure) in the MHD system—into linear constraints by incorporating suitable auxiliary variables. Building on the geometric properties of convex sets, the general GQL framework was recently proposed in <cit.> to study bound-preserving schemes under nonlinear constraints. For the related developments in relativistic MHD, see <cit.>. Very recently, a second-order structure-preserving finite element method was proposed in <cit.> based on convex limiting and a novel splitting. 







In this paper, we develop the LDF OEDG method for the ideal MHD and carry out a PP analysis of the OEDG method with the HLL flux. The contributions, innovations, and significance of this work are outlined as follows:

  
  * For extending the OEDG method to MHD, how to make the OEDG method and the LDF condition compatible is something that must be considered.
        In order to control the oscillations without losing the LDF property of the computed magnetic field, we perform OE steps for Legendre polynomials and LDF polynomials respectively using the LDF DG method as the base scheme, and propose an LDF OEDG scheme.
        The LDF OEDG scheme can inherently balance LDF constraint and oscillation suppression, and is computationally cheap for arbitrary meshes.
  
  * In the framework that proves PP property, it is well known that the decomposition of the cell averages of the numerical solution into a convex combination of the solution values at certain quadrature points is a key ingredient.
       By employing the LDF OEDG scheme, PP limiter, HLL flux and properly discretized Godunov-Powell source term, we provide a provably PP analysis via general convex decomposition <cit.> on Cartesian meshes, and also give derived CFL conditions via some specific convex decompositions, including Zhang-Shu convex decomposition <cit.> and optimal convex decomposition <cit.>. 
  
  * In addition, different from the works in <cit.>, we only consider adding Godunov-Powell source term to the evolution equations of the cell averages, since the PP property of the scheme is only related to the discrete equations satisfied by the cell averages.

Some one- and two-dimensional (2D) MHD examples are utilized to test the accuracy, effectiveness, and robustness of the proposed scheme.
The numerical results show that our scheme can achieve expected convergence order, effectively suppress oscillations, and behaves robust in some tough problems.


The paper is organized as follows.

In section <ref>, we proposed a LDF OEDG method for 2D ideal MHD equations.
In section <ref>, we construct a PP LDF OEDG scheme for 2D ideal MHD equations and prove the PP property of the scheme.
In section <ref>, we show several numerical examples for MHD in one and two dimensions.
Section <ref> gives the concluding remarks for this paper.









§ LOCALLY DIVERGENCE-FREE OEDG METHOD FOR IDEAL MHD


This section presents the LDF OEDG method for the ideal MHD system. Without loss of generality, we focus on the 2D MHD equations, while our method is directly applicable to the one-dimensional case and can be readily extended to the three-dimensional case.





 §.§ Governing equations of 2D MHD system


The 2D conservative MHD system can be written as

    ∂U/∂ t+∂F_1(U)/∂x
        +∂F_2(U)/∂y= 0,

where
U=(ρ, ρu, B, E)^⊤=(ρ, m, B, E)^⊤
is the conservation vector, and the fluxes in the x- and y-directions are respectively given by

    F_1(U)=
      ([                      ρ u_1;  ρu_1^2+(p+B^2/2)-
      B_1^2;           ρ u_1 u_2-B_1B_2;           ρ u_1 u_3-B_1B_3;                          0;              u_1B_2-B_1u_2;              u_1B_3-B_1u_3; (E+p+B^2/2)u_1-
      B_1(B·u);                            ]),  F_2(U)=
      ([                      ρ u_2;           ρ u_2 u_1-B_2B_1;  ρu_2^2+(p+B^2/2)-
      B_2^2;           ρ u_2 u_3-B_2B_3;              u_2B_1-B_2u_1;                          0;              u_2B_3-B_2u_3; (E+p+B^2/2)u_2-
      B_2(B·u);                            ]).

Here, ρ represents the fluid density, m is the momentum vector, and B = (B_1, B_2, B_3) is the magnetic field vector. The variable p stands for thermal pressure, and E is the total energy density, expressed by the ideal state equation as

    E = p/γ-1 + 1/2ρm^2 + 1/2B^2,

where γ is the adiabatic index and · denotes the vector-2 norm.  Additionally, the physical solutions of the MHD equations require a divergence-free constraint on the magnetic field:

    ∇·B = 0,

which implies the absence of magnetic monopoles.

Divide the conservation variables U into two parts:

    R=(ρ, ρ u_1, ρ u_2, ρ u_3, B_3, E)^⊤  and  Q=(B_1, B_2)^⊤,

and correspondingly divide the flux Fi(U) into Fi^R(U) and Fi^Q(U). The MHD system (<ref>) is then split into

    ∂R/∂ t+∂F_1^R(U)/∂x
    	+∂F_2^R(U)/∂y= 0,


    ∂Q/∂ t+∂F_1^Q(U)/∂x
    	+∂F_2^Q(U)/∂y= 0.

The objective of the LDF OEDG method is to find an approximate solution U_h, which comprises R_h and Q_h  obtained through different OEDG discretizations.



Next, we will elaborate on our LDF OEDG method, including the LDF finite element spaces, the semi-discrete LDF DG formulation, and the fully discrete LDF OEDG discretization. It incorporates an LDF OE procedure after each Runge–Kutta stage to eliminate spurious oscillations. The PP property of the LDF OEDG method will be systematically addressed in section <ref>.







 §.§ LDF DG finite element spaces



Let Ω=⋃_ijI_ij with mesh element I_ij=[x_i-1/2, x_i+1/2]×[y_j-1/2, y_j+1/2]   be a partition of the 2D spatial domain Ω.
Denote by Δ x and Δ y the spatial step-sizes in the x- and y-directions, respectively.
We now introduce two finite element DG spaces in which we seek approximations to R and Q, respectively.








The approximation space for R
is defined as follows:

    𝕍_R^k:={v(x) ∈ L^2(Ω):.v(x)|_I_ij∈[ℙ^k(I_ij)]^6  ∀ i,j},

where x=(x, y) denotes the spatial coordinate vector, and ℙ^k(I_ij) is the space of polynomials of total degree less than or equal to k on the element I_ij.
An orthogonal basis of ℙ^k(I_ij) is denoted as {ϕ_ij^(α)(x): |α| ≤ k }, where α = (α_1, α_2) is the multi-index vector with |α| = α_1 + α_2.
For example, we adopt the following (orthogonal) Legendre basis over cell I_ij:

    1,   ξ,   η,   ξ^2-1/3,   ξη,   η^2-1/3,   ⋯

where ξ=2(x-x_i)/Δ x and η=2(y-y_j)/Δ y;
the corresponding multi-index vectors α are

    (0, 0),   (1, 0),   (0, 1),
      (2, 0),   (1, 1),   (0, 2),    ...,

The dimensionality of ℙ^k(I_ij) is D_R^k = (k+1)(k+2)/2. Using the Legendre basis, the approximate solution R_h(x, t)  on element I_ij can be expressed as

    R_h(x, t)=∑_μ=0^k∑_|
    α|=μR_ij^(α)(t)ϕ_ij^(α)(x)   for  x∈ I_ij,

where R_ij^(α)(t) represents a moment of the approximate solution R_h(x, t).


We seek the divergence-free approximation to Q
in the following LDF space <cit.>:

    𝕍_Q^k:={v(x)=
    (v_1(x), v_2(x))^⊤:
    v_1(x), v_2(x) ∈ℙ^k(I_ij),
    .∇·v(x)|_I_ij=0   ∀ i,j
     }.

The local orthogonal basis of 𝕍_Q^k on element I_ij is denoted as {φ_ij^(μ)(x),
μ=0, 1, ⋯, D_Q^k-1},
where D_Q^k=(k+1)(k+4)/2 is the dimensionality of the LDF polynomial space defined on I_ij.
We employ the following local orthogonal LDF basis

    ([ 0; 1;   ]),
      ([ 1; 0;   ]),
      ([ 0; ξ;   ]),
      ([ η; 0;   ]),
      ([  Δ xξ; -Δ yη;       ]),
    
    ([ η^2-1/3;       0;         ]),
      ([       0; ξ^2-1/3;         ]),
      ([ Δ x(ξ^2-1/3);      -2Δ yξη;              ]),
      ([      -2Δ xξη; Δ y(η^2-1/3);              ]),
      ⋯,

where ξ=2(x-x_i)/Δ x and η=2(y-y_j)/Δ y.
With this basis, we express the approximate solution Q_h(x, t) over element I_ij as

    Q_h(x, t)=∑_μ=0^k∑_η=D_Q^μ-1^D_Q^μ-1
    Q_ij^(η)(t)φ_ij^(
    η)(x)   for  x∈ I_ij.

where Q_ij^(η)(t) represents a moment of approximate polynomial solution Q_h(x, t).
Note that the approximate magnetic field Q_h given in the above form naturally satisfies the divergence-free constraint (<ref>) within each element I_ij.






 §.§ Semi-discrete LDF DG formulation


The LDF DG method uses different spatial discretizations for R and Q.

For any test function v(x) ∈ℙ^k(I_ij), the approximate solution R_h(x, t) is expected to satisfy

    d/d t∫_I_ijR_h(x, t) v(x) dx
    =∫_I_ijF_1^R(U_h(x, t)) ∂ v(x)/∂ x dx+∫_I_ijF_2^R(U_h(x, t)) ∂ v(x)/∂ y dx
       -∫_y_j-1/2^y_j+1/2(
    v(x_i+1/2^-, y)F_1^R_, i+1/2(y, t)
    -v(x_i-1/2^+, y)F_1^R_, i-1/2(y, t)
    ) dy 
       -∫_x_i-1/2^x_i+1/2(
    v(x, y_j+1/2^-)F_2^R_, j+1/2(x, t)
    -v(x, y_j-1/2^+)F_2^R_, j-1/2(x, t)
    ) dx,

where the superscripts - and + indicate the left- and right-hand-side limits at the cell interfaces, and



    F_1^R_, i+1/2(y, t):=
      F_1^R(
      U_h(x_i+1/2^-, y, t),
      U_h(x_i+1/2^+, y, t)),


    F_2^R_, j+1/2(x, t):=
      F_2^R(
      U_h(x, y_j+1/2^-, t),
      U_h(x, y_j+1/2^+, t)).

Here, F_1^R and F_2^R denote the numerical fluxes in the x- and y-directions, respectively.
In this work, we use the PP HLL fluxes

    F_ℓ^R(U^-,U^+)=
    𝒱_ℓ^+F_ℓ^R(U^-)-
    𝒱_ℓ^-F_ℓ^R(U^+)+
    𝒱_ℓ^-𝒱_ℓ^+(R^+-R^-)/𝒱_ℓ^+-𝒱_ℓ^-,   ℓ=1, 2,

where 𝒱_ℓ^- and 𝒱_ℓ^+ denote the (properly estimated) non-positive minimum and non-negative maximum wave speeds, respectively,
with ℓ=1 for the x-direction and ℓ=2 for the y-direction.
The estimates of 𝒱_ℓ^- and 𝒱_ℓ^+ are crucial for the PP property; see section <ref> for details. 
Taking the test function v(x) as the (orthogonal) Legendre basis function ϕ_ij^(α)(x) and using the Gauss quadrature of sufficiently high accuracy to approximate the right-hand side of (<ref>), we obtain the following semi-discrete scheme

    (∫_I_ij(ϕ_ij^(α)(x)
    )^2dx)
    dR_ij^(α)(t)/d t
    =Δ x Δ y∑_μ=1^q∑_η=1^qω_μ^Gω_η^GF_1^R(U_h(x_i^(μ), y_j^(η), t)) ∂_xϕ_ij^(α)(x_i^(μ), y_j^(η))
       +Δ x Δ y∑_μ=1^q∑_η=1^qω_μ^Gω_η^GF_2^R(U_h(x_i^(μ), y_j^(η), t)) ∂_yϕ_ij^(α)(x_i^(μ), y_j^(η))
    
       -Δ y ∑_μ=1^qω_μ^G(
    ϕ_ij^(α)(x_i+1/2^-, y_j^(μ))F_1^R_, i+1/2(y_j^(μ), t)
    -ϕ_ij^(α)(x_i-1/2^+, y_j^(μ))F_1^R_, i-1/2(y_j^(μ), t)
    )
       -Δ x∑_μ=1^qω_μ^G(
    ϕ_ij^(α)(x_i^(μ), y_j+1/2^-)F_2^R_, j+1/2(x_i^(μ), t)
    -ϕ_ij^(α)(x_i^(μ), y_j-1/2^+)F_2^R_, j-1/2(x_i^(μ), t)
    ),   ∀
    |α| ≤ k,

where {x_i^(μ)}_μ=1^q and {y_j^(μ)}_μ=1^q
are the q-point Gauss quadrature nodes with 2q-1≥2k in [x_i-1/2, x_i+1/2] and [y_j-1/2, y_j+1/2], respectively, and {ω_μ^G}_μ=1^q are the associated weights.





For any test function v(x) ∈𝕍_Q^k, the approximate solution Q_h(x, t) is expected to satisfy

    d/d t∫_I_ijQ_h(x, t) ·v(x) dx
    =∫_I_ijF_1^Q(U_h(x, t)) ·∂v(x)/∂ x dx+∫_I_ijF_2^Q(U_h(x, t)) ·∂v(x)/∂ y dx
       -∫_y_j-1/2^y_j+1/2(
    v(x_i+1/2^-, y)·F_1^Q_, i+1/2(y, t)
    -v(x_i-1/2^+, y)·F_1^Q_, i-1/2(y, t)
    ) dy 
       -∫_x_i-1/2^x_i+1/2(
    v(x, y_j+1/2^-)·F_2^Q_, j+1/2(x, t)
    -v(x, y_j-1/2^+)·F_2^Q_, j-1/2(x, t)
    ) dx

where

    F_1^Q_, i+1/2(y, t):=
      F_1^Q(
      U_h(x_i+1/2^-, y, t),
      U_h(x_i+1/2^+, y, t)),


    F_2^Q_, j+1/2(x, t):=
      F_2^Q(
      U_h(x, y_j+1/2^-, t),
      U_h(x, y_j+1/2^+, t))

with F_1^Q and F_2^Q denoting the numerical fluxes in the x- and y-directions, respectively.
Here, we use the PP HLL fluxes

    F_ℓ^Q(U^-,U^+)=
    𝒱_ℓ^+F_ℓ^Q(U^-)-
    𝒱_ℓ^-F_ℓ^Q(U^+)+
    𝒱_ℓ^-𝒱_ℓ^+(Q^+-Q^-)/𝒱_ℓ^+-𝒱_ℓ^-,   ℓ=1, 2,

where 𝒱_ℓ^± are the same as those in equation (<ref>). 
Taking the test function v(x) as the orthogonal LDF basis function φ_ij^(m)(x) and using the same Gauss quadrature to approximate the right-hand side of  (<ref>), we can obtain the following semi-discrete scheme

    (∫_I_ij(φ_ij^(m)(x)
    )^2dx)
    dQ_ij^(m)(t)/d t
    =Δ x Δ y∑_μ=1^q∑_η=1^qω_μ^Gω_η^GF_1^Q(U_h(x_i^(μ), y_j^(η), t))·∂_xφ_ij^(m)(x_i^(μ), y_j^(η))
       +Δ x Δ y∑_μ=1^q∑_η=1^qω_μ^Gω_η^GF_2^Q(U_h(x_i^(μ), y_j^(η), t))·∂_yφ_ij^(m)(x_i^(μ), y_j^(η))
    
       -Δ y ∑_μ=1^qω_μ^G(
    φ_ij^(m)(x_i+1/2^-, y_j^(μ))·F_1^Q_, i+1/2(y_j^(μ), t)
    -φ_ij^(m)(x_i-1/2^+, y_j^(μ))·F_1^Q_, i-1/2(y_j^(μ), t)
    )
       -Δ x∑_μ=1^qω_μ^G(
    φ_ij^(m)(x_i^(μ), y_j+1/2^-)·F_2^Q_, j+1/2(x_i^(μ), t)
    -φ_ij^(m)(x_i^(μ), y_j-1/2^+)·F_2^Q_, j-1/2(x_i^(μ), t)
    ), 
                                                           
    m=0, 1, ⋯, D_Q^k-1,

where the Gauss quadrature nodes and weights are the same as those in (<ref>).






 §.§ Fully discrete OEDG schemes



The combination of semi-discrete schemes (<ref>) and (<ref>) can be regarded as
an ordinary differential equation (ODE)
with respect to time t.
For brevity, we rewrite equations (<ref>) and (<ref>) in the following ODE form:

    d/d tU_h=L(U_h),

where U_h denotes the combination of R_h and Q_h.
This ODE system (<ref>) is typically discretized in time by a Runge-Kutta method, for example, the explicit third-order strong-stability-preserving (SSP) Runge-Kutta method:

    {U_h^n, 0   = U_h^n, 
    U_h^n, 1   = U_h^n, 0+Δ t_nL(U_h^n, 0), 
    U_h^n, 2   = 3/4U_h^n, 0+1/4(U_h^n, 1+Δ t_nL(U_h^n, 1)), 
    U_h^n, 3   = 1/3U_h^n, 0+2/3(U_h^n, 2+Δ t_nL(U_h^n, 2)), 
    U_h^n+1   = U_h^n, 3, 
    .

where U_h^n denotes the numerical solution at the nth time level, and Δ t_n is the time step-size.


The LDF Runge–Kutta DG method described above works well in smooth regions, but can produce serious numerical nonphysical oscillations around strong discontinuities.
To address this issue, we propose an LDF OE procedure after each Runge–Kutta stage

    U_σ^n, ℓ= ℱ_OEU_h^n, ℓ,   ℓ=1, 2, 3,

resulting in the fully discrete LDF OEDG method 

    {U_σ^n, 0   = U_σ^n, 
    U_h^n, 1   =
    		U_σ^n, 0+Δ t_nL(U_σ^n, 0),     U_σ^n, 1 = ℱ_OEU_h^n, 1, 
    U_h^n, 2   =3/4U_σ^n, 0+1/4(U_σ^n, 1+Δ t_nL(U_σ^n, 1)),     U_σ^n, 2 = ℱ_OEU_h^n, 2,
    U_h^n, 3    =1/3U_σ^n, 0+2/3(U_σ^n, 2+Δ t_nL(U_σ^n, 2)),     U_σ^n, 3 = ℱ_OEU_h^n, 3,
    U_σ^n+1   = U_σ^n, 3, 
    .

where U_σ^n is the modified numerical solution at the n-th time level, and we define U_σ^0 = U_h^0.

In the following, we will elaborate on the OE procedures ℱ_OE for R_h and Q_h, respectively.  It is important to note that the LDF OE procedure we will introduce is specifically designed to maintain the LDF constraint on the approximate magnetic field.










  §.§.§ Component-wise OE procedure for R_h



The OE procedure applied for R_h is denoted as

    R_σ^n, ℓ= ℱ_OE^R R_h^n, ℓ,   ℓ=1, 2, 3.

The OE operator ℱ_OE^R is defined by (ℱ_OE^R R_h^n, ℓ)(x)=R_σ(
x,  τ)|_τ = Δ t_n with
R_σ(x, τ) ∈𝕍_R^k being the solution to the following damping equations:

    {   d/dτ∫_I_ijR_σ v  dx+∑_m=0^k δ_ij^m(R_h^n, ℓ)∘∫_I_ij(R_σ-P^m-1R_σ) v  dx=0  ∀ v ∈ℙ^k(I_ij), 
       R_σ(x, 0)=R_h^n, ℓ(x),
    .

where
τ is a pseudo-time different from t, “∘" denotes the Hadamard product of two vectors,
P^m is the standard L^2-projection operator into 𝕍_R^m for m≥ 0, and we define P^-1=P^0.
The vector δ_ij^m(R_h^n, ℓ) is defined as

    δ_ij^m(R_h):=
    β_ij^x(σ_i+1/2,j^m(R_h)
    +σ_i-1/2,j^m(R_h))/Δ x+β_ij^y(σ_i,j+1/2^m(R_h)+
    σ_i,j-1/2^m(R_h))/Δ y,

where β_ij^x and β_ij^y denote suitable estimates of the local maximum wave speed in the x- and y-directions, respectively,
and σ_i+1/2,j^m(R_h) and σ_i,j+1/2^m(R_h) correspond to the damping coefficients on the x=x_i+1/2 and y=y_j+1/2 edges, respectively.
In practice, we take
β_ij^x and β_ij^y as the spectral radii of
the Jacobian matrices ∂F_1/∂U( U̅_h^n, ℓ) and ∂F_2/∂U( U̅_h^n, ℓ),
 respectively, and U̅_h^n, ℓ denotes the cell average of
U_h^n, ℓ(x) over I_ij.
In equation (<ref>), the damping coefficients are given by

    σ_i+1/2,j^m(R_h)=(σ_i+1/2,j^m
    (R_h^(1)), σ_i+1/2,j^m(R_h^(2)), ⋯, σ_i+1/2,j^m(R_h^(N_R)))^⊤

and

    σ_i,j+1/2^m(R_h)=(σ_i,j+1/2^m
    (R_h^(1)), σ_i,j+1/2^m(R_h^(2)), ⋯, σ_i,j+1/2^m(R_h^(N_R)))^⊤,

with R_h^(μ) being the μ-th component of R_h, and

    σ_i+1/2,j^m(R_h^(μ))=
      {   0,     if  R_h^(μ)≡avg(R_h^(μ)), 
       (2 m+1) (Δ x)^m/2(2 k-1) m !∑_|α|=m1/Δ y∫_y_j-1/2^y_j+1/2|[[∂^αR_h^(μ)]]_i+1/2,j|dy/R_h^(μ)-avg(R_h^(μ))
           _L^∞(Ω),   otherwise,.
      
    σ_i,j+1/2^m(R_h^(μ))=
      {   0,     if  R_h^(μ)≡avg(R_h^(μ)), 
       (2 m+1) (Δ y)^m/2(2 k-1) m !∑_|α|=m1/Δ x∫_x_i-1/2^x_i+1/2|[[∂^αR_h^(μ)]]_i,j+1/2|dx/R_h^(μ)-avg(R_h^(μ))
           _L^∞(Ω),   otherwise..

In equations (<ref>)-(<ref>),
the vector α = (α_1, α_2) is the
multi-index,
∂^αR_h^(μ) is defined as

    ∂^α R_h^(μ)(x):=
    |α|!/α_1 ! α_2 !∂^|α|/∂ x^α_1∂ y^α_2 R_h^(μ)(x),

and
[[∂^αR_h^(μ)]]_
i+1/2,j and [[∂^αR_h^(μ)]]_i,j+1/2 represent the jump values of ∂^αR_h^(μ) across the x=x_i+1/2 and y=y_j+1/2 edges, respectively.


Note that the coefficient δ_ij^m(R_h^n, ℓ) in the damping equations <ref> only depends on the “initial” value R_σ(x, 0)=R_h^n, ℓ(x).
Consequently, the damping equations <ref> are actually a linear ODE system, and its exact
solution can be analytically given.
We now derive the exact solver of damping equations <ref> for R_σ(x, τ).
Based on the local orthogonal basis {ϕ_ij^(α)(x):  |
α| ≤ k} of ℙ^k(I_ij), the “initial” value R_h^n,ℓ(x) of the  damping equations  (<ref>) can be expanded as

    R_h^n, ℓ(x)=∑_μ=0^k∑_|
    α|=μR_ij^(α)(0)
    ϕ_ij^(α)(x)   for  x∈ I_ij.

Assume that the solution of (<ref>) can be expressed as

    R_σ(x, τ)=∑_μ=0^k∑_|
    α|=μR̂_ij^(α)(τ)
    ϕ_ij^(α)(x)   for  x∈ I_ij.

Note that

    (R_σ-P^m-1R_σ)(x, τ)=∑_μ=max{m, 1}^k ∑_|α|=μR̂_ij^(α)(τ) ϕ_ij^(α)(x).

Substitute (<ref>)-(<ref>) into (<ref>) and take v =ϕ_ij^(α)(x). For |α|=μ≥ 1, we get

    d/dτR̂_ij^(α)(τ) ∫_I_ij(ϕ_ij^(α)(x))^2  dx+∑_m=0^μδ_ij^m(R_h^n, ℓ)∘R̂_ij^(α)(τ)  ∫_I_ij(ϕ_ij^(α)(x))^2
     dx=0,

which implies

    d/dτR̂_ij^(α)(τ) +∑_m=0^μδ_ij^m(R_h^n, ℓ)∘R̂_ij^(α)(τ)  =0,   μ=1,2, ⋯, k,

Solving (<ref>) gives

    R̂_ij^(α)(Δ t_n)=
     exp( -Δ t_n∑_m=0^μδ_ij^m(R_h^n, ℓ) )∘R_ij^(α)(0),   μ=1,2, ⋯, k.

For |α|=0, because

    ∫_I_ij(R_σ-P^-1R_σ)ϕ_ij^(α)(x)
     dx = ∫_I_ij(R_σ-P^0R_σ)ϕ_ij^(α)(x)
     dx≡0,
 we have 
    d/dτR̂_ij^(α)(τ) ∫_I_ij(ϕ_ij^(α)(x))^2  dx=0,

and thus

    R̂_ij^(α)(Δ t_n)=
     R_ij^(α)(0)     α= 0 = (0,0).

Hence, with (<ref>) and (<ref>), one has

    R_σ^n,ℓ(x)
       =(ℱ_OE^R R_h^n, ℓ)(x)=R_σ(x, Δ t_n)
       =R_ij^(0)(0)ϕ_ij^(0)(x)+
    ∑_μ=1^kexp( -Δ t_n∑_m=0^μδ_ij^m(R_h^n, ℓ))∘∑_|α|=μR_ij^(α)(0)ϕ_ij^(α)(x).


The OE procedure (<ref>) implemented after each Runge-Kutta stage modifies the high-order moments of numerical polynomials.
Because this OE procedure does not interfere with DG spatial discretization and RK stage update, it can be easily incorporated to existing DG codes as an independent module.
Only simple addition and multiplication operations for cell averages and derivatives are involved, allowing this OE procedure to be extended to general meshes straightforwardly and efficiently.
Note that the OE procedure is applied directly to the conservative variables without the need for characteristic decomposition, which greatly simplifies the implementation and reduces the computational costs.
It is observed that the OE procedure contains information about wave speed and scales linearly with respect to the solution magnitude, resulting in
being capable of eliminating spurious oscillations effectively for problems across different scales and wave speeds.
Furthermore, as described in <cit.>, the OE method remains stable under the normal CFL condition, and many desirable properties of the original DG method are preserved, including conservation, high-order accuracy, optimal fully-discrete error estimates in theory, and numerical superconvergence for linear scalar conservation laws.







  §.§.§ LDF OE procedure for Q_h



After each Runge–Kutta stage, the LDF OE procedure applied for Q_h is denoted as

    Q_σ^n, ℓ= ℱ_OE^Q Q_h^n, ℓ,   ℓ=1, 2, 3.

The LDF OE operator ℱ_OE^Q is defined by (ℱ_OE^Q Q_h^n, ℓ)(x)=Q_σ(
x, Δ t_n) with
Q_σ(x, τ) ∈𝕍^k_Q (0≤τ≤Δ t_n) being the solution to the following initial value problem:

    {   d/dτ∫_I_ijQ_σ·v dx+∑_m=0^k δ_ij^m(Q_h^n, ℓ) ∫_I_ij(Q_σ-𝒫^m-1Q_σ) ·v dx=0   ∀v∈𝕍_Q^k, 
       Q_σ(x, 0)=Q_h^n, ℓ(x),
    .

where
τ is a pseudo-time different from t,
𝒫^m is the L^2-projection operator into 𝕍_Q^m for m≥ 0, and we define 𝒫^-1=𝒫^0.
The coefficient δ_ij^m(Q_h^n, ℓ) is defined as

    δ_ij^m(Q_h):=
    β_ij^x(σ_i+1/2,j^m(Q_h)
    +σ_i-1/2,j^m(Q_h))/Δ x+β_ij^y(σ_i,j+1/2^m(Q_h)+
    σ_i,j-1/2^m(Q_h))/Δ y,

where β_ij^x and β_ij^y are the same as those in equation (<ref>),
and σ_i+1/2,j^m(Q_h) and σ_i,j+1/2^m(Q_h) denote the damping coefficients on the x=x_i+1/2 and y=y_j+1/2 edges, respectively, given by

    σ_i+1/2,j^m(Q_h)=max_1≤μ≤ 2σ_i+1/2,j^m(Q_h^(μ)),     σ_i,j+1/2^m(Q_h)=max_1≤μ≤ 2σ_i,j+1/2^m(Q_h^(μ)),

with Q_h^(μ) being the μ-th component of Q_h, and

    σ_i+1/2,j^m(Q_h^(μ))=
      {   0,     if  Q_h^(μ)≡avg(Q_h^(μ)), 
       (2 m+1) (Δ x)^m/2(2 k-1) m !∑_|α|=m1/Δ y∫_y_j-1/2^y_j+1/2|
           [[∂^αQ_h^(μ)]]_i+1/2,j|dy/Q_h^(μ)-avg(Q_h^(μ))
           _L^∞(Ω),   otherwise,.
    
    σ_i,j+1/2^m(Q_h^(μ))=
      {   0,     if  Q_h^(μ)≡avg(Q_h^(μ)), 
       (2 m+1) (Δ y)^m/2(2 k-1) m !∑_|α|=m1/Δ x∫_x_i-1/2^x_i+1/2|[[∂^αQ_h^(μ)]]_i,j+1/2|dx/Q_h^(μ)-avg(Q_h^(μ))
           _L^∞(Ω),   otherwise..

In equations (<ref>)-(<ref>),
the vector α = (α_1, α_2) is the
multi-index,
∂^αQ_h^(μ) is defined as

    ∂^α Q_h^(μ)(x):=
    |α|!/α_1 ! α_2 !∂^|α|/∂ x^α_1∂ y^α_2 Q_h^(μ)(x),

and
[[∂^αQ_h^(μ)]]_
i+1/2,j and [[∂^αQ_h^(μ)]]_i,j+1/2 represent the jump values of ∂^αQ_h^(μ) across the x=x_i+1/2 and y=y_j+1/2 edges, respectively.


We now  formulate explicitly the exact solver of the damping equation (<ref>) for Q_σ(x, τ).
Based on the local orthogonal basis  {φ_ij^(μ)(x),
μ=0, 1, ⋯, D_Q^k-1} of 𝕍_Q^k over cell I_ij,
the “initial” value Q_σ(x, 0)=Q_h^n, ℓ(x) of the damping equation (<ref>) can be expressed as

    Q_h^n, ℓ(x)=∑_μ=0^k∑_η=D_Q^μ-1^D_Q^μ-1Q_ij^(η)(0)φ_ij^(η)(x)   for  x∈ I_ij.

Assume that the solution of (<ref>) can be represented as

    Q_σ(x, τ)=∑_μ=0^k∑_η=D_Q^μ-1^D_Q^μ-1Q̂_ij^(η)(τ)φ_ij^(η)(x)   for  x∈ I_ij.

Note that

    (Q_σ-𝒫^m-1Q_σ)(x, τ)=∑_μ=max{m, 1}^k∑_η=D_Q^μ-1^D_Q^μ-1Q̂_ij^(η)(τ)φ_ij^(η)(x)   for  x∈ I_ij.

Substitute (<ref>)-(<ref>) into (<ref>) and take v =φ_ij^(η)(x). For any μ≥ 1 and D_Q^μ-1≤η≤ D_Q^μ-1, one obtains

    d/dτQ̂_ij^(η)(τ)
    ∫_I_ijφ_ij^(η)(x)·φ_ij^(η)(x)  dx+∑_m=0^μδ_ij^m(Q_h^n, ℓ)
    Q̂_ij^(η)(τ)  ∫_I_ijφ_ij^(η)(x)·φ_ij^(η)(x) dx=0,

which can be simplified as

    d/dτQ̂_ij^(η)(τ) +∑_m=0^μδ_ij^m(Q_h^n, ℓ)
    Q̂_ij^(η)(τ)  =0,    D_Q^μ-1≤η≤ D_Q^μ-1.

Solving (<ref>) gives

    Q̂_ij^(η)(Δ t_n)=
     exp( -Δ t_n∑_m=0^μδ_ij^m(Q_h^n, ℓ) )Q_ij^(η)(0),    D_Q^μ-1≤η≤ D_Q^μ-1.

For 0 ≤η≤ D_Q^0 -1 =1,
because

    ∫_I_ij(Q_σ-𝒫^m-1Q_σ)·φ_ij^(η)(x)
     dx≡0   ∀ m ≥ 0,
 we have

    d/dτQ̂_ij^(η)(τ) ∫_I_ijφ_ij^(η)(x)
     ·φ_ij^(η)(x) dx=0,    0 ≤η≤ D_Q^0 -1 =1,

 and thus

    Q̂_ij^(η)(Δ t_n)=
     Q_ij^(η)(0),   η=0, ⋯, D_Q^0-1.

Hence, with (<ref>) and (<ref>), one has

    Q_σ^n,ℓ(x)
       =(ℱ_OEQ_h^n, ℓ)(x)=Q_σ(x, Δ t_n)
       =∑_μ=0^D_Q^0-1Q_ij^(η)(0)φ_ij^(η)(x)+
    ∑_μ=1^kexp( -Δ t_n∑_m=0^μδ_ij^m(Q_h^n,ℓ))
    ∑_η=D_Q^μ-1^D_Q^μ-1Q_ij^(η)(0)φ_ij^(η)(x).



The LDF OE procedure modifies the high-order moments of original LDF magnetic field.
Thanks to the exact divergence-free solution of OE solver in the LDF space,
the modified magnetic field naturally retains the LDF property of the original magnetic field.
Moreover, this LDF OE procedure also inherits many notable advantages mentioned in the component-wise OE procedure, including
non-intrusive to be absorbed into existing LDF DG codes,
straightforward and efficient to be extend on general meshes,
stable under normal CFL condition,
scale-invariant and evolution-invariant in eliminating the spurious
oscillations for problems across various scales and wave speeds,
free of characteristic decomposition,
and retaining desirable conservation, high-order accuracy, and optimal convergence rates of the original DG method.
All these advantages highlight the strengths of LDF OE procedure again.























§ POSITIVITY PRESERVATION FOR LDF OEDG METHOD



For physical solutions, both the density and internal energy (or pressure) must be positive. We define the set of physically admissible states as follows:

    𝒢={U=(ρ, m, B, E)^⊤: ρ>0, ℰ(U):=E-1/2(
    	m^2/ρ + B^2)>0},

where ℰ(U) denotes the internal energy.
When the polynomial degree k ≥ 1,
method does not always ensure that U_σ^n∈𝒢.
Based on the Zhang–Shu framework <cit.>,
to achieve a bound-preserving DG scheme for hyperbolic conservation laws, one must first identify a DG scheme that preserves the updated cell averages within the bounds. Subsequently, a local scaling bound-preserving limiter is used to ensure that the DG solution satisfies these bounds at the points of interest.
According to the discoveries in <cit.>, for the conservative MHD system (<ref>), the PP property of the updated cell averages of the DG solution is closely connected with a set of globally coupled discrete divergence-free conditions, which are, however, incompatible with the standard local scaling limiter. An effective way <cit.> addressing this issue is to consider the following modified MHD system with an additional source term:

    ∂U/∂ t+∂F_1(U)/∂x
    	+∂F_2(U)/∂y=
    	-(∇·B)S(U),

where S(U)=(0, B, u, u·B)^⊤ and the right-hand side of (<ref>) is termed the Godunov–Powell source term <cit.>. As demonstrated in <cit.>, this source term can help to eliminate the effect of divergence error on the PP property.


In this section, we will utilize the optimal convex decomposition approach proposed in <cit.> to derive an efficient, provably PP LDF OEDG method. This method will be built on a PP HLL flux, a properly upwind discretization of the Godunov–Powell source term, and a simple scaling PP limiter. Our method automatically ensures that the updated cell averages of the OEDG solution remain within 𝒢, and the PP limiter enforces the values of the OEDG solution in 𝒢 at any points of interest.






 §.§ Positivity-preserving HLL flux



The HLL fluxes, with ℓ=1 for the x-direction and ℓ=2 for the y-direction,  are given by

    F_ℓ(U^-, U^+)=
    	{   F_ℓ(U^-),     0 ≤𝒱_ℓ,l<𝒱_ℓ,r, 
       𝒱_ℓ,rF_ℓ(U^-)-
    			𝒱_ℓ,lF_ℓ(U^+)+
    			𝒱_ℓ,l𝒱_ℓ,r(U^+-U^-)/𝒱_ℓ,r-𝒱_ℓ,l,      𝒱_ℓ,l<0<𝒱_ℓ,r, 
       F_ℓ(U^+),    𝒱_ℓ,l<𝒱_ℓ,r≤ 0,
    	.

where 𝒱_ℓ,l and 𝒱_ℓ,r are functions of U^- and U^+, denoting the estimated minimum and maximum wave speeds, respectively.



The estimates of {𝒱_ℓ,l} and {𝒱_ℓ,r} are crucial for the PP property <cit.>. 

For any pair of admissible states U and Ũ, define

    α_ℓ,l(U, Ũ):=min{u_ℓ,√(ρ)u_ℓ+√(ρ̃)ũ_ℓ/√(ρ)+√(ρ̃)}-𝒞_ℓ(U)-
    	|B-B̃|/√(ρ)+√(ρ̃),
    	  ℓ=1, 2,


    α_ℓ,r(U, Ũ):=max{u_ℓ,√(ρ)u_ℓ+√(ρ̃)ũ_ℓ/√(ρ)+√(ρ̃)}+𝒞_ℓ(U)+
    	|B-B̃|/√(ρ)+√(ρ̃),
    	  ℓ=1, 2,

and

    α_ℓ,⋆(U, Ũ):=max{|u_ℓ|,|√(ρ)u_ℓ+√(ρ̃)ũ_ℓ/√(ρ)+√(ρ̃)|}+𝒞_ℓ(U)+
    	|B-B̃|/√(ρ)+√(ρ̃),
    	  ℓ=1, 2,

where

    𝒞_ℓ(U):=1/√(2)[𝒞_s^2+
    	B^2/ρ+√((𝒞_s^2+B^2/ρ)^2-4 𝒞_s^2 B_ℓ^2/ρ)]^1/2  with  𝒞_s=√((γ-1)p/2ρ),   ℓ=1, 2.

According to the rigorous PP analysis in <cit.>,  {𝒱_ℓ,l} and {𝒱_ℓ,r} should satisfy

    𝒱_ℓ,l≤α_ℓ,l(U^-, U^+),   𝒱_ℓ,r≥α_ℓ,r(U^+, U^-),   ℓ=1, 2.

We take 𝒱_ℓ,l and 𝒱_ℓ,r as

    𝒱_ℓ,l   =min{α_ℓ,l(U^-, U^+), λ_ℓ,min(U^-), λ_ℓ,min(U^Roe), λ_ℓ,min(U^+)  },   ℓ=1, 2,
    𝒱_ℓ,r   =max{α_ℓ,r(U^+, U^-), λ_ℓ,max(U^-), λ_ℓ,max(U^Roe), λ_ℓ,max(U^+)    },   ℓ=1, 2,

where λ_ℓ,min/max(U) are the minimum and maximum eigenvalues of the Jacobian matrix ∂F_ℓ/∂U, respectively, and  λ_ℓ,min/max(U^Roe) are those of the Roe matrix <cit.>. Let
𝒱_ℓ^-=min{𝒱_ℓ,l, 0} and
𝒱_ℓ^+=max{𝒱_ℓ,r, 0},
then the HLL fluxes in (<ref>) can be rewritten as

    F_ℓ(U^-,U^+)=
    	𝒱_ℓ^+F_ℓ(U^-)-
    		𝒱_ℓ^-F_ℓ(U^+)+
    		𝒱_ℓ^-𝒱_ℓ^+(U^+-U^-)/𝒱_ℓ^+-𝒱_ℓ^-,   ℓ=1, 2.














 §.§ Upwind discretization of Godunov–Powell source term





 the DDF condition is not necessary for PP property due to the application of properly discretized Godunov-Powell source term.
In fact, properly discretized Godunov-Powell source term helps eliminate the effect of the jump of the normal magnetic field on the cell interface.


Since the PP property of the scheme is only related to the discrete equations satisfied by the cell averages, we only consider adding Godunov-Powell source term to the evolution equations of the cell averages.
That is, only the evolution of the cell averages adopts discrete scheme of the symmetrizable MHD equations (<ref>), while the evolutions of the higher-order “moments" still use discrete scheme of the conservative MHD equations, as shown in Section <ref>.









Let U_σ^n(x, y) denote the approximate polynomial solution of degree k at n-th time level after LDF OE procedure, and U̅_ij^n be the cell average over cell I_ij at n-th time level.
For one-order forward Euler method,
the evolution equations satisfied by the cell averages in the symmetrizable MHD equations with the Godunov-Powell source term are given by

    U̅_ij^n+1=U̅_ij^n+Δ t_nL_ij(U_σ^n),

where

    L_ij(U_σ^n)
       :=-1/Δ x∑_μ=1^qω_μ^G[(
    F_1_, i+1/2(y_j^(μ))-
    F_1_, i-1/2(y_j^(μ))).
       +.ℬ_1, i+1/2(y_j^(μ))S(U_σ^n(x_i+1/2^-, y_j^(μ)))-
    ℬ_1, i-1/2(y_j^(μ))S(U_σ^n(x_i-1/2^+, y_j^(μ)))
    ].
    
       -1/Δ y∑_μ=1^qω_μ^G[(
    F_2_, j+1/2(x_i^(μ))-
    F_2_, j-1/2(x_i^(μ))).
       +.ℬ_2, j+1/2(x_i^(μ))S(U_σ^n(x_i^(μ), y_j+1/2^-))-
    ℬ_2, j-1/2(x_i^(μ))S(U_σ^n(x_i^(μ), y_j-1/2^+))
    ].

Here {x_i^(μ)}_μ=1^q and {y_j^(μ)}_μ=1^q
are the q-point Gauss-Legendre quadrature nodes with 2q-1≥2k in [x_i-1/2, x_i+1/2] and [y_j-1/2, y_j+1/2], respectively, and {ω_μ^G}_μ=1^q are the associated weights.
In (<ref>),

    F_1_, i+1/2(y):=
      F_1(
      U_σ^n(x_i+1/2^-, y),
      U_σ^n(x_i+1/2^+, y)),

and

    F_2_, j+1/2(x):=
      F_2(
      U_σ^n(x, y_j+1/2^-),
      U_σ^n(x, y_j+1/2^+)),

where F_1 and F_2 denote HLL numerical fluxes.
In (<ref>), we take

    ℬ_1, i±1/2(y)=
      -𝒱_1,i±1/2^∓(y)/𝒱_1,i±1/2^+(y)-
      𝒱_1,i±1/2^-(y)[(B_1)_σ^n(x_i±1/2^+, y)-
      (B_1)_σ^n(x_i±1/2^-, y)],

and

    ℬ_2, j±1/2(x)=
      -𝒱_2,j±1/2^∓(x)/𝒱_2,j±1/2^+(x)-
      𝒱_2,j±1/2^-(x)[(B_2)_σ^n(x, y_j±1/2^+)-
      (B_2)_σ^n(x, y_j±1/2^-)],

and 
with the definitions of

    𝒱_1,i+1/2^±(y):=𝒱_1^±(
      U_σ^n(x_i+1/2^-, y),
      U_σ^n(x_i+1/2^+, y))

and

    𝒱_2,j+1/2^±(x):=𝒱_2^±(
      U_σ^n(x, y_j+1/2^-),
      U_σ^n(x, y_j+1/2^+)),

where 𝒱_1^± and 𝒱_2^± are the wave speeds in the HLL flux, which are same as that in Section <ref>.
It is worth mentioning that the discretized Godunov-Powell source term in (<ref>) is motivated from the analysis for provably PP property, as shown in the following Theorem 2.

In fact, we can split the discrete divergence into two parts:

    div_I_ijB^n=1/Δ x Δ y∫_I_ij∇·B^n(x) dx=1/Δ x Δ y∮_I_ijB^n(x)·n dl
       =1/Δ x∑_μ=1^qω_μ^G[
    𝒱_1,i+1/2^+(y_j^(μ))
    (B_1)_σ^n(x_i+1/2^-, y_j^(μ))-
    𝒱_1,i+1/2^-(y_j^(μ))
    (B_1)_σ^n(x_i+1/2^+, y_j^(μ))/𝒱_1,i+1/2^+(y_j^(μ))-
    𝒱_1,i+1/2^-(y_j^(μ)).
       +.𝒱_1,i-1/2^+(y_j^(μ))
    (B_1)_σ^n(x_i-1/2^+, y_j^(μ))-
    𝒱_1,i-1/2^-(y_j^(μ))
    (B_1)_σ^n(x_i-1/2^-, y_j^(μ))/𝒱_1,i-1/2^+(y_j^(μ))-
    𝒱_1,i-1/2^-(y_j^(μ))]
       +1/Δ y∑_μ=1^qω_μ^G[
    𝒱_2,j+1/2^+(x_i^(μ))
    (B_2)_σ^n(x_i^(μ), y_j+1/2^-)-
    𝒱_2,j+1/2^-(x_i^(μ))
    (B_2)_σ^n(x_i^(μ), y_j+1/2^+)/𝒱_2,j+1/2^+(x_i^(μ))-
    𝒱_2,j+1/2^-(x_i^(μ)).
       +.𝒱_2,j-1/2^+(x_i^(μ))
    (B_2)_σ^n(x_i^(μ), y_j-1/2^+)-
    𝒱_2,j-1/2^-(x_i^(μ))
    (B_2)_σ^n(x_i^(μ), y_j-1/2^-)/𝒱_2,j-1/2^+(x_i^(μ))-
    𝒱_2,j-1/2^-(x_i^(μ))]
       =div_withinB^n+
    div_facesB^n,

where

    div_withinB^n=1/Δ x∑_μ=1^qω_μ^G[(B_1)_σ^n(x_i+1/2^-,
      y_j^(μ))-
      (B_1)_σ^n(x_i-1/2^+, y_j^(μ))]
       +1/Δ y∑_μ=1^qω_μ^G[(B_2)_σ^n(x_i^(μ), y_j+1/2^-)-
      (B_2)_σ^n(x_i^(μ),y_j-1/2^+)],

and

    div_facesB^n=
      1/Δ x∑_μ=1^qω_μ^G[
    -𝒱_1,i+1/2^-(y_j^(μ))/𝒱_1,i+1/2^+(y_j^(μ))-
    𝒱_1,i+1/2^-(y_j^(μ))((B_1)_σ^n(x_i+1/2^+, y_j^(μ))-
    (B_1)_σ^n(x_i+1/2^-, y_j^(μ)))
    .
       +.𝒱_1,i-1/2^+(y_j^(μ))/𝒱_1,i-1/2^+(y_j^(μ))-
    𝒱_1,i-1/2^-(y_j^(μ))((B_1)_σ^n(x_i-1/2^+, y_j^(μ))-
    (B_1)_σ^n(x_i-1/2^-, y_j^(μ)))]
       +1/Δ y∑_μ=1^qω_μ^G[
    -𝒱_2,j+1/2^-(x_i^(μ))/𝒱_2,j+1/2^+(x_i^(μ))-
    𝒱_2,j+1/2^-(x_i^(μ))((B_2)_σ^n(x_i^(μ),
    y_j+1/2^+)-
    (B_2)_σ^n(x_i^(μ),
    y_j+1/2^-))
    .
       +.𝒱_2,j-1/2^+(x_i^(μ))/𝒱_2,j-1/2^+(x_i^(μ))-
    𝒱_2,j-1/2^-(x_i^(μ))((B_2)_σ^n(x_i^(μ),
    y_j-1/2^+)-
    (B_2)_σ^n(x_i^(μ),
    y_j-1/2^-))].

Note that “div_withinB^n" means the divergence within the
cell I_ij, which becomes zero because of locally divergence-free property,
while “div_facesB^n" involves divergence error across the cell interfaces, which ia consistent with the discretized divergence in Godunov-Powell source term.








 §.§ Convex decomposition of cell averages



In the framework of provably PP high-order schemes, a key ingredient
is the decomposition of the cell averages of the numerical solution into
a convex combination of the solution values at certain quadrature points.
The convex decomposition determines not only the theoretical PP CFL condition of the resulting scheme, but also the points to perform the PP limiter, as will be shown in the Section <ref>.

Let {x_i^(μ)}_μ=1^q and {y_j^(μ)}_μ=1^q
denote the q-point Gauss-Legendre quadrature nodes with 2q-1≥2k in [x_i-1/2, x_i+1/2] and [y_j-1/2, y_j+1/2], respectively, and {ω_μ^G}_μ=1^q be the associated weights.
Denote {x̂_i^(μ)}_μ=1^L and {ŷ_j^(μ)}_μ=1^L
as the L-point Gauss-Lobatto quadrature nodes with 2L-3≥ k in [x_i-1/2, x_i+1/2] and [y_j-1/2, y_j+1/2], respectively, and {ω̂_μ^GL}_μ=1^L are the associated weights with ω̂_1^GL=ω̂_L^GL=1/L(L-1).
Let a_1 and a_2 denote the maximum characteristic speeds
in x- and y-directions.

According to <cit.>, a general feasible convex decomposition for polynomial space ℙ^k(I_ij) on Cartesian mesh can be written as

    U̅_ij=   ∑_μ=1^q ω_μ^G[ω_1^+U_ij(x_i+1/2^-, y_j^(μ))+ω_1^-U_ij(x_i-1/2^+, y_j^(μ))+ω_2^+U_ij(x_i^(μ), y_j+1/2^-)+ω_2^-U_ij(x_i^(μ), y_j-1/2^+)
     ]
       +∑_s=1^S ω_s U_ij(x_ij^(s), y_ij^(s)),

and it simultaneously satisfies the following three conditions:

  
  (i) the convex decomposition holds exactly for all U_ij(x, y)∈ [ℙ^k(I_ij)]^8;
  
  (ii) the weights {ω_1^±, ω_2^±, ω_s} are all positive and their summation equals one;
  
  (iii) the internal node set 𝕊_ij^internal:={(x_ij^(s), y_ij^(s))}_s=1^S⊂ I_ij.

In convex decomposition (<ref>), the boundary node set is given by

    𝕊_ij^boundary:={(x_i+1/2^-, y_j^(μ))}_μ=1^q⋃{(x_i-1/2^+, y_j^(μ))}_μ=1^q⋃{(x_i^(μ),
      y_j+1/2^-, )}_μ=1^q⋃{(x_i^(μ),
      y_j-1/2^+, )}_μ=1^q,

and the cell average of numerical solution is decomposed into a convex combination of solution values at the
quadrature points belonging to 𝕊_ij:=𝕊_ij^boundary⋃𝕊_ij^internal.


Based on the tensor product of the L-point Gauss-Lobatto quadrature with
L=⌈k+3/2⌉ and the q-point Gauss quadrature with
q=k+1,
Zhang and Shu in <cit.> proposed a classic convex decomposition, which is listed as follows.

Zhang-Shu convex decomposition.

    U̅_ij=   ∑_μ=1^q ω_μ^Gω̂_1^GL[κ_1 U_ij(x_i+1/2^-, y_j^(μ))+κ_1
    U_ij(x_i-1/2^+, y_j^(μ)).
       .+κ_2 U_ij(x_i^(μ), y_j+1/2^-, )+κ_2
    U_ij(x_i^(μ), y_j-1/2^+)
     ]
       +∑_s=2^L-1∑_μ=1^qω̂_s^GLω_μ^G[κ_1U_ij(x̂_i^(s), y_j^(μ))+κ_2U_ij(x_i^(μ), ŷ_j^(s))]

with

    κ_1:=a_1/Δ x/a_1/Δ x+a_2/Δ y,   κ_2:=a_2/Δ y/a_1/Δ x+a_2/Δ y.

Cui et al. <cit.> discovered an optimal convex decomposition allowing much fewer nodes and larger time step-size, as shown below for the case of ℙ^2 or ℙ^3 space.

Optimal convex decomposition.

    U̅_ij=   κ_1/2∑_μ=1^q ω_μ^G[ U_ij(x_i+1/2^-, y_j^(μ))+
    U_ij(x_i-1/2^+, y_j^(μ))]
       +κ_2/2∑_μ=1^qω_μ^G[ U_ij(x_i^(μ), y_j+1/2^-, )+
    U_ij(x_i^(μ), y_j-1/2^+)
     ]
       +ω∑_sU_ij(x̂_s, ŷ_s)

with

    {(x̂_s, ŷ_s)}=
    {   (x_i, y_j ±Δ y/2 √(3)√(ϕ_*-ϕ_2/ϕ_*)),     if ϕ_1 ≥ϕ_2, 
       (x_i ±Δ x/2 √(3)√(ϕ_*-ϕ_1/ϕ_*), y_j),     if ϕ_1<ϕ_2,
        .

where

    ϕ_1=a_1/Δ x,   ϕ_2=a_2/Δ y,   ϕ_*=max{ϕ_1, ϕ_2},   ψ=ϕ_1+ϕ_2+2 ϕ_*,   κ_1=ϕ_1/ψ,   κ_2=ϕ_2/ψ,   ω=ϕ_*/ψ .

In Figure <ref>, we show the nodes of optimal convex decomposition and classic convex decompositions for the polynomial space ℙ^2(I_ij), on the square grid, in the case of Δ x/a_1=Δ y/a_2.










 §.§ Positivity-preserving limiter



Let U_ij^n(x) denote the approximate solution polynomial  in I_ij at the n-th time level after the OE procedure, and U̅_ij^n be the cell average of U_ij^n(x) over I_ij.
Following the analysis in <cit.>, the OEDG method with the PP HLL flux and properly discretized Godunov-Powell source term preserves the positivity for the updated cell averages under a CFL condition, if the OEDG solution polynomials satisfy

    U_ij^n(x) ∈𝒢,   ∀x∈𝕊_ij

where 𝕊_ij is all the nodes of convex decomposition (<ref>).


The OE procedure, although enhancing the numerical stability for strong shocks, does not necessarily guarantee the condition (<ref>), which should be enforced by a simple PP limiter. 
For readers' convenience, here we review the PP limiter <cit.> briefly.
The PP limiting procedure mainly consists of two steps.
First, modify the density to enforce the positivity by

    ρ_ij(x)=θ_1(ρ_ij^n(x)-ρ̅_ij^n)+
    	ρ̅_ij^n    with  θ_1=min{1, ρ̅_ij^n-ϵ_1/ρ̅_ij^n-min _x∈𝕊_ijρ_ij^n(x)} .

Then, modify U_ij(x):=
(ρ_ij(x), m_ij^n(x), B_ij^n(x), E_ij^n(x))^⊤ to ensure positive internal energy by

    U_ij^n(x)=θ_2
    	(U_ij(x)-U̅_ij^n)+
    	U̅_ij^n   with  θ_2=min{1, ℰ(U̅_ij^n)-ϵ_2/ℰ(U̅_ij^n)-min _x∈𝕊_ijℰ(U_ij(x))}.

Here ϵ_1=min{10^-13, ρ̅_ij^n} and ϵ_2=min{10^-13, ℰ(U̅_ij^n)} are introduced to avoid rounding error <cit.>.
Note that the modified polynomial solution has the same cell average U̅_ij^n and holds

    U_ij^n(x) ∈𝒢_ϵ=
    	{U: ρ≥ϵ_1, ℰ(U) ≥ϵ_2},   ∀x∈𝕊_ij.

It is worth mentioning that such a limiter does not reduce the high-order accuracy of the DG schemes <cit.>. The PP limiter is performed right after each OE procedure.








 §.§ Algorithmic flow of PP LDF OEDG scheme


Based on the LDF OEDG scheme, properly discretized Godunov-Powell source term and the convex combination of cell averages, we construct the provably PP LDF OEDG scheme as follows.


  
  Step 0 Initialize {U̅_ij^0} and {U_ij^0(x)} by the L^2-projection of initial data onto 𝕍_U^k, see (<ref>)-(<ref>).
      {U̅_ij^0}∈𝒢 due to the convexity of 𝒢.
  
  Step 1 Given admissible cell averages {U̅_ij^n} and approximate polynomial solution U_σ^n(x)
      ∈𝕍_U^k defined by U_ij^n(x), perform PP limiter to modify the polynomials {U_ij^n(x)} such that
      the modified polynomials {U_ij^n(x)} satisfy
      
    U_ij^n(x)∈𝒢,
              ∀x∈𝕊_ij,

      where 𝕊_ij is all the nodes of convex decomposition (<ref>). The PP limiting procedure can refer to Section <ref>.
  
  Step 2 Let U_σ^n(x)
      ∈𝕍_U^k denote the discontinuous piecewise polynomial solution defined by U_ij^n(x). For one-order forward Euler method, update the cell averages by the scheme
      
    U̅_ij^n+1=U̅_ij^n+Δ t_nL_ij(U_σ^n).

      As will be shown in Theorem <ref>, the scheme (<ref>) can ensure U̅_ij^n+1∈𝒢.
  
  Step 3 For one-order forward Euler method, evolve the high-order “moments" of the polynomials {U_ij^n+1(x)} by using the LDF OEDG scheme detailed in Section <ref> to build the discontinuous piecewise polynomial function U_σ^n+1(x)∈𝕍_U^k.
  
  Step 4 If t_n+1<T, assign n← n+1 and go to Step 1; otherwise, output numerical results.





 §.§ Rigorous PP analysis



In order to prove that U̅_ij^n+1∈𝒢 in scheme (<ref>), we first provide several key auxiliary theories, and the rigorous proof of the PP property of the scheme will be given by Theorem <ref> later.


The admissible state set 𝒢 is equivalent to

    𝒢_*={U=(ρ, m, B, E)^⊤: ρ>0, U·n^*+|B^*|^2/2>0, ∀v^*, B^* ∈ℝ^3},

where

    n^*=(|v^*|^2/2,-v^*,
    -B^*, 1)^⊤.



The proof of Lemma <ref> can be found in <cit.>.
For more details about its geometric interpretation, please see the geometric quasilinearization (GQL) approach in <cit.>.



For any U∈𝒢 and any v^*, B^*∈ℝ^3, we have

    S(U) ·n^*=(v-v^*) ·(B-B^*)-v^* ·B^*,


    |√(ρ)(v-v^*) ·(B-B^*)|<U·n^*+|B^*|^2/2.

Furthermore, for any b∈ℝ, it holds

    -b(S(U) ·n^*) ≥ b(v^* ·B^*)-|b|/√(ρ)(U·n^*+|B^*|^2/2).



The proof of Lemma <ref> can be found in <cit.>.



Given admissible states U^(j), 1≤ j≤4, we define

    α̂_1:=𝒞_1(U^(1))+
    Δ y/Δ x+Δ y|B^(1)-B^(2)|/√(ρ^(1))+√(ρ^(2))+
    Δ x/Δ x+Δ y|B^(1)-B^(3)|/√(ρ^(1))+√(ρ^(3))+
    Δ x/Δ x+Δ y|B^(1)-B^(4)|/√(ρ^(1))+√(ρ^(4))+
       +max{(u_1)^(1), Δ y/Δ x+Δ y√(ρ^(1))(u_1)^(1)+√(ρ^(2))(u_1)^(2)/√(ρ^(1))+√(ρ^(2)).
       +Δ x/2(Δ x+Δ y)√(ρ^(1))((u_1)^(1)-(u_2)^(1))+
    √(ρ^(3))((u_1)^(3)-(u_2)^(3))/√(ρ^(1))+√(ρ^(3))
       +.Δ x/2(Δ x+Δ y)√(ρ^(1))((u_1)^(1)+(u_2)^(1))+
    √(ρ^(4))((u_1)^(4)+(u_2)^(4))/√(ρ^(1))+√(ρ^(4))},


    α̂_2:=𝒞_1(U^(2))+
    Δ y/Δ x+Δ y|B^(2)-B^(1)|/√(ρ^(2))+√(ρ^(1))+
    Δ x/Δ x+Δ y|B^(2)-B^(3)|/√(ρ^(2))+√(ρ^(3))+
    Δ x/Δ x+Δ y|B^(2)-B^(4)|/√(ρ^(2))+√(ρ^(4))+
       +max{
    -(u_1)^(2), -Δ y/Δ x+Δ y√(ρ^(2))(u_1)^(2)+√(ρ^(1))(u_1)^(1)/√(ρ^(2))+√(ρ^(1)).
       +-Δ x/2(Δ x+Δ y)√(ρ^(2))((u_1)^(2)+(u_2)^(2))+
    √(ρ^(3))((u_1)^(3)+(u_2)^(3))/√(ρ^(2))+√(ρ^(3))
       +.Δ x/2(Δ x+Δ y)√(ρ^(2))((u_2)^(2)-(u_1)^(2))+
    √(ρ^(4))((u_2)^(4)-(u_1)^(4))/√(ρ^(2))+√(ρ^(4))},


    α̂_3:=𝒞_2(U^(3))+
    Δ x/Δ x+Δ y|B^(3)-B^(4)|/√(ρ^(3))+√(ρ^(4))+
    Δ y/Δ x+Δ y|B^(3)-B^(1)|/√(ρ^(3))+√(ρ^(1))+
    Δ y/Δ x+Δ y|B^(3)-B^(2)|/√(ρ^(3))+√(ρ^(2))+
       +max{(u_2)^(3), Δ x/Δ x+Δ y√(ρ^(3))(u_2)^(3)+√(ρ^(4))(u_2)^(4)/√(ρ^(3))+√(ρ^(4)).
       +Δ y/2(Δ x+Δ y)√(ρ^(3))(-(u_1)^(3)+(u_2)^(3))+
    √(ρ^(1))(-(u_1)^(1)+(u_2)^(1))/√(ρ^(3))+√(ρ^(1))
       +.Δ y/2(Δ x+Δ y)√(ρ^(3))((u_1)^(3)+(u_2)^(3))+
    √(ρ^(2))((u_1)^(2)+(u_2)^(2))/√(ρ^(3))+√(ρ^(2))},


    α̂_4:=𝒞_2(U^(4))+
    Δ x/Δ x+Δ y|B^(4)-B^(3)|/√(ρ^(4))+√(ρ^(3))+
    Δ y/Δ x+Δ y|B^(4)-B^(1)|/√(ρ^(4))+√(ρ^(1))+
    Δ y/Δ x+Δ y|B^(4)-B^(2)|/√(ρ^(4))+√(ρ^(2))+
       +max{
    -(u_2)^(4), -Δ x/Δ x+Δ y√(ρ^(4))(u_2)^(4)+√(ρ^(3))(u_2)^(3)/√(ρ^(4))+√(ρ^(3)).
       +-Δ y/2(Δ x+Δ y)√(ρ^(4))((u_1)^(4)+(u_2)^(4))+
    √(ρ^(1))((u_1)^(1)+(u_2)^(1))/√(ρ^(4))+√(ρ^(1))
       +.Δ y/2(Δ x+Δ y)√(ρ^(4))((u_1)^(4)-(u_2)^(4))+
    √(ρ^(2))((u_1)^(2)-(u_2)^(2))/√(ρ^(4))+√(ρ^(2))}.

Then for any α_j≥α̂_j, the state

    U̅:=1/Δ y(α_1+α_2)+Δ x(α_3+α_4)[
    Δ y(α_1U^(1)-F_1(U^(1)))+Δ y(α_2U^(2)+F_1(U^(2))).
       +.Δ x(α_3U^(3)-F_2(U^(3)))+Δ x(α_4U^(4)+F_2(U^(4)))
    ],

belongs to 𝒢_ρ:={U=(ρ, m, B, E)^⊤: ρ>0}, and satisfies

    U̅·n^*+|B^*|^2/2≥-v^*·B^*(
    Δ y (B_1)^(1)-Δ y (B_1)^(2)+
    Δ x (B_2)^(3)-Δ x (B_2)^(4))/Δ y(α_1+α_2)+Δ x(α_3+α_4),
    ∀v^*,B^*∈ℝ^3.

Furthermore, Δ y(α_1+α_2)+Δ x(α_3+α_4)>0.


 Proof 
 This directly follows from Theorem 1 in <cit.> with N=4, by taking
 
    s_1=s_2=Δ y,   s_3=s_4=Δ x,  ξ^(1)=(1,0),  ξ^(2)=(-1,0),  ξ^(3)=(0,1),  ξ^(4)=(0,-1).





Assume U^-, U^+∈𝒢.
If the approximate wave speeds in the HLL flux satisfy

    𝒱_ℓ,r≥α_r(U^+, U^-),   𝒱_ℓ,l≥α_l(U^-, U^+),   ℓ=1,2,

then

    F_ℓ(U^-, U^+)=
      𝒱_ℓ^-H_ℓ(U^-, U^+)+F_ℓ(U^-)-
      𝒱_ℓ^-U^-


    F_ℓ(U^-, U^+)=
      𝒱_ℓ^+H_ℓ(U^-, U^+)+F_ℓ(U^+)-
      𝒱_ℓ^+U^+

where the intermediate state

    H_ℓ(U^-, U^+):=
      1/𝒱_ℓ^+-𝒱_ℓ^-(
      𝒱_ℓ^+U^+-F_ℓ(U^+)
      -𝒱_ℓ^-U^-+F_ℓ(U^-)
      )

belongs to 𝒢 and satisfies

    H_ℓ(U^-, U^+) ·n^*+|B^*|^2/2≥
      -v^* ·B^*/𝒱_ℓ^+-𝒱_ℓ^-( B_ℓ^+- B_ℓ^-),   ∀v^*, B^* ∈ℝ^3.



Proof 
 This directly follows from Theorem 2 in <cit.> by taking
 ξ=(1,0) for ℓ=1 and ξ=(0,1) for ℓ=2. 

With the above auxiliary theories as the key cornerstone, we now proof the PP property of the scheme (<ref>).
To shorten the notations, we define

    U_i+1/2,j^±, μ:=U_σ^n(x_i+1/2^±, y_j^(μ)),   U_i, j+1/2^μ, ±:=U_σ^n(x_i^(μ), y_j+1/2^±),


    𝒱_1, i+1/2^μ, ±:=
    𝒱_1^±(
    U_σ^n(x_i+1/2^-, y_j^(μ)),
    U_σ^n(x_i+1/2^+, y_j^(μ))
    ),


    𝒱_2, j+1/2^μ, ±:=
    𝒱_2^±(
    U_σ^n(x_i^(μ), y_j+1/2^-),
    U_σ^n(x_i^(μ), y_j+1/2^+)
    ),


    H_1, i+1/2^μ:=
    H_1(
    U_σ^n(x_i+1/2^-, y_j^(μ)),
    U_σ^n(x_i+1/2^+, y_j^(μ))
    ),


    H_2, j+1/2^μ:=
    H_2(
    U_σ^n(x_i^(μ), y_j+1/2^-),
    U_σ^n(x_i^(μ), y_j+1/2^+)
    ),


    α̂_1, i±1/2^μ, ∓=𝒞_1(U_i±1/2,j^∓, μ)+
    Δ y/Δ x+Δ y|B_i±1/2,j^∓, μ - B_i∓1/2,j^±, μ|/√(ρ_i±1/2,j^∓, μ)+√(ρ_i∓1/2,j^±, μ)+
    Δ x/Δ x+Δ y|B_i±1/2,j^∓, μ - B_i,j+1/2^μ,-|/√(ρ_i±1/2,j^∓, μ)+√(ρ_i,j+1/2^μ,-)
       +Δ x/Δ x+Δ y|B_i±1/2,j^∓, μ -B_i,j-1/2^μ,+|/√(ρ_i±1/2,j^∓, μ)+√(ρ_i,j-1/2^μ,+)
       +max{±(u_1)_i±1/2,j^∓,μ,
    ±Δ y/(Δ x+Δ y)√(ρ_i±1/2,j^∓, μ)(u_1)_i±1/2,j^∓,μ+
    √(ρ_i∓1/2,j^±, μ)(u_1)_i∓1/2,j^±,μ/√(ρ_i±1/2,j^∓, μ)+√(ρ_i∓1/2,j^±, μ).
       +Δ x/2(Δ x+Δ y)√(ρ_i±1/2,j^∓, μ)(±(u_1)_i±1/2,j^∓,μ-
    (u_2)_i±1/2,j^∓,μ)+
    √(ρ_i,j+1/2^μ,-)(±(u_1)_i,j+1/2^μ,--
    (u_2)_i,j+1/2^μ,-)/√(ρ_i±1/2,j^∓, μ)+√(ρ_i,j+1/2^μ,-)
       +.
    Δ x/2(Δ x+Δ y)√(ρ_i±1/2,j^∓, μ)(±(u_1)_i±1/2,j^∓,μ+
    (u_2)_i±1/2,j^∓,μ)+
    √(ρ_i,j-1/2^μ,+)(±(u_1)_i,j-1/2^μ,++
    (u_2)_i,j-1/2^μ,+)/√(ρ_i±1/2,j^∓, μ)+√(ρ_i,j-1/2^μ,+)},


    α̂_2, j±1/2^μ, ∓=𝒞_2(U_i,j±1/2^μ,∓)+
    Δ x/Δ x+Δ y|B_i,j±1/2^μ,∓ - B_i,j∓1/2^μ,±|/√(ρ_i,j±1/2^μ,∓)+√(ρ_i,j∓1/2^μ,±)+
    Δ y/Δ x+Δ y|B_i,j±1/2^μ,∓ - B_i+1/2,j^-,μ|/√(ρ_i,j±1/2^μ,∓)+
    √(ρ_i+1/2,j^-,μ)
       +Δ y/Δ x+Δ y|B_i,j±1/2^μ,∓ -B_i-1/2,j^+,μ|/√(ρ_i,j±1/2^μ,∓)+√(ρ_i-1/2,j^+,μ)
       +max{±(u_2)_i,j±1/2^μ,∓,
    ±Δ x/(Δ x+Δ y)√(ρ_i,j±1/2^μ,∓)(u_2)_i,j±1/2^μ,∓+
    √(ρ_i,j∓1/2^μ,±)(u_2)_i,j∓1/2^μ,±/√(ρ_i,j±1/2^μ,∓)+√(ρ_i,j∓1/2^μ,±).
       +Δ y/2(Δ x+Δ y)√(ρ_i,j±1/2^μ,∓)(-  (u_1)_i,j±1/2^μ,∓±(u_2)_i,j±1/2^μ,∓)+
    √(ρ_i+1/2,j^-,μ)(-  (u_1)_i+1/2,j^-,μ±(u_2)_i+1/2,j^-,μ)/√(ρ_i,j±1/2^μ,∓)+√(ρ_i+1/2,j^-,μ)
       +.
    Δ y/2(Δ x+Δ y)√(ρ_i,j±1/2^μ,∓)(  (u_1)_i,j±1/2^μ,∓±(u_2)_i,j±1/2^μ,∓)+
    √(ρ_i-1/2,j^+,μ)(  (u_1)_i-1/2,j^+,μ±(u_2)_i-1/2,j^+,μ)/√(ρ_i,j±1/2^μ,∓)+√(ρ_i-1/2,j^+,μ)}.


The PP property of the scheme (<ref>) is shown as follows.

Let the wave speeds in the HLL flux satisfy the condition (<ref>).
If the polynomial vectors {U_ij^n( x)} are locally divergence-free and satisfy the PP limiting condition (<ref>), then the scheme (<ref>) preserves U̅_ij^n+1∈𝒢 under the CFL-type condition

    Δ t_n≤ c_0min{Δ x ω_1^+/α_1,i+1/2^μ,
      Δ x ω_1^-/α_1,i-1/2^μ,
      Δ y ω_2^+/α_2,j+1/2^μ,
      Δ y ω_2^-/α_2,j-1/2^μ,
       }  ∀ I_ij∈Ω,  1≤μ≤ q,

where c_0 is the maximum allowable CFL number, and

    α_1,i+1/2^μ=α̂_1,i+1/2^μ,--
    𝒱_1,i+1/2^μ,-+
    |ℬ_1, i+1/2(y_j^(μ)) |/√(ρ_i+1/2,j^-,μ),
      α_1,i-1/2^μ=α̂_1,i-1/2^μ,++
    𝒱_1,i-1/2^μ,++
    |ℬ_1, i-1/2(y_j^(μ)) |/√(ρ_i-1/2,j^+,μ),


    α_2,j+1/2^μ=α̂_2,j+1/2^μ,--
    𝒱_2,j+1/2^μ,-+
    |ℬ_2, j+1/2(x_i^(μ)) |/√(ρ_i,j+1/2^μ,-),
      α_2,j-1/2^μ=α̂_2,j-1/2^μ,++
    𝒱_2,j-1/2^μ,++
    |ℬ_2, j-1/2(x_i^(μ)) |/√(ρ_i,j-1/2^μ,+).



We now present the proof of Theorem <ref>.

 Proof 
According to the identities (<ref>)-(<ref>) in Corollary <ref>, the scheme (<ref>) can be written as

    U̅_ij^n+1=U̅_ij^n-   Δ t_n/Δ x∑_μ=1^qω_μ^G[( 𝒱_1,i+1/2^μ,-H_1,i+1/2^μ+F_1(U_i+1/2,j^-,μ)-𝒱_1,i+1/2^μ,-U_i+1/2,j^-,μ).
        -( 𝒱_1,i-1/2^μ,+H_1,i-1/2^μ+F_1(U_i-1/2,j^+,μ)-𝒱_1,i-1/2^μ,+U_i-1/2,j^+,μ)
       .+ℬ_1, i+1/2(y_j^(μ))S(U_i+1/2,j^-,μ)-
    ℬ_1, i-1/2(y_j^(μ))S(U_i-1/2,j^+,μ)]
    
    -   Δ t_n/Δ y∑_μ=1^qω_μ^G[( 𝒱_2,j+1/2^μ,-H_2,j+1/2^μ+F_2(U_i,j+1/2^μ,-)-𝒱_2,j+1/2^μ,-U_i,j+1/2^μ,-).
        -( 𝒱_2,j-1/2^μ,+H_2,j-1/2^μ+F_2(U_i,j-1/2^μ,+)-𝒱_2,j-1/2^μ,+U_i,j-1/2^μ,+)
       .+ℬ_2, j+1/2(x_i^(μ))S(U_i,j+1/2^μ,-)-
    ℬ_2, j-1/2(x_i^(μ))S(U_i,j-1/2^μ,+)].

The scheme (<ref>) above can be reformulated as

    U̅_ij^n+1=U̅_ij^n+
    Ξ_1+Ξ_2+Ξ_3+Ξ_4,

with

    Ξ_1:=   Δ t_n/Δ x∑_μ=1^qω_μ^G[
    (𝒱_1,i+1/2^μ,--α̂_1,i+1/2^μ,-)
    U_i+1/2,j^-,μ-(𝒱_1,i-1/2^μ,++
    α̂_1,i-1/2^μ,+)
    U_i-1/2,j^+,μ]
       +Δ t_n/Δ y∑_μ=1^qω_μ^G[
    (𝒱_2,j+1/2^μ,--α̂_2,j+1/2^μ,-)
    U_i,j+1/2^μ,--(𝒱_2,j-1/2^μ,++
    α̂_2,j-1/2^μ,+)
    U_i,j-1/2^μ,+],


    Ξ_2:=   Δ t_n/Δ x∑_μ=1^qω_μ^G[
    (α̂_1,i+1/2^μ,-U_i+1/2,j^-,μ-
    F_1(U_i+1/2,j^-,μ))+
    (α̂_1,i-1/2^μ,+U_i-1/2,j^+,μ+
    F_1(U_i-1/2,j^+,μ))]
       +Δ t_n/Δ y∑_μ=1^qω_μ^G[
    (α̂_2,j+1/2^μ,-U_i,j+1/2^μ,--
    F_2(U_i,j+1/2^μ,-))+
    (α̂_2,j-1/2^μ,+U_i,j-1/2^μ,++
    F_2(U_i,j-1/2^μ,+))],


    Ξ_3:=   -Δ t_n/Δ x∑_μ=1^qω_μ^G(
    𝒱_1,i+1/2^μ,-H_1,i+1/2^μ-
    𝒱_1,i-1/2^μ,+H_1,i-1/2^μ)
    -Δ t_n/Δ y∑_μ=1^qω_μ^G(
    𝒱_2,j+1/2^μ,-H_2,j+1/2^μ-
    𝒱_2,j-1/2^μ,+H_2,j-1/2^μ),


    Ξ_4:=   -Δ t_n/Δ x∑_μ=1^qω_μ^G(ℬ_1, i+1/2(y_j^(μ))S(U_i+1/2,j^-,μ)-
    ℬ_1, i-1/2(y_j^(μ))S(U_i-1/2,j^+,μ))
       -Δ t_n/Δ y∑_μ=1^qω_μ^G(ℬ_2, j+1/2(x_i^(μ))S(U_i,j+1/2^μ,-)-
    ℬ_2, j-1/2(x_i^(μ))S(U_i,j-1/2^μ,+)).

We define

    U̅_ij^(μ)=1/Δ y(α̂_1,i+1/2^μ,-
    +α̂_1,i-1/2^μ,+)+Δ x(α̂_2,j+1/2^μ,-+α̂_2,j-1/2^μ,+)[Δ y (α̂_1,i+1/2^μ,-U_i+1/2,j^-,μ-
    F_1(U_i+1/2,j^-,μ)).
       +Δ y(α̂_1,i-1/2^μ,+U_i-1/2,j^+,μ+
    F_1(U_i-1/2,j^+,μ))
    +Δ x(α̂_2,j+1/2^μ,-U_i,j+1/2^μ,--
    F_2(U_i,j+1/2^μ,-))
       +.Δ x(α̂_2,j-1/2^μ,+U_i,j-1/2^μ,++
    F_2(U_i,j-1/2^μ,+))],   μ=1, 2, ⋯, q,

then Ξ_2 is equal to

    Ξ_2=Δ t_n/Δ x Δ y∑_μ=1^qω_μ^G(Δ y(α̂_1,i+1/2^μ,-
    +α̂_1,i-1/2^μ,+)+Δ x(α̂_2,j+1/2^μ,-+α̂_2,j-1/2^μ,+)
    )U̅_ij^(μ).


Now, we first show ρ̅_ij^n+1>0.
Thanks to Corollary <ref>, we have, for all 1≤μ≤ q, U̅_ij^(μ)∈𝒢_ρ, and Δ y(α̂_1,i+1/2^μ,-
+α̂_1,i-1/2^μ,+)+Δ x(α̂_2,j+1/2^μ,-+α̂_2,j-1/2^μ,+)>0.
Therefore, Ξ_2∈𝒢_ρ.
According to Corollary <ref>, we have H_1,i+1/2^μ, H_1,i-1/2^μ, H_2,j+1/2^μ and H_2,j-1/2^μ all belong to the 𝒢_ρ,
which means Ξ_3∈𝒢_ρ.
Since the first component of Ξ_4 is zero, we deduce from (<ref>) that

    ρ̅_ij^n+1   >ρ̅_ij^n+Δ t_n/Δ x∑_μ=1^qω_μ^G[
    (𝒱_1,i+1/2^μ,--α̂_1,i+1/2^μ,-)
    ρ_i+1/2,j^-,μ-(𝒱_1,i-1/2^μ,++
    α̂_1,i-1/2^μ,+)
    ρ_i-1/2,j^+,μ]
          +Δ t_n/Δ y∑_μ=1^qω_μ^G[
    (𝒱_2,j+1/2^μ,--α̂_2,j+1/2^μ,-)
    ρ_i,j+1/2^μ,--(𝒱_2,j-1/2^μ,++
    α̂_2,j-1/2^μ,+)
    ρ_i,j-1/2^μ,+]
       =∑_μ=1^q ω_μ^G[ω_1^+ρ_i+1/2,j^-,μ+ω_1^-ρ_i-1/2,j^+,μ+ω_2^+ρ_i,j+1/2^μ,-+ω_2^-ρ_i,j-1/2^μ,+]
    +∑_s=1^S ω_s ρ_σ^n(x_ij^(s), y_ij^(s))
         +Δ t_n/Δ x∑_μ=1^qω_μ^G[
    (𝒱_1,i+1/2^μ,--α̂_1,i+1/2^μ,-)
    ρ_i+1/2,j^-,μ-(𝒱_1,i-1/2^μ,++
    α̂_1,i-1/2^μ,+)
    ρ_i-1/2,j^+,μ]
         +Δ t_n/Δ y∑_μ=1^qω_μ^G[
    (𝒱_2,j+1/2^μ,--α̂_2,j+1/2^μ,-)
    ρ_i,j+1/2^μ,--(𝒱_2,j-1/2^μ,++
    α̂_2,j-1/2^μ,+)
    ρ_i,j-1/2^μ,+],
       ≥∑_μ=1^qω_μ^Gρ_i+1/2,j^-,μ(ω_1^++
    Δ t_n/Δ x(𝒱_1,i+1/2^μ,--α̂_1,i+1/2^μ,-) )
         +∑_μ=1^qω_μ^Gρ_i-1/2,j^+,μ(ω_1^--Δ t_n/Δ x(𝒱_1,i-1/2^μ,++
    α̂_1,i-1/2^μ,+))
         +∑_μ=1^qω_μ^Gρ_i,j+1/2^μ,-(ω_2^++Δ t_n/Δ y(𝒱_2,j+1/2^μ,--
    α̂_2,j+1/2^μ,-))
         +∑_μ=1^qω_μ^Gρ_i,j-1/2^μ,+(ω_2^--Δ t_n/Δ y(𝒱_2,j-1/2^μ,++
    α̂_2,j-1/2^μ,+))
       ≥0,

where we have used the convex decomposition (<ref>) in the above equality, the PP limiting condition (<ref>) in the first inequality, and the CFL-type condition (<ref>) in the last inequality.

we then prove that U̅_ij^n+1·n^*+|B^*|^2/2>0 for any v^*, B^*∈ℝ^3.
It follows from (<ref>) that

    U̅_ij^n+1·n^*+|B^*|^2/2=Π_0+Π_1+Π_2+Π_3+Π_4,

where

    Π_0:=U̅_ij^n·n^*+|B^*|^2/2,


    Π_1:=   Δ t_n/Δ x∑_μ=1^qω_μ^G[
    (𝒱_1,i+1/2^μ,--α̂_1,i+1/2^μ,-)
    (U_i+1/2,j^-,μ·n^*+|B^*|^2/2).
       .-(𝒱_1,i-1/2^μ,++
    α̂_1,i-1/2^μ,+)
    (U_i-1/2,j^+,μ·n^*+|B^*|^2/2)
    ]
       +Δ t_n/Δ y∑_μ=1^qω_μ^G[
    (𝒱_2,j+1/2^μ,--α̂_2,j+1/2^μ,-)
    (U_i,j+1/2^μ,-·n^*+|B^*|^2/2).
       .-(𝒱_2,j-1/2^μ,++
    α̂_2,j-1/2^μ,+)
    (U_i,j-1/2^μ,+·n^*+|B^*|^2/2)
    ],


    Π_2:=Δ t_n/Δ x Δ y∑_μ=1^qω_μ^G(Δ y(α̂_1,i+1/2^μ,-
    +α̂_1,i-1/2^μ,+)+Δ x(α̂_2,j+1/2^μ,-+α̂_2,j-1/2^μ,+)
    )(U̅_ij^(μ)·n^*+|B^*|^2/2),


    Π_3:=   -Δ t_n/Δ x∑_μ=1^qω_μ^G[
    𝒱_1,i+1/2^μ,-(H_1,i+1/2^μ·n^*+|B^*|^2/2)
    -𝒱_1,i-1/2^μ,+(H_1,i-1/2^μ·n^*+|B^*|^2/2)]
       -Δ t_n/Δ y∑_μ=1^qω_μ^G[
    𝒱_2,j+1/2^μ,-(H_2,j+1/2^μ·n^*+|B^*|^2/2)
    -𝒱_2,j-1/2^μ,+(H_2,j-1/2^μ·n^*+|B^*|^2/2)],


    Π_4:=   -Δ t_n/Δ x∑_μ=1^qω_μ^G(ℬ_1, i+1/2(y_j^(μ))S(U_i+1/2,j^-,μ)·n^*-
    ℬ_1, i-1/2(y_j^(μ))S(U_i-1/2,j^+,μ)·n^*
    )
       -Δ t_n/Δ y∑_μ=1^qω_μ^G(ℬ_2, j+1/2(x_i^(μ))S(U_i,j+1/2^μ,-)·n^*-
    ℬ_2, j-1/2(x_i^(μ))S(U_i,j-1/2^μ,+)
    ·n^*).

We now estimate the lower bounds of Π_0, Π_2, Π_3 and Π_4, respectively.
Based on the the convex decomposition (<ref>), we have

    Π_0=   ∑_μ=1^q ω_μ^G[ω_1^+(U_i+1/2,j^-,μ·n^*+|B^*|^2/2)+ω_1^-(U_i-1/2,j^+,μ·n^*+|B^*|^2/2)+ω_2^+(U_i,j+1/2^μ,-·n^*+|B^*|^2/2).
       .+ω_2^-(U_i,j-1/2^μ,+·n^*+|B^*|^2/2)]+∑_s=1^S ω_s (U_σ^n(x_ij^(s), y_ij^(s))·n^*+|B^*|^2/2),
       ≥∑_μ=1^q ω_μ^G[ω_1^+(U_i+1/2,j^-,μ·n^*+|B^*|^2/2)+ω_1^-(U_i-1/2,j^+,μ·n^*+|B^*|^2/2).
         .+ω_2^+(U_i,j+1/2^μ,-·n^*+|B^*|^2/2)
    +ω_2^-(U_i,j-1/2^μ,+·n^*+|B^*|^2/2)],

where the inequality follows from the PP limiting condition (<ref>) and the Lemma <ref>.
According to Corollary <ref>, ∀v^*,B^*∈ℝ^3, we have

    U̅_ij^(μ)·n^*+|B^*|^2/2≥-v^*·B^*(
    Δ y (B_1)_i+1/2,j^-,μ-Δ y (B_1)_i-1/2,j^+,μ+
    Δ x (B_2)_i,j+1/2^μ,--Δ x (B_2)_i,j-1/2^μ,+)/Δ y(α̂_1,i+1/2^μ,-
    +α̂_1,i-1/2^μ,+)+Δ x(α̂_2,j+1/2^μ,-+α̂_2,j-1/2^μ,+).

Therefore, the lower bound of Π_2 can be given by

    Π_2   ≥ -Δ t_n/Δ x Δ y(v^*·B^*)∑_μ=1^qω_μ^G(Δ y (B_1)_i+1/2,j^-,μ-Δ y (B_1)_i-1/2,j^+,μ+
    Δ x (B_2)_i,j+1/2^μ,--Δ x (B_2)_i,j-1/2^μ,+)
       =-Δ t_n(v^*·B^*)div_withinB_σ^n
    =0,

where the last equality uses the LDF property of the magnetic field.
By using the inequality (<ref>) in Corollary <ref>, we give a lower bound of Π_3 as

    Π_3   ≥ -Δ t_n/Δ x∑_μ=1^qω_μ^G[
      𝒱_1,i+1/2^μ,-(-v^*·B^*/𝒱_1,i+1/2^μ,+-𝒱_1,i+1/2^μ,-((B_1)_i+1/2,j^+,μ-
      (B_1)_i+1/2,j^-,μ)).
       .-𝒱_1,i-1/2^μ,+(-v^*·B^*/𝒱_1,i-1/2^μ,+-𝒱_1,i-1/2^μ,-((B_1)_i-1/2,j^+,μ-
      (B_1)_i-1/2,j^-,μ))]
       -Δ t_n/Δ y∑_μ=1^qω_μ^G[
      𝒱_2,j+1/2^μ,-(-v^*·B^*/𝒱_2,j+1/2^μ,+-𝒱_2,j+1/2^μ,-((B_2)_i,j+1/2^μ,+-
      (B_2)_i,j+1/2^μ,-)).
       -.𝒱_2,j-1/2^μ,+(-v^*·B^*/𝒱_2,j-1/2^μ,+-𝒱_2,j-1/2^μ,-((B_2)_i,j-1/2^μ,+-
      (B_2)_i,j-1/2^μ,-))].

Following the inequality (<ref>) in Lemma <ref>, we estimate the lower bound of Π_4 as

    Π_4   ≥Δ t_n/Δ x∑_μ=1^qω_μ^G(
    ℬ_1, i+1/2(y_j^(μ))(v^*·B^*)-
    |ℬ_1, i+1/2(y_j^(μ)) |/√(ρ_i+1/2,j^-,μ)( U_i+1/2,j^-,μ·n^*+
    |B^*|^2/2))
       +Δ t_n/Δ x∑_μ=1^qω_μ^G(
    -ℬ_1, i-1/2(y_j^(μ))(v^*·B^*)-
    |ℬ_1, i-1/2(y_j^(μ)) |/√(ρ_i-1/2,j^+,μ)( U_i-1/2,j^+,μ·n^*+
    |B^*|^2/2))
       +Δ t_n/Δ y∑_μ=1^qω_μ^G(
    ℬ_2, j+1/2(x_i^(μ))(v^*·B^*)-
    |ℬ_2, j+1/2(x_i^(μ)) |/√(ρ_i,j+1/2^μ,-)( U_i,j+1/2^μ,-·n^*+
    |B^*|^2/2))
       +Δ t_n/Δ y∑_μ=1^qω_μ^G(
    -ℬ_2, j-1/2(x_i^(μ))(v^*·B^*)-
    |ℬ_2, j-1/2(x_i^(μ)) |/√(ρ_i,j-1/2^μ,+)( U_i,j-1/2^μ,+·n^*+
    |B^*|^2/2)).

Combining the lower bounds in (<ref>), (<ref>), (<ref>), (<ref>) with (<ref>), we obtain

    U̅_ij^n+1·n^*+
    |B^*|^2/2
       ≥∑_μ=1^qω_μ^G(
    ω_1^++Δ t_n/Δ x(
    𝒱_1,i+1/2^μ,--α̂_1,i+1/2^μ,--
    |ℬ_1, i+1/2(y_j^(μ)) |/√(ρ_i+1/2,j^-,μ)))( U_i+1/2,j^-,μ·n^*+
    |B^*|^2/2)
       +∑_μ=1^qω_μ^G(
    ω_1^--Δ t_n/Δ x(
    𝒱_1,i-1/2^μ,++α̂_1,i-1/2^μ,++
    |ℬ_1, i-1/2(y_j^(μ)) |/√(ρ_i-1/2,j^+,μ)))( U_i-1/2,j^+,μ·n^*+
    |B^*|^2/2)
       +∑_μ=1^qω_μ^G(
    ω_2^++Δ t_n/Δ y(
    𝒱_2,j+1/2^μ,--α̂_2,j+1/2^μ,--
    |ℬ_2, j+1/2(x_i^(μ)) |/√(ρ_i,j+1/2^μ,-)))( U_i,j+1/2^μ,-·n^*+
    |B^*|^2/2)
       +∑_μ=1^qω_μ^G(
    ω_2^--Δ t_n/Δ y(
    𝒱_2,j-1/2^μ,++α̂_2,j-1/2^μ,++
    |ℬ_2, j-1/2(x_i^(μ)) |/√(ρ_i,j-1/2^μ,+)))( U_i,j-1/2^μ,+·n^*+
    |B^*|^2/2)
       >0

where the last inequality uses the CFL condition (<ref>) and PP limiting condition (<ref>).
Therefore, we have

    ρ̅_ij^n+1>0,  U̅_ij^n+1·n^*+
    |B^*|^2/2>0, ∀n^*, B^* ∈ℝ^3,

which implies U̅_ij^n+1∈𝒢 by Lemma <ref>.
The proof is completed. 

As direct consequences of Theorem <ref>, we have the following two corollaries.

(PP via Zhang-Shu convex decomposition).
Let the wave speeds in the HLL flux satisfy the condition (<ref>).
If the polynomial vectors {U_ij^n( x)} are locally divergence-free and satisfy the following PP limiting condition,

    U_ij^n(x)∈𝒢,
      ∀x∈𝕊_ij^Zhang-Shu,

where 𝕊_ij^Zhang-Shu denotes the set of all nodes of Zhang-Shu convex decomposition (<ref>),
then the scheme (<ref>) preserves U̅_ij^n+1∈𝒢 under the CFL-type condition

    Δ t_n≤ c_0min{a_1Δ xΔ y ω̂_1^GL/(a_2Δ x+a_1Δ y)max{α_1,i+1/2^μ, α_1,i-1/2^μ},
      a_2Δ xΔ y ω̂_1^GL/(a_2Δ x+a_1Δ y)max{α_2,j+1/2^μ, α_2,j-1/2^μ}},  
         ∀ I_ij∈Ω,  1≤μ≤ q,

where ω̂_1^GL=1/L(L-1) with L=⌈k+3/2⌉.



(PP via optimal convex decomposition).
Let the wave speeds in the HLL flux satisfy the condition (<ref>).
If the polynomial vectors {U_ij^n( x)} are locally divergence-free and satisfy the following PP limiting condition,

    U_ij^n(x)∈𝒢,
      ∀x∈𝕊_ij^optimal,

where 𝕊_ij^optimal denotes the set of all nodes of optimal convex decomposition (<ref>),
then the scheme (<ref>) preserves U̅_ij^n+1∈𝒢 under the CFL-type condition

    Δ t_n≤ c_0min{a_1Δ x Δ y/max{ 6a_1Δ y+2a_2Δ x,  6a_2Δ x+2a_1Δ y}max{α_1,i+1/2^μ, α_1,i-1/2^μ},.
       .a_2Δ x Δ y/max{ 6a_1Δ y+2a_2Δ x,  6a_2Δ x+2a_1Δ y}max{α_2,j+1/2^μ, α_2,j-1/2^μ}},
         ∀ I_ij∈Ω,  1≤μ≤ q.



In the above discussions, we restrict ourselves to the first-order forward
Euler time discretization.
Since an SSP method is a convex combination of the forward Euler method, our PP
analysis of the proposed scheme remains valid according to the convexity of 𝒢.



§ NUMERICAL TESTS


We present some numerical results of proposed schemes for several 1D and 2D MHD problems, to test the accuracy, the ability of capturing shocks, and the robustness in the extreme medium with strong magnetic field and huge discontinuities.
The tests below are performed on the uniform 1D meshes and 2D rectangular meshes.
We use the HLL flux for all tests.
Without loss of generality, we focus on the polynomials of degree k=2 for DG solutions with the third-order SSP Runge-Kutta time discretization.
The CFL number is set as 0.12.



 §.§ 1D MHD tests


This section tests the proposed OEDG scheme through several 1D MHD examples, including a smooth sine wave problem, three shock tube problems, and the MHD version of the Leblanc problem.



  §.§.§ 1D smooth sine wave problem



We first simulate a smooth MHD problem <cit.> to validate the accuracy of the PP OEDG scheme in 1D.
This problem describes a sine wave propagating with low density, and its exact solution is given by

    (ρ, u, p, B)(x, t)=
    (1+0.99sin(x-t), 1, 0, 0, 1, 0.1, 0, 0),   
    x∈[0, 2π], t≥ 0.

Periodic boundary conditions are used, and γ=1.4.
Table <ref> lists the numerical errors and the corresponding
convergence rates at t=0.1 for the density.
The results show that the 1D ℙ^2-based PP OEDG scheme achieves the expected third-order convergence rate.








  §.§.§ Shock tube problems


Next, we simulate three 1D MHD shock tube problems to examine the shock-capturing and non-oscillatory properties of the PP OEDG method.
The adiabatic index γ is set as 5/3 for the first two shock tube problems and as 2 for the third one.
We use the third-order PP OEDG scheme with 800 uniform cells for all the simulations.
The initial conditions of the first shock tube problem <cit.>

 are  given by

    (ρ, u, B, p)(x, 0)=
    {   (1.08, 1.2, 0.01, 0.5, 2/√(4π), 3.6/√(4π), 2/√(4π), 0.95),     x<0.5,
       (1, 0, 0, 0, 2/√(4π), 4/√(4π), 2/√(4π), 1),     x>0.5,
    .

and the numerical results at time t=0.2 are presented in Figure <ref>.

The second shock tube problem includes a hydrodynamic rarefaction, a switch-on slow shock,
a contact discontinuity, a slow shock, a rotational discontinuity, and a fast rarefaction.
The initial setups are the same as in <cit.>:

    (ρ, u, B, p)(x, 0)=
    {   (1, 0, 0, 0, 0.7, 0, 0, 1),     x<0.5,
       (0.3, 0, 0, 1, 0.7, 1, 0, 0.2),     x>0.5,
    .

and
Figure <ref> shows the numerical results at time t=0.16.
The third shock tube problem was first suggested by Brio and Wu <cit.> with the initial data

    (ρ, u, B, p)(x, 0)=
    {   (1, 0, 0, 0, 0.75, 1, 0, 1),     x<0,
       (0.125, 0, 0, 0, 0.75, -1, 0, 0.1),     x>0,
    .

and the numerical solution at t=0.1 is displayed in Figure <ref>.
From these results, one can see that the PP OEDG method correctly captures
all the discontinuities such as shocks with high resolution and  without any noticeable spurious oscillations.















  §.§.§ Leblanc problem


In order to test the robustness of the 1D PP OEDG scheme, we run a challenging MHD Riemann problem with a strong magnetic field, a large jump of the initial pressure, and extremely low plasma-beta (β=4×10^-8).
This problem
is a variant of the Leblanc problem <cit.> of the gas
dynamics, and was proposed in <cit.>.
The initial setting is as follows:

    (ρ, u, p, B)(x, 0)=
    {   (2, 0, 0, 0, 10^9, 0, 5000, 5000),     x<0,
       (0.001, 0, 0, 0, 1, 0, 5000, 5000),     x>0.
    .

The computational domain is set to [-10, 10], and the adiabatic index γ is taken as 1.4.
Figure <ref> displays the numerical results of density and magnetic pressure at t=0.00003.
As expected,
negative density or pressure does not appear in the cell averages of solution, indicating that the 1D PP OEDG scheme exhibits good stability.









 §.§ 2D MHD tests


This section conducts six 2D MHD tests using ℙ^2–based PP LDF OEDG method, including a smooth sine wave problem, the Orszag–Tang problem, rotor problem, shock cloud interaction, MHD blast wave, and MHD jet flows.




  §.§.§ 2D smooth sine wave problem


To verify the convergence order of the PP LDF OEDG method,
we simulate a 2D smooth sine wave problem <cit.>, which describes a MHD sine wave propagating with low density.
The exact solution of this problem is given by

    (ρ, u, p, B)(x, y, t)=
    (1+0.99sin(x+y-2t), 1, 1, 0, 1, 0.1, 0.1, 0),   
    x∈[0, 2π], y∈[0, 2π], t≥ 0,

and the adiabatic index γ=1.4.
Table <ref> provides the numerical errors and corresponding
convergence rates of density at time t=0.1.
The results show that the ℙ^2–based PP LDF OEDG scheme achieves the third-order accuracy overall.







  §.§.§ Orszag–Tang problem


We now consider a widely-used Orszag–Tang vortex problem <cit.>. For this problem, the vortex starts from a smooth initial field configuration and quickly forms shocks before transitioning to 2D MHD turbulence.
The initial conditions are as follows:

    [ ρ(x, y)=γ^2,    u_1(x, y)=-sin y,   
      u_2(x, y)= sin x,     u_3(x, y)=0,;   p(x, y)=γ,    B_1(x, y)=-sin y,   
      B_2(x, y)=sin(2x),    B_3(x, y)=0, ]

where the adiabatic index γ is taken as 5/3.
The computational domain is set as [0, 2π]×[0, 2π], and the periodic boundary conditions are used on all boundaries.

Figure <ref> shows the density and thermal pressure obtained by third-order PP LDF OEDG scheme at 400×400 resolution.
We can clearly see both sharp shocks and smooth flows, which are
in good agreement with those in <cit.>.
For the Orszag–Tang vortex problem, Jiang and Wu <cit.> previously reported that the negative pressure was occurred at t≈3.9 in their test.
In our implementation for third-order PP LDF OEDG scheme, negative pressure was not observed.
To further verify the convergence of the scheme, Figure <ref> plots the
profiles of density and thermal pressure along the line of y=0.625π at time t=3 under 200×200 and 400×400 resolutions.
In Figure <ref>, the
shock discontinuities for density and thermal pressure are formed
near x=0.4, 0.5, 1.6 and 4.4, which is consistent with the results in <cit.>.
According to the comparison of second-order and fourth-order schemes for pressure profiles
by Yang et al. <cit.>,
the high accuracy scheme results in sharper shock discontinuities at x=0.4 and x=0.5.
In Figure <ref>, it is observed that the shock discontinuities around x=0.4 and x=0.5 with 400×400 resolution are much sharper than those with 200×200 resolution, verifying the convergence of our scheme.











  §.§.§ Rotor problem


We next consider the rotor problem <cit.> to further test the performance of ℙ^2–based PP LDF OEDG method.
In this rotor problem, a dense disk of fluid is placed in a static fluid background, with a gradual velocity tapering layer between the dense disk edge and the static ambient fluid, and the disk will rotate within ambient fluid over time.
The initial setup is given by

    (u_3, B_1, B_2, B_3, p)=(0,2.5 / √(4 π), 0,0,0.5),

and

    (ρ, u_1, u_2)=
    {   (10, -(y-0.5) / r_0, (x-0.5) / r_0)     if  r<r_0, 
       (1+9f, -f(y-0.5) / r, f(x-0.5) / r)     if  r_0<r<r_1, 
       (1, 0, 0)     if  r>r_1,
    .

where r_0=0.1, r_1=0.115, r = √((x-0.5)^2+(y-0.5)^2), and f=(r_1-r)/(r_1-r_0).
The adiabatic index γ is taken as 5/3.
The computational domain is [0, 1]×[0, 1],
and periodic boundary conditions are used.

In Figure <ref>, the contours of density, thermal pressure, Mach number and magnetic pressure are presented at t=0.295 for the third-order PP LDF OEDG method with 400×400 cells.
We can see that the central contours show good conservation of the circular rotation pattern,
and such performance was found by Tóth <cit.> to be challenging for some other MHD numerical schemes.
According to Li and Shu <cit.>, the numerical solution around central part of the Mach number existed some distortion, which was attributed to the divergence error in the magnetic field.
Figure <ref> zooms in the central part of the Mach number at different resolutions.
No distortions are observed at the outskirts
of the central almond-shaped disk region, which suggests that the
divergence error in the magnetic field is well controlled by our LDF OEDG scheme.
To further examine the performance of our scheme, Figure <ref> presents two slices of the magnetic field at time t=0.295 across horizontal and vertical cuts through the centre of the box with different mesh resolutions.
The results show good agreement with previous studies <cit.>.
It was found that in <cit.> higher-order scheme exhibits higher peaks near x=0.35, 0.65 of B_1 and x=0.325, 0.675 of B_2, as well as sharper shock discontinuities around x=0.4, 0.6 of B_1.
In Figure <ref>, the same trend is shown in the high-resolution results compared to the low-resolution results, which reflects the convergence of our scheme.
According to Liu et al. <cit.>,
the overshoots at y=0.05, 0.95 for B_1 and at x=0.38, 0.62 for B_2
were generated by their third-order finite volume scheme.
Figure <ref> shows that no overshoot is observed,
verifying that the PP LDF OEDG scheme can suppress oscillations well.













  §.§.§ Shock cloud interaction


To further validate the shock-capturing and PP performance
of our scheme, we consider shock cloud interaction <cit.>.
This example involves strong MHD shocks interacting with a disruption of dense cloud.
The initial conditions are given by
the left and right sates separated by a discontinuity along the line x=0.6:

    (ρ, 𝐮, p, B)=
    {   (3.86859,0,0,0,167.345,0,2.1826182,-2.1826182),     x<0.6, 
       (1,-11.2536,0,0,1,0,0.56418958,0.56418958),     x>0.6,
    .

and a circular cloud with a radius of 0.15 centered on (0.8, 0.5) is placed on the right side of the discontinuity.
The cloud has the same setups as the ambient fluid except for a higher density 10.
Our computational domain is [0, 1]×[0,1], and adiabatic index γ is 5/3.
For this problem, supersonic inflow condition is employed on the right boundary, and outflow conditions are utilized on the others.



Figure <ref> displays numerical results of density, thermal pressure, velocity, and magnetic pressure at time t=0.06 with 400×400 cells for our third-order PP LDF OEDG scheme.
One can see that the complex flow structures and discontinuities
are well captured with high resolution, and the results show good agreement with those in the literature, e.g., <cit.>.
We found that once the OE procedure is turned off in this test, the code will break down and cause simulation failure, which proves the stability of our scheme.
In addition, we also observed that negative pressure will appear in the cell average of numerical solutions if we remove either the Godunov–Powell source term or the PP limiter. This further confirms the importance of the discretized Godunov–Powell
source term and the PP limiter for preserving positivity.





  §.§.§ MHD blast wave



We now run the blast wave problem to test the ability of our scheme in cases of strong shocks and rarefaction wave.
In the blast wave problem, the over-pressured region, located at the center of strongly magnetized medium with low plasma β, will generate an MHD blast wave, which drives the fast outward magnetic shocks to compress the plasma and magnetic field ahead, and generates rarefaction wave internally.
For this test, we adopted the version of Stone et al. <cit.> with the initial settings as follows:

    [ ρ=1,    u_1=u_2=u_3=0,    B_1=B_2=1 / √(2),    B_3=0,   
     p={   10,     if √((x-1.5)^2+y^2)<0.1, 
       0.1,     otherwise .
    . ]

We use γ=5/3 and fix the boundary conditions.
The simulation is run on the domain [1, 2]×[-0.5, 0.5] until the final time of t=0.2.

In Figure <ref>, we show the density, thermal pressure, velocity, and magnetic pressure for the third-order PP LDF OEDG scheme at time t=0.2 with 400×400 cells.
The results are consistent with those in <cit.>.
One can see that a circular blast wave propagates outward, and a rarefaction wave propagates inward.
The symmetry of the problem is well preserved.
Our scheme captures sharp shocks and remains smooth post-shock regions, which demonstrates the robustness of our method in highly magnetized shock configurations.






  §.§.§ MHD jet flows



The last example considers three challenging MHD jet problems, to further test the performance of our scheme in capturing shocks and maintaining positivity.
MHD jet flows were proposed in Wu and Shu <cit.> by adding a magnetic field in the gas dynamical jet of Balsara <cit.>.
Since there are strong shock waves, shear flows and interface instabilities in high Mach number jet with strong magnetic field and
huge kinetic energy, negative pressure is very likely to be produced in the numerical simulations.

Our computational domain is [-0.5, 0.5]×[0, 1.5], and the adiabatic index γ=1.4.
Initially, we set (ρ, p)=(0.1γ, 1) in the static ambient medium with a magnetic field (0, B_2, 0).
As observed in <cit.>, the larger B_2 takes, the more challenging the test becomes.
On the bottom boundary, a Mach 800 dense jet with (ρ, p, u_1, u_2, u_3)=(γ, 1, 0, 800, 0) is injected along the nozzle {y=0, |x|<0.05}.
We use the fixed inflow condition on the nozzle, and outflow conditions on the other boundary conditions.
For this problem, we set the computational domain as [0, 0.5]×[0, 1.5], i.e., take the right half of the whole domain, and utilize the reflecting boundary condition at x=0.
The domain [0, 0.5]×[0, 1.5] is divided into 200×600 cells.
We first test our scheme with two magnetized cases: (i) B_2=√(200); (ii) B_2=√(2000).
The numerical results of density logarithm log_10ρ, thermal pressure logarithm log_10p, and velocity |u| for two cases are displayed in Figure <ref> and Figure <ref>, respectively.
These numerical results clearly show the time evolution of flow structures, and the results agree with those in <cit.>.
Note that our scheme does a great job of capturing the Mach shock wave at the jet head and the beam/cocoon interface with high resolution.
Specifically, the shear flows near the beam/cocoon interfaces of jets are clearly presented.
In all the simulations, no negative pressure or negative density was observed.
We also run our code on the test of higher Mach 10000 dense jet with stronger magnetized case of B_2=√(20000).
Figure <ref> shows the time resolution of flow structures, including density logarithm log_10ρ, thermal pressure logarithm log_10p, and velocity |u|.
Our scheme exhibits good performance and robustness in such an extreme test. 
We found that if we turn off any of OE procedure, PP limiting procedure, and suitably discretized Godunov-Powell source term, the code will not be able to run jet tests with strong magnetic fields due to negative pressure produced in the cell averages of the numerical solution.
This, again, demonstrates the importance of OE procedure, PP limiting procedure, and suitably discretized Godunov-Powell source term for PP property of the scheme.















§ CONCLUDING REMARKS



In this paper, we extend the OEDG method to ideal MHD equations, and
propose an LDF OEDG scheme by using the LDF DG scheme as the basis scheme.
In the LDF OEDG scheme, the OE procedures are performed on the Legendre basis and LDF basis respectively, as a result, the proposed scheme can effectively control the spurious oscillations while automatically maintaining the LDF property of the computed magnetic field.
By alternately evolving each Runge-Kutta stage and an OE procedure, the implementation of the LDF OEDG scheme is effective and computational cheap on general meshes.
Moreover, based on the PP limiter, HLL flux and proper discretized Godunov-Powell source term added to the evolution equations of the cell averages, we build a framework of PP LDF OEDG scheme on Cartesian meshes, and rigorously analyze the PP property of the proposed scheme.
We derive the PP CFL condition on Cartesian meshes via the general convex decomposition of the cell averages, and also list the PP CFL conditions corresponding to the Zhang-Shu convex decomposition and optimal convex decomposition. 
Using third-order SSP Runge-Kutta time discretization, we test several MHD examples in one and two dimensions. 
The numerical results confirm that the proposed scheme achieves the expected optimal  convergence order and behaves robust even in some tough MHD problems with low density or low pressure.
Our next work is to conduct MHD numerical simulations based on the PP LDF OEDG schemes on unstructured meshes.








§ ACKNOWLEDGMENTS



This work was partially supported by National Natural Science Foundation of China (Grant No. 12171227) and Shenzhen Science and Technology Program (Grant No. RCJC20221008092757098).








