Source-linked AI summary
A Solution Accurate, Efficient and Stable Unsplit Staggered Mesh Scheme for Three Dimensional Magnetohydrodynamics
Dongwook Lee
TL;DR
The paper addresses extending an efficient unsplit MHD scheme to 3D while retaining a broad stability range. It develops reduced and full CTU variants with characteristic-tracing corrections and an upwind-biased CT update. The full scheme runs at CFL numbers close to 1 without the 12 Riemann solves associated with the standard full CTU approach.
Problem
The paper targets an efficient 3D unsplit MHD method that reaches the full CFL stability range without requiring 12 Riemann solves per cell per time step.
Source: §1 · p. 3Method
The authors extend 2D USM to 3D with reduced and full CTU schemes, characteristic-tracing transverse updates, constrained transport, and upwind-biased electric-field averaging.
Source: §5 · pp. 27, 29Results
CFL < 0.5 limits the reduced CTU scheme, while the full CTU scheme operates at CFL numbers close to 1.
Source: §5 · p. 29Takeaways & Limitations
The test suite supports the full scheme as a robust second-order 3D MHD algorithm that preserves the divergence-free constraint near CFL 1.
Source: §5 · p. 29For non-ideal MHD, upwinding in upwind-MEC should be applied only to the induction part of the generalized Ohm’s law.
Source: §2.6.2 · p. 16
Abstract
from arXiv · showhide
In this paper, we extend the unsplit staggered mesh scheme (USM) for 2D magnetohydrodynamics (MHD) (Lee and Deane, 2009) to a full 3D MHD scheme. The scheme is a finite-volume Godunov method consisting of a constrained transport (CT) method and an efficient and accurate single-step, directionally unsplit multidimensional data reconstruction-evolution algorithm, which extends the original 2D corner transport upwind (CTU) method (Colella, 1990). We present two types of data reconstruction-evolution algorithms for 3D: (1) a reduced CTU scheme and (2) a full CTU scheme. The reduced 3D CTU scheme is a variant of a simple 3D extension of the 2D CTU method by Colella (1990) and is considered as a direct extension from the 2D USM scheme. The full 3D CTU scheme is our primary 3D solver which includes all multidimensional cross-derivative terms for stability. The latter method is logically analogous to the 3D unsplit CTU method by Saltzman. The major novelties in our algorithms are twofold. First, we extend the reduced CTU scheme to the full CTU scheme which is able to run with CFL numbers close to unity. Both methods utilize the transverse update technique developed in the 2D USM algorithm to account for transverse fluxes without solving intermediate Riemann problems, which in turn gives cost-effective 3D methods by reducing the total number of Riemann solves. The proposed algorithms are simple and efficient especially when including multidimensional MHD terms that maintain in-plane magnetic field dynamics. Second, we introduce a new CT scheme that makes use of proper upwind information in taking averages of electric fields.
1. Introduction
The paper targets accurate, efficient 3D MHD simulation, where reduced-dimensional models miss some full-3D physics and split schemes have important multidimensional limitations. It extends the 2D USM approach with a 3D CTU method designed for full-CFL stability without 12 Riemann solves per cell per step.
1. Introduction
- Reduced-dimensional simulations cannot capture nonlinear physics specific to full 3D, motivating a 3D MHD solver for magnetized flows and shocks.The introduction cites the 3D energy cascade as an example of physics not represented by reduced systems.Source: §1 · pp. 1–2
- Split MHD schemes can lose expected flow symmetries, fail to preserve in-plane magnetic-field evolution, and produce grid-scale strain-rate artifacts.Split formulations also fail to evolve the normal magnetic field in MHD.Source: §1 · p. 2
- Existing 3D unsplit MHD CTU methods face a stability-cost trade-off: a 6-solve method has CFL < 0.5, while full CFL can require 12 Riemann solves per zone.The cited 6-solve method was chosen for lower complexity, despite its more restrictive CFL bound.Source: §1 · p. 2
- The paper extends 2D USM to 3D, aiming for full-CFL stability without 12 Riemann solves and introducing upwind-biased electric-field averaging for CT.Its design retains a 1D Riemann-solver framework and extends the 2D scheme without modifying the Riemann solvers.Source: §1 · p. 3
2.1. MHD Equations
The section formulates MHD in conservative variables and fluxes, then presents a quasi-linear form in primitive variables for discretizing the coupled equations. The induction equation carries the solenoidal magnetic-field constraint.
2.1. MHD Equations
- The induction equation implies the solenoidal constraint ∇·B = 0.The text defines total pressure as the sum of thermal and magnetic pressure.Source: §2.1 · p. 4
- The conservative state comprises eight MHD variables, with density, momentum, magnetic field, and total energy among them.The section defines thermal pressure from the energy and notes that thermal and magnetic pressure sum to total pressure.Source: §2.1 · p. 4
- The conservative formulation uses x-, y-, and z-direction fluxes, while its quasi-linear form uses primitive variables to discretize the coupled MHD system.The primitive state is listed as density, three velocity components, three magnetic-field components, and pressure.Source: §2.1 · p. 4
2.2. Normal Predictor
The normal predictor constructs second-order interface states from cell-centered primitive variables using characteristic tracing. It treats the normal magnetic field separately so that the divergence-free constraint is maintained.
2.2. Normal Predictor
- The normal predictor treats the normal magnetic field separately while including multidimensional MHD terms and satisfying ∇·B = 0.The text describes the normal-field evolution before the other primitive variables.Source: §2.2 · p. 5
- The hat notation denotes the reduced, usual 1D MHD system, while the bar notation denotes the reassembled full system.For the x-direction, the reduced state and coefficient matrix omit the normal magnetic-field component.Source: §2.2 · pp. 5–6
- MUSCL-Hancock extrapolation builds second-order face states with TVD limiting in characteristic variables, but leaves the normal field unlimited to preserve face continuity.The text notes that characteristic limiting is less prone to spurious oscillations than limiting primitive variables.Source: §2.2 · p. 6
- Characteristic tracing uses directional eigenvalues and eigenvectors with a TVD slope limiter to complete the normal predictor in each direction.The paper describes the x-direction calculation and applies the same approach in the other directions.Source: §2.2 · p. 7
2.3. Transverse Corrector in USM
The transverse corrector adds transverse flux effects to normal predictors using characteristic tracing, reducing intermediate Riemann solves. The reduced 3D CTU scheme omits diagonal corner couplings; the full scheme adds them to reach the full stability range, while CT advances normal magnetic fields subject to continuity and divergence constraints.
2.3. Transverse Corrector in USM
- Characteristic tracing supplies transverse flux-gradient corrections without intermediate Riemann solves, reducing the Riemann-solve count while simplifying multidimensional MHD-term handling.The paper describes three Riemann solves for the hydrodynamic update and three additional solves for the CT magnetic-field update.Source: §2.3.1 · p. 8
- The reduced 3D CTU scheme includes second-order cross derivatives but omits third-order cross derivatives, so it does not fully capture diagonal corner flow and has CFL < 0.5 stability.It extends the 2D CTU correction by adding the third, z-direction transverse correction.Source: §2.3.2 · p. 9
- The transverse corrector uses an upwind-biased slope limiter; using a TVD limiter instead would reduce the scheme’s stability limit.The upwind coupling is the mechanism identified for the full CFL bound in the original 2D CTU scheme.Source: §2.3.2 · p. 9
- The full 3D CTU scheme adds diagonal corner couplings through additional cross-derivative corrections to the reduced scheme’s interface states.These new terms address flows with comparable velocity components in all three directions.Source: §2.3.3 · pp. 9, 11, 12
- The full-scheme corrections reuse eigensystems stored during the normal predictor, leaving upwind differences and dot products as the added calculations.The paper says these additional calculations are much more efficient than those required by the 12-solve CTU scheme.Source: §2.3.1 · p. 8, §2.3.3 · p. 12
- CT advances normal magnetic fields by a half time step using Godunov fluxes, while enforcing face continuity and the divergence-free constraint.The electric fields used for the CT update are constructed from the fluxes.Source: §2.4 · p. 12
2.5. Cell-centered Solution Update from n to n + 1 Time Step
The cell-centered update is followed by a constrained-transport update of staggered magnetic fields. The new upwind-MEC averaging uses directional upwind information to improve stability in biased flows.
2.5. Cell-centered Solution Update from n to n + 1 Time Step
- The unsplit integrator advances cell-centered conservative variables, then corrects cell-centered magnetic fields using face-centered fields; full CTU results use CFL 0.95.The face-centered fields are updated separately through constrained transport.Source: §2.5 · p. 13
2.6.1. The Standard Arithmetic Averaging Approach in MEC: standard-MEC
- Standard-MEC averages four Taylor-expanded electric fields, but simple arithmetic averaging lacks upwinding for directionally biased multidimensional flows.The paper identifies small-angle advection as especially vulnerable to spurious magnetic-field oscillations without upwinding.Source: §2.6.2 · pp. 13–14
2.6.2. Upwind Biased Averaging in MEC: upwind-MEC
- Upwind-MEC uses directional upwind information and high-order terms in corner electric-field averaging to reduce oscillations during small-angle field-loop advection.Its high-order terms add dissipation near discontinuities while remaining consistent for plane-parallel, grid-aligned flows.Source: §2.6.2 · pp. 14, 16
- When local velocities are negligible relative to sound speed, the method switches to standard-MEC; otherwise, upwind-MEC is used for the numerical tests.The empirical switching parameter ε3 is 10^-4.Source: §2.6.2 · pp. 16–17
- For non-ideal MHD, upwinding applies only to the induction term −u×B in the generalized Ohm’s law.The stated scope excludes applying this upwinding to magnetic diffusion, the Hall effect, and the Biermann battery effect.Source: §2.6.2 · p. 16
2.6.3. CT Update from n to n + 1 Time Step
- MEC-constructed electric fields drive the CT update of face-centered magnetic fields, maintaining ∇·B = 0 on the staggered grid.The passage displays the update in the x direction only.Source: §2.6.3 · p. 17
3. Summary
The 3D USM update first constructs normal predictor states, then applies transverse corrections and advances the conservative and magnetic fields. Its full CTU option adds a correction for diagonal upwind information, while the reduced option uses fewer correction steps.
3. Summary
- The normal predictor computes states in each direction and stores eigensystems and characteristic-jump sums for later transverse corrections.The transverse correctors use stored sigma sums to represent transverse flux gradients through characteristic tracing.Source: §3 · p. 17
- The full CTU scheme adds a correction for diagonally moving upwind information at control-volume corners to provide the full stability limit.The reduced scheme instead advances normal fields by a half time step through CT before completing interface-state calculations.Source: §3 · p. 17
- Both CTU variants use three Riemann problems for the CT magnetic-field advance and three final interface solves, for six total.The final interface solves update cell-centered conservative variables.Source: §3 · p. 17
- The algorithm constructs corner electric fields, updates face-centered magnetic fields to the next time step, and averages those fields to correct cell-centered magnetic fields.The unsplit reconstruction-evolution algorithm can also serve as a hydrodynamics solver by omitting magnetic-field steps; that version requires three Riemann solves and no CT update.Source: §3 · pp. 17–18
4. Numerical Results
The numerical results assess the 3D USM scheme across MHD tests, reporting its stability, accuracy, and divergence preservation, as well as shock propagation and pressure behavior in a strong-field blast-wave case.
4. Numerical Results
- The numerical-results suite uses the full 3D CTU method as its default at CFL=0.95 and reports second-order accuracy for smooth flows and magnetic solenoidal constraint preservation to machine round-off.The simulations use Roe-type linearized or HLLD Riemann solvers, with MUSCL-Hancock, PPM, and WENO5 normal predictors.Source: §4 · p. 18
4.1. Field Loop Advection
- The upwind-MEC scheme stabilizes small-angle field-loop advection, preserves the cylindrical shape without oscillations in 3D, and successfully evolves the standard loop to t=1.In the standard test, maximum magnetic pressure rises 48% to 7.41×10^-7 from 5×10^-7 initially, versus 69% growth with standard-MEC; normalized final B3 errors are slightly larger than in [31, 47].Source: §4.1 · pp. 19–20, Figure 3 · p. 20, Figure 4 · p. 21
4.2. Circularly Polarized Alfv´en Wave
- The full 3D CTU scheme is about 20% more computationally efficient than the reduced scheme at CFL=0.95 versus 0.475, with nearly identical errors and second-order accuracy on smooth Alfvén waves.The reported relative CPU cost is CPUrel ≈ 0.8; truncation error depends on CFL, with lower CFL yielding smaller standing-wave error but larger traveling-wave error.Source: §4.2 · pp. 21–22, Figure 5 · p. 22
4.3. Orszag-Tang Problem
- The 3D Orszag-Tang test develops strong discontinuities from smooth initial conditions while preserving flow symmetries and producing density patterns similar to the corresponding 2D problem.The density evolution is shown at t=0.5 and 1.0 on a 128^3 grid.Source: §4.3 · p. 22, Figure 6 · p. 23
4.4. Rotor Problem
- The 3D rotor results closely match the underlying 2D solutions and preserve symmetric oval Mach-number contours near the center at CFL=0.95.The reported contours are at t=0.15 on a 128^3 grid.Source: §4.4 · p. 24, Figure 7 · p. 24
4.5. Cloud & Shock Interaction
- The full 3D CTU and upwind-MEC scheme successfully completes the cloud-shock interaction at CFL=0.95 without numerical instability, capturing the cloud’s disrupted shapes as it enters the planar shock.This succeeds with the Roe solver and van Leer limiter despite the strong shock, high supersonic Mach numbers, wide plasma-beta range, and shear flows.Source: §4.5 · p. 25, Figure 8 · p. 26
4.6. MHD Blast Wave
- With the hybrid Roe-HLLE solver, the zero-field blast produces spherical shocks without grid-aligned shock artifacts, while the intermediate-field case develops elongated structures parallel to the magnetic field.The intermediate field still permits propagation in the y-direction, so the overall shock shape is not radically distorted.Source: §4.6 · pp. 26–27, Figure 9 · p. 27
- In the strongest magnetic-field case, transverse gas displacement is increasingly inhibited and shocks propagate almost entirely parallel to Bx.The explosion becomes highly anisotropic.Source: §4.6 · p. 27
- Weak magneto-sonic waves radiate transverse to the field until total pressure equilibrium is reached in the central region.Source: §4.6 · p. 27
- Pressure remains positive throughout the simulations at high CFL=0.95, demonstrating robustness against the positivity challenge noted by Balsara [7].Balsara [7] identified strong oblique wave propagation as a challenge for maintaining positive pressure.Source: §4.6 · p. 27
5. Conclusion
The 3D USM scheme combines characteristic analysis of normal and transverse MHD fluxes with CTU updates, avoiding intermediate Riemann solves. Its full CTU variant supports CFL numbers close to 1, while upwind-MEC improves field-loop advection by adding dissipation to CT electric-field averaging.
5. Conclusion
- The 3D USM scheme uses characteristic analysis for normal and transverse MHD fluxes, avoiding intermediate Riemann solves and reusing eigensystem evaluations to reduce computational cost.The multidimensional MHD terms are included in both normal and transverse directions.Source: §5 · pp. 27, 29
- The full CTU scheme reaches CFL numbers close to 1, unlike the reduced scheme's limit below 0.5, by including second- and third-order cross-derivative terms.Its relative CPU cost compared with the reduced scheme is less than 1.Source: §5 · p. 29
- Upwind-MEC adds numerical dissipation when averaging electric fields in CT, addressing field-loop distortion in biased-direction advection and maintaining consistency for plane-parallel and grid-aligned flows.The paper contrasts this with existing CT algorithms that can distort field loops in small-angle 2D and 3D advection tests.Source: §5 · p. 29
- The test suite supports a robust, reliable second-order 3D MHD scheme that preserves the divergence-free constraint without observed instability or unphysical-error accumulation at CFL numbers close to 1.The authors report testing on several stringent setups and a wide range of MHD phenomena.Source: §5 · p. 29
- The 3D USM scheme is implemented on uniform and AMR grids and integrated and tested in the official FLASH4 release.Source: §5 · p. 29