High-order conservative finite difference GLM-MHD schemes for cell-centered MHD

1. Introduction

• The study assumes an ideal-gas equation of state for the compressible MHD system.The paper describes the system as compressible MHD in multiple spatial dimensions.Source: §1 · p. 2
• Traditional second-order TVD schemes can be sharp at discontinuities but lose accuracy near local extrema and excessively dissipate smooth flow and contact waves.The stated causes include reduction to first-order near extrema and discretization errors.Source: §1 · p. 2
• High-order finite-volume methods require multidimensional reconstruction and multiple face Riemann problems, whereas finite differences use one-dimensional flux-derivative stencils.Finite-volume methods are better suited to non-uniform grids and adaptive mesh hierarchies.Source: §1 · p. 3
• The proposed cell-centered finite-difference schemes use third- and fifth-order reconstructions and compare their performance on two- and three-dimensional problems.The selected methods include improved WENO and non-polynomial third-order reconstructions, plus WENO-Z and MP5 at fifth order.Source: §1 · p. 3
• GLM divergence cleaning avoids elliptic-cleaning costs and staggered-field complexity while retaining a conservative formulation for mass, momentum, magnetic induction, and energy.Placing primary variables at the same spatial position also eases adding more complex physics.Source: §1 · p. 3

2. The Constrained GLM-MHD Equations

• The GLM-MHD system adds a scalar generalized Lagrange multiplier to the cell-centered conservative MHD equations to control magnetic divergence.The approach modifies Gauss’s and Faraday’s laws of magnetism.Source: §2 · p. 3
• Divergence errors propagate to domain boundaries at finite speed and are damped, while the equations remain hyperbolic and conservative except for the scalar field’s source term.The scalar field ψ satisfies a non-homogeneous equation with a source term.Source: §2 · p. 4
• The GLM modes ±ch are decoupled from the seven ordinary MHD waves and carry jumps in the longitudinal magnetic field and ψ.The GLM waves are assigned the maximum signal speed compatible with the time step.Source: §2 · p. 4
• This decoupling lets the longitudinal field and ψ be solved as a separate 2 × 2 linear system, while standard characteristic decomposition handles the ordinary MHD equations.The solution process separates a 2 × 2 system from the other seven-wave equations.Source: §2 · p. 4
• The source integration is performed analytically, and divergence errors are minimized when αp lies in [0, 1], with smooth flows less sensitive to its value.The minimum grid size is Δh = min(Δx, Δy, Δz).Source: §2 · p. 4

3. Finite Difference schemes

• Cell-centered conservative finite differences evolve point values, reconstruct stable interface fluxes, and use third-order Runge–Kutta time integration, which limits overall accuracy to third order.The global Lax–Friedrichs flux splitting is particularly diffusive, though its extra dissipation becomes less important for higher-order schemes.Source: §3 · pp. 4–5

3.1. Modification for the Constrained GLM-MHD equations

• The GLM-MHD treatment separates the longitudinal magnetic field and multiplier from the seven MHD equations, avoiding a full 9 × 9 characteristic decomposition.The decoupled two-variable system is solved at interfaces before the ordinary 7 × 7 MHD equations are advanced.Source: §3.1 · p. 6

3.2. Third and Fifth-order Accurate Reconstructions

• Four reconstructions are compared: WENO+3 and LimO3 at third order, and WENO-Z and MP5 at fifth order.They include methods using conditional statements (LimO3 and MP5) and methods that avoid them (WENO+3 and WENO-Z).Source: §3.2 · pp. 6–7

3.2.1. Third-Order Improved WENO (WENO+3)

• WENO+3 sets ε = Δx^2 and improves accuracy over the original third-order scheme in smooth regions, while providing essentially non-oscillatory solutions near strong discontinuities and unresolved features.Avoiding loss of accuracy at critical points requires ε = O(Δx^2).Source: §3.2.1 · p. 7

3.2.2. Third-Order Limited reconstruction (LimO3 )

• LimO3 uses a local indicator to switch between limited and unlimited reconstruction, preserving accuracy at smooth extrema while distinguishing them from shallow gradients.Its parameter r controls that distinction; the tests use r = 1.Source: §3.2.2 · pp. 7–8

3.2.3. Fifth-Order Improved WENO: WENO-Z

• WENO-Z uses higher-order information about solution regularity to improve convergence at critical points and reduce dissipation at discontinuities while retaining essentially non-oscillatory behavior.The scheme is an improved fifth-order WENO formulation.Source: §3.2.3 · p. 8

3.2.4. Fifth-order Monotonicity Preserving (MP5)

• MP5 limits a fifth-order interface value to preserve monotonicity near discontinuities and accuracy in smooth regions without first-order degeneration near smooth extrema.Its curvature measure is heuristic and is chosen to reduce room for local extrema to develop.Source: §3.2.4 · p. 9

4. Numerical Tests

• The test problems verify the finite-difference methods, which are implemented in PLUTO to ease comparisons in a common computational framework.The framework is described as flexible.Source: §4 · p. 9
• Unless stated otherwise, Γ = 5/3 and the Courant number is 0.8, 0.4, or 0.3 for one-, two-, or three-dimensional computations, respectively.Source: §4 · p. 10
• Errors are computed with the discrete L1 norm, and magnetic-field divergence is quantified using Eq. (42).The passage introduces the computation of ∇·B for that divergence measure.Source: §4 · p. 10

4.1.1. One Dimensional Propagation

• At 64 zones, fifth-order schemes achieve essentially the accuracy of third-order schemes using 256 zones, and all four schemes meet their expected convergence orders.The expected orders are three for LimO3 and WENO+3, and five for WENO-Z and MP5.Source: §4.1.1 · p. 11

4.1.2. Three Dimensional Oblique Propagation √

• In 3D oblique propagation, all schemes meet their expected orders; fifth-order schemes need roughly one-quarter the resolution of third-order schemes to match accuracy.The fifth-order profiles are nearly indistinguishable from the exact solution by Nx ≳32, versus Nx > 64 for third-order schemes.Source: §4.1.2 · p. 12
• Third- and fifth-order schemes outperform the at-most-second-order CT-PPM algorithm; in the 3D oblique test, their relative CPU costs are 1 : 0.98 : 1.46 : 1.31 for WENO+3, LimO3, WENO-Z, and MP5.The cost ratios are for WENO+3, LimO3, WENO-Z, and MP5, respectively.Source: §4.1.2 · p. 12

4.1. Propagation of Circularly polarized Alfv´en Waves

• Circularly polarized Alfvén waves are exact nonlinear MHD solutions, so their decay directly indicates an algorithm’s numerical viscosity and resistivity.The test uses a standing wave with ρ = 1 and ca = 1.Source: §4.1 · p. 10, §4.1.3 · p. 12

4.1.3. Numerical Dissipation and Long Term Decay in Two Dimensions

• At 120 × 20 mesh resolution, after ≈37 wave transits, fifth-order schemes retain ∼94% of nominal amplitude, versus ∼20% for third-order schemes and ∼3% for second-order TVD with Monotonized Central limiter.Fifth-order schemes preserve the original shape more accurately.Source: §4.1.3 · p. 12

4.2.1. Two-Dimensional Shock tube

• In the 2D rotated Brio-Wu test, all discontinuities are captured and profiles match the reference well; MP5 best describes the normal magnetic field, with ∼0.7% relative peak oscillations.Divergence errors remain bounded, typically ≲10^-2, and tend to saturate for αp ≳0.4 in both 2D and 3D calculations.Source: §4.2.1 · p. 13
• Across 2D shock discontinuities, fifth-order methods use fewer points per discontinuous layer, but their accuracy gain over third-order methods is only a factor of 1.5–2.Interpolation across discontinuities usually degenerates to lower order to suppress spurious oscillations.Source: §4.2.1 · p. 13

4.2.2. Three-Dimensional Shock tube

• For this discontinuous-flow problem, high-order methods provide no particularly significant gain; no prominent over/under-shoots are observed, and normal magnetic-field oscillations are comparable to or smaller than those in.Relative to WENO+3 (= 1), computational costs for this problem are 0.99 : 1.48 : 1.25 for LimO3, WENO-Z, and MP5, respectively.Source: §4.2.2 · p. 14

4.3.1. Two Dimensional Propagation

• At 128^2, fifth-order schemes produce Bx errors ∼4 times smaller than third-order schemes at 256^2, while all schemes quickly converge to their asymptotic orders.The third-order schemes have identical errors and are somewhat more accurate than those reported in.Source: §4.3.1 · p. 14

4.3.2. Three Dimensional Propagation

• The paper introduces a 3D extension of the vortex problem, with the radius defined spherically rather than in the 2D plane.The 3D calculations show excellent agreement with the analytical solution; their L2 results are qualitatively competitive with, though differing configurations and norms complicate direct comparison.Source: §4.3.2 · p. 14
• For fifth-order schemes, 3D divergence errors fall below 10^-8 at resolution 128^3 and decrease rapidly as the mesh thickens.Reported computational-cost ratios for WENO+3, LimO3, WENO-Z, and MP5 are 1 : 0.99 : 1.46 : 1.24.Source: §4.3.2 · p. 14

4.4. Advection of a magnetic field loop

• The loop is transported as a passive scalar at sufficiently large plasma β, making shape preservation a test of dissipation and multidimensional-term balance.In two dimensions, the vector-potential modification near the center removes a singularity that can cause spurious oscillations and erroneous magnetic-energy evaluations.Source: §4.4 · p. 15, §4.4.1 · p. 15

4.4.1. Two Dimensional Propagation

• In two dimensions, all schemes preserve the initial circular shape, but the third-order schemes are substantially more diffusive and MP5 is least diffusive.LimO3 and WENO+3 show no pronounced difference for this problem.Source: §4.4.1 · pp. 15–16
• The two-dimensional divergence error is minimized for fifth-order schemes at αp ≳ 0.4, but for LimO3 and WENO+3 at αp ≲ 0.2.Source: §4.4.1 · p. 16

4.4.2. Three Dimensional Propagation

• The three-dimensional test uses a rotated loop and velocity (1, 1, 2), making its evolution intrinsically three-dimensional despite rotation around one axis.The domain is resolved with 128 × 128 × 256 zones and periodic boundaries.Source: §4.4.2 · p. 16
• All schemes preserve the loop’s shape in three dimensions, while MP5 is least diffusive and LimO3 and WENO+3 have practically coincident dissipation rates.The three-dimensional divergence-error behavior is similar to the two-dimensional case.Source: §4.4.2 · p. 16

• At t = 0.5, WENO+3 and LimO3 show more numerical dissipation than WENO-Z and MP5, while the results agree well with previous investigations.Source: §4.5 · p. 16
• At t = 1, fifth-order schemes, particularly MP5, show a central magnetic island with a high-density spot that is absent from the third-order schemes.The structure may be induced by decreased effective resistivity across the central current sheet.Source: §4.5 · p. 16
• At t = 0.5, divergence errors are comparable with those reported by other investigators and reach their maximum magnitude near discontinuities.Source: §4.5 · p. 16
• The relative computational costs of LimO3, WENO-Z, and MP5 are 1.01, 1.47, and 1.29, respectively, with WENO+3 set to 1.Source: §4.5 · p. 16

• MP5 needs half the grid resolution to match WENO+3 results during the linear growth phase; WENO-Z also gains over LimO3, but less markedly.Magnetic energy grows faster with increased resolution and when moving from third- to fifth-order schemes.Source: §4.6 · p. 17
• Growth rate evolves faster for finer grids and less dissipative schemes, tracking magnetic-field amplification and smaller numerical resistivity.Source: §4.6 · p. 17
• During saturation, fifth-order methods show small-scale structures more clearly, while WENO+3 and LimO3 appear more diffused.The field amplification is eventually limited after t ≳ 8 by tearing-mode instabilities and reconnection events that expel magnetic flux.Source: §4.6 · p. 17
• The relative CPU costs of LimO3, WENO-Z, and MP5 are 0.98, 1.48, and 1.24, respectively, with WENO+3 set to 1.Source: §4.6 · p. 17

• The framework represents primary flow variables at cell centers and uses a generalized Lagrange multiplier for mixed hyperbolic/parabolic propagation and damping of divergence errors.Its one-dimensional reconstruction avoids elliptic cleaning, genuinely multidimensional interpolation, and staggered-mesh complexity.Source: §5 · p. 17
• WENO+3 and LimO3 show third-order accuracy in smooth problems and non-oscillatory transitions at discontinuities.Source: §5 · p. 17
• WENO-Z and MP5 yield high-quality results across selected tests and near-fifth-order accuracy on multidimensional smooth problems, with MP5 slightly better in cost, accuracy, and sharpness at discontinuities.Source: §5 · p. 17
• Fifth-order schemes are ≲50% slower for WENO-Z and ≲30% slower for MP5 than third-order schemes, while three-dimensional problems can yield almost two orders of magnitude in CPU-cost gain.Source: §5 · pp. 17–18
• The finite-difference formulation is competitive with recent finite-volume schemes in accuracy and discontinuity description, and reduces numerical dissipation relative to traditional second-order Godunov-type schemes.Source: §5 · p. 18

Appendix A. Conservative eigenvectors of the GLM-MHD Equations

• The 9 × 9 conservative MHD system is decomposed into left and right eigenvectors using its eigenvalues.The appendix introduces the eigenvectors following partially the notation of.Source: Appendix A · p. 18
• The right eigenvectors are given in matrix form, with Γp/ρ denoting the speed of sound.Source: Appendix A · p. 18
• The appendix prescribes τ, γ1, γ2, and I(f,s)(vi,Bi), and then introduces the left eigenvectors.It defines I(f,s)(vi,Bi) for i = x, y, z.Source: Appendix A · pp. 18–19