AZEuS: AN ADAPTIVE ZONE EULERIAN SCHEME FOR COMPUTATIONAL MAGNETOHYDRODYNAMICS

AZEuS solves the following equations of ideal MHD (with the artificial viscosity and gravity terms included):

$$\begin{eqnarray} \frac{\partial \rho }{\partial t} + \nabla \cdot \left(\rho \vec{v}\right) & = 0; \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} \frac{\partial \vec{s}}{\partial t} + \nabla \cdot \left(\vec{s}\,\vec{v}\right) & =\,\, -\nabla p + (\nabla \times \vec{B})\times \vec{B} -\nabla \cdot \mathsf {Q} - \rho \nabla \phi ;\quad\quad \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} \frac{\partial \vec{B}}{\partial t} + \nabla \times \vec{E} & = 0; \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} \frac{\partial e}{\partial t} + \nabla \cdot (e\vec{v}) & =\,\, -p\nabla \cdot \vec{v} - \mathsf {Q} : \nabla \vec{v}; \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} \frac{\partial e_{T}}{\partial t} + \nabla \cdot (e_{H}\vec{v} + \vec{E}\times \vec{B} + \mathsf {Q}\cdot \vec{v}) & = 0, \end{eqnarray} \tag{ 5 }$$

where ρ is the mass density, $\vec{v}$ is the velocity, $\vec{s} = \rho \vec{v}$ is the momentum density, p is the thermal pressure, $\vec{B}$ is the magnetic induction,⁴ $\mathsf {Q}$ is the von Neumann–Richtmyer artificial viscous stress tensor (von Neumann & Richtmyer 1950; Clarke 2010), ϕ is the gravitational potential and satisfies the Poisson equation (∇²ϕ = 4πGρ), $\vec{E} = - \vec{v}\times \vec{B}$ is the induced electric field, e is the internal energy density, e_H = e + (1/2)ρv² + ρϕ is the hydrodynamical energy density, and e_T = e_H + (1/2)B² is the total energy density. This set of equations (in which Equation (5) can be derived from Equations (1)–(4)) is closed by the ideal gas law, p = (γ − 1)e, where γ is the ratio of specific heats. Figure 1 shows the locations of most of these variables on a fully staggered grid. Other physics terms often found in ZEUS codes such as a second fluid, physical viscosity, radiation, etc., have yet to be implemented.

AZEuS inherits the operator-split methodology of ZEUS-3D, wherein the terms on the right-hand side of Equations (1)–(5) are treated in a source step while those on the left-hand side are accounted for in a separate transport and inductive step. As such, the algorithm is not strictly conservative. However, based as it is on the version of ZEUS-3D described by Clarke (2010), AZEuS can solve either the internal energy equation or the total energy equation, where the latter choice does ensure conservation of total energy to machine round-off error, but at the cost of non-positive-definite thermal pressure. Should positive definite pressures be paramount, the internal energy equation offers a viable option with, in most cases, total energy conserved to within 1% or less (see Clarke 2010 for further discussion).

To accommodate the interpolation schemes, two boundary zones must be specified at the edges of all grids. On a staggered grid, all zone-centered quantities have just the two boundary values, while face-centered quantities have two boundary values plus a value that lies on the face separating the "active zones" from the "boundary zones," henceforth referred to as the "skin" of the grid (Figure 2). As we shall see, skin values for the magnetic field are treated just like active zones, while the momenta on the skin are treated somewhere in between active and boundary values, the difference attributed to the conservation properties of these two quantities.

Zoom In Zoom Out Reset image size

We adopt the following conventions and notation throughout this paper:

• 1.  
  Quantities in coarse and fine zones are denoted with upper and lower cases, respectively: e.g., Q(I, J, K), q(i, j, k). Fluxes for a quantity Q (q) are denoted F_{m, Q} (f_{m, q}), where m = 1, 2, or 3 indicates the component direction.
• 2.  
  If a "coarse" grid or zone is considered to be at level l, its daughter "fine" grid or zone is at level l + 1. The base and coarsest grid, which covers the entire domain, is at level l = 1. The refinement ratio, ν, between level l and l + 1 must be a power of 2 and the same in all directions.
• 3.  
  Grid volumes, areas, lengths, and time steps are ΔV, ΔA_m, Δx_m, and Δt for a coarse grid, and δV, δA_m, δx_m, and δt for a fine grid.
• 4.  
  Indices (i, j, k) correspond to the fine zone at the (left, bottom, back) of a coarse zone with indices (I, J, K). The zone center of a particular fine zone within a coarse zone is designated (i + α, j + β, k + η), where α, β, η = 0, 1, ..., ν − 1. For the 1-face-center of a fine zone, α = 0, 1, ..., ν while β, η = 0, 1, ..., ν − 1, etc.
• 5.  
  Similarly, grid positions for the coarse grid use uppercase indices (e.g., x₁(I)), while grid positions for the fine grid use lowercase indices (e.g., x₁(i)).
• 6.  
  The current time step for the coarse grid is indicated by the uppercase superscript N, while the current fine time step is indicated by the lowercase superscript n. Typically, n = νN. To designate a fine time step within a coarse time step, we use n + τ where 0 ⩽ τ ⩽ ν − 1⁵.
• 7.  
  The region of influence (ROI) of a variable is defined as the area or volume over which that quantity is conserved. For zone-centered scalars on the coarse grid, this is the volume ΔV(I, J, K) (Figure 3(a)). For face-centered but volume-conserved quantities (e.g., S₁), the ROI is the staggered volume (e.g., (1/2)(ΔV(I, J, K) + ΔV(I − 1, J, K))) (Figure 3(b)). Finally, for face-centered but area-conserved quantities (e.g., B₁), the ROI is the area of the face at which the vector component is centered (e.g., ΔA₁(I, J, K)) (Figure 3(c)).

Zoom In Zoom Out Reset image size

Finally, while AZEuS is written in the covariant fashion of ZEUS-2D (Stone & Norman 1992a), our discussion is given in terms of Cartesian-like components with uniform zone sizes within each grid for simplicity. As such, ΔV/δV = ν³, ΔA_m/δA_m = ν², Δx_m/δx_m = ν, and Δt/δt = ν. Some of the modifications necessary for curvilinear coordinates are given in Appendix A.

At the end of a coarse time step, fine and coarse grids are synchronized by overwriting the coarse grid with "better" values from overlying fine grids. Because of the different ROIs on a staggered mesh, the specifics of the overwriting procedure depend on which variable is being overwritten. For zone-centered, volume-conserved quantities (e.g., ρ, e, e_T), the procedure is the same as in BC89:

$$\begin{equation} Q(I,J,K) = \frac{1}{\nu ^3} \sum _{\alpha, \beta, \eta = 0}^{\nu -1}\, q(i+\alpha, j+\beta, k+\eta), \end{equation} \tag{ 6 }$$

where the sum is a triple sum. By inspection, one can confirm that Equation (6) conserves Q (q) locally to within machine round-off error.

For face-centered, volume-conserved quantities such as the momentum density whose ROIs are completely covered by the ROIs of overlying fine zones, we have

$$\begin{equation} S_1(I,J,K) = \frac{1}{\nu ^3} \sum _{\alpha =-\nu /2}^{\nu /2}\, \sum _{\beta, \eta = 0}^{\nu -1}\, {\cal G}(\alpha)\, s_1(i+\alpha, j+\beta, k+\eta), \end{equation} \tag{ 7 }$$

$$\begin{equation} \hbox{where } {\cal G}(\alpha) = \left\lbrace\begin{array}{@{}ll@{}}1/2 & \hbox{if } \alpha = \pm \, \nu /2 \\ 1 & \hbox{otherwise}, \end{array}\right. \end{equation} \tag{ 8 }$$

for the 1-component of the momentum. The factor ${\cal G}(\alpha)$ takes into account that only half of the ROIs of the fine momenta at α = −ν/2, ν/2 cover the ROI of coarse momentum (Figure 3(b)).

Coarse momenta with ROIs partially covered by a fine grid are cospatial with the skin of the overlying fine grid. As skin values of momenta are considered to be boundary values (since one of the fluxes is completely determined from within the boundary), they are not taken to be more reliable than the underlying coarse values (whose fluxes are determined exclusively by interior zones), and thus the coarse values are not overwritten by the fine grid values. Instead, coarse momenta cospatial with a fine grid skin are considered to be adjacent to the fine grid and, as such, are subject to the "flux-correction" step described in the next subsection.

For magnetic field, the conserved quantity is the magnetic flux ($\int \vec{B}\cdot d\vec{A}$). Thus, coarse values of B₁ are overwritten using

$$\begin{equation} B_1(I,J,K) = \frac{1}{\nu ^2} \sum _{\beta, \eta = 0}^{\nu -1}\, b_1(i, j+\beta, k+\eta), \end{equation} \tag{ 9 }$$

where the sum is over all fine ROIs (areas of 1-faces) that cover the coarse ROI. One can easily show that overwritten values of $\vec{B}$ will still satisfy the solenoidal condition—even when combined with values of $\vec{B}$ that are not overwritten—so long as the overlying values of $\vec{b}$ are divergence free and adjacent values of $\vec{B}$ are properly "refluxed" (Section 3.2). In addition, since $\vec{B}$ is an area-conserved quantity, there is never partial coverage of ROIs as there is with the momenta, and thus no analog of ${\cal G}(\alpha)$ (Equation (8)) is necessary for the magnetic field.

A straightforward permutation of indices gives the analogous expressions for the other components of momentum and magnetic field.

A coarse zone adjacent to but not covered by a fine grid shares a face with the fine grid. In order that local mass, momentum, and magnetic flux remain conserved to within machine round-off error, the coarse and fine zones must agree on the fluxes passing across their common face. This is accomplished by keeping track of all coarse fluxes passing across the skin of a contained fine grid, and then subtracting from these the fine fluxes computed during the MHD updates of the fine grids. These "flux corrections" are then subtracted from the coarse zones adjacent to the fine grid during the so-called refluxing step, effectively replacing the coarse fluxes with the fine fluxes along their common face.

For zone-centered, volume-conservative quantities this procedure follows BC89. Thus, for transport in the 1-direction, we have for the flux-corrected quantity, $\widetilde{Q}$,

$$\begin{eqnarray} \widetilde{Q}^{N+1}(I,J,K) &=& Q^{N+1}(I,J,K) - \frac{1}{\Delta V(I,J,K)}\nonumber \\ && \times \bigg [ F_{1,Q}^{N+{\frac{1}{2}}}(I,J,K) - \sum _{\beta,\eta,\tau =0}^{\nu -1} f_{1,q}^{n+\tau +{\frac{1}{2}}}\nonumber\\ && \times (i, j+\beta, k+\eta )\bigg ], \end{eqnarray} \tag{ 10 }$$

where the quantities in square brackets are the flux corrections. For the purpose of illustration, the coarse zone (I, J, K) is taken to be immediately to the right (increasing I) of a fine grid. F^{N + 1/2}_{1, Q} is the time-centered 1-flux of Q (with units $Q {\cal V} \Delta A \Delta t$, where ${\cal V}$ is the coarse velocity)⁶ passing across the 1-face cospatial with the skin of the fine grid, while f^{n + τ + 1/2}_{1, q} are the corresponding fine fluxes (with units qvδAδt, where v is the fine velocity; Figure 4(a)). Note the sum over τ reflecting the fact that there are several fine time steps (typically ν of them) within a single coarse time step.

Zoom In Zoom Out Reset image size

When flux correcting the face-centered, volume-conserved momenta, we must depart from BC89 as the staggered mesh gives rise to four distinct cases that must be dealt with individually. The first case is when both the flux and momentum components are parallel to the fine grid skin normal. Here, the ROI at the boundary is halfway inside the fine grid. For example, consider the 1-flux of S₁ in the ROI straddling the right boundary of a fine grid as shown in Figure 4(b) where the 1-flux correction takes the form

$$\begin{eqnarray} \widetilde{S}_1^{N+1}(I,J,K) &=& S_1^{N+1}(I,J,K) - \frac{1}{\Delta V(I,J,K)}\nonumber\\ && \times \Bigg [ F_{1,S_1}^{N+{\frac{1}{2}}}(I-1,J,K) - \frac{1}{2} \sum _{\beta,\eta,\tau =0}^{\nu -1}\nonumber\\ && \times \bigg (f_{1,s_1}^{n+\tau +{\frac{1}{2}}}\bigg (i-\frac{\nu }{2}-1, j+\beta, k+\eta \bigg)\nonumber\\ && +f_{1,s_1}^{n+\tau +\mbox{$\frac{1}{2}$}}\bigg (i-\frac{\nu }{2}, j+\beta, k+\eta \bigg) \bigg) \Bigg ].\quad\quad \end{eqnarray} \tag{ 11 }$$

Note that in the 1-direction, the coarse and fine momenta pass through the same face, but the fluxes do not. Thus, an average of the fine fluxes at i − ν/2 − 1 and i − ν/2 is needed to properly center the fine fluxes and to ensure local conservation of momentum.

The second case is when the flux is parallel and the momentum component perpendicular to the fine grid skin normal. Here, the adjacent ROI lies entirely outside the fine grid just as a zone-centered quantity. For example, consider the 1-flux of S₂ in the ROI adjacent the right boundary of a fine grid as shown in Figure 4(c), where the flux correction takes the form

$$\begin{eqnarray} \widetilde{S}_2^{N+1}(I,J,K) &=& S_2^{N+1}(I,J,K) - \frac{1}{\Delta V(I,J,K)}\nonumber\\ && \times \Bigg [ F_{1,S_2}^{N+{\frac{1}{2}}}(I,J,K)\, -\,\sum _{\beta =-\nu /2}^{\nu /2}\, \sum _{\eta,\tau =0}^{\nu -1}\, {\cal G}(\beta)\,\nonumber\\ && \times f_{1,s_2}^{n+\tau +{\frac{1}{2}}}(i, j+\beta, k+\eta) \Bigg ]. \end{eqnarray} \tag{ 12 }$$

The factor $\cal {G}(\beta)$ (Equation (8)) ensures that only half of the fine 1-fluxes of the fine ROIs at β = ±ν/2 are included, as a quick glance at Figure 4(c) will verify.

The third case is when the flux is perpendicular and the momentum component parallel to the fine grid skin normal. As with the first case, the ROI at the boundary is halfway inside the fine grid. For example, consider the 1-flux of S₂ in the ROI straddling the lower boundary of a fine grid as shown in Figure 4(d), where the flux correction takes the form

$$\begin{eqnarray} \widetilde{S}_2^{N+1}(I,J,K) &=& S_2^{N+1}(I,J,K) - \frac{1}{\Delta V(I,J,K)} \Bigg \lbrace \frac{\nu -1}{2\nu }\nonumber\\ && \times \Big [ F_{1,S_2}^{N+{\frac{1}{2}}}(I,J,K)-F_{1, S_2}^{N+{\frac{1}{2}}}(I+1,J,K)\Big ] \nonumber\\ && - \sum _{\beta =1}^{\nu /2}\, \sum _{\eta,\tau =0}^{\nu -1}\, {\cal G} (\beta) \Big [ f_{1,s_2}^{n+\tau +{\frac{1}{2}}}(i,j+\beta,k+\eta)\nonumber\\ && - f_{1,s_2}^{n+\tau +{\frac{1}{2}}}(i+\nu,j+\beta,k+\eta) \Big ] \Bigg \rbrace . \end{eqnarray} \tag{ 13 }$$

The fraction (ν − 1)/2ν is the area filling ratio between the fine and coarse fluxes. For example, for ν = 4, (ν − 1)/2ν = 3/8 (Figure 4(d)), and 3/8 of the coarse flux must be replaced with the overlapping fine fluxes. Note that $f_{1,s_2} (i,j,k)$ does not contribute to the fine fluxes that replace part of the coarse flux since it is determined, in part, by values from the boundary region of the fine grid, and thus not taken as reliable enough to replace part of the coarse flux determined exclusively from active zones of the coarse grid.

Equation (13) assumes that both the left and right faces of the coarse ROI need to be flux corrected. Should the fine grid in Figure 4(d) not cover coarse zones (I, J, K) and (I + 1, J, K), flux corrections need only be applied to the left side of the coarse ROI (I, J, K). In this case, the terms $F_{1, S_2}^{N+1/2}(I+1,J,K)$ and $f_{1,s_2}^{n+\tau +1/2}(i+\nu,j+\beta,k+\eta)$ are omitted. Similarly, if the fine grid does not cover coarse zones (I − 1, J, K) and (I, J, K), the terms $F_{1, S_2}^{N+1/2}(I,J,K)$ and $f_{1,s_2}^{n+\tau +1/2}(i,j+\beta,k+\eta)$ are omitted.

The fourth and final case is when both the flux and momentum components are perpendicular to the fine grid skin normal. Here, no flux corrections are required, as can be readily seen by inspection of Figure 4(c).

For the magnetic field components, the operator in the evolution equation (Equation (3)) is the curl, and not the divergence of the hydrodynamical variables. Thus, $\vec{B}$ is a surface-conserved quantity (rather than volume-conserved) and, as such, we make adjustments to the so-called EMFs (defined below) instead of the fluxes; otherwise, the procedure is the same.

Consider the coarse ROI of B₃ immediately to the right of a fine grid (Figure 4(e)). Following Balsara (2001), we have

$$\begin{eqnarray} \widetilde{B}_3^{N+1}(I,J,K) &=& B_3^{N+1}(I,J,K) + \frac{1}{\Delta A_3(I,J,K)}\nonumber\\ && \times \Bigg [ {\cal E}_2^{N+{\frac{1}{2}}}(I,J,K) - \sum _{\beta,\tau =0}^{\nu -1} \varepsilon _2^{n+\tau +{\frac{1}{2}}}\nonumber\\ && \times (i, j+\beta, k) \Bigg ], \end{eqnarray} \tag{ 14 }$$

where ${\cal E}_2^{N+1/2}(I,J,K) = E_2 \Delta x_2 \Delta t$ is the time-centered coarse 2-EMF⁷ ("electro-motive force") located along the 2-edge, and where ε^{n + τ + 1/2}₂(i, j + β, k) is the time-centered fine 2-EMF, located along the same 2-edge as the coarse 2-EMF. Here, the quantities in square brackets are the "EMF corrections," analogous to the flux corrections for the hydrodynamical variables. In ZEUS-3D, the EMFs can be evaluated using various algorithms. In our case, we use the Consistent Method of Characteristics (CMoC) described by Clarke (1996).

Because the magnetic field is an area-conserved quantity, there is no situation where the coarse ROI of a magnetic field component is partially covered by fine ROIs as can happen for the volume-conserved, face-centered momentum. Thus, Equation (14) and its permutations are sufficient to cover all EMF corrections for all field components in all directions. Note, for example, that there are no corrections to be made in the 1-direction for the 1-field, again because of the surface conservative nature of the magnetic field.

Further, induction of magnetic field components penetrating the skin of a fine grid is affected by EMFs computed from quantities taken entirely from within the active portion of the grid. This is in contrast to a momentum component penetrating the skin, half of whose fluxes are computed from boundary values. Thus, we take the skin values of the EMFs and the magnetic components they induce to be just as reliable as those evaluated from within the grid, and it is appropriate to restrict the coarse magnetic field values on the skin of a fine grid with the overlying fine values (e.g., Equation (9)).

Finally, flux-correction equations for coarse zones adjacent to other sides of a fine grid can be obtained by a suitable permutation of indices and subscripts.

In an effort to minimize the errors caused by waves traveling across boundaries between fine and coarse grids, we have introduced higher-order interpolation schemes (e.g., piecewise parabolic interpolation, PPI; Colella & Woodward 1984) to the prolongation step. This improves the results for adaptive refinement and has proven important for static refinement where strong waves are required to cross grid boundaries. By design, these interpolation schemes honor conservation laws and monotonicity.

In all cases and for all variables, we begin by estimating the interface values Q_{L, R}(I) from cubic fits to the coarse grid data (Section 1 of Colella & Woodward 1984) without the monotonization or steepening steps. Now, to fit a parabolic interpolation function, q*₁(x), across the coarse ROI for zone I, we need three constraints. Two come from requiring that q*₁(x) passes through Q_{L, R}(I), while the third comes from requiring that the zone-average value Q(I, J, K) times the zone volume equals the volume integral of q*₁(x). For the original PPI scheme, this final constraint is written as

$$\begin{equation*} Q(I,J,K) \Delta x_1(I) =\, \int _{\Delta x} q_1^*(x)dx. \end{equation*}$$

For our purposes, we need to interpolate fine zone averages from the coarse ones and, as such, the conserved quantity is the Riemann sum and not the integral. Thus, in the 1-direction, our final constraint is

$$\begin{equation} Q(I,J,K) \Delta x_1(I) = \sum ^{\nu -1}_{\alpha =0}\, q_1^*(i+\alpha, j, k) \,\delta x_1(i) \end{equation} \tag{ 15 }$$

$$\begin{equation} \Rightarrow Q(I,J,K) = \frac{1}{\nu } \sum _{\alpha =0}^{\nu -1}\, q_1^*(i+\alpha, j, k) \end{equation} \tag{ 16 }$$

for constant Δx, δx. With this, our parabolic interpolation function is (cf. Equation (1.4) in Colella & Woodward 1984)

$$\begin{equation} q^*_1(i+\alpha) = Q_L(I) + \zeta ( Q_R(I) - Q_L(I) + {\cal H}_1(1 - \zeta) ), \end{equation} \tag{ 17 }$$

where

$$\begin{eqnarray} \zeta &=& \frac{x_1(i+\alpha) - x_1(I)}{\Delta x_1(I)}; \nonumber\\ {\cal H}_1 &=& \frac{1}{f_1(\nu) - f_2(\nu)}\nonumber\\ && \times (Q(I,J,K) - Q_L(I) - f_1(\nu) (Q_R(I) - Q_L(I) ) ); \nonumber\\ f_1(\nu) &=& \frac{1}{2\nu ^2} \sum _{\xi =1}^{\nu }\, (2\xi - 1) = \frac{1}{2};\nonumber\\ f_2(\nu) &=& \frac{1}{4\nu ^3} \sum _{\xi =1}^{\nu }\, (2\xi - 1)^2 = \frac{4\nu ^2-1}{12\nu ^2}. \end{eqnarray} \tag{ 18 }$$

With q*₁(i + α) determined, we set the differences between the fine and coarse zones:

$$\begin{equation} \vspace*{-1pt}\delta q_1 (i+\alpha) = q^*_1(i+\alpha) - Q(I,J,K).\vspace*{-1pt} \end{equation} \tag{ 19 }$$

For zone-centered quantities, we determine the fine-zone differences in the manner of Equation (19) in each of the three directions and set the interpolated fine zone averages to be

$$\begin{eqnarray} q(i+\alpha, j+\beta, k+\eta) &=& Q(I,J,K) + \delta q_1(i+\alpha)\nonumber\\ && + \delta q_2(j+\beta) + \delta q_3 (k+\eta).\quad\quad \end{eqnarray} \tag{ 20 }$$

Given Equation (16) and analogous expressions in the 2- and 3-directions, it is easy to show that this prescription is conservative, that is,

$$\begin{equation} \vspace*{-1pt}Q(I,J,K) = \frac{1}{\nu ^3}\,\sum _{\alpha,\beta,\eta =0}^{\nu -1}q(i+\alpha,j+\beta,k+\eta), \vspace*{-1pt} \end{equation} \tag{ 21 }$$

for uniform grids (ΔV = ν³δV).

For face-centered quantities, we determine the fine-zone differences in the manner of Equation (19) in each of the two orthogonal directions (e.g., in the 2- and 3-directions for S₁) and set the interpolated fine zone averages along each coarse face to be

$$\begin{eqnarray} &&\vspace*{-1pt}s_1(i, j+\beta, k+\eta) = S_1(I,J,K) + \delta s_{1,2}(j+\beta) + \delta s_{1,3} (k+\eta).\nonumber\\ &&\vspace*{-1pt} \end{eqnarray} \tag{ 22 }$$

Exactly the same procedure is used to interpolate B₁ across each coarse 1-face. As with the zone-centered quantities, it is easy to show that this procedure conserves momentum (magnetic) flux on the face:

$$\begin{equation} \vspace*{-1pt} S_1(I,J,K) = \frac{1}{\nu ^2}\,\sum _{\beta,\eta =0}^{\nu -1}s_1(i,j+\beta,k+\eta).\vspace*{-1pt} \end{equation} \tag{ 23 }$$

The next step is to interpolate the face-centered quantities into the interior of the zone, and it is here where the procedure for momentum and magnetic components diverge.

For the volume-conserved momenta, we perform a linear interpolation between the fine momenta on opposing coarse faces. For example, between s₁(i, j + β, k + η) and s₁(i + ν, j + β, k + η) we set

$$\begin{eqnarray} \vspace*{-1pt} s_1(i+\alpha,j+\beta,k+\eta) &=& (1 - \zeta)\, s_1(i,j+\beta,k+\eta)\nonumber\\ && + \zeta \, s_1(i+\nu,j+\beta,k+\eta),\quad \vspace*{-1pt} \end{eqnarray} \tag{ 24 }$$

where α = 1, ..., ν − 1 and ζ = αδx/Δx. This prescription guarantees that over a zone-centered volume, the prolongation of the 1-momentum is conservative. The fact that conservation is over the zone-centered volume and not specifically the 1-momentum ROI is required for Equations (23) and (24) to be consistent with the restriction procedure of Equation (11).

For the area-conserved magnetic field, a simple linear interpolation between opposing coarse faces does not preserve the solenoidal condition. Thus, we turn to a generalized, directionally unsplit version of the algorithm described by Li & Li (2004) in which $\nabla \cdot \vec{b}=0$ so long as $\nabla \cdot \vec{B} = 0$ holds on the underlying coarse grid.

In this approach, we take the coarse and recently interpolated fine values of magnetic field (Equation (22)) and apply the solenoidal condition to determine "intermediate" coarse field values (e.g., B*₁; Figure 5), which are cospatial with fine zone faces. For example, the intermediate values for B₁ are given by

$$\begin{eqnarray} B_1^*(i+\alpha +1) &=& B_1^*(i+\alpha) - \delta x_1(i+\alpha)\nonumber\\ && \times \Bigg (\frac{b_{2,1}^*(i+\alpha,j+\nu,k) - b_{2,1}^*(i+\alpha,j,k)}{\Delta x_2(J)}\nonumber\\ && + \, \frac{b_{3,1}^*(i+\alpha,j,k+\nu) - b_{3,1}^*(i+\alpha,j,k)}{\Delta x_3(K)}\Bigg), \nonumber\\ \end{eqnarray} \tag{ 25 }$$

where

$$\begin{eqnarray} b_{2,1}^*(i+\alpha,j,k) = & \,\, \frac{1}{\nu } \sum _{\eta = 0}^{\nu -1}\, b_2(i+\alpha,j,k+\eta), \nonumber \\ b_{3,1}^*(i+\alpha,j,k) = & \,\, \frac{1}{\nu } \sum _{\beta = 0}^{\nu -1}\, b_3(i+\alpha,j+\beta,k); \nonumber\\ B_1^*(i+0) = & \, B_1(I,J,K), \end{eqnarray} \tag{ 26 }$$

and 0 ⩽ α ⩽ ν − 2. Given B*₁(i + α + 1), and the differences between fine and coarse values at the coarse faces (δb_{1, l}, l = 2, 3; Equation (19)), we then calculate the fine magnetic field at i + α + 1:

$$\begin{eqnarray} b_1(i+\alpha +1,j+\beta,k+\eta) &=& B_1^*(i+\alpha +1)\,\nonumber\\ && + \, \delta b_{1,2}(i+\alpha +1,j+\beta,k)\nonumber \\ && +\, \delta b_{1,3}(i+\alpha +1,j,k+\eta), \quad\quad \end{eqnarray} \tag{ 27 }$$

where

$$\begin{eqnarray} \delta b_{1,2}(i+\alpha,j+\beta,k) &=& (1 - \zeta)\, \delta b_{1,2}(i,j+\beta,k)\nonumber\\ && + \zeta \, \delta b_{1,2}(i+\nu,j+\beta,k),\nonumber \\ \delta b_{1,3}(i+\alpha,j,k+\eta) &=& (1 - \zeta)\, \delta b_{1,3}(i,j,k+\eta)\nonumber\\ && + \zeta \, \delta b_{1,3}(i+\nu,j,k+\eta), \end{eqnarray} \tag{ 28 }$$

and ζ = αδx/Δx.

Zoom In Zoom Out Reset image size

Proceeding incrementally from i + 1 to i + ν − 1 provides all of the values of b₁(i, j, k) between the coarse faces (I, J, K) and (I + 1, J, K). Values for the other magnetic field components are given by a straightforward permutation of indices. In aggregate, these yield a third-order interpolation of fine magnetic field components that preserve the divergence of the underlying coarse magnetic field to machine round-off error.

Finally, a note on the use of vector potentials. Igumenshchev & Narayan (2002) have demonstrated that the vector potential, $\vec{A}$, may be used in the CT algorithm (Evans & Hawley 1988) instead of the magnetic field with results identical to machine round-off error. Thus, one may be tempted to adopt this approach so that prolongation methods simpler than the Li & Li algorithm may be used on $\vec{A}$ as a way to guarantee that the fine grid satisfies the solenoidal condition. Reasons for not taking this approach are given in Appendix B.

Fine grids require prolonged boundary values at the beginning of each fine time step, and thus we also need to perform temporal interpolations on the coarse values. Following BC89, we perform a linear interpolation on the coarse hydrodynamic variables in time and then spatially interpolate them as described above to obtain the necessary boundary information for each fine time step.

For the magnetic field, the fine components penetrating the skins of the fine grid are retained; only those values completely within the boundary region are prolonged from the coarse grid. For this, we require coarse grid EMFs on the fine grid skin and within the boundary region that are time-centered on the current fine time step.

The coarse EMFs coincident with the fine grid skin are replaced with the spatial and temporal sum of the overlying fine EMFs, where the sum over time is necessary since the time step is embedded in our definition (footnote ⁷). Thus and for example, for ${\cal E}_2$ we set

$$\begin{equation} {\cal E}_2^{N+\psi }(I,J,K) = \sum _{\tau =0}^{\tau ^{\prime }}\,\sum _{\beta =0}^{\nu -1}\, \varepsilon _2^{n+\tau +\mbox{$\frac{1}{2}$}} (i,j+\beta,k), \end{equation} \tag{ 29 }$$

where ψ = (2τ' + 1)/2ν, and where 0 ⩽ τ' ⩽ ν − 2 designates the last completed fine time step. This yields values on the fine grid skin that are time-centered between the beginnings of the coarse time step and the current fine time step and take into account the superior quality of the fine grid data.

The coarse EMFs residing entirely within the fine grid boundary region are retained. However, these are time-centered for the coarse grid and, as such, have embedded in them a factor Δt. To render them coeval with fine time step τ' and thus the coarse skin EMFs computed in Equation (29), we must multiply them by (2τ' + 1)/ν = 2ψ to correct the embedded time step factor.

With coarse EMFs on the fine skin and inside the fine boundary region properly time-centered, the coarse magnetic field components are updated to the beginning of the current fine time step using Equation (3), and it is these values that are used in the prolongation methods described in Section 4.1 to create the fine grid magnetic field boundary values.

One might ask why a linear temporal interpolation of the divergence-free coarse magnetic field components is not enough to preserve the solenoidal condition within the boundary region. Indeed, such an approach does guarantee $\nabla \cdot \vec{B}=0$ within the region interpolated, but not in the layer between the interpolated region and its immediate non-interpolated neighbors. The method outlined above based on the coarse EMFs allows interpolations to be performed locally and only where they are needed, while still preserving the solenoidal condition globally to machine round-off error.

We have found it necessary to maintain a certain monotonicity in our prolongations. Failure to do so can lead to negative pressures (even when solving the internal energy equation) and violations of the CFL condition. For example, if the current time step is governed by the sound speed, and the prolongation process leads to an interpolated density less than the surrounding zones, the sound speed in a fine zone could be greater than that which was used in determining the CFL time step, possibly leading to numerical oscillations and loss of stability.

Sequentially stringing together two (or three) 1D PPIs as we do for prolongation of zone-centered quantities (e.g., Equation (20)) can lead to non-monotonic behavior even if each 1D interpolation is separately monotonic. If a PPI-determined value q(i + α, j + β, k + η) is found to lie outside the range set by neighboring coarse zones,

$$\begin{eqnarray*} &&\fl \left.\begin{array}{@{}ll@{}}Q_{\min } = {\rm min}\, (Q(I+\Gamma,J+\Lambda,K+\Upsilon) ); \\ [6pt] Q_{\max } = {\rm max}\, (Q(I+\Gamma,J+\Lambda,K+\Upsilon)),\end{array}\hspace{-6pt}\right\rbrace -1 \le \Gamma, \Lambda, \Upsilon \le +1, \end{eqnarray*}$$

then we "fall back" to piecewise linear interpolations (PLI; van Leer 1977, Appendix A). In rare cases where the PLI-determined value is also non-monotonic (possible only in 3D), we revert to piecewise constant interpolations (PCI, a.k.a. "direct injection" or "donor cell"):

$$\begin{equation*} q(i+\alpha,j+\beta,k+\eta) = Q(I,J,K). \end{equation*}$$

Where PPI yields non-monotonic results, we note that transmission of strong waves across changes in resolution (e.g., static grids) is significantly improved if one first tries PLI rather than falling back directly to PCI or limiting the interpolated value to lie between Q_min and Q_max.

Generally, a fine grid is completely embedded within a coarse grid. The single exception is when both grids share a physical boundary and two items of note must be borne in mind when adapting the physical boundary condition routines in ZEUS-3D to AZEuS.

First, since each grid has only two boundary zones, only part of the coarse boundary region is covered by fine boundary zones. Thus, coarse boundary zones cannot be included in the restriction step, a particular concern when setting magnetic boundary conditions. In AZEuS, we extend the EMF correction scheme of Section 3.2 by retaining three layers of transverse EMF corrections and two layers of longitudinal EMF corrections,⁸ including the skin layer plus two additional layers interior and adjacent to the skin. Immediately before the restriction step, physical boundary conditions are applied to the EMF corrections, which are then used to update the boundary values of the magnetic field components in the coarse grid according to Equation (14). Done in this fashion, there is no risk of introducing monopoles to the coarse boundary zones during the restriction step when two or more grids overlap a physical boundary.

Second, not handled properly, inflow boundary conditions can introduce unexpected violations of conservation laws that can cause unwanted discontinuities in the boundary. In particular, if a boundary variable is to be set according to an analytical function of the coordinates, that variable should be set to the zone average of that function, and not simply to the function value at the location of the variable. While this is good advice for a single-grid application, it is critical for AMR.

For example, suppose the density profile,

$$\begin{equation} \rho (r) = \frac{1}{r^{3/2}}, \end{equation} \tag{ 30 }$$

is to be maintained in cylindrical coordinates along the z = 0 boundary. Let ρ(J) be the density in the coarse zone of dimension (Δz, Δr) centered at r(J) = r, and let ρ(j) and ρ(j + 1) be the densities in the fine zones of dimension (δz, δr) centered at r(j) and r(j + 1) (Figure 6). For a refinement ratio of 2, Δr = 2δr, Δz = 2δz, r(j) = r − (1/2)δr, and r(j + 1) = r + (1/2)δr.

Zoom In Zoom Out Reset image size

If we naively set ρ(J) using Equation (30), the mass of zone J (with volume ΔV = rΔrΔz) is

$$\begin{equation*} \vspace*{-1pt}m(J) = \rho (J) \Delta V(J) = \frac{\Delta r \Delta z}{r^{1/2}}.\vspace*{-1pt} \end{equation*}$$

Similarly, for the fine zones,

$$\begin{equation*} \vspace*{-1pt}m(j) = \frac{\delta r \delta z}{\big (r-\mbox{$\frac{1}{2}$}\delta r\big)^{1/2}}; \qquad m(j+1) = \frac{\delta r \delta z}{\big (r+\mbox{$\frac{1}{2}$}\delta r\big)^{1/2}}.\vspace*{-1pt} \end{equation*}$$

Adding the masses in the four fine zones contained by the coarse zone, we find

$$\begin{equation*} \vspace*{-1pt}2\big [ m(j)+m(j+1)\big ] = m(J)\!\left[ 1 + \frac{3}{32} \!\left(\frac{\delta r}{r}\right)^{\!\!2} + \dots \right ] \,> m(J),\vspace*{-1pt} \end{equation*}$$

and the physical boundary conditions violate mass conservation.

Instead, ρ(J) should be set to the zone average, determined by

$$\begin{equation*} \vspace*{-1pt}\rho (J) = \frac{m(J)}{\Delta V(J)}\vspace*{-1pt} \end{equation*}$$

where the mass function, m(J), is given by

$$\begin{equation} \vspace*{-1pt}m(J) =\, \int _{\Delta V(J)}\!\!\rho (r) dV,\vspace*{-1pt} \end{equation} \tag{ 31 }$$

where ρ(r) is the given density profile (e.g., Equation (30)). Similar expressions apply for m(j) and m(j + 1). Done in this fashion, it is easy to show that

$$\begin{equation*} \vspace*{-1pt}\rho (J)\Delta V(J) = 2 [ \rho (j+1)\delta V(j+1) + \rho (j)\delta V(j) ],\vspace*{-1pt} \end{equation*}$$

and both grids agree on the mass contained within coarse boundary zone J to within machine round-off error. With a little bit of algebra, this can be confirmed analytically for the specific case of Equation (30). For more complicated profiles, a numerical integrator can be employed to perform the necessary integrations in Equation (31).

Unlike Godunov methods, which typically require a single application of boundary conditions at the end of each MHD cycle, the operator-split nature of ZEUS requires boundary conditions to be set several times. For physical boundaries, this is not an issue; physical boundary conditions may be set whenever needed based on data from a single grid. However, adjacent boundaries pose a unique challenge in AZEuS since the ZEUS module is aware only of the grid being updated, and other grids cannot be accessed in the middle of an MHD cycle to update these boundaries.

We have introduced into AZEuS the concept of "self-computing" boundary conditions for adjacent boundaries. In this approach, boundary zones are set at the beginning of each time step using the best information available (either from an adjacent grid of the same resolution or from the prolongation of an underlying coarse grid possibly interpolated in time), and then the full set of operator-split MHD equations are applied to both the boundary and active zones; no adjacent boundary zones are reset to any predetermined quantity inside a single MHD time step. Further, self-computed boundary zones are included in the calculation of the CFL time step, which we find critical (for an operator-split code such as AZEuS) in minimizing transmission errors when waves of any significant amplitude cross adjacent boundaries.

Of course, assumptions about missing data beyond the outermost edge of the grid must be made, and "pollution" from these missing data necessarily propagates inward. For example, consider the left 1-boundary where v₁(i = 1) and v₁(i = 2) represent the boundary values of v₁, v₁(i = 3) the skin value (treated as a boundary value), and v₁(i = 4) the first active value. The pressure gradient at v₁(1) is proportional to p(1) − p(0), yet p(0) is completely unknown (p(1) and p(2) are the only boundary values available). The best we can do for a missing datum such as this is to extrapolate assuming a zero gradient, in which case p(0) = p(1) and no pressure gradient is applied to v₁(1). In this manner, v₁(1) is polluted by the missing datum at the very beginning of the source step.

Additional steps inside a single MHD cycle include the application of artificial viscosity, adiabatic expansion/compression for the internal energy equation, the transport steps, and the induction step. All but the latter contribute to propagating pollution from missing data toward the active portion of the grid. Ideally, one would carry enough boundary zones to prevent pollution from reaching the active grid within a single MHD step, so that the AMR module (with knowledge of all grids) can reset all boundary values to new and unpolluted values before their effects ever reach the active zones. However, for van Leer (1977) interpolation in the transport step and within a single MHD cycle, pollution from missing data reaches the skin and first active zone center on the grid, as well as the first active face and second active zone center if the local velocity points away from the boundary (true regardless of which energy equation is used). Thus, to completely prevent pollution from reaching the active grid, we would have to double the number of boundary zones from two to four.

We have investigated this effect thoroughly and have found no test problem that demonstrates anything but the slightest quantitative effect from pollution by missing data at adjacent boundaries between fine and coarse grids. For adjacent boundaries between grids of like resolution, we have elected to force these grids to overlap (for reasons explained in Section 6) by an amount sufficient to eliminate the problem of missing data pollution altogether.

While most AMR applications rely exclusively on dynamic grids, certain applications, particularly those that exhibit some degree of self-similarity, can benefit enormously from the use of a base of nested, static grids (e.g., Ramsey & Clarke 2011). As such solutions evolve, waves of all types—including strong shocks—must pass sequentially from one grid to another, and it is here where our higher-order prolongation algorithms for adjacent boundaries are critical.

AZEuS has been tested with all 1D shock tube problems from Ryu & Jones (1995, hereafter RJ95) for both static and dynamic grids, and we show the results from two to highlight these features of the code. Both tests use the total energy equation with γ = 5/3, ${\cal C} = 0.75$ (Courant number), and artificial viscosity parameters qcon = 1.0 and qlin = 0.2.

For static grids, Figure 8 shows two solutions for ρ, e_T, v₂, and B₂ from problem (4a) of RJ95 over a domain of x₁ ∈ [ − 0.5, 2.5] (three times larger than RJ95) for a time t = 0.45 (three times longer). The initial discontinuity is placed at x₁ = 0.5, with the left and right states given in the figure caption. The left panels show the domain resolved with 1200 zones and no AMR, while the central panels show the AMR solution with a base resolution of 600 zones and fine static grids with a refinement ratio of 2 placed at x₁ ∈ [ − 0.2, 0.1] and x₁ ∈ [0.7, 1.245] (gray). These locations allow all the physical features, with the exception of the sluggish slow rarefaction, to suffer at least one change in resolution. The right panels show the percent differences between the two solutions.

Zoom In Zoom Out Reset image size

Discounting all zones trapped within a discontinuity (which, even without AMR, are already in error by as much as 100% since discontinuities are supposed to be infinitely sharp), the maximum error one can attribute to the use of static grids in any of the variables is less than 1%. As further evidence of the ability of AZEuS to pass waves of all types across adjacent boundaries, Ramsey & Clarke (2011) find almost no sign of reflection, refraction, or distortion of any type of wave across any of the static grid boundaries in their 2D axisymmetric simulations of protostellar jets in cylindrical coordinates.

For dynamic grids, Figure 9 illustrates the results of test problem (2a) of RJ95. We employ four levels of refinement (sequentially darker shades of gray) above the base grid with a resolution of 40 zones over the domain x₁ ∈ [0.0, 1.0]. The initial discontinuity is located at x₁ = 0.5, with the left and right states given in the figure caption. To track the developing features, refinements are based on three criteria (Khokhlov 1998): contact discontinuities (CD), shocks, and gradients in certain variables, each requiring a threshold value to be set above which a zone is flagged for refinement. Here, we set parameters tolcd = tolshk = 0.1 for CDs and shocks, respectively, and check v₂ for gradients above tolgrad = 0.1. For this problem, the gradient detector is needed to refine both the discontinuity of the initial conditions and rotational discontinuities (x₁ ≃ 0.53 and x₁ ≃ 0.71 in Figure 9), neither of which are detectable as a CD or shock.

Zoom In Zoom Out Reset image size

Additional parameters used to obtain this solution include kcheck = 5 (the number of cycles between successive grid modifications), geffcy = 0.95 (the minimum allowed grid efficiency, defined as the ratio of zones flagged for refinement to the number of zones actually present in a new grid), and ibuff = 2 (the number of buffer zones added around flagged zones). These parameters, in addition to the non-AMR parameters, are summarized in Table 1.

Table 1. Summary and Description of the Significant User-set Parameters in AZEuS

Parameter	Description
${\cal C}$	The Courant number
qcon	Quadratic artificial viscosity parametera
qlin	Linear artificial viscosity parametera
tolcd	Threshold value for the CD detector (range 0–1)
tolshk	Threshold value for the shock detector (range 0–1)
tolgrad	Threshold value(s) for the gradient detector(s) (range 0–1)
kcheck	Number of cycles between successive grid modifications
geffcy	The minimum allowed fractional grid efficiency for creation of new grids (range 0–1)
ibuff	The number of buffer zones added around zones flagged for refinement

Note.^a See Clarke (2010) for a more detailed description of the artificial viscosity parameters.

Download table as:  ASCIITypeset image

7.2.1. MHD Blast

The MHD blast problem of Londrillo & Del Zanna (2000) and Gardiner & Stone (2005) has proven to be a very valuable test of our AMR algorithms and for rooting out problems in the code. It is also a good test of directional biases and the ability of a code to handle the evolution of strong MHD waves. We initialize the problem in the same manner as Clarke (2010) with domain x ∈ [ − 0.5, 0.5], y ∈ [ − 0.5, 0.5], and $(\rho, \vec{v}, B_1, B_2, B_3) = (1, 0, 5\sqrt{2}, 5\sqrt{2}, 0)$ everywhere. Within radius r = 0.125 of the origin, we set the gas pressure to p = 100, and p = 1 elsewhere.

Figure 10 presents our results for the 2D blast problem at t = 0.02 for a base grid of 200² zones and one additional level of refinement with ratio ν = 2. For this test, we have purposely chosen a high grid efficiency of geffcy = 0.92 to strenuously test the ability of AZEuS to handle a large number of grids in a complicated pattern and refine on features that are not preferentially aligned with a coordinate axis. Typically, a lower value for the grid efficiency is used (e.g., geffcy ≃ 0.7) to try to balance minimizing the number of refined zones with the overhead associated with managing an increasing number of small grids.

Zoom In Zoom Out Reset image size

Evidently, AZEuS is able to follow the blast wave closely, regardless of its orientation. Table 2 compares the extrema of the plotted variables between the AMR calculation and uniform grid solutions with 200² and 400² zones. With the exception of the pressure, which exhibits differences between the uniform 400² and AMR solutions of ≲0.5%, the uniform grid results bracket the AMR solution.

Table 2. Extrema for Density, ρ, Gas Pressure, p, and Magnetic Pressure, p_B, in AMR and Uniform Grid Solutions of the 2D MHD Blast Problem at t = 0.02

	AMR		Uniform, 2002		Uniform, 4002
	Min	Max	Min	Max	Min	Max
ρ	0.199	3.36	0.200	3.22	0.189	3.40
p	0.712	32.3	0.771	32.0	0.714	32.1
pB	24.3	75.7	24.9	76.0	23.6	75.6

Download table as:  ASCIITypeset image

For this test, we set tolshk = tolcd = 0.2 and apply tolgrad = 0.2 to e_T. While the gradient detector is useful in refining the initial pressure jump, most (∼99.9%) of the zones flagged for refinement soon thereafter are detected by the CD and shock detectors. Additional parameter settings include kcheck = 10, ibuff = 3, γ = 5/3, ${\cal C} = 0.5$, qcon = 1.0, and qlin = 0.1. All boundaries are set to outflow conditions.

7.2.2. Orszag–Tang MHD Vortex

The 2D vortex problem of Orszag & Tang (1979) has become a standard test for astrophysical MHD codes, and as it has not previously been performed using our version of ZEUS-3D, we present both non-AMR and AMR results here. It is important to note that the same code was used to produce both sets of results: AZEuS is designed to be modular, and by deselecting the "AMR" option at the precompilation step, the code reverts to ZEUS-3D.

For this test, we follow Stone et al. (2008) by initializing a periodic, Cartesian box of size x, y ∈ [0, 1] with initially constant pressure and density, P = 5/12π and ρ = γP = 25/36π, for γ = 5/3. The velocity is initialized to (v_x, v_y, v_z) = (− sin  (2πy), sin  (2πx), 0), and the magnetic field is set through the vector potential A_z = (B₀/4π)cos  (4πx) + (B₀/2π)cos  (2πy), where $B_0 = 1/\sqrt{4\pi }$.

The results for uniform grids with 256² and 512² zones at t = 1/2, as well as the AMR results for a base grid of 128² zones and two levels of refinement (effective resolution of 512² zones), are presented in Figure 11. The bottom right panel shows the distribution of AMR grids at t = 1/2. For clarity, we have only plotted the grids at level 3 (i.e., two levels higher resolution than the base grid). Even then, the filling factor of level l = 2 grids at this time is ≳ 95%. Examining the first three panels of Figure 11 closely, features are noticeably sharper in the 512² solution relative to the 256² results, while the AMR and 512² solutions are indistinguishable.

Zoom In Zoom Out Reset image size

Figure 12 shows slices of the gas pressure as a function of x at t = 1/2 and y = 0.4277, which once again demonstrate that the AMR and 512² uniform grid solutions are virtually identical. Quantitatively, these solutions compare favorably with those from higher-order codes such as ATHENA (Stone et al. 2008), with the discontinuities in the AZEuS solutions being slightly broader.

Zoom In Zoom Out Reset image size

For this problem, we set ${\cal C} = 0.5$, qcon = 1.0, and qlin = 0.1 for both uniform and adaptive grids. For the AMR results, kcheck = 10, geffcy = 0.9, ibuff = 2, and tolshk = 0.2. Neither the CD nor gradient detector was engaged.

7.2.3. Magnetized Accretion Torus

This test is based on the simulations of Hawley (2000) and Mignone et al. (2007) for a magnetized, constant angular momentum torus in axisymmetric spherical (r, ϑ) coordinates and highlights the use of curvilinear coordinates in AZEuS.

The torus structure is described by the equilibrium condition:

$$\begin{equation} \frac{\gamma p}{\left(\gamma - 1\right) \rho } = C - \phi - \frac{1}{2} \frac{l^2_{\rm Kep}}{r^2 \sin ^2 \!\vartheta }, \end{equation} \tag{ 32 }$$

where C is a constant of integration, ϕ = −1/(r − 1) is the pseudo-Newtonian gravitational potential, and l_Kep is the Keplerian angular momentum at the pressure maximum. The pressure p is initially related to the density via the polytropic relation, p = κρ^γ, with γ = 5/3. By specifying the location of the pressure maximum (r_max = 4.7) and the inner edge of the torus (r_min = 3), we can determine the value of C and the (hydrodynamic) structure of the torus.

A poloidal magnetic field is initialized in the torus from the φ-component of the vector potential:

$$\begin{equation*} A_\varphi = \frac{B_0}{\rho _{\rm m}}\, {\rm min}\, (\rho (r,\vartheta) - \rho _{\rm c}, 0), \end{equation*}$$

where B²₀ = 2κρ^γ_m/β_m, β_m and ρ_m are, respectively, the plasma-beta and density at r_max, κ is determined by Equation (32) evaluated at r_max, and ρ_c = ρ_m/2 determines the surface of the last vector equipotential. For the simulations presented here, β_m = 350 and ρ_m = 10. The toroidal velocity (v_φ) in the torus is initialized to the local Keplerian speed, while the poloidal velocity is set to zero everywhere.

Outside the torus, the magnetic field is zero and we initialize a hydrostatic atmosphere with density and temperature contrasts of ρ_atm/ρ_m = 10⁻⁴ and T_atm/T_m = 100, respectively, where T_m is the temperature at r_max.

The domain of the grid is r ∈ [1.5, 20] and ϑ ∈ [0, π/2]. We impose reflecting boundary conditions suitable for a rotation axis at ϑ = 0, reflecting/conducting boundary conditions at the equatorial plane (ϑ = π/2), and outflow boundary conditions at r = 20. At r = r_in = 1.5, we impose "sink" boundary conditions in an attempt to absorb any material reaching the inner boundary. This involves maintaining the density and pressure at their initial values and setting v_r = v_ϑ = v_φ = 0 within the boundary. On the inner skin (r = r_in), v_r is set to the minimum of zero and a linear extrapolation from the grid. Finally, outflow boundary conditions are applied to the EMFs:

$$\begin{eqnarray} {\cal E}_1(r < r_{\rm in}) &=& {\cal E}_1(2r_{\rm in}-r); \nonumber \\ {\cal E}_2(r < r_{\rm in}) &=& 2{\cal E}_2(r_{\rm in}) - {\cal E}_2(2r_{\rm in}-r);\nonumber \\ {\cal E}_3(r < r_{\rm in}) &=& 2{\cal E}_3(r_{\rm in}) - {\cal E}_3(2r_{\rm in}-r). \end{eqnarray} \tag{ 33 }$$

Values of ${\cal E}_2(r_{\rm in})$ and ${\cal E}_3(r_{\rm in})$ are determined by the CMoC algorithm in AZEuS.

Both the uniform grid and AMR solutions are presented here. For the single-grid calculation, we use 592 uniform radial zones and 256 uniform meridional zones. For the AMR solution, we use a base grid of 296 × 128 zones and allow one level of refinement with a refinement ratio ν = 2, giving the same effective resolution as the uniform grid calculation.

For both calculations, ${\cal C}$ = 0.5, qcon = 1.0, and qlin = 0.1. For AMR, kcheck = 50, geffcy = 0.8, ibuff = 2, tolshk = tolcd = 0.2, and tolgrad = 0.2 is applied to B_ϑ. We use the gradient detector to refine on the initial magnetic field configuration (contained entirely within the torus) and the complex field structure that develops later on from the magnetorotational instability (MRI; Balbus & Hawley 1991), since neither is well tracked by the shock and CD detectors.

Following Hawley (2000), we enforce a density floor of 10⁻³ρ(r_in) to prevent the time step from becoming prohibitively small, and we use the internal energy equation to avoid negative pressures. Unlike Hawley (2000), we do not introduce any perturbations on the initial conditions, which slows but does not prevent the onset of MRI.

Figure 13 presents the results of our simulations at t = 300. Evidently and very much unlike the Orszag–Tang vortex, the AMR and uniform grid solutions are visually different, even though they have the same effective resolution. For a pseudo-turbulent, "irreversible" problem such as this, small differences between values in the coarse and fine grids grow exponentially during the simulation and lead to visually different solutions. This is quite unlike the behavior of a "reversible" problem such as the Orszag–Tang vortex, for which small fluctuations grow at worse linearly in time and never manifest as visual differences in the plots. On the other hand, integrated quantities tend to be in better agreement than their detailed distributions. For example, between the uniform grid and AMR solutions, the integrated mass and kinetic energy at t = 300 differ by only 0.9% and 6.5%, respectively. Furthermore, our results are qualitatively similar to Mignone et al. (2007); Figure 13 could easily fit in with their Figure 8 for different Riemann solvers and interpolation schemes.

Zoom In Zoom Out Reset image size

Figure 14 shows the normalized magnetic field divergence, $|\nabla \cdot \vec{B} |\, \Delta x / |\vec{B}|$, at the simulation end time, where the maximum value is 1.34 × 10⁻¹⁵, including both physical and adjacent boundary zones. Indeed, the minimum and maximum values for the divergence throughout the simulation lifetime reach −7.392 × 10⁻¹⁵ and 7.396 × 10⁻¹⁵, respectively, demonstrating that the algorithms employed for prolongation and restriction of the magnetic field in AZEuS are capable of maintaining the solenoidal condition to machine round-off, assuming that the initial conditions are also divergence free. It is also worthwhile to note the absence of grid boundary effects in Figure 14, which one might expect if monopoles existed within the boundaries of the fine grids.

Zoom In Zoom Out Reset image size

The last test is the self-gravitational hydrodynamical collapse calculation first presented by Truelove et al. (1997) and Truelove et al. (1998, hereafter T98). At the time of this writing, self-gravity had not been fully implemented in AZEuS, and so we used the successive over-relaxation (SOR) method (Press et al. 1992) as an easy-to-program, albeit slow, stop-gap measure. Since SOR is incompatible with periodic boundary conditions, we apply inflow boundary conditions instead. This results in some subtle differences from T98, which we discuss below.

A rotating, uniform cloud of mass M = 1 M_☉ and radius R = 5 × 10¹⁶ cm is initialized in the center of a 3D Cartesian box with side length 4R. We use a nearly isothermal equation of state (p = κρ^γ, γ = 1.001), initially uniform rotation Ω = 7.14 × 10⁻¹³ rad s⁻¹ with angular momentum axis in the positive x₃-direction, and energy ratios of

$$\begin{eqnarray} \alpha &=& \frac{5}{2} \left(\frac{3}{4\pi \rho _0 M^2}\right)^{1/3} \frac{c_{s}^2}{G} = 0.16 \quad \hbox{and} \nonumber\\ \beta _\Omega &=& \frac{1}{4\pi } \frac{\Omega ^2}{G\rho _0} = 0.26, \end{eqnarray} \tag{ 34 }$$

where ρ₀ and c_s are, respectively, the initial density and sound speed of the uniform cloud (T98). The remainder of the computational volume is initialized with uniform density and pressure given by ρ_atm = 100ρ_cloud and p_atm = 10p_cloud. Upon the initial uniform density distribution, we apply an azimuthal m = 2 perturbation of the form ρ_cloud(1 + Acos (2φ)), where A = 10% is the perturbation amplitude and φ is the azimuthal angle relative to the cloud center.

Zoom In Zoom Out Reset image size

Following T98, we set the coarsest grid to 32³ zones and designate it as R₈ (meaning that the cloud radius is resolved with eight zones). We immediately add one level of refinement by flagging all zones that have density ρ > ρ_cloud (for a refinement ratio of 4, this gives 32 zones per cloud radius, or R₃₂). Beyond this, we enforce the so-called Truelove criterion:

$$\begin{equation} J = \frac{\Delta x}{\lambda _J} = \Delta x \left(\frac{G\rho }{\pi c_{s}^2}\right)^{1/2} < J_{\rm max}, \end{equation} \tag{ 35 }$$

refining zones that have Jeans numbers J larger than critical value J_max = 0.25. Additional run-time parameters for this test include ${\cal C} = 0.33$, qcon = 1.0, qlin = 0.1, kcheck = 20, geffcy = 0.7, ibuff = 2, and refinement ratio ν = 4.

Figure 15 shows our results at t = 0.59808 = 1.2152t_ff. At this time, there are seven levels of refinement above the base grid (R₁₃₁₀₇₂). Our maximum density at this time is log ρ_max = −9.347399, with density measured in g cm⁻³, an increase over ρ₀ of more than eight orders of magnitude. While our simulation collapses more quickly than T98 and our Figure 15 does not correspond exactly with their Figures 12 and 13, we also reach a somewhat lower maximum density.

The differences between our results and T98 are most likely caused by the differing boundary conditions (ours inflow, theirs periodic). Indeed, T98 allude to this possibility, suggesting that non-periodic boundary conditions could slightly increase the rate of collapse, which we observe.

A.1.1. Conservative Overwrite

When applying the conservative overwriting procedure on a curvilinear grid, the equations of Section 3.1 must be modified to account for the non-constant volumes of zones. For example, Equation (6) for zone-centered quantities becomes

$$\begin{eqnarray} && Q(I,J,K)\, \Delta V_1(I)\,\Delta V_2(J)\,\Delta V_3(K)\,\nonumber\\ && \quad = \sum _{\alpha, \beta, \eta = 0}^{\nu -1} q(i+\alpha, j+\beta, k+\eta)\nonumber\\ && \quad \times \,\delta V_1(i+\alpha)\,\delta V_2(j+\beta)\,\delta V_3(k+\eta), \end{eqnarray} \tag{ A9 }$$

and Equation (7) for face-centered momenta generalizes to

$$\begin{eqnarray} && S_1(I,J,K)\,\Delta V_1(I)\,\Delta V_2(J)\,\Delta V_3(K)\,\nonumber\\ && \quad = \sum _{\alpha =-\nu /2}^{\nu /2} \, \sum _{\beta, \eta = 0}^{\nu -1} {\cal G}^{\prime }(\alpha)\, s_1(i+\alpha, j+\beta, k+\eta)\nonumber\\ && \quad \times \delta V_1(i+\alpha)\,\delta V_2(j+\beta)\,\delta V_3(k+\eta), \end{eqnarray} \tag{ A10 }$$

where

$$\begin{equation} \quad {\cal G}^{\prime }(\alpha) =\left\lbrace \begin{array}{@{}ll}w_{\rm 1,R}(i + \alpha) & \hbox{if } \alpha = -\nu /2; \\ w_{\rm 1,L}(i + \alpha) & \hbox{if } \alpha = +\nu /2; \\ 1 & \hbox{otherwise}, \end{array}\right. \end{equation} \tag{ A11 }$$

and where

$$\begin{eqnarray} w_{\rm 1,L}(i) &=& \frac{1}{\delta V_1(i)} \int _{x_1(i-1/2)}^{x_1(i)} \!dV_1;\nonumber\\ w_{\rm 1,R}(i) &=& \frac{1}{\delta V_1(i)} \int _{x_1(i)}^{x_1(i+1/2)} \!dV_1. \end{eqnarray} \tag{ A12 }$$

While it is still true that the fine zones at α = ±ν/2 are halfway outside the coarse ROI in position space (Figure 3(b)), this is not generally true in volume space. As an example, consider spherical coordinates where the 1-volume differential is dV₁(i) = g₂₁(i)g₃₁(i) dx₁(i) = x²₁(i) dx₁(i). Clearly, dV₁(i) in spherical coordinates is not a linear function of position, and it cannot be assumed that exactly 1/2 of the fine momentum volume is outside the coarse ROI. To correct for this, we calculate the ratio of the half-volume inside the coarse ROI to the actual volume element (Equation (A12)), which is then used as a weighting factor in ${\cal G}^{\prime }(\alpha)$. For reasonable grid parameters, these weighting factors are small corrections (<a few %) relative to the Cartesian factor of 1/2 but nonetheless are included for accuracy.

Suitable expressions for the 2-direction are obtained by a permutation of indices. No corrections are necessary for the 3-direction, since dV₃(k) = dx₃(k) for the curvilinear coordinates discussed here.

For the magnetic field, Equation (9) changes to account for non-uniform areas:

$$\begin{eqnarray} B_1(I,J,K)\,\Delta A_1(I,J,K) &=& \sum _{\beta, \eta = 0}^{\nu -1}\, b_1(i, j+\beta, k+\eta)\,\nonumber\\ && \times \delta A_1(i,j+\beta,k+\eta).\quad\quad \end{eqnarray} \tag{ A13 }$$

A.1.2. Flux Corrections

Fluxes and EMFs as stored by ZEUS-3D already include the appropriate area and length elements, so the required changes for flux corrections in curvilinear coordinates are not substantial. As with the conservative overwrite, the corrections (e.g., Equation (12)) need to be adjusted by replacing occurrences of ${\cal G}$ (Equation (8)) with ${\cal G}^{\prime }$ (Equation (A11)). Similar to Equation (A10), Equation (11) for the flux corrections of momentum components normal to the boundary must also be modified to account for non-constant volumes:

$$\begin{eqnarray} &&\widetilde{S}_1^{N+1}(I,J,K) = S_1^{N+1}(I,J,K) \,-\, \frac{1}{\Delta V(I,J,K)}\nonumber \\ && \quad \times \Bigg \lbrace F_{1,S_1}^{N+{\frac{1}{2}}}(I-1,J,K) - \sum _{\beta,\eta,\tau =0}^{\nu -1} \bigg [ w_{\rm 1,R} \bigg (i-\frac{\nu }{2}\bigg)\,f_{1,s_1}^{n+\tau +{\frac{1}{2}}}\nonumber\\ && \quad \times \bigg (i-\frac{\nu }{2}-1, j+\beta, k+\eta \bigg) + w_{\rm 1,L} \bigg (i-\frac{\nu }{2}\bigg)\,f_{1,s_1}^{n+\tau +{\frac{1}{2}}}\nonumber\\ && \quad \times \bigg (i-\frac{\nu }{2}, j+\beta, k+\eta \bigg) \bigg ] \Bigg \rbrace . \end{eqnarray} \tag{ A14 }$$

By the same argument, Equation (13) for coarse momenta with flux components parallel to a grid boundary also needs to be modified, as the original factor of $(\nu - 1) / 2\nu \equiv {\cal R}$ is not correct for curvilinear coordinates:

$$\begin{eqnarray} && \widetilde{S}_2^{N+1}(I,J,K) = S_2^{N+1}(I,J,K)- \frac{1}{\Delta V(I,J,K)}\nonumber\\ && \quad \times \Bigg \lbrace {\cal L}(J,\nu)\bigg [ F_{1,S_2}^{N+{\frac{1}{2}}}(I,J,K)-F_{1, S_2}^{N+{\frac{1}{2}}}(I+1,J,K)\bigg ] \nonumber\\ && \quad - \sum _{\beta =1}^{\nu /2}\, \sum _{\eta,\tau =0}^{\nu -1}\, {\cal G}^{\prime }(\beta)\bigg [ f_{1,s_2}^{n+\tau +{\frac{1}{2}}}(i,j+\beta,k+\eta)\nonumber\\ && \quad - f_{1,s_2}^{n+\tau +{\frac{1}{2}}}(i+\nu,j+\beta,k+\eta) \bigg] \Bigg \rbrace . \end{eqnarray} \tag{ A15 }$$

Here, ${\cal L}(J,\nu) = W_{2,L}(J,\nu)$ if S₂(I, J, K) is at the upper edge of the fine grid (high j), and ${\cal L}(J,\nu) = W_{2,R}(J,\nu)$ if S₂(I, J, K) is at the lower edge of the fine grid (low j), with

$$\begin{equation*} W_{2,L}(J,\nu) = \frac{1}{\Delta V_2(J)} \int _{x_2(J-1/2)}^{x_2(J-1/2+{\cal R})} \!dV_2; \end{equation*}$$

$$\begin{equation*} W_{2,R}(J,\nu) = \frac{1}{\Delta V_2(J)} \int _{x_2(J+1/2-{\cal R})}^{x_2(J+1/2)} \!dV_2. \end{equation*}$$

Like the function ${\cal G}^{\prime }$, ${\cal L}$ accounts for the portion of the coarse zone flux in volume space that is being replaced by fine fluxes. As before, expressions for the 2-direction are obtained by a suitable permutation of indices.

The changes required for prolongation on curvilinear grids deal entirely with making the interpolations conservative for curvilinear volumes and areas. In the case of 1D PPI, Equation (16) summarizing the conservation constraint is modified to

$$\begin{equation} Q(I,J,K) \Delta V_1(I) = \sum ^{\nu -1}_{\alpha =0}\, q(i+\alpha, j, k) \delta V_1(i), \end{equation} \tag{ A16 }$$

which then affects the rest of the scheme (Equations (17) and (18)):

$$\begin{equation} q^*_1(i+\alpha) = Q_L(I) + \zeta [ Q_R(I) - Q_L(I) + {\cal H}_1^{\prime } (1 - \zeta ) ] , \end{equation} \tag{ A17 }$$

where

$$\begin{eqnarray*} \zeta &=& \frac{x_1(i+\alpha) - x_1(I)}{\Delta x_1(I)};\nonumber \\ {\cal H}_1^{\prime } &=& \frac{1}{f_1^{\prime }(\nu,i) - f_2^{\prime }(\nu,i)}\lbrace \Delta V_1(I) [Q(I,J,K) - Q_L(I)]\nonumber\\ && - f_1^{\prime }(\nu,i) [Q_R(I) - Q_L(I) ] \rbrace; \nonumber \\ f_1^{\prime }(\nu,i) &=& \frac{1}{2\nu } \sum _{\xi =1}^{\nu } (2\xi - 1)\delta V_1(i+\xi -1); \nonumber \\ f_2^{\prime }(\nu,i) &=& \frac{1}{4\nu ^2} \sum _{\xi =1}^{\nu } (2\xi - 1)^2\delta V_1(i+\xi -1). \nonumber \end{eqnarray*}$$

As for PLI, the original 1D scheme in the 1-direction takes the form (van Leer 1977)

$$\begin{equation} q^*_1(i+\alpha) = Q(I) + (2\zeta -1)\,\Delta Q^{\prime }, \end{equation} \tag{ A18 }$$

where

$$\begin{equation*} \Delta Q^{\prime } = \left\lbrace {\begin{array}{@{}ll}\frac{\Delta Q_{\rm R}\, \Delta Q_{\rm L}}{\Delta Q_{\rm R} + \Delta Q_{\rm L}} & {\hbox{if } \Delta Q_{\rm R}\, \Delta Q_{\rm L} > 0}; \\ [12pt]{0,} & {\hbox{otherwise}},\end{array}}\right. \end{equation*}$$

and where

$$\begin{equation*} \Delta Q_{\rm R} = Q(I+1)-Q(I); \quad \Delta Q_{\rm L} = Q(I)-Q(I-1). \end{equation*}$$

Equation (A18) will not generally conserve mass, momentum, or magnetic flux on a curvilinear grid. To correct this, we relax the constraint that the linear interpolation profile must pass through Q(I) at the zone center. This releases one degree of freedom that can then be used with Equation (A16), resulting in

$$\begin{eqnarray} &&q^*_1(i+\alpha) = Q(I) + \Delta Q^{\prime }\left[ 2\zeta -1 - \frac{1}{\Delta V_1(I)}\,f_3^{\prime }(\nu,i)\right] ,\nonumber\\ \end{eqnarray} \tag{ A19 }$$

where

$$\begin{equation*} f_3^{\prime }(\nu,i) = \frac{1}{\nu }\,\sum _{\xi =1}^{\nu } \, [(2\xi - 1) - \nu ]\,\delta V_1(i+\xi -1).\nonumber \end{equation*}$$

The additional term in Equation (A19) can be viewed either as a shift in the zone-centered intercept or as choosing a modified value of ζ corresponding to the volume-centered rather than the spatial-centered coordinate.

Unfortunately, this "shift" can push q*₁(i + α) beyond neighboring values of Q(I) at the edges of the interpolation profile, resulting in a loss of monotonicity. However, since ζ = 0 or 1 at the left or right side of the interpolation profile, we can write

$$\begin{eqnarray} \Delta Q_{\rm L}^{\prime } &=& Q(I) - q^*_1(i+\alpha) = \Delta Q^{\prime } (1 + {\cal K}_1^{\prime }) \nonumber\\ {\rm or}\quad \Delta Q_{\rm R}^{\prime } &=& q^*_1(i+\alpha) - Q(I) = \Delta Q^{\prime } (1 - {\cal K}_1^{\prime }), \quad\quad \end{eqnarray} \tag{ A20 }$$

where ${\cal K}_1^{\prime } = f_3^{\prime }(\nu,i)\, / \, \Delta V_1(I)$. If we limit the slope of the interpolation profile to

$$\begin{equation} \Delta Q^{\prime \prime } = \left\lbrace\hspace{-2pt} \begin{array}{@{}ll}\,\hbox{sign}\, \big(\Delta Q^{\prime }\big)\, |\Delta Q_{\rm L}|\, / \, (1 + {\cal K}_1) & \hbox{if}\ |\Delta Q_{\rm L}^{\prime }| > |\Delta Q_{\rm L}| \\ \,\hbox{sign}\,\big (\Delta Q^{\prime }\big)\, |\Delta Q_{\rm R}|\, / \, (1 - {\cal K}_1) & \hbox{if}\ |\Delta Q_{\rm R}^{\prime }| > |\Delta Q_{\rm R}| \\ \phantom{x}\Delta Q^{\prime } & \hbox{otherwise}, \end{array}\right. \end{equation} \tag{ A21 }$$

and, because of the properties of the PLI slope ΔQ', the first two cases of Equation (A21) will never occur simultaneously. Thus, the equation for conservative, monotonic, 1D PLI in curvilinear coordinates is

$$\begin{equation} q^*_1(i+\alpha) = Q(I) + \Delta Q^{\prime \prime } [(2\zeta -1) - {\cal K}_1]. \end{equation} \tag{ A22 }$$

The previous discussion applies to interpolation of momenta in the directions perpendicular to the component normal (e.g., s₁ in the 2- and 3-directions). The covariant procedure for linear interpolation in between coarse faces (e.g., s₁ in the 1-direction; Equation (24)) is as follows:

$$\begin{eqnarray} && s_1(i+\alpha,j+\beta,k+\eta)\,\delta V_1(i+\alpha) = (1 - \zeta)\, s_1(i,j+\beta,k+\eta)\,\nonumber \\ && \quad \times \delta V_1(i)+ \zeta \, s_1(i+\nu,j+\beta,k+\eta)\,\delta V_1(i+\nu). \end{eqnarray} \tag{ A23 }$$

Finally, the generalized Li & Li (2004) algorithm is adapted to curvilinear coordinates simply by replacing $\vec{B}$ (and $\vec{b}$) in Equations (25)–(28) with the magnetic "flux functions":

$$\begin{equation} \vec{\Psi } = \big (\Psi _1, \Psi _2, \Psi _3 \big) = \left(g_{21}\,g_{31}\,\frac{\Delta V_2}{\Delta x_2} B_1,\, g_{31}\,g_{32}\, B_2,\, g_{21}\, B_3 \right). \end{equation} \tag{ A24 }$$

With this modification, the solenoidal condition can be written in "Cartesian-like" form, regardless of the geometry:

$$\begin{equation*} \nabla \cdot \vec{\Psi } = \frac{\partial \Psi _1}{\partial x_1} + \frac{\partial \Psi _2}{\partial x_2} + \frac{\partial \Psi _3}{\partial x_3} = 0, \end{equation*}$$

and the prolongation of magnetic field in curvilinear coordinates proceeds exactly as described in Section 4.1.