THE ATHENA ASTROPHYSICAL MAGNETOHYDRODYNAMICS CODE IN CYLINDRICAL GEOMETRY

Expanding the derivative operators in cylindrical coordinates in Equation (1a), we have the continuity equation in conservative variable form:

$$\begin{equation} \partial _t \rho + {\frac{1}{R}}\partial _R (R \,\rho v_R) + {\frac{1}{R}}\partial _\phi (\rho v_\phi) + \partial _z (\rho v_z) = 0. \end{equation} \tag{ 10 }$$

For the momentum Equation (1b), we have, in terms of the symmetric momentum tensor M (Equation (5)):

$$\begin{equation} \partial _t(\rho v_R) + {\frac{1}{R}}\partial _R (R M_{RR}) + {\frac{1}{R}}\partial _\phi M_{\phi R} + \partial _z M_{zR} = {\frac{1}{R}}M_{\phi \phi } \end{equation} \tag{ 11a }$$

$$\begin{equation} \partial _t(\rho v_\phi) + {\frac{1}{R}}\partial _R (R M_{R\phi }) + {\frac{1}{R}}\partial _\phi M_{\phi \phi } + \partial _z M_{z\phi } = -{\frac{1}{R}}M_{\phi R} \end{equation} \tag{ 11b }$$

$$\begin{equation} \partial _t(\rho v_z) + {\frac{1}{R}}\partial _R (R M_{Rz}) + {\frac{1}{R}}\partial _\phi M_{\phi z} + \partial _z M_{zz} = 0. \end{equation} \tag{ 11c }$$

Mignone et al. (2007) note that the symmetric character of M allows further simplification of Equation (11b) in cylindrical coordinates, since

$$\begin{equation} (\nabla \cdot \mathbf {M})_\phi = {\frac{1}{R^2}} \partial _R (R^2 \,M_{R\phi }) + {\frac{1}{R}}\partial _\phi M_{\phi \phi } + \partial _z M_{z\phi }. \end{equation} \tag{ 12 }$$

This leads to the so-called angular momentum conserving form of the ϕ-momentum equation:

$$\begin{eqnarray} &&\partial _t(\rho \,R \,v_\phi) + {\frac{1}{R}}\partial _R (R^2 M_{R\phi }) + {\frac{1}{R}}\partial _\phi (R M_{\phi \phi }) + \partial _z (R M_{z\phi }) = 0.\nonumber\\ \end{eqnarray} \tag{ 13 }$$

Note that in this form, the conserved quantity is angular momentum and there is no source term. However, since the original Cartesian version of Athena makes use of linear momenta in the flux calculations and since we do not wish to alter those calculations, we rewrite this equation once more to obtain the ϕ-momentum equation that we shall use:

$$\begin{equation} \partial _t(\rho \,v_\phi) + {\frac{1}{R^2}} \partial _R (R^2 \,M_{R\phi }) + {\frac{1}{R}}\partial _\phi M_{\phi \phi } + \partial _z M_{z\phi } = 0. \end{equation} \tag{ 14 }$$

This leaves the term M_ϕϕ/R = (ρv²_ϕ − B²_ϕ + P*)/R appearing in the R-momentum Equation (11a) as the only geometric source term in the cylindrical coordinate expansion of Equation (1b). In practice, this source term is often (partially) balanced by the radial component of the gravitational source term

$$\begin{equation} -\rho \,\boldsymbol{\nabla } \Phi \end{equation} \tag{ 15 }$$

for many astrophysical applications.

For the total energy Equation (1c) in conservative form, we have

$$\begin{eqnarray} &&\partial _t E + {\frac{1}{R}}\partial _R [ R \,((E + P^*) v_R - B_R \,(\boldsymbol{B} \cdot \boldsymbol{v})) ] + {\frac{1}{R}}\partial _\phi ((E + P^*) v_\phi \nonumber \\ &&\quad - B_\phi (\boldsymbol{B} \cdot \boldsymbol{v})) + \partial _z ((E + P^*) v_z - B_z (\boldsymbol{B} \cdot \boldsymbol{v})) = 0. \end{eqnarray} \tag{ 16 }$$

The gravitational source term for the energy equation is

$$\begin{equation} -\rho \,\boldsymbol{v} \cdot \boldsymbol{\nabla } \Phi. \end{equation} \tag{ 17 }$$

Finally, for the induction Equation (1d), we have in terms of the antisymmetric induction tensor, J (Equation (6)):

$$\begin{equation} \partial _t B_R + {\frac{1}{R}}\partial _\phi J_{\phi R} + \partial _z J_{zR} = 0 \end{equation} \tag{ 18a }$$

$$\begin{equation} \partial _t B_\phi + {\frac{1}{R}}\partial _R (R \,J_{R\phi }) + \partial _z J_{z\phi } = {\frac{J_{R \phi }}{R}} \end{equation} \tag{ 18b }$$

$$\begin{equation} \partial _t B_z + {\frac{1}{R}}\partial _R (R \,J_{Rz}) + {\frac{1}{R}}\partial _\phi J_{\phi z} = 0. \end{equation} \tag{ 18c }$$

It is important for the preservation of the divergence constraint that we avoid source terms in the magnetic fluxes. Thus, we rewrite the ϕ-induction Equation (18b) as

$$\begin{equation} \partial _t \left({\frac{B_\phi }{R}}\right) + {\frac{1}{R}}\partial _R \left(R \,{\frac{J_{R\phi }}{R}} \right) + \partial _z \left({\frac{J_{z\phi }}{R}}\right) = 0. \end{equation} \tag{ 19 }$$

Note that in this form, the conserved quantity is B_ϕ/R and there is no source term. Equation (19) has also modified the fluxes of B_ϕ, but since Athena does not actually use flux differences to evolve the magnetic fields (see Section 8), this is not a serious problem. We will use a reduced form equivalent to Equation (19) for B_ϕ in the reconstruction step (see Section 4):

$$\begin{equation} \partial _t B_\phi + \partial _R J_{R\phi } + \partial _z J_{z\phi } = 0. \end{equation} \tag{ 20 }$$

In this form, we use the fluxes J_Rϕ and J_zϕ originally appearing in Equation (18b), not the modified fluxes J_Rϕ/R and J_zϕ/R appearing in Equation (19).

Expanding the derivative operators in cylindrical coordinates in Equation (1a) and projecting in the R-direction, we have for the continuity equation in primitive variable form:

$$\begin{equation} \partial _t \rho + \rho \partial _R v_R + v_R \partial _R \rho = -{\frac{1}{R}}\rho v_R. \end{equation} \tag{ 24 }$$

The left-hand side of Equation (24) contains all the terms from Equation (21), and the term on the right-hand side is the first component of the geometric source term vector, $\boldsymbol{s}_{\rm geom}$. Furthermore, if we make the substitution R ↦ x and ignore the source term, we recover the x-projection of the continuity equation in Cartesian coordinates.

For the momentum equation, we begin with the conservative form of Equation (1b) and use the continuity equation and divergence-free constraint to eliminate terms and obtain

$$\begin{equation} \rho \,\partial _t \boldsymbol{v} + \rho (\boldsymbol{v} \cdot \nabla)\boldsymbol{v} - (\boldsymbol{B} \cdot \nabla)\boldsymbol{B} + \nabla P^* = 0. \end{equation} \tag{ 25 }$$

By explicitly enforcing $\nabla \cdot \boldsymbol{B}=0$ here, we ensure that any numerical error in the divergence of the magnetic field cannot influence the evolution of momentum during the reconstruction step.

Next, we divide through by ρ, substitute P* = P + B²/2, project in the R-direction, expand the partials, and move the source terms to the right-hand side to obtain

$$\begin{eqnarray} &&\partial _t v_R + v_R \partial _R v_R + {\frac{1}{\rho }} \partial _R P + {\frac{1}{\rho }} B_\phi \partial _R B_\phi + {\frac{1}{\rho }} B_z \partial _R B_z\nonumber\\ &&\quad = {\frac{1}{R}} \left(v_\phi ^2 - {\frac{1}{\rho }} B_\phi ^2\right), \end{eqnarray} \tag{ 26a }$$

$$\begin{equation} \partial _t v_\phi + v_R \partial _R v_\phi - {\frac{1}{\rho }} B_R \partial _R B_\phi = -{\frac{1}{R}}(v_\phi v_R - {\frac{1}{\rho }} B_\phi B_R), \end{equation} \tag{ 26b }$$

$$\begin{equation} \partial _t v_z + v_R \partial _R v_z - {\frac{1}{\rho }} B_R \partial _R B_z = 0. \end{equation} \tag{ 26c }$$

Recall that the ϕ-momentum Equation (13) can be expressed in angular momentum conserving form and thus avoid a geometric source term. However, we must include the source term on the right-hand side of Equation (26b) in primitive variable form in order to preserve the specific structure of the coefficient matrix, A, on the left-hand side of Equation (21). Finally, the gravity source terms in the momentum equation are given by the components of $-\boldsymbol{\nabla } \Phi$ in cylindrical coordinates.

We begin with the internal energy equation in coordinate-free form:

$$\begin{equation} \partial _t P + \boldsymbol{v} \cdot \nabla P + \gamma P \,\nabla \cdot \boldsymbol{v} = 0. \end{equation} \tag{ 27 }$$

Then, projecting the equations in the R-direction, expanding the partials, and moving the source term to the right-hand side, we obtain

$$\begin{equation} \partial _t P + v_R \partial _R P + \gamma P \partial _R v_R = -{\frac{1}{R}}\gamma P v_R. \end{equation} \tag{ 28 }$$

In primitive form, there is no gravity source term in the energy equation.

For the induction Equation (1d), also written as $\partial _t \boldsymbol{B} - \nabla \times (\boldsymbol{v} \times \boldsymbol{B}) = 0$, we begin with the components given in Equations (18a), (18c), and (20). Moving terms proportional to R⁻¹∂_R(RB_R), R⁻¹∂_ϕB_ϕ, and ∂_zB_z to the right-hand side, we obtain

$$\begin{eqnarray} &&\partial _t B_R + {\frac{1}{R}}\partial _\phi (v_\phi B_R) - B_\phi {\frac{1}{R}}\partial _\phi v_R + \partial _z (v_z B_R) - B_z \partial _z v_R\nonumber\\ &&\quad = v_R \left[ {\frac{1}{R}}\partial _\phi B_\phi + \partial _z B_z \right], \end{eqnarray} \tag{ 29a }$$

$$\begin{eqnarray} &&\partial _t B_\phi + \partial _R (v_R B_\phi) - B_R \partial _R v_\phi + \partial _z (v_z B_\phi) - B_z \partial _z v_\phi\nonumber\\ &&\quad = v_\phi \left[ {\frac{1}{R}}\partial _R (R B_R) + \partial _z B_z \right] - {\frac{1}{R}}v_\phi B_R, \end{eqnarray} \tag{ 29b }$$

$$\begin{eqnarray} &&\partial _t B_z + {\frac{1}{R}}\partial _R [ R (v_R B_z)] - B_R \partial _R v_z + {\frac{1}{R}}\partial _\phi (v_\phi B_z) - B_\phi {\frac{1}{R}}\partial _\phi v_z\nonumber\\ &&\quad = v_z \left[{\frac{1}{R}}\partial _R (R B_R) + {\frac{1}{R}}\partial _\phi B_\phi \right]. \end{eqnarray} \tag{ 29c }$$

In two dimensions, the divergence-free constraint in cylindrical coordinates with ∂_z ≡ 0 implies that the right-hand side of Equation (29c) is identically zero. Applying the divergence constraint in Equation (29c), projecting in the R-direction, and expanding the R-partials on the left-hand side, we obtain

$$\begin{equation} \partial _t B_R = 0, \end{equation} \tag{ 30a }$$

$$\begin{eqnarray} &&\partial _t B_\phi + B_\phi \partial _R v_R - B_R \partial _R v_\phi + v_R \partial _R B_\phi = v_\phi {\frac{1}{R}}\partial _R (R B_R)\nonumber\\ &&\quad - {\frac{1}{R}}v_\phi B_R, \end{eqnarray} \tag{ 30b }$$

$$\begin{equation} \partial _t B_z + B_z \partial _R v_R - B_R \partial _R v_z + v_R \partial _R B_z = -{\frac{1}{R}}v_R B_z. \end{equation} \tag{ 30c }$$

Note that the left-hand sides of Equations (30b) and (30c) exactly match the Cartesian form (GS05, Equation (30)) if we make the substitution R ↦ x. Also, the divergence term on the right-hand side of Equation (30b) matches the divergence term on the right-hand side of the Cartesian form, except that it appears in cylindrical coordinate form. However, the additional source terms −v_ϕB_R/R in Equation (30b) and −v_RB_z/R in Equation (30c), which vanish as R → ∞, are curvature-related terms that are unique to the system in cylindrical coordinates.

In three dimensions, the cancellation of the $\nabla \cdot \boldsymbol{B}$ terms on the right-hand side of Equation (29c) is no longer possible. Instead, GS08 introduce an algorithm that adds a limited amount of the MHD source terms from the transverse directions to the source terms in each splitting direction. This is done in such a way that the overall induction equation is not altered and so that the sum of the MHD source terms is minimized. These constraints take the form of minmod limiter functions that reduce to the underlying two-dimensional algorithm in the limit of two-dimensional, grid-aligned problems (see GS08, Section 3). The three-dimensional algorithm in cylindrical coordinates yields the system

$$\begin{eqnarray} &&\partial _t B_R + \left\lbrace {\frac{1}{R}}\partial _\phi (v_\phi B_R - B_\phi v_R) - v_R L_{R\phi } (\partial _z B_z) \right\rbrace \nonumber \\ &&\quad+ \bigg\lbrace \partial_z(v_z B_R - B_z v_R) - v_R L_{Rz} \bigg({\frac{1}{R}}\partial _\phi B_\phi \bigg) \bigg\rbrace = 0, \end{eqnarray} \tag{ 31a }$$

$$\begin{eqnarray} &&\partial _t B_\phi + \lbrace \partial _R(v_R B_\phi - B_R v_\phi) - v_\phi L_{\phi R} (\partial _z B_z) \rbrace\nonumber \\ &&\quad+ \left\lbrace \partial _z(v_z B_\phi - B_z v_\phi) - v_\phi L_{\phi z} \left[{\frac{1}{R}}\partial _R (R B_R)\right] \right\rbrace = 0,\qquad \end{eqnarray} \tag{ 31b }$$

$$\begin{eqnarray} &&\partial _t B_z + \left\lbrace {\frac{1}{R}}\partial _R [ R(v_R B_z - B_R v_z) ] - v_z L_{zR} \left({\frac{1}{R}}\partial _\phi B_\phi\right) \right\rbrace\nonumber \\ &&\quad + \left\lbrace {\frac{1}{R}}\partial _\phi (v_\phi B_z - B_\phi v_z) - v_z L_{z \phi } \left[{\frac{1}{R}}\partial _R (R B_R)\right] \right\rbrace = 0,\nonumber\\ \end{eqnarray} \tag{ 31c }$$

where

$$\begin{equation} L_{R\phi } (\partial _z B_z) \equiv {\rm minmod} \left(-{\frac{1}{R}}\partial _\phi B_\phi, \;\partial _z B_z\right), \end{equation} \tag{ 32a }$$

$$\begin{equation} L_{\phi R} (\partial _z B_z) \equiv {\rm minmod} \left(-{\frac{1}{R}}\partial _R (R B_R), \;\partial _z B_z\right), \end{equation} \tag{ 32b }$$

$$\begin{equation} L_{zR} \left({\frac{1}{R}}\partial _\phi B_\phi\right) \equiv {\rm minmod} \left(-{\frac{1}{R}}\partial _R (R B_R), \;{\frac{1}{R}}\partial _\phi B_\phi\right). \end{equation} \tag{ 32c }$$

Note that the limiter L_ij is only applied to the equation for B_i projected in the j-direction. We require that

$$\begin{equation} L_{Rz} \left({\frac{1}{R}}\partial _\phi B_\phi\right) = -L_{R\phi } (\partial _z B_z), \end{equation} \tag{ 33a }$$

$$\begin{equation} L_{\phi z} \left[{\frac{1}{R}}\partial _R (R B_R)\right] = -L_{\phi R} (\partial _z B_z), \end{equation} \tag{ 33b }$$

$$\begin{equation} L_{z \phi } \left[{\frac{1}{R}}\partial _R (R B_R)\right] = -L_{zR} \left({\frac{1}{R}}\partial _\phi B_\phi\right), \end{equation} \tag{ 33c }$$

so that the limiters cancel pairwise when summed over all projections. The limiters defined in Equations (32) and (33) are the same as in the three-dimensional Cartesian formulae (GS08, Equations (11) and (15)), with ∂_xB_x ↦ R⁻¹∂_R(RB_R) and ∂_yB_y ↦ R⁻¹∂_ϕB_ϕ.

Projecting Equations (31) in the R-direction, expanding the partials (except for those appearing in the limiter functions), moving all of the source terms to the right-hand side, and using the properties of the minmod function together with the $\nabla \cdot \boldsymbol{B}=0$ constraint, we obtain

$$\begin{equation} \partial _t B_R = 0, \end{equation} \tag{ 34a }$$

$$\begin{eqnarray} &&\partial _t B_\phi + B_\phi \partial _R v_R + v_R \partial _R B_\phi - B_R \partial _R v_\phi\nonumber \\ &&\quad= v_\phi \;{\rm minmod} \left[{\frac{1}{R}}\partial _R (R B_R), \;-{\frac{1}{R}}\partial _\phi B_\phi \right] - {\frac{1}{R}}v_\phi B_R,\qquad \end{eqnarray} \tag{ 34b }$$

$$\begin{eqnarray} &&\partial _t B_z + B_z \partial _R v_R + v_R \partial _R B_z - B_R \partial _R v_z\nonumber \\ &&\quad = v_z \;{\rm minmod} \left[{\frac{1}{R}}\partial _R (R B_R), \;-\partial _z B_z\right] - {\frac{1}{R}}v_R B_z.\quad \end{eqnarray} \tag{ 34c }$$

Note that for the two-dimensional case with ∂_z ≡ 0, the $\nabla \cdot \boldsymbol{B}$ constraint implies that the arguments of the minmod function in Equation (34b) are equal and that the minmod function in Equation (34c) evaluates to zero, so that we recover the two-dimensional system in Equations (30). The minmod terms on the right-hand side of Equations (34b) and (34c) are analogous to the corresponding terms in Cartesian coordinates derived in GS08, and the remaining terms are geometric.

In summary, for the primitive variable equations in cylindrical coordinates, the MHD source term vectors are given by Equations (18) and (19) of GS08 via the substitutions ∂_xB_x ↦ R⁻¹∂_R(RB_R) and ∂_yB_y ↦ R⁻¹∂_ϕB_ϕ, and the geometric source term vector is given by

$$\begin{equation} \boldsymbol{s}_{\rm geom} \equiv \dsty \left[\begin{array}{@{}c@{}} -{\dsty\frac{1}{R}}\rho v_R \\[8pt] \dsty{\frac{1}{R}}\big(v_\phi ^2 - {\frac{1}{\rho }}B_\phi ^2\big) \\[8pt] \dsty-{\frac{1}{R}}\big(v_\phi v_R - {\frac{1}{\rho }}B_\phi B_R\big) \\[8pt] \dsty0 \\[4pt] \dsty-{\frac{1}{R}}\gamma P v_R \\[8pt] \dsty-{\frac{1}{R}}v_\phi B_R \\[8pt] \dsty-{\frac{1}{R}}v_R B_z \end{array} \right]. \end{equation} \tag{ 35 }$$

Since the geometric source terms arise directly from the scale factors in the R-partials, we associate the geometric source term $\boldsymbol{s}_{\rm geom}$ exclusively with the R-direction. Note that $\Vert \boldsymbol{s}_{\rm geom} \Vert \rightarrow 0$ in the limit of vanishing curvature, i.e., as R → ∞.

We emphasize that the geometric source terms in Equation (35) are used only in obtaining the L/R states, not for the final FV update. Finally, the gravity source terms for the L/R states are given by the cylindrical coordinate components of $-\boldsymbol{\nabla } \Phi$ in the momentum equation, and there is no gravity source term in the energy equation.

The piecewise linear method (PLM) of reconstruction approximates each primitive variable by defining in the ith zone,

$$\begin{equation} a(s) \equiv a_{L,i} + s \,\Delta a_i \equiv a_{R,i} - (1-s) \,\Delta a_i, \end{equation} \tag{ 38 }$$

where Δa_i ≡ a_R,i − a_L,i represents the difference of some quantity a over the zone, and a_L,i and a_R,i are the values of a at the left and right interfaces of the zone, respectively. Thus, to specify a(s) completely, we need only to define Δa_i and a_L,i for each zone as functions of the volume averages, a_i.

5.1.1. PLM in Cartesian Coordinates

From the consistency requirement in Equation (36) with Cartesian coordinates,

$$\begin{equation} a_i = \frac{1}{\Delta x} \int _{x_{i-1/2}}^{x_{i+1/2}} a(x) \,dx = \int _0^1 a(s) \,ds. \end{equation} \tag{ 39 }$$

Substituting Equation (38) into Equation (39) and integrating, we obtain

$$\begin{equation} a_i = a_{L,i} + {\frac{1}{2}}\Delta a_i = a_{R,i} - {\frac{1}{2}}\Delta a_i, \end{equation} \tag{ 40 }$$

from which

$$\begin{equation} a_{L,i} = a_i - {\frac{1}{2}}\Delta a_i, \end{equation} \tag{ 41a }$$

$$\begin{equation} a_{R,i} = a_i + {\frac{1}{2}}\Delta a_i. \end{equation} \tag{ 41b }$$

The usual Cartesian formulae for the differences over the zone are

$$\begin{equation} \Delta a_i \equiv \left\lbrace \begin{array}{@{}ll} \dsty{\frac{1}{2}}\left(a_{i+1} - a_{i-1} \right), &{\rm centered} \\[3pt] \dsty\left(a_{i+1} - a_i \right), &{\rm forward} \\[3pt] \dsty\left(a_i - a_{i-1} \right), &{\rm backward}\end{array} \right.. \end{equation} \tag{ 42 }$$

It is clear from Equation (42) that constant and linear profiles are reconstructed exactly. In order to make the reconstruction TVD, these differences (or "slopes," as they are commonly called) are limited using a slight variant of the monotonized central-difference (MC) limiter as described by Leveque (2002), and flattened to avoid the introduction of new extrema. Leveque argues that limiting should be performed in characteristic variables so that the accuracy of the reconstruction for smooth wave families is not adversely affected by limiting in other non-smooth wave families. This requires a special description for systems of conservation laws including a bounded linear transformation from primitive to characteristic variables, with inverse transformation from characteristic variables back to primitive. We do not go into detail here, but note that the inverse transformation is not guaranteed to be monotonicity-preserving, hence an additional monotonization is performed on the resulting primitive variable differences. In the previous implementation, this was done in a non-conservative manner, and we have since implemented a related scheme which preserves the consistency requirement of Equation (36).

Finally, by applying the monotonized, limited differences Δa_i in Equations (40) and (41), for smooth profiles we obtain second-order accurate approximations to the values of a across each zone at time tⁿ, except possibly at local extrema.

5.1.2. PLM in Cylindrical Coordinates

From the consistency requirement of Equation (36) with cylindrical coordinates,

$$\begin{eqnarray} a_i &=& \frac{1}{R_i \,\Delta R} \int _{R_{i-1/2}}^{R_{i+1/2}} a(R) RdR\nonumber\\ &=& \frac{1}{R_i} \int _0^1 a(s) (R_{i-1/2} + s \,\Delta R) \, ds. \end{eqnarray} \tag{ 43 }$$

Substituting Equation (38) into Equation (43) and integrating, we obtain

$$\begin{equation} a_i = a_{L,i} + {\frac{1}{2}}\Delta a_i \left(1 + \gamma _i \right) = a_{R,i} - {\frac{1}{2}}\Delta a_i \left(1 - \gamma _i \right), \end{equation} \tag{ 44 }$$

where

$$\begin{equation} \gamma _i \equiv \frac{\Delta R}{6 \,R_i} \end{equation} \tag{ 45 }$$

is a correction factor for curvature. Note that γ_i → 0 for fixed ΔR as R → ∞, or for fixed R as ΔR → 0, and from the formulae in Equation (44), we recover the Cartesian formulae in Equation (40). Solving for a_L,i and a_R,i, we obtain

$$\begin{equation} a_{L,i} = a_i - {\frac{1}{2}}\Delta a_i \left(1 + \gamma _i \right), \end{equation} \tag{ 46a }$$

$$\begin{equation} a_{R,i} = a_i + {\frac{1}{2}}\Delta a_i \left(1 - \gamma _i \right). \end{equation} \tag{ 46b }$$

Next, we wish to define the differences Δa_i such that constant and linear profiles are reconstructed exactly. Assuming a linear profile, say a(R) = CR with some constant slope C, we enforce the consistency relation of Equation (43) to obtain

$$\begin{equation} a_i \equiv \left< a \right>_i = \left< C R \right>_i = C \left< R \right>_i, \end{equation} \tag{ 47 }$$

where

$$\begin{equation} \left< R \right>_i = R_i + \frac{(\Delta R)^2}{12\,R_i} \end{equation} \tag{ 48 }$$

is the R-coordinate of the volume centroid of the ith zone. To obtain an exact reconstruction for a linear profile, we require that Δa_i = C ΔR. First, we consider the Cartesian formula for a centered difference to obtain

$$\begin{eqnarray} {\frac{1}{2}}(a_{i+1} - a_{i-1}) &=& {\frac{1}{2}}\left(C \left<R_{i+1}\right> - C\left<R_{i-1}\right> \right)\nonumber\\ &=& C \,\Delta R \,\left(1 - \frac{(\Delta R)^2}{12 \,R_{i+1} \,R_{i-1}} \right). \end{eqnarray} \tag{ 49 }$$

If we divide the Cartesian centered-difference formula in Equation (42) by the term in parentheses from Equation (49), we obtain the desired difference, C ΔR, exactly. This—along with similar calculations for the forward- and backward-difference slopes—suggests the following definitions:

$$\begin{eqnarray} \Delta a_i \equiv \left\lbrace \begin{array}{@{}ll} {{\dsty\frac{1}{2}}\left(a_{i+1} - a_{i-1} \right)}/{\left(1 - \frac{(\Delta R)^2}{12 \,R_{i+1} \,R_{i-1}} \right)}, &{\rm centered} \\[3pt] \dsty{\left(a_{i+1} - a_i \right)}/{\left(1 - \frac{(\Delta R)^2}{12 \,R_{i+1} \,R_i} \right)}, &{\rm forward} \\[3pt] \dsty{\left(a_i - a_{i-1} \right)}/{\left(1 - \frac{(\Delta R)^2}{12 \,R_i \,R_{i-1}} \right)}, &{\rm backward}\end{array} \right.\!\!\!.\nonumber\\ \end{eqnarray} \tag{ 50 }$$

It is clear from Equation (50) that constant profiles yield Δa_i = 0, hence these are also reconstructed exactly. In order to make the reconstruction TVD, these differences are limited and monotonized as in the Cartesian case. Finally, by applying the resulting differences to Equations (38) and (46), for smooth profiles we obtain second-order accurate approximations to the values of a across each zone at time tⁿ, except possibly at local extrema.

The piecewise parabolic method (PPM) of reconstruction approximates the profile of each primitive variable in the ith zone as

$$\begin{eqnarray} a(s) &\equiv& a_{L,i} + s \,\Delta a_i + s(1-s)a_{6,i}\nonumber\\ &\equiv& a_{R,i} - (1-s) \,\Delta a_i + s(1-s)a_{6,i}, \end{eqnarray} \tag{ 51 }$$

where, as for linear reconstruction, Δa_i ≡ a_R,i − a_L,i represents the average difference of some quantity a over the zone, and a_L,i and a_R,i are the values of a at the left and right interfaces of the zone, respectively. The term a_6,i is the so-called parabolic coefficient (see CW, Equation (1.5)). To specify a(s) completely, we need only to define Δa_i, a_L,i, and a_6,i for each zone as functions of the set of volume averages, {a_i}.

5.2.1. PPM Reconstruction in Cylindrical Geometry

To satisfy the consistency requirement of Equation (36) in cylindrical coordinates, we substitute Equation (51) into Equation (43) and integrate, then solve for a_6,i in terms of the quantities a_i and Δa_i:

$$\begin{equation} a_{6,i} \equiv 6 \left[ a_i - a_{L,i} - {\frac{1}{2}}\Delta a_i \left(1 + \gamma _i \right) \right]. \end{equation} \tag{ 52 }$$

Equation (52) is equivalent to the definition of the parabolic coefficient appearing in Equation (18) of Blondin & Lufkin (1993, hereafter BL). Recall that γ_i → 0 for fixed ΔR as R → ∞ or for fixed R as ΔR → 0 (see Equation (45)), and in these limits we recover the Cartesian version of the parabolic coefficient, which is the same as Equation (52) except with γ_i = 0 (see CW, Equation (1.5)).

We wish to define the values of a_L,i, a_R,i, and Δa_i ≡ a_R,i − a_L,i such that constant, linear, and parabolic profiles are reconstructed exactly, or at least to second order. By constructing a quartic polynomial from the {a_i}, one can show (see BL, Section 3) that the cylindrical formulae can be obtained from the corresponding Cartesian formulae (see CW, Equation (1.6)) by making the canonical substitution a_l ↦ a_l R_l. This yields

$$\begin{eqnarray} &&a_{L,i} \,R_{i-1/2} = {\frac{1}{2}}\left(a_i \,R_i + a_{i-1} \,R_{i-1} \right) - {\frac{1}{6}} \left(\delta a_i \,R_i - \delta a_{i-1} \,R_{i-1} \right),\nonumber\\ \end{eqnarray}\vspace*{-12pt} \tag{ 53a }$$

$$\begin{eqnarray} &&a_{R,i} \,R_{i+1/2} = {\frac{1}{2}}\left(a_{i+1} \,R_{i+1} + a_i \,R_i \right) - {\frac{1}{6}} \left(\delta a_{i+1} \,R_{i+1} - \delta a_i \,R_i \right).\nonumber\\ \end{eqnarray} \tag{ 53b }$$

Here, the centered, forward, and backward differences in zone i are

$$\begin{equation} \delta a_i \equiv \left\lbrace \begin{array}{@{}ll} {\dsty\frac{1}{2}}\left(a_{i+1} \,R_{i+1} - a_{i-1} \,R_{i-1} \right)/R_i, &{\rm centered} \\[6pt] \dsty\left(a_{i+1} \,R_{i+1} - a_{i} \,R_{i} \right)/R_{i+1/2}, &{\rm forward} \\[6pt] \dsty\left(a_{i} \,R_{i} - a_{i-1} \,R_{i-1} \right)/R_{i-1/2}, &{\rm backward}\end{array} \right.\!\!\!, \end{equation} \tag{ 54a }$$

and similarly for the i + 1 and i − 1 zones. Note that the corresponding Cartesian formulae are given by taking R_i → 1 in Equations (53) and (54).

First, assuming a constant profile, a(R) = C, and taking the volume average across the ith zone, we have that a_i = C. Using centered differences, it follows that δa_i = C ΔR/R_i, hence from Equations (53), we see that a_L,i = C as desired. For forward- and backward differences, a constant profile yields a_L,i = C[1 + O((ΔR)³)].

Next, assuming a linear profile, a(R) = CR, and volume averaging, we have a_i given using Equations (47) and (48) above. Using either centered, forward, or backward differences, it follows that δa_i = 2C ΔR, hence from Equations (53), we see that a_L,i = CR_i−1/2 and a_R,i = CR_i+1/2, as desired.

Finally, assuming a parabolic profile, a(R) = CR², it is straightforward to show that the volume average of R across the ith zone is

$$\begin{equation} a_i = C \left< R^2 \right>_i = C \left(R_i^2 + \frac{(\Delta R)^2}{4} \right). \end{equation} \tag{ 55 }$$

Using centered differences, it follows that δa_i = C ΔR[3R_i +5(ΔR)²/(4R_i)], hence from Equations (53), we see that a_L,i = CR²_i−1/2, as desired. For forward- or backward differences, it can be shown that a_L,i = CR²_i−1/2[1 + O(ΔR/R_i−1/2)³]. Thus, we conclude that parabolic profiles can be recovered up to the required order.

As in Section 5.1, the slopes (meaning the centered, forward, and backward difference δa_i's) are monotonized in characteristic form using a slight variant of the MC limiter, as in Leveque (2002). Then, after computing the parabolic interpolant, the slopes are re-monotonized to ensure that the interpolation introduces no new extrema. Following BL, the values of a_L,i and a_R,i are reset to

$$\begin{equation} a^*_{L,i} = a^*_{R,i} = a_i \end{equation} \tag{ 56 }$$

whenever a_i is a local extremum with respect to a_L,i and a_R,i, or to

$$\begin{equation} a^*_{L,i} = \frac{6a_i - a_{R,i}\,(4 + 3\gamma _i)}{2 - 3\gamma _i}, \end{equation} \tag{ 57a }$$

$$\begin{equation} a^*_{R,i} = \frac{6a_i - a_{L,i}\,(4 - 3\gamma _i)}{2 + 3\gamma _i}, \end{equation} \tag{ 57b }$$

whenever they are close enough to a_i so that the parabolic interpolation function introduces new extrema. The test for this case,

$$\begin{equation} | a_{R,i} - a_{L,i} | \ge | a_{6,i} |, \end{equation} \tag{ 58 }$$

is geometry independent. As a result, for smooth profiles we obtain second-order accurate approximations to the values of a across each zone at time tⁿ, except possibly at local extrema. By taking γ_i = 0, we recover the Cartesian versions of Equations (57) (see CW, Equation (1.10)).

Here, we describe the evaluation of the L/R states at time t^n+1/2 based on a piecewise linear reconstruction of the underlying profile at time tⁿ, defined by Equations (38), (46), and (50) of Section 5.1. Substituting a(s) into Equations (63) and integrating, we obtain

$$\begin{equation} f^a_{L,i+1/2}(\chi _R) = a_{R,i} - {\frac{1}{2}}\chi _R \,\Delta a_i \,(1-\beta _{R,i}(\chi _R)) \end{equation} \tag{ 64a }$$

on the right side of the zone (left of the interface) and

$$\begin{equation} f^a_{R,i-1/2}(\chi _L) = a_{L,i} + {\frac{1}{2}}\chi _L \,\Delta a_i \,(1+\beta _{L,i}(\chi _L)) \end{equation} \tag{ 64b }$$

at the left of the zone (right of the interface). Here, a_L,i and a_R,i are the values of a at the left and right interfaces of the ith zone, respectively, Δa_i is the monotonized difference of a across the zone from Equation (50), and we have defined the functions

$$\begin{eqnarray} \beta _{R,i}(\chi _R) &\equiv & \frac{\chi _R \,\Delta R}{6(R_{i+1/2}-{\frac{1}{2}}\chi _R \,\Delta R) }, \end{eqnarray} \tag{ 65a }$$

$$\begin{eqnarray} \beta _{L,i}(\chi _L) &\equiv & \frac{\chi _L \,\Delta R}{6(R_{i-1/2}+{\frac{1}{2}}\chi _L \,\Delta R) }, \end{eqnarray} \tag{ 65b }$$

as additional correction factors due to the curvature of the zone. Note that β_R,i(χ_R), β_L,i(χ_L) → 0 as R → ∞, i.e., in the limit of vanishing curvature, in which case Equations (64) reduce to the Cartesian formulae. Note further that in averaging a(s) over the whole zone, i.e., taking χ_L/R = 1, we have β_R,i(1) = β_L,i(1) = γ_i (see Equation (45)), and from Equations (64), we recover the averages of Equation (44).

Here, we describe the evaluation of the L/R states at time t^n+1/2 based on a piecewise parabolic reconstruction of the underlying profile at time tⁿ, defined by Equations (51), (52), (53), and (54) of Section 5.2. Substituting a(s) into Equations (63) and integrating, we obtain expressions analogous to Equations (64) for the time average of the right(left)-moving waves over the domain of dependence of the left(right) interface state at R_i+1/2(R_i−1/2) upon these waves:

$$\begin{eqnarray} f^a_{L,i+1/2}(\chi _R) &=& a_{R,i} - {\frac{1}{2}}\chi _R \left[ \Delta a_i - \left(1 - {\frac{2}{3}}\chi _R \right) a_{6,i} \right] \nonumber \\ && + {\frac{1}{2}}\chi _R [ \Delta a_i - \left(1 - \chi _R \right) a_{6,i} ] \beta _{R,i}(\chi _R),\qquad \end{eqnarray} \tag{ 66a }$$

$$\begin{eqnarray} f^a_{R,i-1/2}(\chi _L) &=& a_{L,i} + {\frac{1}{2}}\chi _L \left[ \Delta a_i + \left(1 - {\frac{2}{3}}\chi _L \right) a_{6,i} \right] \nonumber \\ && + {\frac{1}{2}}\chi _L [ \Delta a_i + \left(1 - \chi _L \right) a_{6,i} ] \beta _{L,i}(\chi _L).\qquad \end{eqnarray} \tag{ 66b }$$

The functions β_R,i(χ_R) and β_L,i(χ_L) defined in Equations (65) are the correction factors due to the curvature of the zone. The PLM result in Equations (64) corresponds to setting a_6,i = 0. Note that Equations (66) also reduce to the Cartesian formulae (see CW, Equation (1.12)) when β → 0 as R → ∞. Note further that in averaging a(s) over the whole zone, i.e. taking χ_L/R = 1, we have β_R,i(1) = β_L,i(1) = γ_i, and the results in Equations (66) are consistent with Equation (52).

Once time-averaged L/R states have been obtained in the characteristic variables (as in Equations (64) or (66)), we convert back to the primitive variables using R, the matrix consisting of the right-eigenvectors of A. Since L transforms from primitive to characteristic variables via $\boldsymbol{a} = \mathbf {L}\boldsymbol{w}$, and since RL = I, where I is the standard identity matrix, $\mathbf {R}\boldsymbol{a}=\boldsymbol{w}$ accomplishes the inverse transformation.

Defining a^ν,n+1/2_L,i+1/2 ≡ f^a_L,i+1/2(χ^ν_R) and a^ν,n+1/2_R,i−1/2 ≡ f^a_R,i−1/2(χ^ν_L) for each characteristic variable, we obtain the total time-averaged L/R states in primitive variable form by summing the projections of these contributions onto the right eigenspace:

$$\begin{equation} \boldsymbol{w}^{n+1/2}_{L/R,i \pm 1/2} = \mathbf {R} \,\boldsymbol{a}^{n+1/2}_{L/R,i \pm 1/2} = \sum _{\nu = 1}^{M} a^{\nu,n+1/2}_{L/R,i \pm 1/2} \boldsymbol{r}^\nu. \end{equation} \tag{ 67 }$$

For example, in the Cartesian case (Stone et al. 2008, Equation (42)), using Equation (64a) or (66a) with β = 0, we have

$$\begin{eqnarray} \boldsymbol{w}^{n+1/2}_{L,i+1/2} &=& \boldsymbol{w}_{i}^n + \left(\frac{1}{2} - \frac{\lambda ^M \,\Delta t}{2\Delta x} \right) \left(\Delta \boldsymbol{w} \right)_i - \frac{\Delta t}{2\Delta x}\nonumber\\ &&\times\, \sum _\nu (\lambda ^\nu - \lambda ^M) [ \boldsymbol{l}^\nu \cdot (\Delta \boldsymbol{w})_i] \boldsymbol{r}^\nu. \end{eqnarray} \tag{ 68 }$$

As written in Equation (67), the inverse transformation includes a sum over all waves. However, following CW and Colella, contributions to $\boldsymbol{w}$ on a given interface from waves that propagate in the opposite direction may be discarded to yield a more robust solution for strongly nonlinear problems, i.e., the waves are upwinded in the appropriate direction. Thus, for left-moving waves with λ < 0, we may set χ_L/R = 0 in Equation (64a) to obtain f^a_L,i+1/2(0) = a_R,i. Similarly for right-moving waves with λ>0, we may set χ_L/R = 0 in Equation (64b) to obtain f^a_R,i−1/2(0) = a_L,i.

Stone et al. (2008) have noted that this upwinding destroys the formal second-order convergence for smooth flows. In this case, the one-dimensional L/R states are accurate to the desired order if and only if we account for all waves, including those propagating toward the interfaces from outside of the zone. With this approach, the sum in Equation (68) includes all λ^ν. This means that for a given zone we integrate an extrapolation of the local reconstruction profile over a domain of dependence that lies outside the zone. However, for non-smooth flows, we have found that this can lead to significant errors near discontinuities, since we may end up extrapolating the local reconstruction beyond a point of discontinuity into a region where it is no longer a good approximation of the profile. Thus, for flows that may contain discontinuities, we follow CW and restrict to integration only over characteristics that propagate toward the interfaces from the interior of the zone, as previously described. However, we do not use the reference states for waves propagating away from the interface. Thus, we compute the states using

$$\begin{eqnarray} &&\boldsymbol{w}^{n+1/2}_{L,i+1/2} = \boldsymbol{w}_{i}^n + \frac{1}{2} \left(\Delta \boldsymbol{w} \right)_i - \frac{\Delta t}{2\Delta x} \sum _{\nu :\;\lambda ^\nu > 0} \lambda ^\nu \,[ \boldsymbol{l}^\nu \cdot (\Delta \boldsymbol{w})_i] \boldsymbol{r}^\nu,\nonumber\\ \end{eqnarray}\vspace*{-12pt} \tag{ 69a }$$

$$\begin{eqnarray} &&\boldsymbol{w}^{n+1/2}_{R,i-1/2} = \boldsymbol{w}_{i}^n - \frac{1}{2} \left(\Delta \boldsymbol{w} \right)_i - \frac{\Delta t}{2\Delta x} \sum _{\nu :\;\lambda ^\nu < 0} \lambda ^\nu \,[ \boldsymbol{l}^\nu \cdot (\Delta \boldsymbol{w})_i] \boldsymbol{r}^\nu,\nonumber\\ \end{eqnarray} \tag{ 69b }$$

for Cartesian PLM, which we refer to as "upwind-only" integration. For PLM in cylindrical coordinates, factors (1 − β_R,i(ξ^ν_R)) and (1 + β_L,i(ξ^ν_L)) are included in the sums in Equations (69a) and (69b), respectively, using Equations (64a) and (64b). Equations (66a) and (66b) are used to obtain analogous expressions for PPM in cylindrical coordinates.

To see this, we rewrite the induction equation as

$$\begin{equation} \partial _t \boldsymbol{B} + \nabla \times {\mathbf {\mathcal {E}}}= 0, \end{equation} \tag{ 78 }$$

where ${\mathbf {\mathcal {E}}}= -\boldsymbol{v} \times \boldsymbol{B}$ is the electric field in ideal MHD (the electromotive force (EMF)). Then, integrating over the oriented bounding surface of grid cell (i, j, k) and applying Stokes' theorem, we find that

$$\begin{eqnarray} B^{n+1}_{R; \;i\pm 1/2,j,k} &=& B^{n}_{R; \;i\pm 1/2,j,k} - \frac{\Delta t}{R_{i\pm 1/2} \,\Delta \phi }\nonumber\\ &&\times\, \big({\mathcal {E}}^{n+ 1/2}_{z; \;i\pm 1/2,j+ 1/2,k} - {\mathcal {E}}^{n+ 1/2}_{z; \;i\pm 1/2,j- 1/2,k} \big), \nonumber \\ &&+ \frac{\Delta t}{\Delta z} \big({\mathcal {E}}^{n+ 1/2}_{\phi ; \;i\pm 1/2,j,k+ 1/2} - {\mathcal {E}}^{n+ 1/2}_{\phi ; \;i\pm 1/2,j,k- 1/2} \big)\nonumber\\ \end{eqnarray} \tag{ 79a }$$

$$\begin{eqnarray} B^{n+1}_{\phi ; \;i,j\pm 1/2,k} &=& B^{n}_{\phi ; \;i,j\pm 1/2,k} + \frac{\Delta t}{\Delta R}\nonumber\\ &&\times \big({\mathcal {E}}^{n+ 1/2}_{z; \;i+ 1/2,j\pm 1/2,k} - {\mathcal {E}}^{n+ 1/2}_{z; \;i- 1/2,j\pm 1/2.k} \big) \nonumber \\ &&- \frac{\Delta t}{\Delta z} \big({\mathcal {E}}^{n+ 1/2}_{R; \;i,j\pm 1/2,k+ 1/2} - {\mathcal {E}}^{n+ 1/2}_{R; \;i,j\pm 1/2,k- 1/2} \big),\nonumber\\ \end{eqnarray}\vspace*{-18pt} \tag{ 79b }$$

$$\begin{eqnarray} B^{n+1}_{z; \;i,j,k\pm 1/2} &=& B^{n}_{z; \;i,j,k\pm 1/2} - \frac{\Delta t}{R_i \,\Delta R}\nonumber\\ &&\times\, \big({\mathcal {E}}^{n+ 1/2}_{\phi ; \;i+ 1/2,j,k\pm 1/2} - {\mathcal {E}}^{n+ 1/2}_{\phi ; \;i- 1/2,j,k\pm 1/2} \big) \nonumber \\ &&+ \frac{\Delta t}{R_i \,\Delta \phi } \big({\mathcal {E}}^{n+ 1/2}_{R; \;i,j+ 1/2,k\pm 1/2} - {\mathcal {E}}^{n+ 1/2}_{R; \;i,j- 1/2,k\pm 1/2}\big),\nonumber\\ \end{eqnarray} \tag{ 79c }$$

where

$$\begin{eqnarray} B^{n}_{R; \;i\pm 1/2,j,k} &\equiv & \frac{1}{R_{i\pm 1/2} \,\Delta \phi \,\Delta z} \int _{z_{k- 1/2}}^{z_{k+ 1/2}}\!\!\! \int _{\phi _{j- 1/2}}^{\phi _{j+ 1/2}} B_R(R_{i\pm 1/2},\phi,z,t^n)\nonumber\\ &&\times \,R_{i\pm 1/2} \,d\phi \,dz, \end{eqnarray}\vspace*{-10pt} \tag{ 80a }$$

$$\begin{eqnarray} B^{n}_{\phi ; \;i,j\pm 1/2,k} &\equiv & \frac{1}{\Delta R \,\Delta z} \int _{z_{k- 1/2}}^{z_{k+ 1/2}}\!\!\! \int _{R_{i- 1/2}}^{R_{i+ 1/2}} B_\phi (R,\phi _{j\pm 1/2},z,t^n) \,dR \,dz,\nonumber\\ \end{eqnarray}\vspace*{-10pt} \tag{ 80b }$$

$$\begin{eqnarray} &&B^{n}_{z; \;i,j,k\pm 1/2} \equiv \frac{1}{R_i \,\Delta R \,\Delta \phi }\nonumber\\ &&\quad\times\, \int _{\phi _{j- 1/2}}^{\phi _{j+ 1/2}}\!\!\! \int _{R_{i- 1/2}}^{R_{i+ 1/2}} B_z(R,\phi,z_{k\pm 1/2},t^n) \,R \,dR \,d\phi, \end{eqnarray} \tag{ 80c }$$

are the interface area-averaged components of the magnetic field normal to each surface (i.e., the magnetic flux per unit area) and

$$\begin{eqnarray} &&{\mathcal {E}}^{n+ 1/2}_{R; \;i,j\pm 1/2,k\pm 1/2} \equiv \frac{1}{\Delta t \,\Delta R}\nonumber\\ &&\quad\times\, \int _{t^n}^{t^{n+1}}\!\!\!\! \int _{R_{i- 1/2}}^{R_{i+ 1/2}} {\mathcal {E}}_R(R,\phi _{j\pm 1/2},z_{k\pm 1/2},t) \,dR \,dt, \end{eqnarray} \tag{ 81a }$$

$$\begin{eqnarray} &&{\mathcal {E}}^{n+ 1/2}_{\phi ; \;i\pm 1/2,j,k\pm 1/2} \equiv \frac{1}{\Delta t \,R_{i\pm 1/2} \,\Delta \phi }\nonumber\\ &&\quad\times\, \int _{t^n}^{t^{n+1}}\!\!\!\! \int _{\phi _{i- 1/2}}^{\phi _{i+ 1/2}} {\mathcal {E}}_\phi (R_{i\pm 1/2},\phi,z_{k\pm 1/2},t) \,R_{i\pm 1/2} \,d\phi \,dt,\nonumber\\ \end{eqnarray} \tag{ 81b }$$

$$\begin{eqnarray} &&{\mathcal {E}}^{n+ 1/2}_{z; \;i\pm 1/2,j\pm 1/2,k} {\equiv} \frac{1}{\Delta t \,\Delta z} \int _{t^n}^{t^{n+1}}\!\!\!\! \int _{z_{i- 1/2}}^{z_{i+ 1/2}} {\mathcal {E}}_z(R_{i\pm 1/2},\phi _{j\pm 1/2},z,t) \,dz \,dt\nonumber\\ \end{eqnarray} \tag{ 81c }$$

are the corner-centered EMFs averaged over the edges bounding each surface. The EMFs in Equations (81) are approximated to some desired order of accuracy and the surface-averaged field components are evolved using Equations (79).

The interface-centered, area-averaged magnetic field components in Equations (80) comprise the fundamental representation of the magnetic field in Athena. However, one often needs to refer to the cell-centered, volume-averaged magnetic field components as well. Therefore, we adopt the averages

$$\begin{equation} B^{R,n}_{i,j,k} \equiv \frac{1}{2R_i} \big(R_{i- 1/2} B^{n}_{R; \;i- 1/2,j,k} + R_{i+ 1/2} B^{n}_{R; \;i+ 1/2,j,k} \big),\quad \end{equation} \tag{ 82a }$$

$$\begin{equation} B^{\phi,n}_{i,j,k} \equiv \frac{1}{2} \big(B^{n}_{\phi ; \;i,j- 1/2,k} + B^{n}_{\phi ; \;i,j+ 1/2,k} \big), \end{equation} \tag{ 82b }$$

$$\begin{equation} B^{z,n}_{i,j,k} \equiv \frac{1}{2} \big(B^{n}_{z; \;i,j,k- 1/2} + B^{n}_{z; \;i,j,k+ 1/2} \big). \end{equation} \tag{ 82c }$$

Note the use of an R-weighted average of the B_R interface values in Equation (82a); it is straightforward to show that this is the appropriate second-order accurate average in cylindrical coordinates. Equations (82) imply consistency relations between the Godunov fluxes computed by the Riemann solver (the fluxes of the volume-averaged magnetic field components) and the corner-centered EMFs (the fluxes of the area-averaged magnetic field components). These relations define how they are computed from each other (see GS05, for details).

The primary modification to the CTU+CT algorithm described in GS05 for cylindrical coordinates concerns the calculation of the upwinded, corner-centered EMF component ${\mathcal {E}}_z$. As we shall demonstrate, we must combine spatial gradients of different curvature to form ${\mathcal {E}}_z$. However, to form ${\mathcal {E}}_R$ or ${\mathcal {E}}_\phi$, we combine spatial gradients of the same curvature, hence the effect of that curvature is canceled out and no subsequent modification is necessary.

For example, to compute ${\mathcal {E}}_{z; \;i-1/2,j-1/2}$, we estimate $(\partial _\phi {\mathcal {E}}_z)_{i-1/2,j-3/4}$ and use a centered-difference scheme to calculate one estimate

$$\begin{eqnarray} &&{\mathcal {E}}_{z; \;i-1/2,j-1/2} = {\mathcal {E}}_{z; \;i-1/2,j-1} + {\frac{R_{i-1/2} \Delta \phi }{2}} (\partial _\phi {\mathcal {E}}_z)_{i-1/2,j-3/4}.\nonumber\\ \end{eqnarray} \tag{ 83 }$$

In the same manner, we integrate ${\mathcal {E}}_z$ to the corner from each of the remaining adjacent interface centers and take the arithmetic average:

$$\begin{eqnarray} &&{\mathcal {E}}_{z; \;i-1/2,j-1/2} = {\frac{1}{4}}({\mathcal {E}}_{z; \;i-1/2,j-1} + {\mathcal {E}}_{z; \;i-1/2,j}\nonumber \\ &&\quad + {\mathcal {E}}_{z; \;i-1,j-1/2} + {\mathcal {E}}_{z; \;i,j-1/2}) \nonumber \\ &&\quad+ {\frac{R_{i-1/2} \,\Delta \phi }{8}} \big[ (\partial _\phi {\mathcal {E}}_z)_{i-1/2,j-3/4} - (\partial _\phi {\mathcal {E}}_z)_{i-1/2,j-1/4} \big] \nonumber \\ &&\quad+ {\frac{\Delta R}{8}} \big[ (\partial _R {\mathcal {E}}_z)_{i-3/4,j-1/2} - \left(\partial _R {\mathcal {E}}_z \right)_{i-1/4,j-1/2} \big]. \end{eqnarray} \tag{ 84 }$$

To ensure stability, for $(\partial _\phi {\mathcal {E}}_z)_{i-1/2,j-3/4}$, we use the upwinding scheme (GS05, Equation (50)) based on the sign of the mass flux at the center of each interface:

$$\begin{eqnarray} &&\dsty(\partial _\phi {\mathcal {E}}_z)_{i-1/2,j-3/4}\nonumber\\ &&\dsty\quad = \left\lbrace \begin{array}{@{}ll} (\partial _\phi {\mathcal {E}}_z)_{i-1,j-3/4}, & v_{R; \;i-1/2,j-1}>0 \\ \dsty(\partial _\phi {\mathcal {E}}_z)_{i,j-3/4}, & v_{R; \;i-1/2,j-1}<0 \\ \dsty{\frac{1}{2}}\big[ (\partial _\phi {\mathcal {E}}_z)_{i-1,j-3/4} + (\partial _\phi {\mathcal {E}}_z)_{i,j-3/4} \big], & {\rm otherwise} \end{array} \right..\nonumber\\ \end{eqnarray} \tag{ 85 }$$

The formulae for the remaining gradients are analogous.

To obtain the estimates of $\partial _\phi {\mathcal {E}}_z$ needed in Equation (85), we use a centered-difference scheme based on the cell-centered EMFs, computed using volume averages of ρ, $\rho \boldsymbol{v}$, and $\boldsymbol{B}$, and on the interface-centered EMFs, which come directly from the fluxes:

$$\begin{eqnarray} (\partial _\phi {\mathcal {E}}_z)_{i,j-3/4} &=& {\frac{2}{R_{i} \,\Delta \phi }} ({\mathcal {E}}_{z; \;i,j-1/2} - {\mathcal {E}}_{z; \;i,j-1}), \end{eqnarray} \tag{ 86a }$$

$$\begin{eqnarray} &&(\partial _\phi {\mathcal {E}}_z)_{i-1,j-3/4} = {\frac{2}{R_{i-1} \,\Delta \phi }} ({\mathcal {E}}_{z; \;i-1,j-1/2} - {\mathcal {E}}_{z; \;i-1,j-1}).\nonumber\\ \end{eqnarray} \tag{ 86b }$$

Note that this scheme has no dependence on the particular type of Riemann solver used to calculate the fluxes. Furthermore, note the different radial scale factors appearing in Equations (86) resulting from the combination of ϕ-gradients at different radii; the factors of Δϕ cancel the factor of Δϕ from Equation (84), but the radial scale factors themselves do not cancel and must be inserted into the algorithm. On the other hand, when ϕ-gradients are combined at the same radius, as is the case for the corner integration of ${\mathcal {E}}_{R}$, both the Δϕ and the corresponding radial scale factors cancel, hence no modification is required.

In this problem, we investigate various steady equilibria in order to evaluate the code's ability to balance forces. While there are many possible tests to choose from, we give here two representative examples demonstrating simple magnetohydrostatic equilibria.

First, we consider the axisymmetric magnetic field

$$\begin{equation} \boldsymbol{B} = \frac{B_0}{R} \hat{\phi }, \end{equation} \tag{ 87 }$$

for which the outward magnetic pressure and inward tension forces sum to zero. Note that the magnetic field given by Equation (87) satisfies the $\nabla \cdot \boldsymbol{B} = 0$ constraint. We use B₀ = 1 and set the velocity $\boldsymbol{v}=0$, mass density ρ = 1, and gas pressure P = 1. Figure 1 shows the convergence of the L₂ norm of the L₁ error vector (rms error) for the solution at time t = 10, defined as

$$\begin{equation} \delta \boldsymbol{q} = \frac{1}{N} \sum _i \big| \boldsymbol{q}_i - \boldsymbol{q}_i^0 \big|, \end{equation} \tag{ 88 }$$

where $\boldsymbol{q}_i^0$ is the initial solution. For reference, we plot a line of slope −2 (dashed) alongside the error to demonstrate that the convergence is second order in 1/N. These data were computed using the HLLD fluxes and third-order reconstruction; the results were similar for all combinations of Roe or HLLD fluxes, second- or third-order reconstruction, and one-dimensional, two-dimensional, or three-dimensional integrators. However, because the two-dimensional and three-dimensional algorithms differ significantly from the one-dimensional version, especially in their inclusion of transverse flux gradients and CT updates, we also present the results of the same test using these integrators on grids which are essentially one dimensional, but contain a few grid cells in each transverse direction considered. By symmetry, it is clear that any number of cells may be used in the transverse directions, but since the grid cell volumes change with R in multidimensions, the CFL condition will determine the timestep based on grid cell volume as well as the maximum signal speed, so the absolute errors should not be compared between, say, the one-dimensional and two-dimensional algorithms. Only the order of convergence of each individual algorithm is meaningful. Additionally, we have performed tests with the outward acceleration v²_ϕ/R from a solid-body rotation profile v_ϕ = Ω₀R balanced by the gradient of the static gravitational potential Φ = (Ω₀R)²/2, as well as with constant v_z ≠ 0, and find similar results in all cases.

Zoom In Zoom Out Reset image size

Second, we consider the non-axisymmetric magnetic field

$$\begin{equation} \boldsymbol{B} = \frac{B_0 \,\cos (\psi)}{R}\hat{R}, \end{equation} \tag{ 89 }$$

which when combined with the gas pressure

$$\begin{equation} P = P_0 + \frac{B_0^2 \,[1 + \sin ^2(\psi)]}{2 R^2} \end{equation} \tag{ 90 }$$

yields zero net force in the ϕ-direction. The combination of static gravitational potential

$$\begin{equation} \Phi _1 = - \frac{B_0^2}{2\rho _0 R^2} \end{equation} \tag{ 91 }$$

and mass density

$$\begin{equation} \rho = \rho _0 [1 + \sin ^2(\psi)] \end{equation} \tag{ 92 }$$

together balance the gradient of the gas pressure in the R-direction. Here, we use the angular coordinate ψ = 2π(ϕ − ϕ_min)/(ϕ_max − ϕ_min) so that B_R, P, and ρ are all periodic in the ϕ-domain. Note that the magnetic field given by Equation (89) satisfies the $\nabla \cdot \boldsymbol{B}=0$ constraint. We use B₀ = 1, P₀ = 1, ρ₀ = 1, and once again use the solid-body rotation profile v_ϕ = Ω₀R balanced by the gradient of the static gravitational potential Φ₂ = (Ω₀R)²/2, where Ω₀ = π/4. Note that the total potential is given by Φ = Φ₁ + Φ₂.

Figure 2 shows the convergence of the rms error for the solution at time t = 10. For reference, we once again plot a line of slope −2 (dashed) alongside the error to demonstrate that the convergence is second order. These data were computed using the HLLD fluxes, third-order reconstruction, and the two-dimensional integrator; the results were similar for all combinations of the Roe or HLLD fluxes, second- or third-order reconstruction, and two-dimensional or three-dimensional integrators. We use a computational domain of size N-by-N with R ∈ [1, 2], ϕ ∈ [0, π/4], and z = 0. We use periodic boundary conditions in the ϕ- and z-directions, and Dirichlet boundary conditions in the R-direction. Similar results were obtained using a Neumann boundary condition in the R-direction.

Zoom In Zoom Out Reset image size

In this problem, we investigate the stability of rotating disks evolved with our code, using the two-dimensional integrator. Given a differential rotation profile, Ω(R), Rayleigh's criterion for stability to axisymmetric, infinitesimal disturbances is that specific angular momentum increase outward:

$$\begin{equation} \partial _R [ (R^2 \Omega (R))^2 ] > 0. \end{equation} \tag{ 93 }$$

While it is possible that systems satisfying Rayleigh's criterion are still subject to growth of finite-amplitude non-axisymmetric disturbances, laboratory measurements at Reynolds number up to 2 × 10⁶ have found that Couette flows violating Rayleigh's criterion show large angular momentum transport associated with turbulence, while those satisfying Rayleigh's criterion do not (Ji et al. 2006).

For our test, we consider power-law rotational profiles of the form Ω(R) ∝ R^−q, where q is a constant, the so-called "shear parameter." Rayleigh's criterion in a differentially rotating system of this form predicts stability when q < 2, and instability for q>2. For example, Keplerian rotation with q = 1.5 is predicted to be stable (for unmagnetized flows).

Using a constant background density and pressure, we set

$$\begin{equation} v_\phi (R) = R \,\Omega (R) = \Omega _0 R^{1-q} \end{equation} \tag{ 94 }$$

and set the gravitational potential so that rotational equilibrium is achieved. Next, we perturb v_ϕ to

$$\begin{equation} \tilde{v}_\phi = v_\phi + \delta v_\phi, \end{equation} \tag{ 95 }$$

where δv_ϕ is a random variable uniformly distributed in [ − , ], and is small, typically on the order of 10⁻⁴. The same initial perturbation is used for each value of q considered. We use a grid of 200 × 400 cells over the domain [3, 7] × [0, π/2], which is chosen so that $R_{\rm{avg}} \,\Delta \phi \sim 2 \Delta R$. We set ρ₀ = 200, P₀ = 1, and use an adiabatic index of γ = 5/3 which gives c_s ≈ 0.09. With Ω₀ = 2π, $v_{\phi,\rm{min}} \approx 0.9$, which puts the Mach numbers in a range of approximately 10–20 over the domain. Thus, the flow is rotationally dominated.

As a diagnostic of instability, we compute a scaled mean perturbed angular momentum flux,

$$\begin{equation} \frac{\langle R \,\rho v_R \,\delta v_\phi \rangle } {\langle R P \rangle } = \frac{\int\!\!\int R \,\rho v_R (v_\phi -R\Omega) RdR\,d\phi }{\int\!\! \int RP RdR\,d\phi }. \end{equation} \tag{ 96 }$$

For stable flows, this will remain on the order of the initial perturbation, but for unstable flows, it will diverge exponentially. Figure 3 shows the values of the dimensionless angular momentum flux as a function of time for t ∈ [0, 300] for various values of the shear parameter near the marginal stability limit of q = 2. Consistent with Rayleigh's criterion, the flows with q < 2 remain stable, and those with q>2 go unstable. For the q = 2.05 case, the instability reaches saturation more quickly (around t = 90) and the mass flies off the grid, but the characteristic exponential growth is observed before this point. This test demonstrates the code's accurate conservation of angular momentum near the boundary of rotational stability.

Zoom In Zoom Out Reset image size

Additionally, we investigate the long-term stability of the rotation profiles given by the shear parameters q = 1, typical of galactic disk systems, and q = 1.5, typical of Keplerian systems. For this test we use the unperturbed equilibrium solutions as initial data. As a diagnostic of the error, we compute the cumulative mean of the dimensionless background angular momentum flux (proportional to the radial accretion rate):

$$\begin{equation} \frac{\langle \langle R \,\rho v_R \,(\Omega R) \rangle \rangle }{\langle \langle R P \rangle \rangle } = \frac{\int\!\! \int\!\! \int R \,\rho v_R (\Omega R) \,R\,dR\,d\phi \,dt}{\int\!\! \int\!\! \int RP \,R\,dR\,d\phi \,dt}. \end{equation} \tag{ 97 }$$

Figure 4 shows the dimensionless cumulative mean as a function of time for t ∈ [0, 300]. We use the same computational domain and background state as above. Since v_ϕ = ΩR ∝ R^1−q is constant for q = 1, and constant profiles are reconstructed exactly in our code, we observe relatively small errors for this level of discretization, indicating that angular momentum is conserved very well for systems of astrophysical interest.

Zoom In Zoom Out Reset image size

In this problem, we investigate a strong two-dimensional shock using the HLLC solver. We use the parameter set of GS08 and compare the outputs of our cylindrical code with those from the Cartesian version of Athena. For the Cartesian version we use the domain (x, y) ∈ [ − 0.5, 0.5] × [ − 0.75, 0.75], and for the cylindrical version we use the domain (R, ϕ) ∈ [1, 2] × [ − 0.5, 0.5] so that the physical domain spans an arc length of $R_{\rm{mid}}(\phi _{\rm{max}}-\phi _{\rm{min}})=1.5$ at $R_{\rm{mid}}=1.5$, giving a roughly similar domain sizes. The initial conditions consist of a circular region of hot gas with radius R = 0.1 and pressure P = 10 in an ambient medium of uniform pressure P₀ = 0.1 and density ρ₀ = 1. We use a computational grid of 200 × 300 cells, third-order reconstruction, and the HLLC fluxes with upwind-only integration for the L/R states in each version of the test. Contour plots of the density, pressure and specific kinetic energy densities of the evolved state at time t = 0.2 are shown in Figure 5 using the cylindrical (left column) or Cartesian (right column) version of Athena. For additional comparison, one-dimensional plots of these variables along a horizontal line through the center of the blast are shown in Figure 6, demonstrating excellent agreement between the Cartesian and cylindrical versions of Athena. Note in the Cartesian version that symmetry is perfectly preserved by the integrator, which is most easily seen in the grid noise in the interior of the shell in Figure 5. Symmetry is also preserved rather well in the cylindrical version although in this case the grid is non-uniform.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this problem, we investigate a steady, axisymmetric, rotating hydrodynamic wind as a test of angular momentum transport across a sonic transition. We adopt a Newtonian gravitational potential of the form Φ_g = −GM/R. The constants of motion are given by

$$\begin{equation} K = P \rho ^{-\gamma }, \end{equation} \tag{ 98 }$$

$$\begin{equation} \dot{M} = R \rho v_R, \end{equation} \tag{ 99 }$$

$$\begin{equation} J = R v_\phi. \end{equation} \tag{ 100 }$$

This flow must satisfy the Bernoulli equation, $\mathcal {B} = {\rm constant}$, along streamlines for

$$\begin{equation} \mathcal {B} \equiv {\frac{1}{2}}v^2 + h + \Phi _g, \end{equation} \tag{ 101 }$$

where

$$\begin{equation} h \equiv \int \frac{dP}{\rho } = \frac{c_s^2}{\gamma -1} \end{equation} \tag{ 102 }$$

is the specific enthalpy of the gas.

We scale density and pressure to their values at infinity, ρ_∞ ≡ ρ(∞) and P_∞ ≡ P(∞), and the radial coordinate to some finite fiducial value, R_B. In terms of α ≡ ρ/ρ_∞ and χ ≡ R/R_B, the constant entropy parameter is K = c²_∞/(γ ρ^γ−1_∞), where c²_∞ ≡ γP_∞/ρ_∞ is the square of the sound speed at infinity. At any radius, the local sound speed, radial velocity, and specific enthalpy satisfy

$$\begin{equation} c_s^2 = c_\infty ^2 \alpha ^{\gamma -1}, \end{equation} \tag{ 103 }$$

$$\begin{equation} v_R^2 = \mathcal {M}_R^2 c_\infty ^2 \alpha ^{\gamma -1}, \end{equation} \tag{ 104 }$$

$$\begin{equation} h = \frac{1}{\gamma -1} c_\infty ^2 \alpha ^{\gamma -1}, \end{equation} \tag{ 105 }$$

where $\mathcal {M}_R \equiv v_R / c_s$ is the radial Mach number.

We define the dimensionless radial mass flux by

$$\begin{equation} \lambda \equiv \frac{\dot{M}}{R_B \rho _\infty c_\infty } = \chi \mathcal {M}_R \alpha ^{(\gamma +1)/2} \end{equation} \tag{ 106 }$$

and the dimensionless angular momentum by

$$\begin{equation} \omega \equiv \frac{J}{R_B c_\infty } = \frac{\chi v_\phi }{c_\infty }. \end{equation} \tag{ 107 }$$

Solving Equation (106) for α and introducing β ≡ 2(γ − 1)/(γ + 1), we find that $\alpha ^{\gamma -1} = [\lambda /(\chi \mathcal {M}_R)]^\beta$. Finally, taking R_B to be the Bondi radius, R_B ≡ GM/c²_∞, we obtain the dimensionless Bernoulli equation

$$\begin{eqnarray} &&\left(\frac{1}{2} \mathcal {M}_R^{2-\beta } + \frac{1}{\gamma -1} \mathcal {M}_R^{-\beta } \right) \lambda ^{\beta }\nonumber\\ &&\quad = \left[ \frac{\tilde{\mathcal {B}}}{\gamma -1} \chi ^{\beta } + \chi ^{\beta -1} - \frac{\omega ^2}{2} \chi ^{\beta -2} \right], \end{eqnarray} \tag{ 108 }$$

where $\tilde{\mathcal {B}} \equiv \mathcal {B}/h_\infty$ is the dimensionless Bernoulli constant, and h_∞ ≡ c²_∞/(γ − 1) is the specific enthalpy at infinity.

Figure 7 shows the contours of λ for various values of ω, using an adiabatic index of γ = 5/3. For the ω = 0 case, we recover a cylindrical version of Parker's spherically symmetric wind (see, e.g., Spitzer 1978). The bold lines represent the transonic solutions (wind and accretion) passing through the X-type saddle points. These solutions can be found by first writing Equation (108) as $\mathcal {F}(\mathcal {M}_R) f(\lambda) = \mathcal {G}(\chi)$ and requiring that $\mathcal {F}^{\prime }=\mathcal {G}^{\prime }=0$. The first constraint implies that $\mathcal {M}_R = 1$, i.e., the saddle point is the sonic point. The second constraint yields a quadratic in χ, and for ω ∈ (0, ω_max), there are two distinct, positive solutions,

$$\begin{equation} \chi _\pm = \frac{3-\gamma \pm \sqrt{(3-\gamma)^2-16\tilde{\mathcal {B}}\omega ^2}}{4\tilde{\mathcal {B}}}, \end{equation} \tag{ 109 }$$

where $\omega _{\rm max} \equiv (3-\gamma)/(4\tilde{\mathcal {B}}^{1/2})$. It can be shown that χ₋ gives an O-type critical point and χ₊ gives the desired transonic critical point. For ω>ω_max, no transonic solutions exist. Note further that no transonic solutions exist for χ < χ_min, where χ_min represents the point at which the transonic wind solution for which $\mathcal {M}_R \rightarrow \infty$ as χ → ∞ joins the transonic accretion solution for which $\mathcal {M}_R \rightarrow 0$ as χ → ∞ (see, e.g., the lower-left panel of Figure 7). We define χ_min to be the smallest value of χ ⩾ 0 for which $\mathcal {G}(\chi) \ge 0$, which is given by

$$\begin{equation} \chi _{\rm min} = \frac{-(\gamma -1)+\sqrt{(\gamma -1)^2+2\tilde{\mathcal {B}}\omega ^2}}{2\tilde{\mathcal {B}}}. \end{equation} \tag{ 110 }$$

Finally, the critical value of the radial mass flux, λ_c, is defined by Equation (108) with χ = χ₊ and $\mathcal {M}_R = 1$, which is given by

$$\begin{equation} \lambda _{\rm c} = \big[ \chi _+^{\beta -2} (\chi _+ - \omega ^2) \big]^{1/\beta }. \end{equation} \tag{ 111 }$$

Zoom In Zoom Out Reset image size

For our code test, we use γ = 5/3 (i.e. β = 1/2), ω = 0.3, and $\tilde{\mathcal {B}}=1$, which give the transonic solution shown in Figure 7(c). The critical point occurs at χ₊ ≈ 0.479, with λ_c ≈ 1.377. We solve the problem on the domain χ ∈ [χ₋, 2], where χ₋ ≈ 0.188, using bisection with a tolerance of = 10⁻¹⁰ to evaluate $\mathcal {M}_R$ at each χ from Equation (108). Once $\mathcal {M}_R$ is known, v_R, ρ, and P follow algebraically. We choose units such that GM = c_∞ = 1, which yields R_B = 1. We fix the solution at the inner and outer boundaries and evolve the initial solution long enough for it to settle into equilibrium.

Figure 8 shows the convergence of the L₂ norm of the L₁ error vector for the solution at time t = 5.0. These data were computed using the Roe fluxes, second-order reconstruction, and the one dimension; the results were similar for all combinations of Roe or HLLC fluxes and second- or third-order reconstruction. The test was also performed using the two-dimensional and three-dimensional integrator with a few grid cells in the transverse directions.

Zoom In Zoom Out Reset image size

This test clearly demonstrates the code's ability to maintain smooth, steady hydrodynamic flows in both subsonic and supersonic regimes, as well as its ability to conserve angular momentum to second order in cylindrical geometry.

In this problem, we investigate the advection of a weak field loop in two-dimensional and three-dimensional cylindrical coordinates, analogous to the Cartesian test appearing in GS05 and GS08. The main difference with our test is that we advect the field loop in the ϕ-direction only as opposed to a more general advection oblique to the grid. We use the computational domain (R, ϕ) ∈ [1, 2] × [ − 2/3, 2/3], which has the same total area as the Cartesian version of the test. We use periodic boundary conditions in ϕ and fixed boundary conditions in R. We use uniform initial density ρ₀ = 1 and pressure P₀ = 1 with a solid-body rotation profile of v_ϕ = Ω₀R, where we set Ω₀ = 4/3 so that the field loop is advected once across the grid by t = 1. The initial z-component of the magnetic field is 0, and the R- and ϕ-components are set using the z-component of the magnetic vector potential

$$\begin{equation} A_z \equiv \left\lbrace \begin{array}{@{}ll} A_0 (a_0-r) & {\rm for }\; r \le a_0, \\ 0 & {\rm for }\; r > a_0 \end{array} \right., \end{equation} \tag{ 112 }$$

where a₀ is the radius of the field loop and $r \equiv \sqrt{R^2 + R_0^2 - 2 R R_0 \cos (\phi -\phi _0)}$ is the distance from the center of the loop (R₀, ϕ₀). We use A₀ = 10⁻³ and a₀ = 0.3 so that inside the field loop P/B = 10⁶ and the field loop should be advected passively.

Figure 9 shows the magnetic energy density B²/2 and magnetic field lines at times t = 0 and t = 2 for the two-dimensional problem. The field lines are the contours of A_z which can be readily computed since the field is planar and the CTU+CT algorithm preserves the $\nabla \cdot \boldsymbol{B} = 0$ condition. Note that the circular shape of the field lines is nicely preserved.

Zoom In Zoom Out Reset image size

Figure 10 shows the time evolution of the volume-averaged magnetic energy density. The dissipation is well-described by a power law of the form 〈B²/2〉 = A(1 − (t/τ)^α) (for t ≪ τ), with A = 7.02 × 10⁻⁸, τ = 1.46 × 10⁴, and α = 0.342, and with a residual error of 0.0580. Note that the overall dissipation in this problem is less than that of the Cartesian version since the advection is only in one direction; GS05 found τ = 1.06 × 10⁴ and α = 0.291 for a similar fit.

Zoom In Zoom Out Reset image size

In this problem, we investigate two-dimensional and three-dimensional MHD shocks in a strongly magnetized medium with low plasma-β, denoted β_p ≡ 2P/B. We run two problems, one with B₀ = 1 and β_p = 0.2 using the parameter set of GS08, and another with B₀ = 10 and β_p = 0.02 using the parameter set of Londrillo & Del Zanna (2000). In each case, we compare the outputs of the cylindrical and Cartesian versions of Athena.

The moderate B-field, β_p = 0.2 case uses HLLD fluxes and the same setup as the hydrodynamic blast described in Section 10.3, but with a uniform background magnetic field of strength B₀ = 1 oriented at a 45° angle to the positive $\hat{x}$- or $\hat{R}$-axis. For the two-dimensional case, contour plots of the density, pressure, specific kinetic energy, and magnetic energy are shown in Figure 11 based on the cylindrical (left column) or Cartesian (right column) versions of Athena. Figure 12 shows a plot of these variables along a horizontal line through the center of the blast. Note that in the Cartesian version, the background field is uniformly inclined to the grid, but in the cylindrical version, the angle the background field makes with the grid changes as a function of ϕ; nonetheless, a high degree of symmetry is observed in the solution. We have also tested a three-dimensional analog of this problem on a 200³ grid with z ∈ [ − 0.5, 0.5]. Again, plots of selected variables along a horizontal line through the center of the blast (in the z = 0 plane) in Figure 13 show agreement between the Cartesian and cylindrical versions of Athena.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The strong B-field, β_p = 0.02 case uses the domain (x, y) ∈ [ − 0.5, 0.5] × [ − 0.5, 0.5] for the Cartesian version, and for the cylindrical version, uses the domain (R, ϕ) ∈ [1, 2] × [2/3, 2/3], giving a roughly similar domain size in each case. The initial conditions consist of a circular region of hot gas with radius R₀ = 0.125 and pressure P = 100 in an ambient medium of uniform pressure P₀ = 1 and density ρ₀ = 1. There is uniform background magnetic field of strength B₀ = 10 oriented parallel to the x-axis. We use a computational grid of 200² cells, third-order reconstruction and the HLLD fluxes with upwind-only integration for the L/R states. The density, pressure, specific kinetic energy, and magnetic energy along a horizontal line through the center of the blast at time t = 0.02 are shown in Figure 14, with comparison to the Cartesian version of Athena. We have also conducted a three-dimensional analog of this test on a 200³ grid with z ∈ [ − 0.5, 0.5]. Figure 15 shows a comparison of contour plots with the Cartesian version of Athena, and Figure 16 shows a comparison along a horizontal line through the center of the blast (in the z = 0 plane).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this problem, we investigate a cylindrical version of the Weber–Davis wind solution as described in Sakurai (1985). We assume a steady, axisymmetric, two-dimensional MHD flow with planar magnetic field, and a gravitational potential Φ_g = −GM/R. The constants of motion are K = Pρ^−γ, $\dot{M} = R\rho v_R$, f = RB_R, β = B_R/(ρv_R), as well as

$$\begin{eqnarray} \Omega &=& {\frac{1}{R}}\left(v_\phi - \frac{B_\phi }{\beta \rho } \right), \end{eqnarray} \tag{ 113a }$$

$$\begin{eqnarray} J &=& R (v_\phi - \beta B_\phi) = R^2 \Omega + R u_\phi (1 - \beta ^2 \rho), \end{eqnarray} \tag{ 113b }$$

$$\begin{eqnarray} \mathcal {B} &=& {\frac{1}{2}}\big(u_R^2 + u_\phi ^2\big) + h + \Phi _g - {\frac{1}{2}}(\Omega R)^2. \end{eqnarray} \tag{ 113c }$$

Here, $\boldsymbol{u} \equiv (v_R, \,v_\phi -R\Omega,\,0)$ is the velocity in a frame rotating at angular velocity Ω, and in this frame $\boldsymbol{B} = \beta \rho \boldsymbol{u}$, so that the fluid feels no force from the magnetic field. Note that the Bernoulli parameter in the rotating frame includes a centrifugal potential contribution. The Alfvén Mach number in the rotating frame is given by $\mathcal {M}_A \equiv u/c_A = 1/\sqrt{\beta ^2 \rho }$. Let R_A and ρ_A denote the radius and density, respectively, at the Alfvén Mach point, i.e., where $\mathcal {M}_A = 1$. Then $\beta = 1/\sqrt{\rho _A}$ and J = R²_AΩ.

Letting x ≡ R/R_A and $y \equiv \rho /\rho _A = \rho \beta ^2 = \mathcal {M}_A^{-2}$ denote the scaled radius and density, respectively, the Bernoulli parameter is

$$\begin{eqnarray} \mathcal {B} &=& \frac{GM}{R_A} \left[ \frac{\eta }{2 x^2 y^2} + \frac{\omega }{2} \left(\left(\frac{1/x - x}{1-y} \right)^2 - x^2 \right)\right.\nonumber\\ &&\left. +\, \frac{\theta }{\gamma - 1} y^{\gamma -1} - \frac{1}{x} \right], \end{eqnarray} \tag{ 114 }$$

where we have defined the scaled parameters

$$\begin{eqnarray} \eta &\equiv & \frac{\dot{M}^2}{R_A \rho _A^2 GM}, \end{eqnarray} \tag{ 115a }$$

$$\begin{eqnarray} \theta &\equiv & \frac{\gamma K \rho _A^{\gamma -1} R_A}{GM}, \end{eqnarray} \tag{ 115b }$$

$$\begin{equation} \omega \equiv \frac{R_A^3 \Omega ^2}{GM}. \end{equation} \tag{ 115c }$$

Letting $\tilde{\mathcal {B}} \equiv \mathcal {B}/(GM/R_A)$, we have

$$\begin{equation} x \frac{\partial \tilde{\mathcal {B}}}{\partial x} = -\frac{\eta }{x^2 y^2} + \omega \left(\frac{(x^2 - 1/x^2)}{(1-y)^2} - x^2 \right) + \frac{1}{x}, \end{equation} \tag{ 116a }$$

$$\begin{equation} y \frac{\partial \tilde{\mathcal {B}}}{\partial y} = -\frac{\eta }{x^2 y^2} + \omega y \frac{(1/x - x)^2}{(1-y)^3} + \theta y^{\gamma -1}. \end{equation} \tag{ 116b }$$

To find wind solutions, we solve the Bernoulli Equation (114) under the constraint that the solution be locally flat at the slow- and fast-magnetosonic points, i.e., $\partial \tilde{\mathcal {B}}/\partial x = \partial \tilde{\mathcal {B}}/\partial y = 0$ at (x_s,  y_s) and (x_f,  y_f). We will specify θ and ω, and let η, $\tilde{\mathcal {B}}$, x_s, y_s, x_f, and y_f vary. Note that this becomes a system of six equations in six unknowns, so if a solution exists, it must be unique. Following Sakurai (1985), we use the parameters γ = 1.2, θ = 1.5 and ω = 0.3 as for the spherically symmetric solar wind model of Weber and Davis. The numerical solution is given in Table 1.

Table 1. Weber–Davis Wind Parameters

Parameter	Value
γ	1.2
θ	1.5
ω	0.3
η	2.3609
$\tilde{\mathcal {B}}$	7.8745
xs	0.5243
ys	2.4986
xf	1.6383
yf	0.5374

Download table as:  ASCIITypeset image

We solve the problem on the domain x ∈ [0.4, 1.8], so that the wind solution passes through all three critical points. We choose units such that GM = 1 and fix the initial solution at the inner and outer boundaries, evolving it long enough for equilibrium to develop. Figure 17 shows the convergence of the rms error in the L₁-norm of the solution at t = 5.0 compared to the initial solution. These data were computed using the Roe fluxes, second-order reconstruction, and the one-dimensional integrator; the results were similar for all combinations of Roe/HLLD fluxes and second-/third-order reconstruction. Because the two-dimensional and three-dimensional algorithms differ significantly from the one-dimensional version, especially in their treatment of magnetic fields, we present the results of the same test using these integrators on grids which are essentially one-dimensional, but contain a few grid cells in each transverse direction considered.

Zoom In Zoom Out Reset image size

Evidently, the algorithm yields second-order convergence for smooth MHD flows, and is able to maintain a steady magnetized solution, also conserving angular momentum as it is exchanged between the fluid and magnetic field.