CONSTRAINED-TRANSPORT MAGNETOHYDRODYNAMICS WITH ADAPTIVE MESH REFINEMENT IN CHARM

The algorithmic operations are carried out on a discrete representation of the continuous D-dimensional space given by the cubic lattice, ${{\boldsymbol{i}}}\equiv (i_0, \ldots, i_{{\mathbf {D}}- 1}) \in \mathbb {Z}^{\mathbf {D}}$. The computational domain, referred to as a grid Γ, is a bound subset of $\mathbb {Z}^{\mathbf {D}}$ and provides a finite-volume discretization of the continuous space into a collection of control volumes, faces, and edges. Each control volume is identified by an index ${{\boldsymbol{i}}}\equiv (i_0, \ldots, i_{{\mathbf {D}}- 1}) \in \Gamma$ and corresponds to a region of space,

$$\begin{equation} V_{{\boldsymbol{i}}} = [{\boldsymbol{i}}h, ({\boldsymbol{i}}+ {\boldsymbol{v}})h], \end{equation} \tag{ 5 }$$

where h is the mesh spacing and ${\boldsymbol{v}}\equiv (1,\ldots,1)$ is the vector whose components are all equal to one. The face-centered discretization of space based on the same control volumes is $\Gamma _{\rm f}^{{{\boldsymbol{e}}^{d}}} = \lbrace {\boldsymbol{i}}\pm \frac{1}{2}{{\boldsymbol{e}}^{d}}: {\boldsymbol{i}}\in \Gamma \rbrace$, where ${{\boldsymbol{e}}^{d}}$ is the unit vector in the d direction. $\Gamma _{\rm f}^{{{\boldsymbol{e}}^{d}}}$ indexes the cell faces of Γ normal to ${{\boldsymbol{e}}^{d}}$ representing the areas

$$\begin{equation} A_{{\boldsymbol{i}}+ \frac{1}{2}{{\boldsymbol{e}}^{d}}} = [({\boldsymbol{i}}+ {{\boldsymbol{e}}^{d}})h, ({\boldsymbol{i}}+ {\boldsymbol{v}})h], \quad {\boldsymbol{i}}+ \frac{1}{2}{{\boldsymbol{e}}^{d}}\in \Gamma _{\rm f}^{{{\boldsymbol{e}}^{d}}}. \end{equation} \tag{ 6 }$$

Finally, the edge-centered discretization of space is $\Gamma _{\rm e}^{{{\boldsymbol{e}}^{d}}} = \lbrace {\boldsymbol{i}}\pm \frac{1}{2}{{\boldsymbol{e}}^{d_1}}\pm \frac{1}{2}{{\boldsymbol{e}}^{d_2}}: {\boldsymbol{i}}\in \Gamma,\,d\ne d_1\ne d_2\rbrace$. $\Gamma _{\rm e}^{{{\boldsymbol{e}}^{d}}}$ indexes the edges of the cells in Γ aligned with ${{\boldsymbol{e}}^{d}}$ representing the lengths

$$\begin{equation} L_{{\boldsymbol{i}}+ \frac{1}{2}({{\boldsymbol{e}}^{d_1}}+{{\boldsymbol{e}}^{d_2}})} = [{\boldsymbol{i}}h+\frac{1}{2}({{\boldsymbol{e}}^{d_1}}+{{\boldsymbol{e}}^{d_2}}), ({\boldsymbol{i}}+ {\boldsymbol{v}})h],\; {\boldsymbol{i}}+ \frac{1}{2}({{\boldsymbol{e}}^{d_1}}+{{\boldsymbol{e}}^{d_2}}) \in \Gamma _{\rm e}^{{{\boldsymbol{e}}^{d}}}. \end{equation} \tag{ 7 }$$

Figure 1 illustrates a control volume with the various types of discretization.

Zoom In Zoom Out Reset image size

Given the above discretization, we define cell-centered discrete variables on Γ as

$$\begin{equation*} \phi : \Gamma \rightarrow {\mathbb {R}}^m, \end{equation*}$$

and denote by $\phi _{\boldsymbol{i}}\in {\mathbb {R}}^m$ the value of ϕ at cell ${\boldsymbol{i}}\in \Gamma$. Similarly, we define face-centered vector fields on $\Gamma _{\rm f}^{{{\boldsymbol{e}}^{d}}}$ as

$$\begin{equation*} {\vec{F}} = (F_0,\dots,F_{{\mathbf {D}}-1}) \hbox{, } F_d : \Gamma _{\rm f}^{{\boldsymbol{e}}^d} \rightarrow {\mathbb {R}}^m, \end{equation*}$$

and denote by $F_{d,{\boldsymbol{i}}+ \frac{1}{2}{{\boldsymbol{e}}^{d}}} \in {\mathbb {R}}^m$ the value of F_d at ${\boldsymbol{i}}+ \frac{1}{2}{{\boldsymbol{e}}^{d}}\in \Gamma _{\rm f}^{{{\boldsymbol{e}}^{d}}}$, and also define edge-centered vector fields on $\Gamma _{\rm e}^{{{\boldsymbol{e}}^{d}}}$ as

$$\begin{equation*} {\vec{E}} = (E_0,\dots, E_{{\mathbf {D}}-1}) \hbox{, } E_d : \Gamma _{\rm e}^{{\boldsymbol{e}}^d} \rightarrow {\mathbb {R}}^m \end{equation*}$$

and denote by $E_{d,{\boldsymbol{i}}+ \frac{1}{2}({{\boldsymbol{e}}^{d_1}}+{{\boldsymbol{e}}^{d_2}})} \in {\mathbb {R}}^m$ the value of E_d at ${\boldsymbol{i}}+ \frac{1}{2}({{\boldsymbol{e}}^{d_1}}+{{\boldsymbol{e}}^{d_2}}) \in \Gamma _{\rm e}^{{{\boldsymbol{e}}^{d}}}.$ Finally, for face-centered fields we introduce the discretized divergence operator

$$\begin{equation} (\vec{D}\cdot \vec{F})_{\boldsymbol{i}}= \frac{1}{h} \sum ^{{\mathbf {D}}-1}_{d=0} \big(F_{d,{\boldsymbol{i}}+\frac{1}{2}{\boldsymbol{e}}^d} - F_{d,{\boldsymbol{i}}- \frac{1}{2}{\boldsymbol{e}}^d}\big),\quad {\boldsymbol{i}}\in \Gamma, \end{equation} \tag{ 8 }$$

and for edge-centered fields the discretized curl operator

$$\begin{eqnarray} (\vec{D}\times \vec{E})_{d,{\boldsymbol{i}}+\frac{1}{2}{{\boldsymbol{e}}^{d}}} &=& \frac{1}{h}\! \sum _{d_1, d_2}\! \varepsilon _{dd_1d_2} \big(E_{d_1,{\boldsymbol{i}}+\frac{1}{2}{\boldsymbol{e}}^d+\frac{1}{2}{{\boldsymbol{e}}^{d_2}}} - E_{d_1,{\boldsymbol{i}}+\frac{1}{2}{\boldsymbol{e}}^d-\frac{1}{2}{{\boldsymbol{e}}^{d_2}}}\big),\nonumber\\ &&{\boldsymbol{i}}+ \frac{1}{2}{{\boldsymbol{e}}^{d}}\in \Gamma _{\rm f}^{{{\boldsymbol{e}}^{d}}}. \end{eqnarray} \tag{ 9 }$$

Time is also discretized into a number of finite intervals of variable size. In particular, we write t^{n + 1} = tⁿ + Δtⁿ, where t indicates the solution time, n indicates the integration step, and Δtⁿ indicates the time-step interval.

The coupled system of Equations (1)–(4), with the allowance for a non-zero source term S(U), can be cast in the following compact form:

$$\begin{equation} \frac{\partial U}{\partial t} + {\bf \nabla \cdot F} = S, \end{equation} \tag{ 10 }$$

where the conservative variables, U, and the associated fluxes along the direction d, F_d, are defined as

$$\begin{eqnarray} U &=& {\left(\begin{array}{@{}c@{}}\rho \\ \rho u_0 \\ \vdots \\ \rho u_{{\mathbf {D}}-1} \\ \rho e \\ B_0 \\ \vdots \\ B_{{\mathbf {D}}-1} \end{array}\right)}, \nonumber\\ F_d(U) &=& {\left(\begin{array}{@{}c@{}}\rho u_d \\ \rho u_0 u_d + p \,\delta _{d0} -B_0B_d \\ \vdots \\ \rho u_{{\mathbf {D}}-1} u_d + p\,\delta _{d{\mathbf {D}}-1} -B_{{\mathbf {D}}-1}B_d \\ u_d\left(\rho e+p\right)-B_d{\bf B}\,{\bm \cdot}\, {\bf u} \\ u_0B_d-B_0u_d \\ \vdots \\ u_{{\mathbf {D}}-1}B_d-B_{{\mathbf {D}}-1}u_d \end{array}\right)}. \end{eqnarray} \tag{ 11 }$$

Following Godunov's approach and its higher-order extensions, we can conveniently formulate a numerical integration scheme based on the conservative properties of the system (10). In this approach one follows the evolution of cell-centered volume-averaged conservative variables, defined as

$$\begin{equation} U^n_{\boldsymbol{i}}= \frac{1}{V_{\boldsymbol{i}}} \int _{V_{\boldsymbol{i}}} U(t^n,{\boldsymbol{x}})\,{\it dV}. \end{equation} \tag{ 12 }$$

The evolution equation is obtained upon suitable manipulation of the original continuous Equation (10), and reads (Leveque 1998)

$$\begin{equation} U_{\boldsymbol{i}}^{n+1} = U_{\boldsymbol{i}}^n- \frac{\Delta t}{\Delta x}\, \sum ^{{\mathbf {D}}-1}_{d=0} \left(F^{n+\frac{1}{2}}_{d,{\boldsymbol{i}}+\frac{1}{2}{\boldsymbol{e}}^d} - F^{n+\frac{1}{2}}_{d,{\boldsymbol{i}}- \frac{1}{2}{\boldsymbol{e}}^d}\right) + S^{n+\frac{1}{2}}, \end{equation} \tag{ 13 }$$

where the face-centered time-averaged fluxes along the direction d are defined as

$$\begin{equation} F^{n+\frac{1}{2}}_{d,{\boldsymbol{i}}+\frac{1}{2}{\boldsymbol{e}}^d} = \frac{1}{\Delta t A_{{\boldsymbol{i}}+\frac{1}{2}{\boldsymbol{e}}^d}} \int _{t^n}^{t^{n+1}}dt\int _{A_{{\boldsymbol{i}}+\frac{1}{2}{\boldsymbol{e}}^d}} F_d(t,{\boldsymbol{x}}) \,{\it dA}. \end{equation} \tag{ 14 }$$

If the source term is non-stiff, we can obtain a second-order estimate for it by the simple time average, $S^{n+\frac{1}{2}}\simeq \frac{1}{2}(S^n+S^{n+1})$.

Note that given the flux along a direction d₁, $F_{d_1,{\boldsymbol{i}}+\frac{1}{2}{{\boldsymbol{e}}^{d_1}}}$, the component corresponding to the magnetic field along the direction d₂ ≠ d₁, defines a face-centered electric field along d according to

$$\begin{equation} E_{d,{\boldsymbol{i}}+\frac{1}{2}{{\boldsymbol{e}}^{d_1}}}=-\varepsilon _{dd_1d_2}\,F_{d_1,{\boldsymbol{i}}+\frac{1}{2}{{\boldsymbol{e}}^{d_1}}} ^{B_{d_2}}(U), \end{equation} \tag{ 15 }$$

where we use subscripts to indicate directions and centering, and superscripts for components. Indeed, Godunov's scheme also updates in time the cell-centered magnetic field variables. As already pointed out, however, the updated magnetic field in general does not remain solenoidal. Therefore, we adopt instead a CT discretization strategy in which the primary description of the magnetic field we will be using face-centered area-averaged variables defined as

$$\begin{equation} B^n_{d,{\boldsymbol{i}}+\frac{1}{2}{\boldsymbol{e}}^d} = \frac{1}{A_{{\boldsymbol{i}}+\frac{1}{2}{\boldsymbol{e}}^d}} \int _{A_{{\boldsymbol{i}}+\frac{1}{2}{\boldsymbol{e}}^d}} B_d(t^n,{\boldsymbol{x}}) \,{\it dA}. \end{equation} \tag{ 16 }$$

Note that an estimate of the cell-centered magnetic field variables is still needed in order to construct the fluxes in Equation (14). We will return to this point shortly. The evolution equation for the face-centered magnetic field variables is obtained again from a manipulation of Faraday's law and reads

$$\begin{eqnarray} B_{d,{\boldsymbol{i}}+\frac{1}{2}{\boldsymbol{e}}^d}^{n+1} &=& B_{d,{\boldsymbol{i}}+\frac{1}{2}{\boldsymbol{e}}^d}^n -\varepsilon _{dd_1d_2} \frac{\Delta t}{\Delta x}\nonumber\\ &&\times\big[ \big(E^{n+\frac{1}{2}}_{d_1,{\boldsymbol{i}}+\frac{1}{2}{\boldsymbol{e}}^d+\frac{1}{2}{{\boldsymbol{e}}^{d_2}}} - E^{n+\frac{1}{2}}_{d_1,{\boldsymbol{i}}+\frac{1}{2}{\boldsymbol{e}}^d-\frac{1}{2}{{\boldsymbol{e}}^{d_2}}}\big)\nonumber \\ && -\big(E^{n+\frac{1}{2}}_{d_2,{\boldsymbol{i}}+\frac{1}{2}{\boldsymbol{e}}^d+\frac{1}{2}{{\boldsymbol{e}}^{d_1}}} - E^{n+\frac{1}{2}}_{d_2,{\boldsymbol{i}}+\frac{1}{2}{\boldsymbol{e}}^d-\frac{1}{2}{{\boldsymbol{e}}^{d_1}}}\big)\big] \nonumber \\ &&+ S_{B_d,{\boldsymbol{i}}+\frac{1}{2}{\boldsymbol{e}}^d}^{n+\frac{1}{2}} \end{eqnarray} \tag{ 17 }$$

with d ≠ d₁ ≠ d₂, 0 ⩽ d, d₁, d₂ < D. The edge-centered time-averaged electric field is formally defined as

$$\begin{equation} E^{n+\frac{1}{2}}_{d,{\boldsymbol{i}}+\frac{1}{2}{{\boldsymbol{e}}^{d_1}}+\frac{1}{2}{{\boldsymbol{e}}^{d_2}}} = \frac{1}{\Delta t\,L_{{\boldsymbol{i}}+ \frac{1}{2}({{\boldsymbol{e}}^{d_1}}+{{\boldsymbol{e}}^{d_2}})}} \int _{t^n}^{t^{n+1}}\!\!dt\int _{L_{{\boldsymbol{i}}+ \frac{1}{2}({{\boldsymbol{e}}^{d_1}}+{{\boldsymbol{e}}^{d_2}})}}\!\!\!\! E_d(t,{\boldsymbol{x}}) \,{\it dL}, \end{equation} \tag{ 18 }$$

and, as above, the time-centered estimate of the non-stiff source term for the face-centered magnetic field components is obtained by simple arithmetic averaging. The cell-centered magnetic field variables are then defined in terms of the primary face-centered values using the following second-order accurate reconstruction scheme (Ryu et al. 1998; Balsara 2001; Gardiner & Stone 2005):

$$\begin{equation} B_{d,{\boldsymbol{i}}} = \frac{1}{2}\big(B_{d,{\boldsymbol{i}}-\frac{1}{2}{\boldsymbol{e}}^d}+B_{d,{\boldsymbol{i}}+\frac{1}{2}{\boldsymbol{e}}^d} \big). \end{equation} \tag{ 19 }$$

An important part of CT schemes concerns the calculation of the time-averaged, edge-centered electric fields entering the time update (Equation (17)). For example, one can employ some type of bilinear interpolation. Dai & Woodward (1998) interpolate in space and time the cell-centered velocity and magnetic field at time tⁿ with the solution from the Godunov scheme at time t^{n + 1} to obtain time- and corner-centered electric fields. On the other hand, Ryu et al. (1998) and Balsara (2001) take advantage of Equation (15) and interpolate the face-centered variables that allow reconstruction of electric fields at cell edges. The interpolation scheme proposed in Ryu et al. (1998) has the property that for plane-parallel configurations their multidimensional scheme reduces to the 1D scheme. This same property is shared by the upwind scheme proposed in Gardiner & Stone (2005) which is adopted here and is described in more detail in Section 2.2.6. This property is important because it guarantees self-consistency between (1) the electric fields used for the time update of face-centered magnetic variables, (2) the MHD fluxes used for the time update of the cell-centered variables, and (3) the synchronization step of the magnetic variables given in Equation (19) (Gardiner & Stone 2005). In the most general case, the electric field is obtained by solving a 2D Riemann problem (Londrillo & Del Zanna 2000; Londrillo & del Zanna 2004; Teyssier et al. 2006; Fromang et al. 2006), of which the upwind scheme mentioned above is a special case.

Note that by applying the divergence operator in Equation (8) to the face-centered magnetic field evolved according to Equation (17), one finds that the divergence of the field does not change in time. So the magnetic field remains solenoidal if initially so. This is the property of the CT scheme.

The accuracy and stability of the numerical solution depend principally on how the time-averaged fluxes and electric fields entering Equations (13) and (17), respectively, are computed. In the following sections, we shall describe the algorithmic details characterizing such calculation.

The first part of the algorithm consists in reconstructing the numerical solution within the spatial domain of the mesh. This requires estimating gradients and eventually extrapolating state variables in time and space from cell centers to cell faces with the use of characteristic tracing. Such operations are usually done in primitive space, where the variables are defined as

$$\begin{equation} W=(\rho,u_0,\dots,u_{{\mathbf {D}}-1},p_g,B_0,\dots,B_{{\mathbf {D}}-1})^T, \end{equation} \tag{ 21 }$$

as it simplifies the characteristic analysis of the system. The evolution of the system in terms of primitive variables is given by a system of equations in non-conservative form,

$$\begin{equation} \frac{\partial W}{\partial t} + \sum ^{{\mathbf {D}}-1}_{d=0} A_d(W) \frac{\partial W_d}{\partial x_d} = S^W, \end{equation} \tag{ 22 }$$

where S^W = ∇_UW · S and A_d = ∇_UW · ∇_UF_d · ∇_WU, and ∇_UW and ∇_UW are the Jacobian of the transformation from primitive to conserved variables and vice versa. Given their importance for the construction of a sound numerical scheme, the properties of the operators A_d have been studied in great detail in the literature (e.g., Brio & Wu 1988; Roe & Balsara 1996). Here it suffices to make the following observations. In 1D, d = 0, the parallel component of the magnetic field is constant, B₀ = const., and A₀ effectively becomes a 7 × 7 matrix. In general the operator A₀ is characterized by seven eigenvalues, and associated left and right eigenvectors, with values u, u ± c_s, u ± c_A, and u ± c_f obeying the hierarchy: c_s ⩽ c_A ⩽ c_f. The first eigenvalue listed above corresponds to the usual entropy wave, and the other six to the three pairs of MHD waves (slow magnetosonic, Alfvén, and fast magnetosonic) propagating downstream (+) or upstream (−) in the flow, respectively. Because up to five of the eigenvalues may actually coincide, the system is not strictly hyperbolic, so care must be taken to avoid singularities with expressions involving the eigenvectors (Brio & Wu 1988; Roe & Balsara 1996). Finally, in more than 1D, because B_d is not affected by gradients in the d direction, it has been customary to use the 1D analysis when formulating the predictor step for higher-order Godunov-like MHD algorithms. However, that leads to neglect of terms ∝∂B_d/∂x_d, which are not necessarily null in multidimensions. In fact, it is important to include these terms to avoid degrading the solution accuracy, as recently pointed out by Gardiner & Stone (2005).

After the computing primitive from conservative variables, $W^n_{\boldsymbol{i}}=W^n_{\boldsymbol{i}}(U^n_{\boldsymbol{i}})$ and the synchronization of cell-centered and face-centered magnetic fields according to Equation (19), we proceed to the calculation of the slopes along each direction d as follows. First, central and side slopes are estimated, respectively, as

$$\begin{eqnarray} \hspace{2pc}\delta _{d,0} W_{\boldsymbol{i}}= & \frac{1}{2}\big(W^n_{{\boldsymbol{i}}+ {\boldsymbol{e}}^d} - W^n_{{\boldsymbol{i}}- {\boldsymbol{e}}^d}\big), \end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray} \delta _{d,-} W_{\boldsymbol{i}}= & W^n_{\boldsymbol{i}}- W^n_{{\boldsymbol{i}}- {\boldsymbol{e}}^d}, \end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray} \delta _{d,+} W_{\boldsymbol{i}}= & W^n_{{\boldsymbol{i}}+ {\boldsymbol{e}}^d} - W^n_{\boldsymbol{i}}, \end{eqnarray} \tag{ 25 }$$

and then limiting is applied component-wise either in primitive or in characteristic space. We use van Leer's (also known as Monotonized Central) limiter defined as

$$\begin{equation*} \delta ^{vL}(\delta _0,\delta _-,\delta _+) = \left\lbrace \begin{array}{@{}l@{\;\;}l@{}}\mbox{sgn}(\delta _0)\; \min (|\delta _0|,2|\delta _-|,2|\delta _+|) & \mbox{if $\delta _-\delta _+ > 0$} \\\ms 0 & {\rm otherwise}. \end{array}\right. \end{equation*}$$

In the case of primitive limiting the limiter is applied directly to each component k of the primitive slopes (23)–(25), i.e.,

$$\begin{equation} \delta _d W_{\boldsymbol{i}}^k=\delta ^{vL}\big(\delta _{d,0}W_{\boldsymbol{i}}^k,\delta _{d,-}W_{\boldsymbol{i}}^k,\delta _{d,+} W_{\boldsymbol{i}}^k\big). \end{equation} \tag{ 26 }$$

In the case of characteristic limiting the limiter is applied to the components of the primitive slopes in characteristic space, namely,

$$\begin{eqnarray} \hspace{-6.5pc}\delta _d W_{\boldsymbol{i}}&=& \sum _k \alpha ^kr^k, \end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray} \hspace{-3.5pc}\alpha ^k&=&\delta ^{vL}\big(\alpha ^k_0,\alpha ^k_-,\alpha ^k_+\big), \end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray} \alpha ^k_\# &=& l^k \cdot \delta _{d,\#} W_{\boldsymbol{i}},\quad \#=0,+,-, \end{eqnarray} \tag{ 29 }$$

where $l^k= l^k(W^n_{\boldsymbol{i}})$ and $r^k = r^k(W^n_{\boldsymbol{i}})$ are the left and right eigenvectors of the operator $A_{{\boldsymbol{i}},d}$.

Next we extrapolate the primitive variables from cell centers to face centers along each coordinate direction, d, by taking into account the time-averaged effect due to the slopes just computed. This is done most conveniently by using the 1D versions of the MHD equations in primitive form. There is no evolution in the d-component of the magnetic field due to derivatives along the d-direction, so the normal components of the magnetic field on cell faces are straightforwardly provided by the face-centered component of the magnetic field, without need for even geometrical extrapolation. On the other hand, as pointed out by Gardiner & Stone (2005) in multidimensional MHD the reconstruction step must include terms proportional to ∂B_d/∂x_d terms (no summation over repeated indices is implied here). These terms arise from the requirement to balance the divergence terms in more than 1D and their neglect can cause serious degradation of the numerical solution. However, here the source terms are added only during the normal predictor step and are omitted during the corrections associated with transverse flux gradients. Our tests suggest that neither the accuracy nor the robustness of the code is affected by this choice. The reconstruction of the primitive variables onto cell faces can be done with various degrees of accuracy. Typically, one uses a piecewise-linear (PLM) or piecewise-parabolic (PPM) reconstruction scheme (Colella & Woodward 1984). Although we mostly use a PPM algorithm for the tests presented here, for simplicity we illustrate the case of PLM reconstruction. So in this case the extrapolation of the primitive variables in space and time from cell centers to faces along the direction d takes the form

$$\begin{eqnarray} &&{\left(\begin{array}{@{}c@{}}\hat{W}^{n}_{{\boldsymbol{i}}, \pm, d} \\ B^{n}_{d,{\boldsymbol{i}},\pm,d} \end{array}\right)} &=& \; {\left(\begin{array}{@{}c@{}}\hat{W}^n_{\boldsymbol{i}}\\ B^{n}_{d,{\boldsymbol{i}}+{{\boldsymbol{e}}^{d}}} \end{array}\right)} +\frac{1}{2}\left[ {\left(\begin{array}{@{}cc@{}}\pm I & 0 \\ 0 & 0 \end{array}\right)}\right.\nonumber\\ &&\left.-\,\frac{\Delta t}{\Delta x} {\left(\begin{array}{@{}c@{}}\hat{A}_{{\boldsymbol{i}},d} \\ 0 \end{array}\right)} P_\pm (\delta _d \hat{W}^n_{\boldsymbol{i}}) \right] +\frac{\Delta t}{2} {\left(\begin{array}{@{}c@{}}S_{d,{\boldsymbol{i}}}^{n,\rm MHD} \\\ms 0 \end{array}\right)},\qquad \end{eqnarray} \tag{ 30 }$$

where we have explicitly separated out the reconstruction of the normal component of the magnetic field (the hat symbol in the notation indicates that the components corresponding to the normal magnetic field are omitted). In addition, we have used the projector operator defined as

$$\begin{equation} P_\pm (W) = \sum _{\pm \lambda _k > 0} (l_k \cdot W) r_k, \end{equation} \tag{ 31 }$$

where λ_k are eigenvalues of $A_{{\boldsymbol{i}},d}$. This projector operator filters out the components of the gradients that propagate away from the cell interface. However, when a Riemann solver of the Harten–Lax–van Leer (HLL) family (Harten et al. 1983) is employed, in order to obtain second-order accuracy the filter is switched off and both in the PPM and PLM cases, and the summation is carried over all waves, irrespective of their sign (this is further discussed in Section 2.2.5). Finally, S_{d, MHD} represents the MHD source term required in multidimensional MHD which we implement in the form (Crockett et al. 2005)

$$\begin{equation} S_{d,{\boldsymbol{i}}}^{n,\rm MHD}= {\left(\begin{array}{@{}c@{}}0\\ \frac{B_{0}}{\rho } \\ \vdots \\ \frac{B_{{\mathbf {D}}-1}}{\rho } \\ {\bf B}\,{\bm \cdot}\, {\bm u}\\ u_{d_1} \\ u_{d_2} \end{array}\right)}_{\boldsymbol{i}}^n \left(\frac{\partial B_d}{\partial x_d}\right)_{\boldsymbol{i}}^n, \end{equation} \tag{ 32 }$$

where

$$\begin{equation} \left(\frac{\partial B_d}{\partial x_d}\right)_{\boldsymbol{i}}^n = \frac{B_{d,{\boldsymbol{i}}+\frac{1}{2}e^d}^n-B_{d,{\boldsymbol{i}}-\frac{1}{2}e^d}^n}{\Delta x}, \end{equation} \tag{ 33 }$$

and, as usual, d ≠ d₁ ≠ d₂ and, in this particular case, 0 ⩽ d₁ < d₂ < D. The normal predictor step is completed by the final corrections for a non-stiff source term according to

$$\begin{equation} W^{n}_{{\boldsymbol{i}}, \pm, d}= W^{n}_{{\boldsymbol{i}}, \pm, d} +\frac{\Delta t}{2} S^{n}_{\boldsymbol{i}}. \end{equation} \tag{ 34 }$$

After the primitive variables have been extrapolated to cell faces, we add corrections due to gradients parallel to the cell faces. We find it convenient at this point to convert back to a conservative representation, thus we compute $U^{n}_{{\boldsymbol{i}},\pm,d}=U^{n}_{{\boldsymbol{i}},\pm,d}(W^{n}_{{\boldsymbol{i}},\pm,d})$. Following the CTU scheme the corrections are expressed in terms of transverse flux gradients, with fluxes obtained from a Riemann solver, ${\cal R}$. In accord with the CT scheme, however, for the update of magnetic field variables we use gradients of edge-centered electric fields suitably interpolated in space and time from their face- and cell-centered values. Both the interpolation procedure, ${\cal I}_E$, and the Riemann solver, ${\cal R}$, will be specified at the end of this section. It suffices here to say that the time centering and interpolation accuracy of the interpolated edge-centered electric fields is consistent with that of the MHD fluxes returned by the Riemann solver. Note that no additional "multidimensional MHD source term," as required by a straightforward 3D extension of Gardiner & Stone (2005), appears in this part of the algorithm, which considerably simplifies it.

The steps involved in the modified CTU update can then be summarized as follows.

• 1.  
  Use a Riemann solver, ${\cal R}(U^{\rm Left},U^{\rm Right},d)$, to obtain a first estimate of the fluxes across cell faces along each direction, d

  $$\begin{equation} F^{{\rm 1D}}_{d,{\boldsymbol{i}}+ \frac{1}{2}{{\boldsymbol{e}}^{d}}} = {\cal R} \big(U^n_{{\boldsymbol{i}}, +, d}, U^n_{{\boldsymbol{i}}+ {{\boldsymbol{e}}^{d}}, -, d},d\big), \end{equation} \tag{ 35 }$$

  where the "1D" notation loosely indicates face-centered quantities corrected for the one-dimensional gradients during the normal predictor step.
• 2.  
  Use the newly obtained fluxes with the primitive solution at time tⁿ to interpolate the electric fields from cell faces and cell centers to cell edges

  $$\begin{equation} E^{{\rm 1D}}_{d,{\boldsymbol{i}}+\frac{1}{2}{{\boldsymbol{e}}^{d_1}}+ \frac{1}{2}{{\boldsymbol{e}}^{d_2}}} = {\cal I}_E\big(F^{{\rm 1D}}_{d_1,* + \frac{1}{2}{{\boldsymbol{e}}^{d_1}}}, F^{{\rm 1D}}_{d_2,*+\frac{1}{2}{{\boldsymbol{e}}^{d_2}}},W^{n}_{*}\big), \end{equation} \tag{ 36 }$$

  where the "*" symbol indicates that the interpolation requires values of the arguments at various cell centers and faces.
• 3.  
  As part of the CTU prescription to obtain (1, 1, 1) diagonal coupling, apply corrections to density, momentum, and energy components of $U_{{\boldsymbol{i}}, \pm, d_1}$, due to one set of transverse flux derivatives along d₂. For each face d₁, there will be D  −  1 such corrected states, one for each direction perpendicular to d₁. We indicate these states with a "2D" notation and write

  $$\begin{equation} U^{{\rm 2D}} _{{\boldsymbol{i}}, \pm, d_1, d_2} = U^n_{{\boldsymbol{i}}, \pm, d_1} - \frac{\Delta t}{3\Delta x} \big(F^{{\rm 1D}}_{d_2,{\boldsymbol{i}}+ \frac{1}{2}{{\boldsymbol{e}}^{d_2}}} - F^{{\rm 1D}}_{d_2,{\boldsymbol{i}}- \frac{1}{2}{{\boldsymbol{e}}^{d_2}}}\big), \end{equation} \tag{ 37 }$$

  where the 1D notation loosely indicates face-centered quantities corrected for the one-dimensional gradients during the normal predictor step.
• 4.  
  Likewise, use the CT scheme to correct magnetic field components affected by the same set of transverse flux derivatives. There are D − 1 magnetic field components of $U^n_{{\boldsymbol{i}}, \pm, d_1}$ on the face d₁ which are affected by the transverse flux along d₂. First, the component along d₁, which we indicate with $U_{B_{d_1}}$, is corrected as

  $$\begin{eqnarray} &&U_{B_{d_1},{\boldsymbol{i}}\pm,d_1,d_2}^{{\rm 2D}} = U_{B_{d_1},{\boldsymbol{i}}\pm,d_1}^{n} -\varepsilon _{d_1d_3d_2} \frac{\Delta t}{3\Delta x}\nonumber\\ &&\qquad\times\big(E^{{\rm 1D}}_{d_3,{\boldsymbol{i}}\pm \frac{1}{2}{{\boldsymbol{e}}^{d_1}}+\frac{1}{2}{{\boldsymbol{e}}^{d_2}}} - E^{{\rm 1D}}_{d_3,{\boldsymbol{i}}\pm \frac{1}{2}{{\boldsymbol{e}}^{d_1}}-\frac{1}{2}{{\boldsymbol{e}}^{d_2}}}\big). \end{eqnarray} \tag{ 38 }$$

  Second, the magnetic field component, $U_{B_{d_3}}$, along the remaining direction d₃ ≠ d₂ ≠ d₁. This component is cell-centered with respect to d₃, so its correction is given by the average of the relevant face-centered contributions at ${\boldsymbol{i}}+\frac{1}{2}{\boldsymbol{e}}^{d_3}$ and ${\boldsymbol{i}}-\frac{1}{2}{\boldsymbol{e}}^{d_3}$ (Gardiner & Stone 2008), namely,

  $$\begin{eqnarray} && U_{B_{d_3},{\boldsymbol{i}}\pm,d_1,d_2}^{{\rm 2D}} = U_{B_{d_3},{\boldsymbol{i}}\pm,d_1}^{n} -\varepsilon _{d_1d_3d_2} \frac{\Delta t}{6\Delta x}\nonumber\\ &&\qquad\times\big[\big(E^{{\rm 1D}}_{d_1,{\boldsymbol{i}}+\frac{1}{2}{\boldsymbol{e}}^{d_3}+\frac{1}{2}{{\boldsymbol{e}}^{d_2}}} - E^{{\rm 1D}}_{d_1,{\boldsymbol{i}}+\frac{1}{2}{\boldsymbol{e}}^{d_3}-\frac{1}{2}{{\boldsymbol{e}}^{d_2}}}\big)\nonumber \\ &&\qquad+\big(E^{{\rm 1D}}_{d_1,{\boldsymbol{i}}-\frac{1}{2}{\boldsymbol{e}}^{d_3}+\frac{1}{2}{{\boldsymbol{e}}^{d_2}}} - E^{{\rm 1D}}_{d_1,{\boldsymbol{i}}-\frac{1}{2}{\boldsymbol{e}}^{d_3}-\frac{1}{2}{{\boldsymbol{e}}^{d_2}}}\big) \big]. \end{eqnarray} \tag{ 39 }$$
• 5.  
  Next use a Riemann solver to obtain fluxes for each pair of states corrected for transverse fluxes. This provides D − 1 fluxes per cell face:

  $$\begin{eqnarray} F_{d_1,{\boldsymbol{i}}+ \frac{1}{2}{\boldsymbol{e}}^{d_1}, d_2}^{{\rm 2D}} &=& {\cal R} \big(U_{{\boldsymbol{i}}, +, d_1, d_2}^{{\rm 2D}}, U_{{\boldsymbol{i}}+ {\boldsymbol{e}}^{d_1}, -, d_1, d_2}^{{\rm 2D}}, d_1\big),\nonumber\\ d_1 &\ne & d_2,\quad 0 \le d_1,\quad d_2 < {\mathbf {D}}. \end{eqnarray} \tag{ 40 }$$
• 6.  
  Obtain new interpolated values of the electric field from the averages of the above computed fluxes

  $$\begin{equation} E_{d,{\boldsymbol{i}}+\frac{1}{2}{{\boldsymbol{e}}^{d_1}}+ \frac{1}{2}{{\boldsymbol{e}}^{d_2}}}^{{\rm 2D}} = {\cal I}_E\big(\tilde{F}_{d_1,* + \frac{1}{2}{{\boldsymbol{e}}^{d_1}}}^{{\rm 2D}}, \tilde{F}_{d_2,*+\frac{1}{2}{{\boldsymbol{e}}^{d_2}}}^{{\rm 2D}},W^{n}_{*}\big), \end{equation} \tag{ 41 }$$

  where

  $$\begin{eqnarray} \tilde{F}_{d_1,{\boldsymbol{i}}+ \frac{1}{2}{{\boldsymbol{e}}^{d_1}}}^{{\rm 2D}} &= &\frac{1}{2}\big(F_{d_1,{\boldsymbol{i}}+ \frac{1}{2}{{\boldsymbol{e}}^{d_1}},d_2}^{{\rm 2D}} +F_{d_1,{\boldsymbol{i}}+ \frac{1}{2}{{\boldsymbol{e}}^{d_1}},d_3}^{{\rm 2D}}\big), \end{eqnarray} \tag{ 42 }$$

  $$\begin{eqnarray} \tilde{F}_{d_2,{\boldsymbol{i}}+ \frac{1}{2}{{\boldsymbol{e}}^{d_2}}}^{{\rm 2D}} & = &\frac{1}{2}\big(F_{d_2,{\boldsymbol{i}}+ \frac{1}{2}{{\boldsymbol{e}}^{d_2}},d_1}^{{\rm 2D}} +F_{d_2,{\boldsymbol{i}}+ \frac{1}{2}{{\boldsymbol{e}}^{d_2}},d_3}^{{\rm 2D}}\big). \end{eqnarray} \tag{ 43 }$$
• 7.  
  Compute final corrections to the density, momentum, and energy components of $U_{{\boldsymbol{i}}, \pm, d}$ due to transverse fluxes using the above Riemann solutions, and obtain sought second-order accurate time-centered cell-interface values according to

  $$\begin{eqnarray} U^{n+\frac{1}{2}}_{{\boldsymbol{i}}, \pm, d} &=& \; U_{{\boldsymbol{i}}, \pm, d} \; - \; \frac{\Delta t}{2\Delta x} \big(F_{d_1,{\boldsymbol{i}}+ \frac{1}{2}{\boldsymbol{e}}^{d_1}, d_2}^{{\rm 2D}} - F_{d_1,{\boldsymbol{i}}- \frac{1}{2}{\boldsymbol{e}}^{d_1}, d_2}^{{\rm 2D}}\big)\nonumber \\ && - \frac{\Delta t}{2\Delta x} \big(F_{d_2,{\boldsymbol{i}}+ \frac{1}{2}{\boldsymbol{e}}^{d_2}, d_1}^{{\rm 2D}} - F_{d_2,{\boldsymbol{i}}- \frac{1}{2}{\boldsymbol{e}}^{d_2}, d_1}^{{\rm 2D}}\big),\nonumber\\ && d \ne d_1 \ne d_2,\quad 0 \le d, d_1,\quad d_2 < {\mathbf {D}}. \end{eqnarray} \tag{ 44 }$$
• 8.  
  Likewise, use the CT scheme to apply final corrections to magnetic field components. In analogy to Equation (38) in step 4 we write

  $$\begin{eqnarray} &&U_{B_d,{\boldsymbol{i}},\pm,d}^{n+\frac{1}{2}} = U_{B_d,{\boldsymbol{i}},\pm,d}^{n} -\varepsilon _{dd_1d_2} \frac{\Delta t}{2\Delta x}\nonumber\\ &&\qquad\times\big[\big(E_{d_1,{\boldsymbol{i}}\pm \frac{1}{2}{{\boldsymbol{e}}^{d}}+\frac{1}{2}{{\boldsymbol{e}}^{d_2}}}^{{\rm 2D}} - E_{d_1,{\boldsymbol{i}}\pm \frac{1}{2}{{\boldsymbol{e}}^{d}}-\frac{1}{2}{{\boldsymbol{e}}^{d_2}}}^{{\rm 2D}} \big)\nonumber \\ &&\qquad-\big(E_{d_2,{\boldsymbol{i}}\pm \frac{1}{2}{{\boldsymbol{e}}^{d}}+\frac{1}{2}{{\boldsymbol{e}}^{d_1}}}^{{\rm 2D}} - E_{d_2,{\boldsymbol{i}}\pm \frac{1}{2}{{\boldsymbol{e}}^{d}}-\frac{1}{2}{{\boldsymbol{e}}^{d_1}}}^{{\rm 2D}} \big)\big], \end{eqnarray} \tag{ 45 }$$

  and, in analogy to Equation (39), for the magnetic field components parallel to the cell face we write

  $$\begin{eqnarray} &&U_{B_{d_1},{\boldsymbol{i}}\pm,d}^{n+\frac{1}{2}} = U_{B_{d_1},{\boldsymbol{i}}\pm,d}^{n} -\varepsilon _{dd_1d_2} \frac{\Delta t}{2\Delta x}\nonumber\\ &&\qquad\times\big[\big(E_{d,{\boldsymbol{i}}+\frac{1}{2}{\boldsymbol{e}}^{d_1}+\frac{1}{2}{{\boldsymbol{e}}^{d_2}}}^{{\rm 2D}} - E_{d,{\boldsymbol{i}}+\frac{1}{2}{\boldsymbol{e}}^{d_1}-\frac{1}{2}{{\boldsymbol{e}}^{d_2}}}^{{\rm 2D}} \big)\nonumber \\ &&\qquad+\big(E_{d,{\boldsymbol{i}}-\frac{1}{2}{\boldsymbol{e}}^{d_1}+\frac{1}{2}{{\boldsymbol{e}}^{d_2}}}^{{\rm 2D}} - E_{d,{\boldsymbol{i}}-\frac{1}{2}{\boldsymbol{e}}^{d_1}-\frac{1}{2}{{\boldsymbol{e}}^{d_2}}}^{{\rm 2D}} \big)\nonumber \\ &&\qquad-\big(E_{d_2,{\boldsymbol{i}}+\frac{1}{2}{\boldsymbol{e}}^{d_1}+\frac{1}{2}{{\boldsymbol{e}}^{d}}}^{{\rm 2D}} - E_{d_2,{\boldsymbol{i}}+\frac{1}{2}{\boldsymbol{e}}^{d_1}-\frac{1}{2}{{\boldsymbol{e}}^{d}}}^{{\rm 2D}} \big)\nonumber \\ &&\qquad- \big(E_{d_2,{\boldsymbol{i}}-\frac{1}{2}{\boldsymbol{e}}^{d_1}+\frac{1}{2}{{\boldsymbol{e}}^{d}}}^{{\rm 2D}} - E_{d_2,{\boldsymbol{i}}-\frac{1}{2}{\boldsymbol{e}}^{d_1}-\frac{1}{2}{{\boldsymbol{e}}^{d}}}^{{\rm 2D}} \big) \big]. \end{eqnarray} \tag{ 46 }$$
• 9.  
  Compute the final second-order estimate of time-averaged fluxes

  $$\begin{equation} F^{n+\frac{1}{2}}_{d,{\boldsymbol{i}}+ \frac{1}{2}{\boldsymbol{e}}^d} = {\cal R} \big(U^{n+\frac{1}{2}}_{{\boldsymbol{i}}, +, d}, U^{n+\frac{1}{2}}_{{\boldsymbol{i}}+ {\boldsymbol{e}}^d, -, d}, d\big). \end{equation} \tag{ 47 }$$
• 10.  
  Compute the final estimate of the electric field using the above fluxes and an updated value for the cell-centered conservative variables using the averaged fluxes in Equation (42), namely,

  $$\begin{equation} E^{n+\frac{1}{2}}_{d,{\boldsymbol{i}}+\frac{1}{2}{{\boldsymbol{e}}^{d_1}}+ \frac{1}{2}{{\boldsymbol{e}}^{d_2}}} = {\cal I}_E\big(F^{n+\frac{1}{2}}_{d_1,* + \frac{1}{2}{{\boldsymbol{e}}^{d_1}}}, F^{n+\frac{1}{2}}_{d_2,*+\frac{1}{2}{{\boldsymbol{e}}^{d_2}}},\tilde{W}^{{\rm 2D}} _{*}\big), \end{equation} \tag{ 48 }$$

  where

  $$\begin{equation*} \tilde{W}_{{\boldsymbol{i}}}^{{\rm 2D}} = W_{{\boldsymbol{i}}}^n - \frac{1}{2}\frac{\Delta t}{\Delta x}\; \nabla _UW\cdot \sum \limits ^{{\mathbf {D}}- 1}_{d=0}\big(\tilde{F}_{d,{\boldsymbol{i}}+ \frac{1}{2}{\boldsymbol{e}}^d}^{{\rm 2D}} - \tilde{F}_{d,{\boldsymbol{i}}- \frac{1}{2}{\boldsymbol{e}}^d}^{{\rm 2D}}\big), \end{equation*}$$

  and the operator ∇_UW symbolizes the transformation from conservative to primitive variables.
• 11.  
  Finally update cell-centered conservative variables using Equation (13) with the fluxes in Equation (47) and the face-centered magnetic field variables using Equation (17) with the electric field in Equation (48) and synchronize the magnetic variables using Equation (19).

Note that Gardiner & Stone (2008) use the 2D fluxes, F^2D, and electric fields, E^2D, to construct the time-centered states at cell faces. This simplifies considerably the complexity of their algorithm with respect to the original CTU scheme, requiring only 6 Riemann solves instead of 12. The associated CFL stability condition, however, requires a time-step size half the usual value, so that for pure fluid application the total computational cost remains rather unaffected.

For the purpose of this paper, we have implemented the HLLD solver recently developed by Miyoshi & Kusano (2005). This is an extended version of the original HLL solver for the MHD, which includes the entropy, Alfvén, and fast magnetosonic waves. The solver appears quite accurate and robust and relatively inexpensive. Unlike traditional solvers (e.g., exact, linear, Roe's type, etc.), HLL-type solvers directly compute the hyperbolic fluxes, not the Riemann state, at cell interfaces. Furthermore, the flux solution is built upon the full spatially reconstructed state to the left and right sides of the cell interface. These states must also include the characteristic components that are not transmitted through the cell interface. For this reason, in order for the fluxes returned by HLL-type solvers to be second-order accurate, the projector P is modified in such a way that the summation is carried over all waves, irrespective of their sign. The solver is extensively documented in the original paper and its description will not be repeated here.

In this section, we describe the scheme used to interpolate the face-centered electric fields returned by the Riemann solver onto cell edges. In Section 2.2.4, this was indicated with the notation, ${\cal I}_E$. It is important for the stability of the overall algorithm to choose this interpolation scheme in such a way that there is consistency between the time update of the face-centered magnetic field, the cell-centered variables, and the synchronization step of the magnetic variables, Equation (19). For this purpose, we have adopted the upwind scheme described in Gardiner & Stone (2005, 2008). This will be described next, for the case of 3D. In this case each cell edge is shared by four adjacent faces from which the electric field can be interpolated. Hence, these four interpolated values will be arithmetically averaged. The scheme for ${\cal I}_E()$ uses three arguments,

$$\begin{equation*} E_{d,{\boldsymbol{i}}+\frac{1}{2}{{\boldsymbol{e}}^{d_1}}+\frac{1}{2}{{\boldsymbol{e}}^{d_2}}} ={\cal I}_E\big(F_{d_1,*+\frac{1}{2}{{\boldsymbol{e}}^{d_1}}},F_{d_2,*+\frac{1}{2}{{\boldsymbol{e}}^{d_2}}},W_*\big), \end{equation*}$$

where the "*" symbol indicates that the interpolation requires values of the arguments at different cells and faces. The face-centered fluxes define the face-centered electric fields according to Equation (15), and the cell-centered primitive state variable, $W_{\boldsymbol{i}}$, is used to define a cell-centered electric field according to the usual formula

$$\begin{equation*} E_{d,{\boldsymbol{i}}}=-\varepsilon _{dd_1d_2}\big(u_{d_1,{\boldsymbol{i}}}B_{d_2,{\boldsymbol{i}}}-u_{d_2,{\boldsymbol{i}}}B_{d_1,{\boldsymbol{i}}}\big), \end{equation*}$$

with the velocity and magnetic field variables given by the components of $W_{\boldsymbol{i}}$. The cell-centered electric field is used together with face-centered electric fields (15) to define the following transverse quasi-cell-centered gradients

$$\begin{equation} \left(\frac{\delta E_d}{\delta x_{d_1}}\right)_{{\boldsymbol{i}}+\frac{1}{4}{{\boldsymbol{e}}^{d_1}}} = 2\frac{E_{d,{\boldsymbol{i}}+\frac{1}{2}{{\boldsymbol{e}}^{d_1}}}-E_{d,{\boldsymbol{i}}}}{\Delta x}. \end{equation} \tag{ 49 }$$

In turn, these gradients are used to define quasi-face-centered gradients of the electric field by the upwind scheme (Gardiner & Stone 2005)

$$\begin{eqnarray*} &&\left(\frac{\delta E_d}{\delta x_{d_2}}\right)_{{\boldsymbol{i}}+\frac{1}{2}{{\boldsymbol{e}}^{d_1}}+\frac{1}{4}{{\boldsymbol{e}}^{d_2}}} \nonumber\\ &&=\left\lbrace \begin{array}{@{}l@{\quad }l@{}}\left(\frac{\delta E_d}{\delta x_{d_2}}\right)_{{\boldsymbol{i}}+\frac{1}{4}{{\boldsymbol{e}}^{d_2}}} & \mbox{if $u_{{\boldsymbol{i}}+\frac{1}{2}{{\boldsymbol{e}}^{d_2}}}> 0$} \\ \left(\frac{\delta E_d}{\delta x_{d_2}}\right)_{{\boldsymbol{i}}+\frac{3}{4}{{\boldsymbol{e}}^{d_2}}} & \mbox{if $u_{{\boldsymbol{i}}+\frac{1}{2}{{\boldsymbol{e}}^{d_2}}}< 0$} \\ \frac{1}{2}\left[\left(\frac{\delta E_d}{\delta x_{d_2}}\right)_{{\boldsymbol{i}}+\frac{1}{4}{{\boldsymbol{e}}^{d_2}}}+\left(\frac{\delta E_d}{\delta x_{d_2}}\right)_{{\boldsymbol{i}}+\frac{3}{4}{{\boldsymbol{e}}^{d_2}}}\right] & {\rm otherwise}, \end{array}\right.\qquad \end{eqnarray*}$$

where the sign of the velocity at the cell interface is given by the sign of the mass flux returned by the Riemann solver. So, the upwinding character of the interpolation scheme is based on the contact mode alone, which may raise stability concerns. However, no signs of unstable behavior appeared during the testing of the code. A more accurate edge interpolation scheme can be obtained at the higher cost of a 2D Riemann solver (Londrillo & Del Zanna 2000; Fromang et al. 2006). With the above definitions, the interpolated edge-centered electric field is defined as

$$\begin{eqnarray} &&E_{d,{\boldsymbol{i}}+\frac{1}{2}{{\boldsymbol{e}}^{d_1}}+\frac{1}{2}{{\boldsymbol{e}}^{d_2}}} = \frac{1}{4} \big(E_{d,{\boldsymbol{i}}+\frac{1}{2}{{\boldsymbol{e}}^{d_1}}} +E_{d,{\boldsymbol{i}}+{{\boldsymbol{e}}^{d_2}}+\frac{1}{2}{{\boldsymbol{e}}^{d_1}}} +E_{d,{\boldsymbol{i}}+\frac{1}{2}{{\boldsymbol{e}}^{d_2}}}\nonumber\\ &&\qquad+E_{d,{\boldsymbol{i}}+{{\boldsymbol{e}}^{d_1}}+\frac{1}{2}{{\boldsymbol{e}}^{d_2}}}\big)\nonumber \\ &&\qquad+\frac{\Delta x}{8} \left[ \left(\frac{\delta E_d}{\delta x_{d_2}}\right)_{{\boldsymbol{i}}+\frac{1}{2}{{\boldsymbol{e}}^{d_1}}+\frac{1}{4}{{\boldsymbol{e}}^{d_2}}} - \left(\frac{\delta E_d}{\delta x_{d_2}}\right)_{{\boldsymbol{i}}+\frac{1}{2}{{\boldsymbol{e}}^{d_1}}+\frac{3}{4}{{\boldsymbol{e}}^{d_2}}}\right]\nonumber \\ &&\qquad+\frac{\Delta x}{8}\left[ \left(\frac{\delta E_d}{\delta x_{d_1}}\right)_{{\boldsymbol{i}}+\frac{1}{2}{{\boldsymbol{e}}^{d_2}}+\frac{1}{4}{{\boldsymbol{e}}^{d_1}}} - \left(\frac{\delta E_d}{\delta x_{d_1}}\right)_{{\boldsymbol{i}}+\frac{1}{2}{{\boldsymbol{e}}^{d_2}}+\frac{3}{4}{{\boldsymbol{e}}^{d_1}}} \right].\qquad \end{eqnarray} \tag{ 50 }$$

The test based on the propagation of the Alfvén wave is characterized by the following unperturbed state and superposed perturbation (Crockett et al. 2005):

$$\begin{equation} W= {\left(\begin{array}{@{}c@{}}\rho _0 \\ u_x \\ u_y \\ u_z \\ P_0 \\ B_0/\sqrt{2} \\ B_0/\sqrt{2} \\ 0 \\ \end{array}\right)}, \quad \delta W = {\left(\begin{array}{@{}c@{}}0 \\ 0 \\ 0 \\ -c_A\\ 0\\ 0\\ 0\\ B_0\\ \end{array}\right)} \delta \sin ({\bf k}\, {\bm \cdot}\, {\bf x}), \end{equation} \tag{ 54 }$$

where ρ₀ = P₀ = B₀ = 1 and u_x = u_y = u_z = 0, $c_A=B_0/\sqrt{\rho _0}=1$, is the Alfvén speed, and k and x are the wave vector and position vector, respectively. The convergence rates for the perturbed quantities are summarized in Table 2. The velocity and magnetic field converge with second-order accuracy. The error on the unperturbed variables (not reported) is much smaller and converges at the same rate as the perturbed variables until dominated by the machine roundoff error. As already pointed out, tests with different orientation of $\bf k$ show the same convergence rates.

For fast and slow magnetosonic waves, the unperturbed state and the superposed perturbations read (Crockett et al. 2005)

$$\begin{eqnarray} W &=& {\left(\begin{array}{@{}c@{}}\rho _0 \\ u_x \\ u_y \\ u_z \\ P_0 \\ B_0 \hat{b}_x\\ B_0 \hat{b}_y\\ 0 \\ \end{array}\right)},\nonumber\\ \delta W &=& {\left(\begin{array}{@{}c@{}}\rho _0 \\ (\sqrt{2}c^2_w\hat{b}_y -c_g^2 \hat{k}_y)/c_w\\ -(\sqrt{2}c^2_w\hat{b}_x-c_g^2\hat{k}_x)/c_w\\ 0 \\ \rho _0c_g^2\\ -\sqrt{2}B_0(c_w^2-c_g^2)\hat{k}_y/c_A^2\\ \sqrt{2}B_0(c_w^2-c_g^2)\hat{k}_x/c_A^2\\ 0\\ \end{array}\right)} \delta \sin ({\bf k}\, {\bm \cdot}\, {\bf x}), \end{eqnarray} \tag{ 55 }$$

where $\hat{b}$ is the unit vector along the magnetic field vector and is at π/4 radians with respect to the unit vector $\hat{k}\equiv {\bf k}/k$, $c_g=\sqrt{\gamma P_0/\rho _0}$ is the gas sound speed, c_w is the speed of fast or slow MHD waves, and all other symbols take the meaning and values as in the previous section.

The convergence rates are reported in Tables 3 and 4 for fast and slow waves, respectively. As with Alfvén waves, the errors in the perturbed variables converge with second-order accuracy, while the error in the unperturbed variables are much smaller and converge at the same rate until they are affected by machine precision.

We begin with the Riemann problem presented in Brio & Wu (1988), listed at the top of Table 5. The x-component of the magnetic field is B_x = 0.75 and the adiabatic index γ = 2 as in the original paper. Following common practice, the problem is solved on a grid with 800 points along the x-axis and the solution is computed until time t = 0.1. The results are shown in Figure 2. All the solution features are well reproduced, including, from left to right, a fast rarefaction followed by a slow compound wave, moving to the left, and a contact discontinuity, slow shock and fast rarefaction moving to the right (Brio & Wu 1988). As usual with shock-capturing methods, shocks are resolved with a couple of zones throughout the duration of the calculation, while the contact discontinuities, captured here within a few zones, tend to spread out over time. There are some oscillations in the x-component of the velocity field. These are not present if we adopt a PLM reconstruction scheme, but worsen if we switch from a characteristic to a primitive limiting scheme.

Zoom In Zoom Out Reset image size

The second Riemann problem listed in Table 5 is taken from Dai & Woodward (1998). In this case, the x-component of the magnetic field is $B_x=4/\sqrt{4\pi }$ and the adiabatic index γ = 5/3. The problem is solved on a grid with 512 points along the x-axis and the solution is computed until time t = 0.2. The results are shown in Figure 3. The initial conditions for this problem expose the full eigenstructure of the MHD system, as they produce three pairs of MHD waves traveling in opposite directions with respect to the initial discontinuity, in addition to the contact wave. The waves include fast and slow shocks responsible, among others, for the jumps in the pressure and velocity fields, the contact wave which appears in the density field alone, and the rotational discontinuity which affects the magnetic field components alone. As in the previous case, while all discontinuities are well reproduced, shocks are the sharpest features captured with about two zones.

Zoom In Zoom Out Reset image size

The last 1D Riemann problem listed in Table 5 is similar to the ones presented in Einfeldt et al. (1991). Its purpose is to test the method/code robustness in the case of flows in which the energy is dominated by the kinetic component and unphysical states with negative density or internal energy can arise. This problem is solved using a grid with 256 grid points. The solution at time t = 0.1 is shown in Figure 4. The results in Figure 4 show a stable solution which correctly reproduces the two rarefaction waves propagating away from the grid midpoint. A very similar test has also been performed by Stone et al. (2008), with similar results, and Miyoshi & Kusano (2005). The latter authors actually use a faster expansion velocity, namely $u_{x,\frac{L}{R}}=\mp 3$, in their initial conditions and show that when used with a (first-order) Godunov's method the numerical solution remains stable. While we are able to reproduce their result we find that at high resolution some spurious oscillations can be generated when a higher-order Godunov's method is used. This may suggest the need for artificial viscosity in the case of highly supersonic expansions with higher-order Godunov's methods.

Zoom In Zoom Out Reset image size

We have repeated the problem in Section 3.2.2 with the initial discontinuity inclined with respect to the grid such that its normal has components $\vec{n}=(2,1)/\sqrt{5}$. This problem tests the ability of the code to reproduce 1D solutions when they are not aligned with the grid. Besides the numerical tests, there are a few other complications related with the numerical set-up of this problem. First, the boundary conditions need revision. To avoid the implementation of "shift-periodic" boundary conditions (Tóth 2000), we use a domain with size $(2L,L)\,2/\sqrt{5}$. This allows us to accommodate two adjacent identical problems along the direction, $\vec{n}=(2,1)/\sqrt{5}$, and apply periodic boundary conditions to the domain.⁴ In addition, since the magnetic field rotates across the discontinuity, the initialization is non-trivial and, unless a special precaution is taken, e.g., by deriving the magnetic field from a vector potential, the initial magnetic field will not be divergence-free (this is not an issue in the special case in which $\vec{n}$ is along the diagonal and there is symmetry among the coordinate axis). In our case no such precaution is taken and we just remap the initial conditions onto the rotated grid. As a result there is a jump in the normal component of the magnetic field across the discontinuity. This causes some minor artifacts with respect to the 1D solution. This is acceptable since we are interested in making sure that the structure of the 1D solution is reproduced with fidelity and the waves propagate at the correct speed.

In order to solve the problem, we have covered the domain with a grid of 1144 × 572 cells. This corresponds to 512 cells along the direction perpendicular to $\vec{n}$, which is equivalent to the resolution used in Section 3.2.2. The results are shown in Figure 5 where we plot the values of the solution along the first row of the computational domain, starting from $\frac{1}{2}L+\delta _L$ to $\frac{3}{2}L+\delta _L$. The starting point is shifted by δ_L ∼ 0.2 to the left, to hide the perturbation of W_left due to the interaction with the state to its left (which is W_right). The vectorial components (u, v, w) are the equivalent of (x, y, z) in the rotated system.

Zoom In Zoom Out Reset image size

Figure 5 shows that the 1D solution is correctly recovered when the plane of the discontinuity is inclined with respect to the grid. We notice some oscillation at the fast magnetosonic shocks, which to some extent are also present in Figure 7 of Gardiner & Stone (2008). These are probably associated with the oscillations in the normal component of the magnetic field, also reported in the last panel of Figure 5. These can perhaps be suppressed with more aggressive limiters (Londrillo & del Zanna 2004). There is also a spurious feature, ahead of the left-moving rotational discontinuity, which is probably due to the non-solenoidal character of the initial magnetic field. However, none of these features affect either the jump conditions or the wave speeds, as attested to by the very good correspondence between the solution in Figure 5 and the 1D counterpart in Figure 3.

The first test is the compressible Orszag–Tang vortex problem (Orszag & Tang 1979). We solve this problem on a computational domain of size L = 1 with periodic boundary conditions and a rectangular grid of 200 × 200 cells. We use an adiabatic index γ = 5/3. We present results obtained using the PPM reconstruction scheme and primitive limiting, but very similar results are produced using characteristic limiting. The initial conditions in terms of the primitive variables are as follows:

$$\begin{eqnarray} W &=&[\rho _0,-u_0\sin (2\pi y),u_0 \sin (2\pi x),0,P_0,\nonumber\\ &&-B_0\sin (2\pi y), B_0 \sin (4\pi x),0]^T, \end{eqnarray} \tag{ 57 }$$

where ρ₀ = 25/36π, P₀ = 5/12π, and u₀ and B₀ are defined in terms of the sonic Mach number, ${\cal M}=1$, and plasma beta, β = 10/3, respectively. Although the initial conditions are smooth, eventually the flow develops a complex structure with sharp features and discontinuities. In Figure 6 we show the contours of the numerical solution for density, gas pressure, kinetic energy, and magnetic pressure at time t = 0.5. 1D cuts along the line y = 0.4277 for the thermal pressure (top) and magnetic pressure (bottom) are also shown in Figure 7. The code maintains the symmetry of the solution with respect to the central point. In addition, Figure 7 shows that discontinuous features are captured within a few zones. Finally, the solution can be compared with similar plots at the same solution time produced by other authors (e.g., Ryu et al. 1998; Tóth 2000; Stone et al. 2008), from which it appears that the code produces the correct flow structures.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Another 2D problem commonly used as a test for multidimensional MHD codes is the rotor problem described in Balsara & Spicer (1999). It consists of a rotating disk of dense material, threaded by a magnetic field initially directed along the x-axis, and embedded in a tenuous ambient medium at rest. As the rotor spins, it winds up the magnetic field lines, generating Alfvén waves which propagate into the ambient medium. The problem is solved on a computational domain of size L = 1 with periodic boundary conditions and covered with a rectangular grid of 400 × 400 cells. The adiabatic index is γ = 1.4. We use the PPM reconstruction scheme and primitive limiting, although the use of characteristic limiting produces very similar results. The setup of the initial conditions corresponds to the first rotor problem discussed in Tóth (2000), i.e.,

$$\begin{equation} W [x,t=0] = \left\lbrace \begin{array}{@{}lll}W_{{\rm disk}} & \mbox{if} & r< r_{{\rm disk}},\\\ms W_{{\rm amb}} & \mbox{if} & r> r_{{\rm amb}}, \end{array} \right. \end{equation} \tag{ 58 }$$

where $W_{{\rm disk}}\,{=}\,[\rho _{{\rm disk}},-\frac{u_0}{r_0}(y-0.5),\frac{u_0}{r_0}(x-0.5),0, P_0,B_0,0,0]^T$, W_amb = (ρ₀, 0, 0, 0, P₀, B₀, 0, 0)^T, $r\equiv \sqrt{x^2+y^2}$, r_disk defines the disk radius, and r_amb delimits the ambient medium. In the transition region between the ambient medium and the rotor the density and velocity fields are interpolated according to ρ = (ρ_disk − ρ₀)f(r) + ρ₀, u_x = −f(r)u₀(y − 0.5)/r, and u_y = f(r)u₀(x − 0.5)/r, with f(r) = (r_amb − r)/(r_amb − r_disk), whereas density and pressure remain uniform. The results for this test are presented in Figure 8, which effectively reproduces Figure 18 of Tóth (2000). The numerical solution appears very well behaved and the code seems to pass this test as well.

Zoom In Zoom Out Reset image size

In this section, we test the ability of the code to follow the advection of a loop of weak magnetic field frozen in a background flow. This problem is non-trivial for conservative schemes and Gardiner & Stone (2005, 2008) emphasize that spurious results will be produced as a result of improper account of the multidimensional MHD source terms entering the predictor step, as discussed in Section 2.2.3. The initial conditions are detailed in Gardiner & Stone (2008). The loop is basically a tube of magnetic flux frozen in a medium with unit density and pressure, ρ = p_gas = 1 and uniform advection velocity, u_loop, to be specified below.

We carry out the test both in a 2D or a 3D geometry. For the 2D case, we align the loop axis with the x₂ coordinate axis of the computational domain. The vector potential from which the magnetic field is initialized is then given by, A = (0, 0, A₂), with

$$\begin{equation} A_2= \left\lbrace \begin{array}{@{}lll}B_0(R-r) & \mbox{if} & r\le R\\\ms 0 & \mbox{if} & r>R, \end{array} \right. \end{equation} \tag{ 59 }$$

where B₀ = 10⁻³, $r=\sqrt{\vphantom{A^A}\smash{\hbox{${x_0^2+x_1^2}$}}}$, and R = 0.3 is the radius of the tube. The computational domain itself consists of a rectangular box of dimensions (2N, N) with periodic boundaries, and the loop advection velocity is u_loop = (2, 1).

For the 3D configuration, the axis of the tube is tilted around the x₁-axis by $\theta =\arctan 2$ radians (clockwise) and the tube is advected with velocity u_loop = (1, 1, 2). The domain is still periodic but has dimensions (N, N, 2N). The vector potential is still defined as in Equation (59) but with respect to the new coordinates:

$$\begin{eqnarray} \nonumber x^\prime _0 &=& x_0\cos \theta +x_2\sin \theta, \\ x^\prime _1 &=& x_1, \\ x^\prime _2 &=& {-}x_0\sin \theta +x_2\cos \theta. \nonumber \end{eqnarray} \tag{ 60 }$$

The results of this test are first illustrated in Figures 9 and 10 which show the magnetic pressure at t = 2 for the 2D and 3D cases, respectively. In the former case N = 64, while in the latter case N = 128. The results are comparable to other MHD implementations, in particular Gardiner & Stone (2005, 2008). As pointed out by those authors, the magnetic pressure suffers dissipation mostly close to the loop center (where the curl of $\vec{B}$ is singular) and boundary. However, both in the 2D and 3D cases, the loop retains to a good extent its initial symmetry and energy. This latter point is further illustrated in Figure 11, which shows the evolution of the normalized magnetic pressure up to t = 2 for three different resolutions. One important aspect of the 3D version of this test is to check the ability of the scheme to keep the magnetic field component along the loop axis, B₃, close to zero. This is important here because we have not strictly followed the recommendation of Gardiner & Stone (2008) when implementing the multidimensional MHD source terms in the predictor step. The evolution of B₃ is reported in Figure 12, again for three different values of the resolution, up to t = 2. Due to the simpler form of the employed multidimensional MHD source terms, 〈|B₃|〉 is larger than found in Gardiner & Stone (2008) at the same time (t = 1), but only slightly so and B₃ remains negligible compared to the total magnetic field.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our first test involving AMR will be concerned with the advection of a loop of weak magnetic field, the same problem that was already addressed in Section 3.3.3. Here we focus on the 3D version of the problem, although analogous results are obtained in 2D. Thus, we use exactly the same initial conditions as for the 3D problem in Section 3.3.3. The domain is again periodic with dimensions (N, N, 2N), N = 64. We allow for one level of refinement with refinement ratio n_ref = 2. We refine regions where the magnetic field energy density is above a threshold B_threshold = 10⁻⁴ in code units. Effectively the magnetic loop is always refined (the same result is obtained if one refines regions where the magnetic field has a normalized gradient above a minimal threshold). Thus, we expect to obtain the same results as for the case in which N = 128. In Figures 11 and 12, we plot the time evolution of the normalized magnetic pressure and of the B₃ magnetic field component using open circles. In both cases, the open circles sit right on top of the solid line which indeed corresponds to the uniform box with AMR equivalent resolution, N = 128.

As a second example of an AMR–MHD application we compute the interaction of a cloud with a strong shock wave. This problem is fully 3D and it also employs the full CTU scheme. It mostly tests the code robustness in conditions of high Mach number shocks, high magnetic to thermal pressure ratio, strong shear flows, all of which are recurrent in astrophysical applications.

This process of shock–cloud interaction is common in the interstellar medium where shocks produced by supernova explosions interact with the surrounding multiphase medium (Klein et al. 1994; Mac Low et al. 1994). Related processes, characterized by similar hydrodynamic structures, are the supersonic motion of an over-dense cloud through a thin magnetized medium (Schiano et al. 1995; Jones et al. 1996; Vietri et al. 1997; Miniati et al. 1999b; Gregori et al. 1999, 2000), or supersonic clouds collisions (Lattanzio et al. 1985; Klein et al. 1995; Miniati et al. 1997, 1999a).

The initial conditions for our problem are as follows: the background gas has unit density and thermal pressure, and is at rest; a cloud with the same pressure but 10 times higher density, moves through the thin gas with a velocity v_c = −3.47871373 along the x direction. A plane-parallel shock with Mach number ${\cal M}=10$ propagates along the same axis but in the opposite direction to the cloud. The initial magnetic field is uniform, of unit strength and aligned with the x-axis. The computational box has dimensions [0, 1] × [0, 0.5] × [0, 0.5]. The boundary conditions correspond to supersonic inflow and outflow for the lower and upper boundaries of the x-axis, and are periodic otherwise. The calculation is carried out in 3D on a base grid of 64 × 32 × 32 cells. Two additional levels of refinement, with refinement ratio 4, are generated dynamically in regions where the normalized undivided density gradient |Δρ|/ρ > 0.1, and/or in the presence of shocks according to the criteria |ΔP|/P > 0.1 and ∇ · v < 0.

Figure 14 shows from top to bottom the solution for the density, gas pressure, magnetic field magnitude, and plasma beta parameter (β = P_gas/P_B), at time t = 0.021251, corresponding to 160 cycles on the finest level. The main features, discussed at length in the above papers, are correctly reproduced. The plane shock front moving from the left crushes the cloud. As the cloud moves to the left, it creates a bow shock in front of it, where the pressure has its highest value. The cloud motion also generates a low-pressure region at its rear, where the magnetic pressure dominates the gas pressure and the beta plasma is much less then unity. In addition, as the fluid flows past the cloud, the magnetic field lines entrained in the cloud body fold on themselves creating a current sheet (Jones et al. 1996).

Zoom In Zoom Out Reset image size

Note that the maximum value of the normalized divergence of the magnetic field |Δx∇ · B|/|B|, is a few × 10⁻¹³. This is completely negligible with respect to the solution value and demonstrates that our implementation of the above operators for the coarse–fine magnetic field interpolation and refluxing operations is correct.