MHD

Introduction

Castro implements a constrained transport (CT) corner transport upwind (CTU) ideal MHD scheme based on the work of Miniati & Martin [7]. MHD is enabled by compiling with USE_MHD = TRUE. This replaces the pure hydrodynamics solver, but uses the same driver as the CTU hydrodynamics solver. This means that all of the source terms supported by the hydrodynamics solver are also supported by MHD.

Note

The MHD solver is still under development and should not be used for science simulations. See the issues in the github repo.

The MHD solver supports 3-d only.

Currently the MHD solver is single-level only. AMR support is forthcoming.

Equations and Data Structures

The ideal MHD equations we solve appear as:

\begin{array}{r} \begin{aligned} \frac{\partial ρ X_{k}}{\partial t} + \frac{\partial}{\partial x_{j}} (ρ U_{j} X_{k}) & = ρ {\dot{ω}}_{k} \\ \frac{\partial ρ U_{j}}{\partial t} + \frac{\partial}{\partial x_{j}} (ρ U_{i} U_{j} + p δ_{i j} - B_{j} B_{i}) & = S_{ρ U_{j}} \\ \frac{\partial ρ E}{\partial t} + \frac{\partial}{\partial x_{j}} [U_{j} (ρ E + p) - B_{j} B_{i} u_{i}] & = S_{ρ E} + ρ H_{nuc} \\ \frac{\partial B_{i}}{\partial t} & = - \frac{\partial}{\partial x_{j}} (U_{j} B_{i} - B_{j} U_{i}) \end{aligned} \end{array}

where

p = p_{g} + \frac{1}{2} B^{2}

and the $S$ sources represent the hydrodynamical sources and the remainder are reaction sources. The MHD solver uses the same time-advancement driver as the hydrodynamic CTU driver (Strang+CTU Evolution), using Strang splitting for the reactions (by default).

Note

Note: we are following the convention in the MHD community of setting the permeability to unity. This is why the magnetic pressure has the form $\frac{1}{2} B^{2}$ instead of $\frac{1}{8 π} B^{2}$ . In effect, we are carrying our magnetic field as $B^{'} = B / \sqrt{4 π}$ .

If you wish to express the magnetic field in Gauss, then you will need to multiply the value Castro carries by $\sqrt{4 π}$ .

The constrained transport algorithm algebraically ensures that $\nabla \cdot B = 0$ . Throughout the algorithm, the magnetic fields are face-centered, the electric fields are edge-centered, and the conserved state is cell-centered (or averages), unless specified otherwise. We use the following indexing conventions:

U(i,j,k) is the cell-centered state, $U_{i, j, k}$

qleft(i,j,k) is the interface value of a state, for example, in the x-direction this would be $q_{i - 1 / 2, j, k}$

Bx(i,j,k) is the face-centered x-component of the magnetic field, $B_{x, i - 1 / 2, j, k}$

Ex(i,j,k) is the edge-centered x-component of the electric field, $E_{x, i, j - 1 / 2, k - 1 / 2}$

Problem Initialization

There is an additional initialization routine for MHD, problem_initialize_mhd_data(), that is used to initialize the face-centered magnetic field components. This is done separately from the main conserved fluid state.

The conserved fluid state is initialized in problem_initialize_state_data() just as with pure hydrodynamics problems. Note that you do not need to include the magnetic energy contribution to the total energy density, UEDEN. After this initialization, the driver handles the addition of the magnetic contribution.

Hydrodynamics Update

We use piecewise linear or piecewise parabolic reconstruction with a characteristic projection and the full 12-Riemann solve corner transport upwind method. The HLLD Riemann solver is used.

Within the solver, we use the same indexing into the primitive state as defined in Primitive State Vector Integer Keys, with the additions of QMAGX, QMAGY, and QMAGZ for the cell-centered magnetic field components and QPTOT for the total pressure (gas + magnetic).

Just like with pure hydrodynamics, the reconstruction type is controlled by castro.ppm_type, with 0 selecting piecewise linear and 1 selecting piecewise parabolic.

For the piecewise linear method, the slope limiting is controlled by the same plm_iorder runtime parameter as for hydrodynamics. Additionally, we have an option mhd_limit_characteristic that controls whether you want to do the slope limiting on the characteristic variables (the default) or the primitive variables.

Electric Update

Coupled to the hydrodynamics update is the electric field update. Here we update the components of E using the contact upwind scheme first proposed in [40]. The updated electric field then gives the magnetic field via Faraday’s law and the discretization ensures that $\nabla \cdot B = 0$ .