Systematic construction of upwind constrained transport schemes for MHD
Introduction
Magnetohydrodynamics (MHD) is the basic modelization framework to treat plasmas at the macroscopic level, that is neglecting kinetic effects and as a single fluid, an approximation commonly used for applications to laboratory, space and astrophysical plasmas. Magnetic fields and currents created by the moving charges play a fundamental role in the dynamics of the fluid, which is considered to be locally neutral, and the set of hydrodynamical (Euler) equations must be supplemented by the magnetic contributions to the global energy and momentum, and by a specific prescription for the evolution of the magnetic field itself, the so-called induction equation, that is Faraday's law combined to a constitutive relation between the current and the electric field (Ohm's law).
In extending the same numerical techniques employed for the Euler equations to multi-dimensional MHD, a major challenge dwells in preserving the divergence-free constraint of the magnetic field which inherently follows from the curl-type character of Faraday's law. This is especially true for Godunov-type shock-capturing schemes based on the properties of the hyperbolic set of conservation laws, let us call them standard upwind procedures, where spatial partial derivatives do not commute when discontinuities are present and spurious effects (numerical magnetic monopoles) arise.
The research in this field, which is crucial for building accurate and robust finite-volume (FV) or finite-difference (FD) shock-capturing numerical codes for computational MHD, started more than twenty years ago and several different methods have been proposed. Here, for sake of conciseness, we simply refer the reader to the paper by Tóth (2000) [1] for a comprehensive discussion and comparison of the early schemes. Summarizing in brief, among the proposed solution methods are schemes based on the cleaning of the numerical monopoles, either by solving an elliptic (Poisson) equation and removing the monopole contribution from the updated fields [2], or by adding specific source terms and an additional evolution equation for the divergence of B itself to preserve the hyperbolic character of the MHD set [3], [4], [5]. Other methods evolve in time the vector potential [6], [7] or use different alternatives, such as for the so-called flux distribution schemes [8].
Conversely, a radically different strategy is adopted in the so-called constrained transport (CT) methods. These schemes all rely on the curl-type nature of the induction equation and on its discretization based on Stokes' theorem (rather than on Gauss' one as needed for the equations with the divergence operator), as first realized in pioneering works where the evolution equation for the magnetic field alone was solved [9], [10], [11], [12]. The CT method was later extended to the full MHD system of hyperbolic equations in the context of Godunov-type schemes [13], [14], [15]. In FV-CT schemes, magnetic field components are stored as surface integrals at cell interfaces as primary variables to be evolved via the induction equation, while the corresponding fluxes are line-averaged electric field (namely electromotive force, emf) components located at cell edges, to recover the discretized version of Stokes' theorem. By doing so, the solenoidal constraint can be preserved exactly during time evolution.
A major difficulty of the CT formalism is the computation of upwind-stable emf components located at zone-edges [15]. This can be achieved either by properly averaging the interface fluxes computed when solving the 1D Riemann problems at zone interfaces, or by using genuine (but much more complex) 2D Riemann solvers computed directly at cell edges [16], [17], [18], [19]. In the latter case the dissipative part of the multidimensional emf can be shown to behave as a proper resistive term for the induction equation [20]. As far as the former (simpler) case is concerned, in the original work by [15] the emf was obtained as the arithmetic of the four upwind fluxes nearest to the zone edge. It was then recognized (e.g. [21]) that this approach has insufficient numerical dissipation and it does not reduce to the plane-parallel algorithm for grid-align flow. Gardiner & Stone (2005) [21] suggested that this issue could be solved by doubling the dissipation and introduced a recipe to construct a stable and non-oscillatory upwind emf with optimal numerical dissipation based on the direction of the contact mode. This approach (here referred to as the CT-Contact method) is, however, mainly supported by empirical results as there is no formal justification that the emf derivative should obey such selection rule. In addition, the method can be at most 2nd-order accurate thus making the generalization to higher-order methods not feasible.
A rigorous approach to this problem for both FV and FD Godunov-type schemes for computational MHD was originally proposed by Londrillo and Del Zanna (2004) [22] with their upwind constrained transport (UCT) method. According to the UCT methodology, the continuity property of the magnetic field components at cell interfaces (which follows from the solenoidal constraint) is considered as a built-in condition in a numerical scheme, enabling face-centered fields to be evolved as primary variables. At the same time, staggered magnetic field components enter as single-state variables in the fluid fluxes at the corresponding cell interfaces, and as two-state reconstructed values at cell edges in the four-state emf for the induction equation. Time-splitting techniques should be avoided as they prevent exact cancellation of terms at the numerical level. The emf components constructed using information from the four neighboring upwind states must also automatically reduce to the correct 1D numerical fluxes for plane parallel flows and discontinuities aligned with the grid directions. According to the authors, these are the necessary conditions to preserve the divergence-free condition and to avoid the occurrence of numerical monopoles that may arise while computing the divergence of fluid-like fluxes numerically. In the original work of [22], a second-order FV scheme based on Roe-type MHD solver (UCT-Roe) and a high-order scheme based on characteristic-free reconstruction and a two-wave HLL approximate Riemann solver (UCT-HLL) were proposed.
The UCT-HLL scheme was further simplified by Del Zanna et al. (2007) [23] in the context of general relativistic MHD, and recipes were given to build a FD UCT schemes of arbitrary order of accuracy by testing several reconstruction methods. High-order FD-CT methods were also recently proposed by [24] who, instead, constructed the emf by simply doubling the amount of numerical dissipation. While FD approaches are based on a point value representation of primary variables and avoid multi-dimensional reconstructions, FV-CT schemes of higher than 2nd-order accuracy are much more arduous to construct albeit they are likely to increase robustness, see the review by Balsara [25] and references therein. Efforts in this direction were taken by Balsara (2009) [26], [27] in the context of ADER-WENO FV schemes who designed genuinely third- and fourth-order spatially accurate numerical methods for MHD. More recently, fourth-order FV schemes using high-order quadrature rules evaluated on cell interfaces have been proposed by Felker & Stone (2018) [28] and Verma (2019) [29]. Here the construction of the higher-order emf follows the general guidelines of the UCT-HLL (or Lax-Friedrichs) approach introduced by [23] and later resumed for truly 2D Riemann problem by [16].
The goal of the present work is to systematically construct UCT schemes for classical MHD, using a variety of 1D Riemann solvers avoiding the full spectral resolution in characteristic waves, and providing the correct averaging recipes to build the four-state emf fluxes at zone edges, extending the scheme by [23] to less dissipative solvers like HLLD [30], where the Riemann fan is split into four intermediate state to include also the contact and Alfvénic contributions, other than simply the fast magnetosonic ones. This novel UCT scheme is tested by performing several multi-dimensional numerical tests, and comparison is also made with other CT popular schemes based on emf averaging, including the simple arithmetic averaging method, those based on doubling the diffusive contribution of 1D fluxes, as in [21], the above mentioned UCT-HLL one, and a novel UCT version of the GFORCE scheme by [31].
Our paper is structured as follows. In §2 we introduce the basic CT discretization and general notations while in §3 we review basic averaging CT schemes. The UCT method and the original Roe and HLL schemes are discussed in §4, while the new UCT-based composition formulae are presented in §5. Numerical benchmarks are introduced in §6 and a final summary is reported in §7.
Section snippets
Notations and general formalism for CT schemes
The ideal MHD equations are characterized by two coupled sub-systems, one for the time evolution of the set of conservative hydro-like flow variables (mass, momentum and energy densities denoted, respectively, with ρ, and ): and an induction equation for the evolution of the magnetic field where the curl operator appears instead of a divergence and where the electric field is not a primary variable, depending on the
CT schemes based on emf averaging
The set of conserved variables may be extended to include the zone-centered representation of magnetic field components as well, usually provided by simple spatial average from the two neighboring faces along the relevant direction at the beginning of any timestep. The solution to the full Riemann problem (8 equations and variables for 3D MHD) thus provides point-value upwind fluxes for the zone-centered magnetic field as well. Indeed, indicating with a square bracket the flux component we
The upwind constrained transport (UCT) method: original framework
The CT discretization scheme outlined so far allows to preserve exactly the solenoidal condition for the magnetic field but this is not enough to avoid the presence of spurious magnetic monopole terms when computing the magnetic forces, if fluxes are calculated by using reconstructed values of the magnetic field components as for the other fluid variables. As first realized in [33] and systematically demonstrated in [22], the only way to properly take into account the specific smoothness
The UCT method: a novel generalized composition formula for Jacobian-free Riemann solvers
We now present a general formalism for constructing emf averaging schemes with in-built upwind dissipation properties and using component-wise Riemann solver, that is, not directly employing characteristic information. Our starting point is the definition of the inter-cell numerical flux function for which we assume the approximate form given by Eq. (12). We shall also assume that the dissipative terms of the induction system can be expressed as linear combinations of the left and right
Numerical benchmarks
In what follows we compare different emf-averaging schemes in terms of accuracy, robustness and dissipation properties. Our selection includes:
- •
Thew arithmetic averaging, given by Eq. (19);
- •
the CT-Contact emf averaging, given by Eq. (20);
- •
the CT-Flux emf, given by Eq. (23);
- •
the UCT-HLL scheme following our new composition formula, Eq. (33) with coefficients given by Eq. (32). Extensive numerical testing (including several other additional tests not shown here) has demonstrated our formulation of
Summary
The systematic construction and comparison of averaging schemes to evaluate the electromotive force (emf) at zone edges in constrained transport MHD has been the subject of this work. The upwind constrained transport (UCT) formalism, originally developed by [22], has been reconsidered under a more general perspective where the edge-averaged electric field can be constructed using the information available from 1D face-centered, component-wise Riemann solvers. This approach offers enhanced
Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Acknowledgements
The authors wish to thank P. Londrillo and G. Bodo for useful discussions. We also acknowledge the OCCAM supercomuting facility available at the Competence Centre for Scientific Computing at the University of Torino.
References (62)
The constraint in shock-capturing magnetohydrodynamics codes
J. Comput. Phys.
(2000)- et al.
The effect of nonzero on the numerical solution of the magnetohydrodynamic equations
J. Comput. Phys.
(1980) - et al.
A solution-adaptive upwind scheme for ideal magnetohydrodynamics
J. Comput. Phys.
(1999) - et al.
Divergence correction techniques for Maxwell solvers based on a hyperbolic model
J. Comput. Phys.
(2000) - et al.
Hyperbolic divergence cleaning for the MHD equations
J. Comput. Phys.
(2002) - et al.
An unstaggered constrained transport method for the 3D ideal magnetohydrodynamic equations
J. Comput. Phys.
(2011) Flux-corrected transport techniques for multidimensional compressible magnetohydrodynamics
J. Comput. Phys.
(1991)- et al.
A staggered mesh algorithm using high order Godunov fluxes to ensure solenoidal magnetic fields in magnetohydrodynamic simulations
J. Comput. Phys.
(1999) Multidimensional HLLE Riemann solver: application to Euler and magnetohydrodynamic flows
J. Comput. Phys.
(2010)A two-dimensional HLLC Riemann solver for conservation laws: application to Euler and magnetohydrodynamic flows
J. Comput. Phys.
(2012)
Multidimensional Riemann problem with self-similar internal structure. Part I - application to hyperbolic conservation laws on structured meshes
J. Comput. Phys.
Multidimensional Riemann problem with self-similar internal structure. Part II - application to hyperbolic conservation laws on unstructured meshes
J. Comput. Phys.
Multidimensional Riemann problem with self-similar internal structure. Part III - a multidimensional analogue of the HLLI Riemann solver for conservative hyperbolic systems
J. Comput. Phys.
An unsplit Godunov method for ideal MHD via constrained transport
J. Comput. Phys.
On the divergence-free condition in Godunov-type schemes for ideal magnetohydrodynamics: the upwind constrained transport method
J. Comput. Phys.
Divergence-free reconstruction of magnetic fields and WENO schemes for magnetohydrodynamics
J. Comput. Phys.
Efficient, high accuracy ADER-WENO schemes for hydrodynamics and divergence-free magnetohydrodynamics
J. Comput. Phys.
A fourth-order accurate finite volume method for ideal MHD via upwind constrained transport
J. Comput. Phys.
A multi-state HLL approximate Riemann solver for ideal magnetohydrodynamics
J. Comput. Phys.
MUSTA fluxes for systems of conservation laws
J. Comput. Phys.
An HLLC Riemann solver for magneto-hydrodynamics
J. Comput. Phys.
Von Neumann stability analysis of globally divergence-free RKDG schemes for the induction equation using multidimensional Riemann solvers
J. Comput. Phys.
A new efficient formulation of the HLLEM Riemann solver for general conservative and non-conservative hyperbolic systems
J. Comput. Phys.
A simple and accurate Riemann solver for isothermal MHD
J. Comput. Phys.
Accurate monotonicity-preserving schemes with Runge Kutta time stepping
J. Comput. Phys.
Roe matrices for ideal MHD and systematic construction of roe matrices for systems of conservation laws
J. Comput. Phys.
An unsplit Godunov method for ideal MHD via constrained transport in three dimensions
J. Comput. Phys.
A second-order unsplit Godunov scheme for cell-centered MHD: the CTU-GLM scheme
J. Comput. Phys.
An unsplit staggered mesh scheme for multidimensional magnetohydrodynamics
J. Comput. Phys.
High-order conservative finite difference GLM-MHD schemes for cell-centered MHD
J. Comput. Phys.
A robust numerical scheme for highly compressible magnetohydrodynamics: nonlinear stability, implementation and tests
J. Comput. Phys.
Cited by (24)
Pseudo-spectral solver versus grid-based solver: A quantitative accuracy test using GMHD3D and PLUTO4.4
2024, Computers and FluidsA 4<sup>th</sup>-order accurate finite volume method for ideal classical and special relativistic MHD based on pointwise reconstructions
2024, Journal of Computational PhysicsA polarization study of jets interacting with turbulent magnetic fields
2023, Monthly Notices of the Royal Astronomical Society