Systematic construction of upwind constrained transport schemes for MHD

https://doi.org/10.1016/j.jcp.2020.109748Get rights and content

Highlights

  • With this work we assess the validity of novel UCT schemes for the solution of the MHD equations in multiple dimensions.

  • The proposed schemes improve on the original formulation of Londrillo & Del Zanna (2004) [22] by providing more general guidelines allowing different Riemann solvers to be incorporated.

  • This approach offers enhanced flexibility allowing new upwind techniques to be incorporated in CT-MHD schemes at the modest cost of storing transverse velocity, weight coefficients for the fluxes and diffusion terms for the magnetic field.

Abstract

The constrained transport (CT) method reflects the state of the art numerical technique for preserving the divergence-free condition of magnetic field to machine accuracy in multi-dimensional MHD simulations performed with Godunov-type, or upwind, conservative codes. The evolution of the different magnetic field components, located at zone interfaces using a staggered representation, is achieved by calculating the electric field components at cell edges, in a way that has to be consistent with the Riemann solver used for the update of cell-centered fluid quantities at interfaces. Albeit several approaches have been undertaken, the purpose of this work is, on the one hand, to compare existing methods in terms of robustness and accuracy and, on the other, to extend the upwind constrained transport (UCT) method by Londrillo & Del Zanna (2004) [22] and Del Zanna et al. (2007) [23] for the systematic construction of new averaging schemes. In particular, we propose a general formula for the upwind fluxes of the induction equation which simply involves the information available from the base Riemann solver employed for the fluid part, provided it does not require full spectral decomposition, and 1D reconstructions of velocity and magnetic field components from nearby intercell faces to cell edges. Our results are presented here in the context of second-order schemes for classical MHD, but they can be easily generalized to higher than second order schemes, either based on finite volumes or finite differences, and to other physical systems retaining the same structure of the equations, such as that of relativistic or general relativistic 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 B 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 ρ, ρv and E):Ut+F=0,whereU=(ρρvE),F=(ρvρvvBB+Ipt(E+pt)v(vB)B), and an induction equation for the evolution of the magnetic fieldBt+×E=0, where the curl operator appears instead of a divergence and where the electric field E=v×B is not a primary variable, depending on the

CT schemes based on emf averaging

The set of conserved variables Uc 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)

  • D.S. Balsara

    Multidimensional Riemann problem with self-similar internal structure. Part I - application to hyperbolic conservation laws on structured meshes

    J. Comput. Phys.

    (2014)
  • D.S. Balsara et al.

    Multidimensional Riemann problem with self-similar internal structure. Part II - application to hyperbolic conservation laws on unstructured meshes

    J. Comput. Phys.

    (2015)
  • D.S. Balsara et al.

    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.

    (2017)
  • T.A. Gardiner et al.

    An unsplit Godunov method for ideal MHD via constrained transport

    J. Comput. Phys.

    (2005)
  • P. Londrillo et al.

    On the divergence-free condition in Godunov-type schemes for ideal magnetohydrodynamics: the upwind constrained transport method

    J. Comput. Phys.

    (2004)
  • D.S. Balsara

    Divergence-free reconstruction of magnetic fields and WENO schemes for magnetohydrodynamics

    J. Comput. Phys.

    (2009)
  • D.S. Balsara et al.

    Efficient, high accuracy ADER-WENO schemes for hydrodynamics and divergence-free magnetohydrodynamics

    J. Comput. Phys.

    (2009)
  • K.G. Felker et al.

    A fourth-order accurate finite volume method for ideal MHD via upwind constrained transport

    J. Comput. Phys.

    (2018)
  • T. Miyoshi et al.

    A multi-state HLL approximate Riemann solver for ideal magnetohydrodynamics

    J. Comput. Phys.

    (2005)
  • E.F. Toro et al.

    MUSTA fluxes for systems of conservation laws

    J. Comput. Phys.

    (2006)
  • S. Li

    An HLLC Riemann solver for magneto-hydrodynamics

    J. Comput. Phys.

    (2005)
  • D.S. Balsara et al.

    Von Neumann stability analysis of globally divergence-free RKDG schemes for the induction equation using multidimensional Riemann solvers

    J. Comput. Phys.

    (2017)
  • M. Dumbser et al.

    A new efficient formulation of the HLLEM Riemann solver for general conservative and non-conservative hyperbolic systems

    J. Comput. Phys.

    (2016)
  • A. Mignone

    A simple and accurate Riemann solver for isothermal MHD

    J. Comput. Phys.

    (2007)
  • A. Suresh et al.

    Accurate monotonicity-preserving schemes with Runge Kutta time stepping

    J. Comput. Phys.

    (1997)
  • P. Cargo et al.

    Roe matrices for ideal MHD and systematic construction of roe matrices for systems of conservation laws

    J. Comput. Phys.

    (1997)
  • T.A. Gardiner et al.

    An unsplit Godunov method for ideal MHD via constrained transport in three dimensions

    J. Comput. Phys.

    (2008)
  • A. Mignone et al.

    A second-order unsplit Godunov scheme for cell-centered MHD: the CTU-GLM scheme

    J. Comput. Phys.

    (2010)
  • D. Lee et al.

    An unsplit staggered mesh scheme for multidimensional magnetohydrodynamics

    J. Comput. Phys.

    (2009)
  • A. Mignone et al.

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

    J. Comput. Phys.

    (2010)
  • K. Waagan et al.

    A robust numerical scheme for highly compressible magnetohydrodynamics: nonlinear stability, implementation and tests

    J. Comput. Phys.

    (2011)
  • Cited by (24)

    View all citing articles on Scopus
    View full text