THE PLUTO CODE FOR ADAPTIVE MESH COMPUTATIONS IN ASTROPHYSICAL FLUID DYNAMICS

We consider a Newtonian fluid with density ρ, velocity v = (v_x, v_y, v_z), and magnetic induction B = (B_x, B_y, B_z) and write the single-fluid MHD equations as

$$\begin{eqnarray} \frac{\partial \rho }{\partial t} + \nabla \cdot \left(\rho {\bf v}\right) &=& 0 \,, \nonumber\\ \displaystyle \frac{\partial {\bf (\rho {\bf v})}}{\partial t} + \nabla \cdot [ \rho {\bf v}{\bf v} - {\bf B}{\bf B}] + \nabla p_t&=& \nabla \cdot \mathsf {$\tau$ } + \rho {\bf g}, \nonumber\\&&[-3mm] \\ \displaystyle \frac{\partial {\cal E}}{\partial t} + \nabla \cdot [ ({\cal E} + p_t){\bf v} - ({\bf v}\cdot {\bf B}){\bf B}]& =& \nabla \cdot \mathsfpi_{\cal E} - \Lambda + \rho {\bf v}\cdot {\bf g}\,, \nonumber\\ \displaystyle \frac{\partial {\bf B}}{\partial t} - \nabla \times \left({\bf v}\times {\bf B}\right) &=& -\nabla \times \left(\mathsf {$\eta$ }{\bf J}\right)\,, \nonumber \end{eqnarray} \tag{ 1 }$$

where p_t = p + B²/2 is the total (thermal+magnetic) pressure, ${\cal E}$ is the total energy density, g is the gravitational acceleration term, and Λ accounts for optically thin radiative losses or heating. Divergence terms on the right-hand side account for dissipative physical processes and are described in detail in Section 2.1.1. Proper closure is provided by choosing an equation of state (EoS) which, for an ideal gas, allows us to write the total energy density as

$$\begin{equation} {\cal E} = \frac{p}{\Gamma - 1} + \frac{1}{2}\rho {\bf v}^2 + \frac{1}{2}{\bf B}^2{,} \end{equation} \tag{ 2 }$$

with Γ being the specific heat ratio. Alternatively, by adopting a barotropic or an isothermal EoS, the energy equation can be discarded and one simply has, respectively, p = p(ρ) or p = c²_sρ (where c_s is the constant speed of sound).

Chemical species and passive scalars are advected with the fluid and are described in terms of their number fraction X_α, where α = 1, ..., N_ions label the particular ion. They obey non-homogeneous transport equations of the form

$$\begin{equation} \frac{\partial (\rho X_\alpha)}{\partial t} + \nabla \cdot \left(\rho X_\alpha {\bf v}\right) = \rho S_\alpha , \end{equation} \tag{ 3 }$$

where the source term S_α describes the coupling between different chemical elements inside the reaction network (see, for instance, Teşileanu et al. 2008).

2.1.1. Non-ideal Effects

Non-ideal effects due to dissipative processes are described by the differential operators included on the right-hand side of Equation (1). Viscous stresses may be included through the viscosity tensor $\mathsf {$\tau$ }$ defined by

$$\begin{equation} \mathsf {$\tau$ } = \rho \nu \left[ \nabla {\bf v} + (\nabla {\bf v})^{\rm T} - \frac{2}{3} \mathsf{I}\nabla \cdot {\bf v}\right]\,, \end{equation} \tag{ 4 }$$

where ν is the kinematic viscosity and $\mathsf{I}$ is the identity matrix. Similarly, magnetic resistivity is accounted for by prescribing the resistive η tensor (diagonal). Dissipative terms contribute to the net energy balance through the additional flux $\mathsfpi_{\cal{E}}$ appearing on the right-hand side of Equation (1):

$$\begin{equation} \mathsfpi_{\cal E} = {\bf F}_c + {\bf v}\cdot \mathsf {$\tau$ } - \mathsf {$\eta$ }\cdot {\bf J}\times {\bf B}\,, \end{equation} \tag{ 5 }$$

where the different terms give the energy flux contributions due to, respectively, thermal conductivity, viscous stresses, and magnetic resistivity.

The thermal conduction flux F_c smoothly varies between classical and saturated regimes and reads

$$\begin{equation} {\bf F}_c = \frac{q}{|{\bf F}_{\rm class}| + q} {\bf F}_{\rm class}, \end{equation} \tag{ 6 }$$

where q = 5ϕρc³_iso is the magnitude of the saturated flux (Cowie & McKee 1977), ϕ is a parameter of order unity accounting for uncertainties in the estimate of q, c_iso is the isothermal speed of sound, and

$$\begin{equation} {\bf F}_{\rm class} = \kappa _\parallel \hat{{\bf b}}(\hat{{\bf b}}\cdot \nabla T) + \kappa _\perp [\nabla T - \hat{{\bf b}}(\hat{{\bf b}}\cdot \nabla T)] \end{equation} \tag{ 7 }$$

is the classical heat flux with conductivity coefficients κ_∥ and κ_⊥ along and across the magnetic field lines, respectively (Orlando et al. 2008). Indeed, the presence of a partially ordered magnetic field introduces a large anisotropic behavior by channeling the heat flux along the field lines while suppressing it in the transverse direction (here $\hat{{\bf b}} = {\bf B}/|{\bf B}|$ is a unit vector along the field line). We point out that, in the classical limit q → ∞, thermal conduction is described by a purely parabolic operator and flux discretization follows standard FD. In the saturated limit (|∇T| → ∞), on the other hand, the equation becomes hyperbolic and thus an upwind discretization of the flux is more appropriate (Balsara et al. 2008). This is discussed in more detail in Appendix A.

A (special) relativistic extension of the previous equations requires the solution of energy momentum and number density conservation. Written in divergence form we have

$$\begin{eqnarray} \frac{\partial (\rho \gamma)}{\partial t} + \nabla \cdot \left(\rho \gamma {\bf v}\right) &=& 0 \,, \nonumber\\ \frac{\partial {\bf m}}{\partial t} + \nabla \cdot [ w\gamma ^2{\bf v}{\bf v} - {\bf B}{\bf B} - {\bf E}{\bf E}] + \nabla p_t &=&0 \,, \nonumber\\&&[-2mm] \\[-2mm] \frac{\partial {\bf B}}{\partial t} - \nabla \times \left({\bf v}\times {\bf B}\right) &=& 0 \,, \nonumber\\ \displaystyle \frac{\partial \cal E}{\partial t} + \nabla \cdot \left({\bf m}-\rho \gamma {\bf v}\right) &=& 0 \,,\nonumber \end{eqnarray} \tag{ 8 }$$

where ρ is the rest-mass density, γ is the Lorentz factor, velocities are given in units of the speed of light (c = 1), and the fluid momentum m accounts for matter and electromagnetic terms: m = wγ²v + E × B, where E = −v × B is the electric field and w is the gas enthalpy. The total pressure and energy include thermal and magnetic contributions and can be written as

$$\begin{equation} p_t = p + \frac{{\bf B}^2 + {\bf E}^2}{2} \,,\quad {\cal E} = w\gamma ^2 - p +\frac{{\bf B}^2 + {\bf E}^2}{2} - \rho \gamma \,. \end{equation} \tag{ 9 }$$

Finally, the gas enthalpy w is related to ρ and p via an EoS, which can be either the ideal gas law,

$$\begin{equation} w = \rho + \frac{\Gamma p}{\Gamma -1} \,, \end{equation} \tag{ 10 }$$

or the Taub-Mathews (TM, Mathews 1971) EoS

$$\begin{equation} w = \frac{5}{2}p + \sqrt{\frac{9}{4}p^2 + \rho ^2} \,, \end{equation} \tag{ 11 }$$

which provides an analytic approximation of the Synge relativistic perfect gas (Mignone & McKinney 2007).

A relativistic formulation of the dissipative terms will not be presented here and will be discussed elsewhere.

In the following, we adopt an orthonormal system of coordinates specified by the unit vectors $\hat{{\bf e}}_d$ (d is used to label the direction, e.g., d = {x, y, z} in Cartesian coordinates) and conveniently assume that conserved variables ${\bf U}\,{=}\,(\rho ,\rho {\bf v},{\cal E}, {\bf B}, \rho X_\alpha)$—for the MHD equations—and ${\bf U}\,{=}\,(\rho \gamma ,{\bf m},{\cal E}, {\bf B})$—for RMHD—satisfy the following hyperbolic/parabolic partial differential equations

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

where $\mathsf{F}$ and $\mathsfpi$ are, respectively, the hyperbolic and parabolic flux tensors. The source term S_p is a point-local source term which accounts for body forces (such as gravity), cooling, chemical reactions, and the source term for the scalar multiplier (see Equation (14) below). We note that equations containing curl or gradient operators can always be cast in this form by suitable vector identities. For instance, the projection of ∇ × E in the coordinate direction given by the unit vector $\hat{{\bf e}}_d$ can be rewritten as

$$\begin{equation} (\nabla \times {\bf E})\cdot \hat{{\bf e}}_d \equiv \nabla \cdot ({\bf E}\times \hat{{\bf e}}_d) + {\bf E}\cdot (\nabla \times \hat{{\bf e}}_d)\,, \end{equation} \tag{ 13 }$$

where the second term on the right-hand side should be included as an additional source term in Equation (12) whenever different from zero (e.g., in cylindrical geometry). Similarly, one can rewrite the gradient operator as $\nabla p = \nabla \cdot (\mathsf {I}p)$.

Several algorithms employed in PLUTO are best implemented in terms of primitive variables, V  =  (ρ, v, B, p). In the following, we shall assume a one-to-one mapping between the two sets of variables, provided by appropriate conversion functions, that is, V  =  V(U) and U = U(V).

The system of conservation laws is advanced in time using the CTU (Colella 1990) method recently described by MT. Here, we outline the algorithm in a more concise manner and extend its applicability in the presence of parabolic (diffusion) terms and higher order reconstruction. Although the algorithms illustrated here are explicit in time, we will also consider the description of more sophisticated and effective approaches for the treatment of parabolic operators in a forthcoming paper.

We shall assume hereafter an equally-spaced grid with computational cells centered in (x_i, y_j, z_k) having size Δx  ×  Δy  ×  Δz. For the sake of exposition, we omit the integer subscripts i, j, or k when referring to the cell center and only keep the half-increment notation in denoting face values, e.g., F_{d, ±} ≡ F_{y, i, j + 1/2, k} when d = y. Following Mignone et al. (2007), we use $\Delta {\cal V}$ and A_{d, ±} to denote, respectively, the cell volume and areas of the lower and upper interfaces orthogonal to $\hat{{\bf e}}_d$.

An explicit second-order accurate discretization of Equation (12), based on a time-centered flux computation, reads

$$\begin{equation} \bar{{\bf U}}^{n+1} = \bar{{\bf U}}^n + \Delta t^n\sum _d\left({\bf{\cal L}}_{H,d}^{n+{1}/{2}} + {\bf{\cal L}}_{P,d}^{n+{1}/{2}} \right)\,, \end{equation} \tag{ 15 }$$

where $\bar{{\bf U}}$ is the volume-averaged array of conserved values inside the cell (i, j, k), Δtⁿ is the explicit time step, whereas ${\bf{\cal L}}_{H,d}$ and ${\bf{\cal L}}_{P,d}$ are the increment operators corresponding to the hyperbolic and parabolic flux terms, respectively,

$$\begin{eqnarray} {\bf{\cal L}}_{H,d}^{n+{1}/{2}} &=& -\frac{ A_{d,+}{\bf F}^{n+{1}/{2}}_{d,+} - A_{d,-}{\bf F}^{n+{1}/{2}}_{d,-}}{\Delta {\cal V}_d} + \hat{{\bf S}}_d^{n+{1}/{2}} \,,\quad \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray} {\bf{\cal L}}_{P,d}^{n+{1}/{2}} &=& \frac{ A_{d,+}\boldsymbol{\Pi }^{n+{1}/{2}}_{d,+} - A_{d,-}\boldsymbol{\Pi }^{n+{1}/{2}}_{d,-}}{\Delta {\cal V}_d} \,. \end{eqnarray} \tag{ 17 }$$

In the previous expression, F_{d, ±} and $\boldsymbol{\Pi }_{d,\pm }$ are, respectively, right (+) and left (−) face- and time-centered approximations to the hyperbolic and parabolic flux components in the direction of $\hat{{\bf e}}_d$. The source term $\hat{{\bf S}}_d$ represents the directional contribution to the total source vector $\sum \hat{{\bf S}}_d\equiv \hat{{\bf S}}_{\rm body}+\hat{{\bf S}}_{\rm geo}$ including body forces and geometrical terms implicitly arising when differentiating the tensor flux on a curvilinear grid. Cooling, chemical reaction terms, and the source term in Equation (14) (namely, (c²_h/c_p²)ψ) are treated separately in an operator-split fashion.

The computation of F_{d, ±} requires solving, at cell interfaces, a Riemann problem between time-centered adjacent discontinuous states, i.e.,

$$\begin{eqnarray} {\bf F}^{n+{1}/{2}}_{x,i+{1}/{2}} &=& {\cal R}\left({\bf U}^{n+{1}/{2}}_{i,+}, {\bf U}^{n+{1}/{2}}_{i+1,-}\right) \,, \nonumber\\ {\bf F}^{n+{1}/{2}}_{y,j+{1}/{2}} &=& {\cal R}\left({\bf U}^{n+{1}/{2}}_{j,+}, {\bf U}^{n+{1}/{2}}_{j+1,-}\right) \,, \\ {\bf F}^{n+{1}/{2}}_{z,k+{1}/{2}}&=& {\cal R}\left({\bf U}^{n+{1}/{2}}_{k,+}, {\bf U}^{n+{1}/{2}}_{k+1,-}\right) \,,\nonumber \end{eqnarray} \tag{ 18 }$$

where ${\cal R}(\cdot ,\cdot)$ is the numerical flux resulting from the solution of the Riemann problem. With PLUTO, different Riemann solvers may be chosen at runtime depending on the selected physical module (see Paper I): Rusanov (Lax-Friedrichs), HLL, HLLC (contact), and HLLD are common to both classical and relativistic MHD modules while the Roe solver is available for hydro and MHD. The input states in Equation (18) are obtained by first carrying an evolution step in the normal direction, followed by correcting the resulting values with a transverse flux gradient. This yields the corner-coupled states which, in the absence of parabolic terms (i.e., when ${\mathsfpi}=\mathsf{0}$), are constructed exactly as illustrated by MT. For this reason, we will not repeat it here.

When $\mathsfpi\,{\ne}\, 0$, on the other hand, we adopt a slightly different formulation that does not require any change in the computational stencil. At constant x-faces, for instance, we modify the corner-coupled states to

$$\begin{equation} {\bf U}^{n+\frac{1}{2}}_{i,\pm } = {\bf U}^{*}_{i,\pm } + \frac{\Delta t^n}{2}\left[\sum _{d\ne x}{\bf{\cal L}}_{H,d}^n + \sum _d {\bf{\cal L}}_{P,d}^n\right]\,, \end{equation} \tag{ 19 }$$

where, using Godunov's first-order method, we compute, for example,

$$\begin{equation} {\bf F}^n_{y,j+{1}/{2}} = {\cal R}\big(\bar{{\bf U}}^n_j, \bar{{\bf U}}^n_{j+1}\big)\,, \end{equation} \tag{ 20 }$$

i.e., by solving a Riemann problem between cell-centered states. States at constant y- and z-faces are constructed in a similar manner. The normal predictors U*_± can be obtained either in characteristic or primitive variables as outlined in Sections 3.2 and 3.3, respectively. Parabolic (dissipative) terms are discretized in a flux-conservative form using standard central FD approximations to the derivatives. For any term in the form Π = g(U)∂_xf(U), for instance, we evaluate the right interface flux appearing in Equation (17) with the second-order accurate expression

$$\begin{equation} \Pi _{x,+} \approx g\left(\frac{{\bf U}_i + {\bf U}_{i+1}}{2}\right) \frac{\displaystyle f({\bf U}_{i+1})-f({\bf U}_{i})}{\Delta x}\,, \end{equation} \tag{ 21 }$$

and similarly for the other directions (a similar approach is used by Tóth et al. 2008). We take the solution available at the cell center at t = tⁿ in Equation (19) and at the half time step n + 1/2 in Equation (15). For this latter update, time- and cell-centered conserved quantities may be readily obtained as

$$\begin{equation} {\bf U}^{n+\frac{1}{2}} = \bar{{\bf U}}^n + \frac{\Delta t^n}{2}\sum _d\big({\bf{\cal L}}_{H,d}^n + {\bf{\cal L}}_{P,d}^n\big)\,. \end{equation} \tag{ 22 }$$

The algorithm requires a total of six solutions to the Riemann problem per zone per step and it is stable under the Courant–Friedrichs–Levy (CFL) condition C_a ⩽ 1 (in two dimensions) or C_a ⩽ 1/2 (in three dimensions), where

$$\begin{equation} C_a = \Delta t^{n+1}\max _{ijk}\left[\max _d\left(\frac{\lambda ^{\max }_d}{\Delta x_d}\right) + \sum _d\frac{2{\cal D}^{\max }_d}{\Delta x_d^2}\right] \end{equation} \tag{ 23 }$$

is the CFL number, λ^max_d and ${\cal D}^{\max }_d$ are the (local) largest signal speed and diffusion coefficient in the direction given by d = {x, y, z}, respectively. Equation (23) is used to retrieve the time step Δt^{n + 1} for the next time level if no cooling or reaction terms are present. Otherwise we further limit the time step so that the relative change in pressure and chemical species remains below a certain threshold _c:

$$\begin{equation} \Delta t^{n+1} \rightarrow \min \left[\Delta t^{n+1}, \frac{\epsilon _c\Delta t^n}{\max (|\delta p/p|,|\delta X_\kappa |)}\right], \end{equation} \tag{ 24 }$$

where δp/p and δX_κ are the maximum fractional variations during the source step (Teşileanu et al. 2008). We note that the time step limitation given by Equation (24) does not depend on the mesh size and can be estimated on unrefined cells only. This allows to take full advantage of the adaptivity in time as explained in Section 4.2.

The computation of the normal predictor states can be carried out in characteristic variables by projecting the vector V of primitive variables onto the left eigenvectors l^κ_i ≡ l^κ(V_i) of the primitive system. Specializing to the x-direction,

$$\begin{equation} w^\kappa _{i,l} = {\bf l}_i^\kappa \cdot {\bf V}_{i+l}\,,\quad l = -S,\ldots ,S{,} \end{equation} \tag{ 25 }$$

where κ = 1, ..., N_wave labels the characteristic wave with speed λ^κ and the projection extends to all neighboring zones required by the interpolation stencil of width 2S + 1. The employment of characteristic fields rather than primitive variables requires the additional computational cost associated with the full spectral decomposition of the primitive form of the equations. Nevertheless, it has shown to produce better-behaved solutions for highly nonlinear problems, notably for higher-order methods.

For each zone i and characteristic field κ, we first interpolate w^κ_{i, l} to obtain w^κ_{i, ±}, that is, the rightmost (+) and leftmost (−) interface values from within the cell. The interpolation can be carried out using either fourth-, third- or second-order piecewise reconstruction as outlined later in this section.

Extrapolation in time to tⁿ + Δtⁿ/2 is then carried out by using an upwind selection rule that discards waves not reaching a given interface in Δt/2. The result of this construction, omitting the κ index for the sake of exposition, reads

$$\begin{equation} w^{*}_{i,+} = w^{\rm ref}_{i,+} + \beta _{+}\left\lbrace w_{i,+} - \frac{\nu }{2}[\delta w_i + \delta ^2w_i(3 - 2\nu)] - w^{\rm ref}_{i,+}\right\rbrace{,} \end{equation} \tag{ 26 }$$

$$\begin{equation} w^{*}_{i,-} = w^{\rm ref}_{i,-} + \beta _{-}\left\lbrace w_{i,-} - \frac{\nu }{2} [\delta w_i - \delta ^2w_i(3 + 2\nu)] - w^{\rm ref}_{i,-}\right\rbrace{,} \end{equation} \tag{ 27 }$$

where ν ≡ λΔt/Δx is the Courant number of the κth wave, β_± = (1 ± sign(ν))/2, whereas δw_i and δ²w_i are defined by

$$\begin{equation} \delta w_i = w_{i,+} - w_{i,-}\,,\qquad \delta ^2 w_i = w_{i,+} - 2w_i + w_{i,-}\,. \end{equation} \tag{ 28 }$$

The choice of the reference state w^ref_{i, ±} is somewhat arbitrary and one can simply set w^ref_{i, ±} = w_{i, 0} (Rider et al. 2007), which has been found to work well for flows containing strong discontinuities. Alternatively, one can use the original prescription (Colella & Woodward 1984)

$$\begin{equation} w^{\rm ref}_{i,+} = w_{i,+} - \frac{\nu _{\max }}{2} [\delta w_i + \delta ^2w_i\left(3 - 2\nu _{\max }\right)]\,, \end{equation} \tag{ 29 }$$

$$\begin{equation} w^{\rm ref}_{i,-} = w_{i,-} - \frac{\nu _{\min }}{2} [\delta w_i - \delta ^2w_i\left(3 + 2\nu _{\min }\right)]\,, \end{equation} \tag{ 30 }$$

where ν_max = max (0, max _κ(ν_κ)) and ν_min = min (0, min _κ(ν_κ)) are chosen so as to minimize the size of the term susceptible to characteristic limiting (Colella & Woodward 1984; Miller & Colella 2002; Mignone et al. 2005). However, we note that, in the presence of smooth flows, both choices may reduce the formal second-order accuracy of the scheme since, in the limit of small Δt, contributions carried by waves not reaching a zone edge are not included when reconstructing the corresponding interface value (Stone et al. 2008). In these situations, a better choice is to construct the normal predictors without introducing a specific reference state or, equivalently, by assigning w^ref_{i, ±} = w_{i, ±} in Equations (26) and (27).

The time-centered interface values obtained in characteristic variables through Equations (26) and (27) are finally used as coefficients in the right-eigenvector expansion to recover the primitive variables:

$$\begin{equation} {\bf V}^{*}_{i,\pm } = \sum _\kappa w^{\kappa ,*}_{i,\pm } {\bf r}^\kappa _{i} + \frac{\Delta t^n}{2}\left({\bf S}_{g,i}^n + {\bf S}_{B_x,i}^n\frac{\Delta B_x}{\Delta x} + {\bf S}_{\psi ,i}^n\frac{\Delta \psi }{\Delta x}\right)\,, \end{equation} \tag{ 31 }$$

where r^κ is the right eigenvector associated with the κth wave, S_g is a source term accounting for body forces and geometrical factors, while ${\bf S}_{\rm B_x}$ and S_ψ arise from calculating the interface states in primitive variables rather than conservative ones and are essential for the accuracy of the scheme in multiple dimensions. See MT for a detailed discussion on the implementation of these terms.

As mentioned, the construction of the left and right interface values can be carried out using different interpolation techniques. Although some of the available options have already been presented in the original paper (Mignone et al. 2007), here we briefly outline the implementation details for three selected schemes providing (respectively) fourth-, third-, and second-order spatially accurate interface values in the limit of vanishing time step. Throughout this section we will make frequent usage of the undivided differences of characteristic variables (Equation (25)) such as

$$\begin{equation} \Delta w_{i,+{1}/{2}} = w_{i,+1} - w_{i,0}\,,\quad \Delta w_{i,-{1}/{2}} = w_{i,0} - w_{i,-1}\, \end{equation} \tag{ 32 }$$

for the ith zone.

Piecewise Parabolic Method. The original PPM reconstruction by Colella & Woodward (1984; see also Miller & Colella 2002) can be directly applied to characteristic variables giving the following fourth-order limited interface values:

$$\begin{equation} w_{i,\pm } = \frac{w_{i,0} + w_{i,\pm 1}}{2} \mp \frac{\overline{\Delta w}_{i,\pm 1} - \overline{\Delta w}_{i,0}}{6}\,, \end{equation} \tag{ 33 }$$

where slope limiting is used to ensure that w_{i, ±} are bounded between w_i and w_{i ± 1}:

$$\begin{equation} \overline{\Delta w}_{i,0} = {\rm mm}\left[ \frac{\Delta w_{i,+{1}/{2}}+\Delta w_{i,-{1}/{2}}}{2}, 2{\rm mm}\big(\Delta w_{i,-{1}/{2}}, \Delta w_{i,+{1}/{2}}\big)\right] \end{equation} \tag{ 34 }$$

and

$$\begin{equation} {\rm mm}(a,b) = \frac{{\rm {sign}}(a) + {\rm {sign}}(b)}{2} \min (|a|, |b|) \, \end{equation} \tag{ 35 }$$

is the MinMod function. The original interface values defined by Equation (33) must then be corrected to avoid the appearance of local extrema. By defining δ_± = w_{i, ±} − w_{i, 0}, we further apply the following parabolic limiter

$$\begin{equation} \delta _\pm = \left\lbrace \begin{array}{@{}ll}0 & \quad \hbox{if}\quad \delta _+\delta _- > 0 \\ -2\delta _\mp & \quad \hbox{if}\quad |\delta _\pm | \ge 2|\delta _\mp |\,, \end{array}\right. \end{equation} \tag{ 36 }$$

where the first condition flattens the distribution when w_{i, 0} is a local maximum or minimum whereas the second condition prevents the appearance of an extremum in the parabolic profile. Once limiting has been applied, the final interface values are obtained as

$$\begin{equation} w_{i,\pm } = w_{i,0} + \delta _\pm \,. \end{equation} \tag{ 37 }$$

In some circumstances, we have found that further application of the parabolic limiter (Equation (36)) to primitive variables may reduce oscillations.

Third-order-improved WENO. As an alternative to the popular TVD limiters, the third-order-improved WENO reconstruction proposed by Yamaleev & Carpenter (2009; see also Mignone et al. 2010) may be used. The interpolation still employs a three-point stencil but provides a piecewise parabolic profile that preserves the accuracy at smooth extrema, thus avoiding the well-known clipping of classical second-order TVD limiters. Left and right states are recovered by a convex combination of second-order interpolants into a weighted average of the order of three. The nonlinear weights are adjusted by the local smoothness of the solution so that essentially zero weights are given to non-smooth stencils while optimal weights are prescribed in smooth regions. In compact notation

$$\begin{eqnarray} w_{i,+} &=& \displaystyle w_{i,0} + \frac{a_+\Delta w_{i,+{1}/{2}} + {1}/{2}a_-\Delta w_{i,-{1}/{2}}}{2a_++a_-}\,,\nonumber\\&&[-2mm] \\[-2mm] w_{i,-} &=& \displaystyle w_{i,0} - \frac{a_-\Delta w_{i,-{1}/{2}} + {1}/{2}a_+\Delta w_{i,+{1}/{2}}}{2a_-+a_+}\,, \nonumber \end{eqnarray} \tag{ 38 }$$

where

$$\begin{equation} a_\pm = 1 + \frac{\big(\Delta w_{i,+{1}/{2}} - \Delta w_{i,-{1}/{2}}\big)^2}{\Delta x^2 + \Delta w_{i,\pm /{1}{2}}^2}\,. \end{equation} \tag{ 39 }$$

As one can see, an attractive feature of WENO reconstruction consists in completely avoiding the usage of conditional statements. The improved WENO scheme has enhanced accuracy with respect to the traditional third-order scheme of Jiang & Shu (1996) in regions where the solution is smooth and provides oscillation-free profiles near strong discontinuities.

Linear reconstruction. Second-order TVD limiting is provided by

$$\begin{equation} w_{i,\pm } = w_{i,0} \pm \frac{\overline{\Delta w}_{i,0}}{2}\,, \end{equation} \tag{ 40 }$$

where $\overline{\Delta w}_{i,0}$ is a standard limiter function such as the monotonized-central (MC) limiter (Equation (34)). Other, less steep forms of limiting are the harmonic mean (van Leer 1974):

$$\begin{equation} \overline{\Delta w}_{i,0} = \left\lbrace \begin{array}{@{}ll}\displaystyle \frac{2\Delta w_{i,+{1}/{2}}\Delta w_{i,-{1}/{2}}}{\Delta w_{i,+{1}/{2}}+\Delta w_{i,-{1}/{2}}} \ & \quad {\rm if}\quad \Delta w_{i,+{1}/{2}}\Delta w_{i,-{1}/{2}} > 0 \,, \\ 0 & \quad {\rm otherwise}\,; \end{array}\right. \end{equation} \tag{ 41 }$$

the Van Albada limiter (van Albada et al. 1982):

$${\fontsize{7.8}{9.8}\selectfont{\begin{equation} \overline{\Delta w}_{i,0} = \left\lbrace \begin{array}{@{}ll}\displaystyle \frac{\Delta w_{i,+{1}/{2}}\Delta w_{i,-{1}/{2}}(\Delta w_{i,+{1}/{2}}+\Delta w_{i,-{1}/{2}})}{\Delta w_{i,+{1}/{2}}^2+\Delta w_{i,-{1}/{2}}^2} & {\rm if}\,\, \Delta w_{i,+{1}/{2}}\Delta w_{i,-{1}/{2}} > 0 \\ 0 & {\rm otherwise} \end{array}\right.\,{,} \end{equation}}} \tag{ 42 }$$

or the MinMod limiter, Equation (35).

Linear reconstruction may also be locally employed in place of a higher-order method whenever a strong shock is detected. This is part of a built-in hybrid mechanism that selectively identifies zones within a strong shock in order to introduce additional dissipation by simultaneously switching to the HLL Riemann solver. Even if the occurrence of such situations is usually limited to very few grid zones, this fail-safe mechanism does not sacrifice the second-order accuracy of the scheme and has been found to noticeably improve the robustness of the algorithm avoiding the occurrence of unphysical states. This is described in more detail in Appendix B. However, in order to assess the robustness and limits of applicability of the present algorithm, it will not be employed for the test problems presented here unless otherwise stated.

The normal predictor states can also be directly computed in primitive variables using a simpler formulation that avoids the characteristic projection step. In this case, one-dimensional (1D) left and right states are obtained at t^{n + 1/2} using

$$\begin{equation} {\bf V}^*_{i,\pm } = {\bf V}_i^n + \beta _\pm \left[ {\bf V}({\bf U}_i^{*})-{\bf V}_i^n \pm \frac{\Delta {\bf V}_i}{2}\right]\,, \end{equation} \tag{ 43 }$$

where β₊ = (1 + sgn(λ_max))/2 and β₋ = (1 − sgn(λ_min))/2 may be used to introduce a weak form of upwind limiting, in a similar fashion to Section 3.2. Time-centered conservative variables U*_i follow from a simple conservative MUSCL-Hancock step:

$$\begin{equation} {\bf U}^*_i = {\bf U}_i^n - \frac{\Delta t^n}{2\Delta {\cal V}_x} \left[A_{x,+}{\bf F}\left({\bf V}^n_{i,+}\right) - A_{x,-}{\bf F}\left({\bf V}^n_{i,-}\right)\right] + \frac{\Delta t^n}{2}\hat{{\bf S}}_d, \end{equation} \tag{ 44 }$$

where A_{x, ±} and $\Delta {\cal V_x}$ are the area and volume elements in the x direction, respectively, and Vⁿ_{i, ±} are obtained using linear reconstruction (Equation (40)) of the primitive variables. This approach offers an ease of implementation over the characteristic tracing step as it does not require the eigenvector decomposition of the equations or their primitive form. Notice that, since this step is performed in conservative variables, the multidimensional terms proportional to ∂B_x/∂x do not need to be included.

In PLUTO-CHOMBO, zones are tagged for refinement whenever a prescribed function χ(U) of the conserved variables and of its derivatives exceeds a prescribed threshold, i.e., χ(U) > χ_r. Generally speaking, the refinement criterion may be problem dependent, thus requiring the user to provide an appropriate definition of χ(U). The default choice adopts a criterion based on the second derivative error norm (Löhner 1987), where

$$\begin{equation} \chi ({\bf U}) = \sqrt{ \frac{\sum _d \big|\Delta _{d,+{1}/{2}}\sigma - \Delta _{d,-{1}/{2}}\sigma \big|^2}{\sum _d \big(\big|\Delta _{d,+{1}/{2}}\sigma \big| + \big|\Delta _{d,-{1}/{2}}\sigma \big| +\epsilon \sigma _{d,{\rm ref}}\big)^2 }, } \end{equation} \tag{ 50 }$$

where σ ≡ σ(U) is a function of the conserved variables, $\Delta _{d,\pm \frac{1}{2}}\sigma$ are the undivided forward and backward differences in the direction d, e.g., $\Delta _{x,\pm \frac{1}{2}}\sigma = \pm (\sigma _{i\pm 1} - \sigma _i)$. The last term appearing in the denominator, σ_{d, ref}, prevents regions of small ripples from being refined (Fryxell et al. 2000) and it is defined by

$$\begin{equation} \sigma _{x,{\rm ref}} = |\sigma _{i+1}| + 2|\sigma _i| + |\sigma _{i-1}|. \end{equation} \tag{ 51 }$$

Similar expressions hold when d = y or d = z. In the computations reported in this paper we use = 0.01 as the default value.

In the usual AMR strategy, grids belonging to level ℓ are advanced in time by a sequence of steps with typical size

$$\begin{equation} \Delta t^{\ell ,n+1} = \frac{\Delta t^{0}_{\min }}{(2)^\ell }\,, \end{equation} \tag{ 52 }$$

where we assume, for simplicity, a grid jump of 2. Here, Δt⁰_min is chosen by collecting and re-scaling to the base grid the time steps from all levels available from the previous integration step n:

$$\begin{equation} \Delta t^{0}_{\min } = \min _{0\le \ell \le \ell _{\max }} [(2)^\ell \Delta t^{\ell ,n}]{,} \end{equation} \tag{ 53 }$$

where Δt^{ℓ, n} is computed using Equation (24). However, this procedure may become inefficient in the presence of source terms whose timescale does not depend on the grid size. As an illustrative example, consider a strong radiative shock propagating through a static cold medium. In the optically thin limit, radiative losses are assumed to be local functions of the state vector, but they do not involve spatial derivatives. If the fastest timescale in the problem is dictated by the cooling process, the time step estimate should then become approximately the same on all levels, Δt ≈ Δt^ℓ_rad ≈ Δt⁰_rad, regardless of the mesh size. However from the previous equations, one can see that finer levels with ℓ > 0 will advance with a time step (2)^ℓ smaller than required by the single grid estimate. Equation (52) is nevertheless essential for proper synchronization between nested levels.

Simple considerations show that this deficiency may be cured by treating split and leaf cells differently. Split zones in a given level ℓ are, in fact, overwritten during the projection step using the more accurate solution computed on children cells belonging to level ℓ + 1. Thus, accurate integration of the source term is not important for these cells and could even be skipped. From these considerations, one may as well evaluate the source-term-related time step on leaf cells only, where the accuracy and stability of the computed solution is essential. This simple expedient should speed up the computations by a factor of approximately (2)^ℓ, thus allowing us to take full advantage of the refinement offered by the AMR algorithm without the time step restriction. Besides, this should not alter or degrade the solution computed during this single hierarchical integration step as long as the projection step precedes the regrid process.

The proposed modification is expected to be particularly efficient in those problems where radiative losses are stronger in proximity of steep gradients.

Both PLUTO and PLUTO-CHOMBO support calculations in parallel computing environments through the MPI library. Since the time evolution of the AMR hierarchy is performed recursively, from lower to upper levels, each level of refinement is parallelized independently by distributing its boxes to the set of processors. The computation on a single box has no internal parallelization. Boxes are assigned to processors by balancing the computational load on each processor. Currently, the workload of a single box is estimated by the number of grid points of the box, considering that the integration requires approximately the same amount of flops per grid point. This is not strictly true in some specific case, e.g., in the presence of optically thin radiative losses, and a strategy to improve the load balance in such situations is currently under development. On the basis of the box workloads, CHOMBO's load balancer uses the Kernighan–Lin algorithm for solving knapsack problems.

CHOMBO weak scaling performance has been thoroughly benchmarked: defining the initial setup by spatially replicating a 3D hydrodynamical problem proportionally to the number of CPUs employed, the execution time stays constant with excellent approximation (see Van Straalen et al. 2009 for more details).

To test the parallel performance of PLUTO-CHOMBO in real applications, we performed a number of strong scaling tests by computing the execution time as a function of the number of processors for a given setup. While this is a usual benchmark for static grid calculations, in the case of AMR computations this diagnostic is strongly problem dependent and difficult to interpret.

In order to find some practical rule to improve the scaling of AMR calculations, we investigated the dependency of the parallel performance on some parameters characterizing the adaptive grid structure: the maximum box size allowed and the number of levels of refinement employed using different refinement ratios. As a general rule, the parallel performance deteriorates when the number of blocks per level becomes comparable to the number of processors or, alternatively, when the ideal workload per processor (i.e., the number of grid cells of a level divided by the number of CPUs) becomes comparable to the maximum box size. As we will show, decreasing the maximum possible block size can sensibly increase the number of boxes of a level and therefore improves the parallel performance. On the other hand, using less refinement levels with larger refinement ratios to achieve the same maximum resolution can lower the execution time, reducing the parallel communication volume and avoiding the integration of intermediate levels.

In Sections 5 and 6, we present several parallel scaling tests and their dependence on the aforementioned grid parameters.

The shock tube test problem is a common benchmark for an accurate description of both continuous and discontinuous flow features. In the following, we consider 1D and 3D configurations of standard shock tubes proposed by Torrilhon (2003) and Ryu & Jones (1995).

5.1.1. One-dimensional Shock Tube

Following Torrilhon (2003), we address the capability of the AMR scheme to handle and refine discontinuous features as well as to correctly resolve the non-uniqueness issue of MHD Riemann problems in FV schemes (Torrilhon 2003; Torrilhon & Balsara 2004, and references therein). Left and right states are given by

$${\fontsize{9.7}{11.9}\selectfont{\begin{equation} \left\lbrace \begin{array}{@{}lr}{\bf V}_L = \displaystyle \left(1,0,0,0,1,1,0,1\right)^{\rm T} & \mathrm{for}\quad x_1 < 0 \,, \\ {\bf V}_R = \displaystyle \left(0.2,0,0,0,1,\cos (\alpha),\sin (\alpha),0.2\right)^{\rm T} & \mathrm{for}\quad x_1 > 0 \,, \end{array}\right. \end{equation}}} \tag{ 54 }$$

where V = (ρ, v_x, v_y, v_z, B_x, B_y, B_z, p) is the vector of primitive variables. As discussed by Torrilhon (2003), for a wide range of initial conditions MHD Riemann problems have unique solutions, consisting of Lax shocks, contact discontinuities, or rarefaction waves. Nevertheless, certain sets of initial values exist that can result in non-unique solutions. When this occurs, along with the regular solution arises one that allows for irregular MHD waves, for example, a compound wave. A special case where the latter appears is when the initial transverse velocity components are zero and the magnetic field vectors on either side of the interface are anti-parallel. Such a case was noted by Brio & Wu (1988) and can be reproduced by simply choosing α = π in our initial condition.

A 1D non-unique solution is calculated using a static grid with 400 zones, x ∈ [ − 1, 1.5]. Left and right boundaries are set to outflow and the evolution stops at time t = 0.4, before the fast waves reach the borders. The resulting density profile is shown in Figure 6. The solid lines denote the two admissible exact solutions: the regular (red) and the one containing a compound wave (black), the latter situated at x ∼ − 0.24. It is clear that the solution obtained with the Godunov-type code is the one with the compound wave (symbols).

Zoom In Zoom Out Reset image size

The crucial problematic of this test occurs when α is close to but not exactly equal to π. Torrilhon (2003) has proven that regardless of scheme, the numerical solution will erroneously tend to converge to an "irregular" one similar to α = π (pseudoconvergence), even if the initial conditions should have a unique, regular solution. This pathology can be cured either with high-order schemes (Torrilhon & Balsara 2004) or with a dramatic increase in resolution on the region of interest, proving AMR to be a useful tool. To demonstrate this we choose α = 3 for the field's twist.

Starting from a coarse grid of 512 computational zones, we vary the number of refinement levels, with a consecutive jump ratio of 2. The 10 level run (solid red line) incorporates also a single jump ratio of 4 between the sixth and seventh refinement levels, reaching a maximum equivalent resolution of 1,048,576 zones (see Figure 7). The refinement criterion is set upon the variable σ = (B²_x + B_y² + B²_z)/ρ using Equation (50) with a threshold χ_r = 0.03, whereas integration is performed using PPM with a Roe Riemann solver and a Courant number C_a = 0.9. As resolution increases the compound wave disentangles and the solution converges to the expected regular form (Torrilhon 2003; Fromang et al. 2006). In Table 2, we compare the CPU running time of the AMR runs versus static uniform grid computations at the same effective resolution. With 5 and 10 levels of refinement (effective resolutions 16,384 and 1,048,576 zones, respectively) the AMR approach is ∼32 and ∼1978 times faster than the uniform mesh computation, respectively.

Zoom In Zoom Out Reset image size

Table 1. Systems of Coordinates Adopted in PLUTO-CHOMBO

Label	Cartesian	Cylindrical
x1	x	r
x2	y	z
x3	z	/
V1	x	r2/2
V2	y	z
V3	z	/
A1, +	ΔyΔz	r+Δz
A2, +	ΔzΔx	rΔr
A3, +	ΔxΔy	/
${\cal V}$	ΔxΔyΔz	rΔrΔz

Download table as:  ASCIITypeset image

Table 2. CPU Running Time for the 1D MHD Shock Tube Using Both Static and AMR Computations

Static Run		AMR Run			Gain
Nx	Time	Level	Ref Ratio	Time
	(s)			(s)
512	0.5	0	2	0.5	1
1024	1.9	1	2	1.1	1.7
2048	7.5	2	2	2.3	3.3
4096	31.6	3	2	4.5	7.0
8192	138.0	4	2	8.7	15.9
16384	546.1	5	2	16.9	32.3
1048576	2.237 × 106a	10	2 (4)	1131.1	1977.8

Notes. The first and second columns give the number of points N_x and corresponding CPU for the static grid run (no AMR). The third, fourth, and fifth columns give, respectively, the number of levels, the refinement ratio, and CPU time for the AMR run at the equivalent resolution. The last row refers to the solid red line of Figure 7, where a jump ratio of four was introduced between levels six and seven to reach an equivalent of ∼10⁶ grid points. The last column shows the corresponding gain factor calculated as the ratio between static and AMR execution time. ^aCPU time has been inferred from ideal scaling.

Download table as:  ASCIITypeset image

5.1.2. Three-dimensional Shock Tube

The second Riemann problem was proposed by Ryu & Jones (1995) and later considered also by Tóth (2000), Balsara & Shu (2000), Gardiner & Stone (2008), Mignone & Tzeferacos (2010), and Mignone et al. (2010). An initial discontinuity is described in terms of primitive variables as

$$\begin{equation} \left\lbrace \begin{array}{@{}ll}{\bf V}_L = \displaystyle \left(1.08,1.2,0.01,0.5,\frac{2}{\sqrt{4\pi }},\frac{3.6}{\sqrt{4\pi }} \frac{2}{\sqrt{4\pi }},0.95\right)^{\rm T} &\\ \ms \qquad \mathrm{for}\quad x_1 < 0 &\\ {\bf V}_R = \displaystyle \left(1,0,0,0,\frac{2}{\sqrt{4\pi }},\frac{4}{\sqrt{4\pi }}, \frac{2}{\sqrt{4\pi }},1\right)^{\rm T} &\\ \ms \qquad \mathrm{for}\quad x_1 > 0 &\end{array}\right., \end{equation} \tag{ 55 }$$

where V = (ρ, v₁, v₂, v₃, B₁, B₂, B₃, p) is the vector of primitive variables. The subscript "1" gives the direction perpendicular to the initial surface of discontinuity whereas "2" and "3" correspond to the transverse directions. We first obtain a 1D solution on the domain x ∈ [ − 0.75, 0.75] using 6144 grid points, stopping the computations at t = 0.2.

In order to test the ability of the AMR scheme to maintain the translational invariance and properly refine the flow discontinuities, the shock tube is rotated in a 3D Cartesian domain. The coarse level of the computational domain consists of [384 × 4 × 4] zones and spans [ − 0.75, 0.75] in the x-direction while y, z ∈ [0, 0.015625]. The rotation angles, γ around the y-axis and α around the z-axis, are chosen so that the planar symmetry is satisfied by an integer shift of cells (n_x, n_y, n_z). The rotation matrix can be found in Mignone et al. (2010). By choosing tan α = −r₁/r₂ and tan β = tan γ/cos α = r₁/r₃, one can show (Gardiner & Stone 2008) that the three integer shifts n_x, n_y, n_z must obey

$$\begin{equation} n_x - n_y \frac{r_1}{r_2} + n_z \frac{r_1}{r_3} = 0 \,, \end{equation} \tag{ 56 }$$

where cubed cells have been assumed and (r₁, r₂, r₃) = (1, 2, 4) will be used. Computations stop at t = 0.2cos αcos γ, once again before the fast waves reach the boundaries. We employ four refinement levels with consecutive jumps of 2, corresponding to an equivalent resolution of 6144 × 64 × 64 zones. The refinement criterion is based on the normalized second derivative of (|B_x| + |B_y|)ρ with a threshold value χ_r = 0.1. Integration is done with PPM reconstruction, a Roe Riemann solver, and a Courant number of C_a = 0.4.

The primitive variable profiles (symbols) are displayed in Figure 8 along the x-direction,⁵ together with the 1D reference solution in the x ∈ [ − 0.25, 0.55] region. In agreement with the solution of Ryu & Jones (1995), the wave pattern produced consists of a contact discontinuity that separates two fast shocks, two slow shocks, and a pair of rotational discontinuities. A 3D close-up of the top-hat feature in the density profile is shown in Figure 9, along with AMR levels and mesh. The discontinuities are captured correctly, and the AMR grid structure respects the plane symmetry. Our results favorably compare with those of Gardiner & Stone (2008), Mignone & Tzeferacos (2010), Mignone et al. (2010), and previous similar 2D configurations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The AMR computation took approximately 3 hr and 53 minutes on two 2.26 GHz quad-core Intel Xeon processors (eight cores in total). For the sake of comparison, we repeated the same computation on a uniform mesh of 768 × 8 × 8 zones (1/8 of the effective resolution) with the static grid version of PLUTO employing ≈79 s. Thus, extrapolating from ideal scaling, the computational cost of the fixed grid calculation is expected to increase (at least) by a factor 2¹² giving an overall gain of the AMR over the uniform grid approach of ∼23.

The next test problem considers the 2D advection of a magnetic field loop. This test, proposed by Gardiner & Stone (2005), aims to benchmark the scheme's dissipative properties and the correct discretization balance of multidimensional terms through monitoring the preservation of the initial circular shape of the loop.

As in Gardiner & Stone (2005) and Fromang et al. (2006), we define the computational domain by x ∈ [ − 1, 1] and y ∈ [ − 0.5, 0.5] discretized on a coarse grid of 64 × 32 grid cells. In the initial condition, both density and pressure are uniform and equal to 1, while the velocity of the flow is given by ${\bf v} = V_0\cos \alpha \hat{{\bf e}}_x + V_0\sin \alpha \hat{{\bf e}}_y$ with $V_0 = \sqrt{5}$, $\sin \alpha = 1/\sqrt{5}$, and $\cos \alpha = 2/\sqrt{5}$. The magnetic field is then defined through the magnetic vector potential as

$$\begin{equation} A_z = \left\lbrace \begin{array}{@{}ll}A_0(R-r) & \hbox{if} \quad r \le R \,, \\ 0 & \hbox{if} \quad r > R \,, \end{array} \right. \end{equation} \tag{ 57 }$$

where A₀ = 10⁻³, R = 0.3, and $r=\sqrt{x^2 + y^2}$. The simulation evolves until t = 2, when the loop has thus performed two crossings through the periodic boundaries. The test is repeated with two, three, and four levels of refinement (jump ratio of 2), resulting to equivalent resolutions of [256 × 128], [512 × 256], and [1024 × 512], respectively. Refinement is triggered whenever the second derivative error norm of (B²_x + B_y²) × 10⁶, computed via Equation (50), exceeds the threshold χ_r = 0.1. The integration is carried out utilizing WENO reconstruction and the Roe Riemann solver, with C_a = 0.4.

The temporal evolution of magnetic energy density is seen in Figure 10, with four levels of refinement. As the field loop is transported inside the computational domain, the grid structure changes to follow the evolution and the field lines retain the initial circular form. An efficient way to quantitatively measure the diffusive properties of the scheme is to monitor the dissipation of the magnetic energy. In the top panel of Figure 11, we plot the normalized mean magnetic energy as a function of time. By increasing the levels of refinement the dissipation of magnetic energy decreases, with 〈B²〉 ranging from ∼94% to ∼98% of the initial value. In order to quantify the computational gain of the AMR scheme, we repeat the simulations with a uniform grid resolved onto as many points as the equivalent resolution, without changing the employed numerical method. The speed-up is reported in the bottom panel of Figure 11.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Magnetic reconnection refers to a topological rearrangement of magnetic field lines with opposite directions, accompanied with a conversion of magnetic energy into kinetic and thermal energy of the plasma. This is believed to be the basic mechanism behind energy release during solar flares. The first solution to the problem was given independently by Sweet (1958) and Parker (1957), treating it as a 2D boundary layer problem in the laminar limit.

According to the Sweet–Parker model, the magnetic field's convective inflow is balanced by ohmic diffusion. Along with the assumption of continuity, this yields a relation between reconnection and plasma parameters. If L and δ are the boundary layer's half length and width, respectively, we can write the reconnection rate ${\cal E}$ as

$$\begin{equation} {\cal E}\equiv \frac{u_{{\rm in}}}{u_{{\rm out}}}\sim \frac{\delta }{L}\sim \frac{1}{\sqrt{S}}. \end{equation} \tag{ 58 }$$

With u_in and u_out we denote the inflow and outflow speeds, into and out of the boundary layer, respectively. The Lundquist number for the boundary layer is defined as S = u_AL/η, with u_A being the Alfvén velocity directly upstream of the layer, η the magnetic resistivity, and L the layer's half-length. This dependency of the reconnection rate with the square root of magnetic resistivity is called the Sweet–Parker scaling and has been verified both numerically (Biskamp 1986; Uzdensky & Kulsrud 2000) and experimentally (Ji et al. 1998).

Following the guidelines of the Geospace Environment Modeling (GEM) Magnetic Reconnection Challenge (Birn et al. 2001), the computational domain is a 2D Cartesian box, with x ∈ [ − L_x/2, L_x/2] and y ∈ [ − L_y/2, L_y/2] where we choose L_x = 25.6 and L_y = 12.8. The initial condition consists of a Harris current sheet: the magnetic field configuration is described by B_x(y) = B₀tanh (y/λ) whereas the flow's density is ρ = ρ₀sech²(y/λ) + ρ_∞, where λ = 0.5, ρ₀ = 1, and ρ_∞ = 0.2. The flow's thermal pressure is deduced assuming equilibrium with magnetic pressure, P = B²₀/2 = 0.5. The initial magnetic field components are perturbed via

$$\begin{equation} \begin{array}{ll}dB_x = - \Psi _0(\pi /L_y)\sin (\pi y/L_y)\cos (2\pi x/L_x)\\[5pt] dB_y = \Psi _0(2\pi /L_x)\sin (2\pi x/L_x)\cos (\pi y /L_y). \end{array} \end{equation} \tag{ 59 }$$

where Ψ₀ = 0.1. The coarse grid consists of [64  ×  32] points, and additional levels of refinement are triggered using the following criterion based on the current density:

$$\begin{equation} \frac{|\Delta _xB_y - \Delta _yB_x|}{|\Delta _xB_y|+|\Delta _yB_x|+\sqrt{\rho }/\xi } > \frac{\chi _{\min } - \chi _{\max }}{1 + (1-\xi)^2} + \chi _{\max }, \end{equation} \tag{ 60 }$$

where Δ_xB_y and Δ_yB_x are the undivided central differences of B_y and B_x in the x- and y-direction, respectively, and ξ = Δx⁰/Δx^ℓ ⩾ 1 is the ratio of grid spacings between the base (0) and current level (l). The threshold values χ_min = 0.2 and χ_max = 0.37 are chosen in such a way that refinement becomes increasingly harder for higher levels. We perform test cases using either five levels of refinement with a consecutive jump ratio of 2, or three levels with a jump ratio of 2:4:4 reaching, in both cases, an equivalent resolution of 2048 × 1024 mesh points. Boundaries are periodic in the x-direction, whereas perfectly conducting boundary walls are set at y = ±L_y/2. We follow the computations until t = 100 using PPM reconstruction with a Roe Riemann solver and a Courant number C_a = 0.8. In Figure 12, we display the temporal evolution of the current density for a case where the uniform resistivity is set to η = 8 × 10⁻³. A reconnection layer is created in the center of the domain, which predisposes resistive reconnection (Biskamp 1986). In agreement with Reynolds et al. (2006) the maximum value of the current density decreases with time (Figure 4 of that study). As seen in the first two panels (t = 0, 50), the refinement criterion is adequate to capture correctly both the boundary layer and the borders of the magnetic island structure. In the rightmost panel (t = 100), we also draw sample magnetic field lines to better visualize the reconnection region.

Zoom In Zoom Out Reset image size

In order to compare our numerical results with theory, we repeated the computation, varying the value of the magnetic resistivity η. For small values of resistivity (large Lundquist numbers S), the boundary layer is elongated and presents large aspect ratios of A ≡ L/δ. Biskamp (1986) reports that for A beyond the critical value of A_crit ≃ 100 the boundary layer is tearing unstable, limiting the resistivity range in which the Sweet–Parker reconnection model operates. Since in the Sun's corona S can reach values of ∼10¹⁴ ≫ S_max, secondary island formation must be taken into account (Cassak & Drake 2009). In this context, ensuring that we respect Biskamp's stability criterion, we calculate the temporal average of δ/L, analogous to the reconnection rate ${\cal E}$, for various η and reproduce the Sweet–Parker scaling (Figure 13). The boundary layer's half -width (δ) and half-length (L) are estimated from the e-folding distance of the peak of the electric current, while the AMR scheme allows us to economically resolve the layer's thickness with enough grid points.

Zoom In Zoom Out Reset image size

Parallel performance for this problem is shown in Figure 14, where we plot the speed-up $S=T_1/T_{N_{\rm CPU}}$ as a function of the number of processors N_CPU for the three- and five-level AMR computations as well as for fixed uniform grid runs carried out at the equivalent resolution of 2048 × 1024 zones. Here, T₁ is the same reference constant for all calculations and is equal to the (inferred) running time of the single processor static mesh computation while $T_{N_{\rm CPU}}$ is the execution time measured with N_CPU processors. The scaling reveals an efficiency (defined as S/N_CPU) larger than 0.8 for less than 256 processors with the three- and five-level computations being, respectively, eight to nine and four to five times faster than the fixed grid approach. The number of blocks on the finest level is maximum at the end of integration and is slightly larger for the three-level run (1058 versus 835). This result indicates that using fewer levels of refinement with larger grid ratios can be more efficient than introducing more levels with consecutive jumps of 2, most likely because of the reduced integration cost due to the missing intermediate levels and the decreased overhead associated with coarse-fine level communication and grid generation process. Efficiency quickly drops when the number of CPU tends to become, within a factor of between two and three, comparable to the number of blocks.

Zoom In Zoom Out Reset image size

The current sheet problem proposed by Gardiner & Stone (2005) and later considered by Fromang et al. (2006) in the AMR context is particularly sensitive to numerical diffusion. The test problem follows the evolution of two current sheets, initialized through a discontinuous magnetic field configuration. Driven solely by numerical resistivity, reconnection processes take place, making the resulting solution highly susceptible to grid resolution.

The initial condition is discretized onto a Cartesian 2D grid x, y ∈ [0, 2], with 64 × 64 zones at the coarse level. The fluid has uniform density ρ = 1 and thermal pressure P = 0.1. Its bulk flow velocity v is set to zero, allowing only for a small perturbation in v_x = v₀ sin (π y), where v₀ = 0.1. The initial magnetic field has only one non-vanishing component in the vertical direction,

$$\begin{equation} B_y = \left\lbrace \begin{array}{@{}ll}-B_0 & \hbox{if} \quad |x-1| \le 0.5{,} \\ \,\,\,\,B_0 & \hbox{otherwise}{,} \end{array} \right. \end{equation} \tag{ 61 }$$

where B₀ = 1, resulting in a magnetically dominated configuration. Boundaries are periodic and the integration terminates at t = 4. We activate refinement whenever the maximum between the two error norms (given by Equation (50)) computed with the specific internal energy and the y-component of the magnetic field exceeds the threshold value χ_r = 0.2. In order to filter out noise between magnetic island we set = 0.05. We allow for four levels of refinement and carry out integration with the Roe Riemann solver, the improved WENO reconstruction, and a CFL number of 0.6.

In Figure 15, we show temporal snapshots of pressure profiles for t = 0.5, 1.0, 1.5, 2, 3.5, and 4.0, along with sample magnetic field lines. Since no resistivity has been specified, the elongated current sheets are prone to tearing instability as numerical resistivity is small with respect to Biskamp's criterion (Biskamp 1986). Secondary islands, "plasmoids," promptly form and propagate parallel to the field in the y-direction. As reconnection occurs, the field is dissipated and its magnetic energy is transformed into thermal energy, driving Alfvén and compressional waves which further seed reconnection (Gardiner & Stone 2005; Lee & Deane 2009). Due to the dependency of numerical resistivity on the field topology, reconnection events are most probable at the nodal points of v_x. At the later stages of evolution, the plasmoids eventually merge in proximity of the anti-nodes of the transverse speed, forming four larger islands, in agreement with the results presented in Gardiner & Stone (2005). A close-up on the bottom left island is shown in Figure 16, at the end of integration. The refinement criterion on the second derivative of thermal pressure, as suggested by Fromang et al. (2006), efficiently captures the island features of the solution.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this example, we consider the dynamical evolution of two fluids of different densities initially in hydrostatic equilibrium. If the lighter fluid is supporting the heavier fluid against gravity, the configuration is known to be subject to the Rayleigh–Taylor instability. Our computational domain is the box spanned by x, z ∈ [ − 1/2, 1/2], y ∈ [ − 3/2, 1/2] with gravity pointing in the negative y-direction, i.e., g = (0, −1, 0). The two fluids are separated by an interface initially lying in the xz-plane at y = 0, with the heavy fluid ρ = ρ_h at y > 0 and ρ = ρ_l at y < 0. In the current example we employ ρ_h = 4, ρ_l = 1 and specify the pressure from the hydrostatic equilibrium condition:

$$\begin{equation} p(y) = \frac{100}{\gamma } - \rho y \,, \end{equation} \tag{ 62 }$$

so that the sound crossing time in the light fluid is 0.1. The seed for instability is set by perturbing the velocity field at the center of the interface:

$$\begin{equation} v_y = - \frac{\exp [-(5r)^2]}{10\cosh (10y)^2}, \end{equation} \tag{ 63 }$$

where $r=\sqrt{x^2+z^2}$. We assume a constant magnetic field B = (B_x, 0, 0) parallel to the interface and oriented in the x-direction with different field strengths, B_x = 0 (hydro limit), B_x = 0.2B_c (moderate field), and B_x = 0.6B_c (strong field). Here, $B_c = \sqrt{(\rho _h - \rho _l)gL}$ is the critical field value above which instabilities are suppressed (Stone & Gardiner 2007). We use the PPM method with the Roe Riemann solver and a Courant number C_a = 0.45. The base grid has 16 × 32 × 16 cells and we perform two sets of computations using (1) four refinement levels with a grid spacing ratio of 2 and (2) two refinement levels with a grid jump of 4, in both cases achieving the same effective resolution (256 × 512 × 256). Refinement is triggered using the second derivative error norm of density and a threshold value χ_ref = 0.5. Periodic boundary conditions are imposed at the x- and z-boundaries while fixed boundaries are set at y = −3/2 and y = 1/2.

Results for the unmagnetized, weakly, and strongly magnetized cases are shown at t = 3 in Figure 17. In all cases, we observe the development of a central mushroom-shaped finger. Secondary instabilities due to Kelvin–Helmholtz modes develop on the side and small-scale structures are gradually suppressed in the direction of the field as the strength is increased. Indeed, as pointed out by Stone & Gardiner (2007), the presence of a uniform magnetic field has the effects of reducing the growth rate of modes parallel to it, although the interface still remains Rayleigh–Taylor unstable in the perpendicular direction due to interchange modes. As a net effect, the evolution becomes increasingly anisotropic as the field strengthens.

Zoom In Zoom Out Reset image size

Parallel performance is plotted, for the moderate field case, in Figure 18 where we measure the speed-up factors of the four- and two-level computations versus the number of CPUs. Here, the speed-up is defined by $S=T_1/T_{N_{\rm CPU}}$, where T₁ is the inferred running time relative to the four-level computation on a single processor. The two-level case shows improved performance over the four-level calculation both in terms of CPU cost as well as parallel efficiency S/N_CPU. Indeed the relative gain between the two cases approaches a factor of two for an increasing number of CPUs. Likewise, the efficiency of the fewer level case remains above ≳ 0.8 up to 512 processors and drops to ∼0.61 (versus ∼0.44 for the four-level case) at the largest number of employed cores (1024), which is less than the number of blocks on the finest grid.

Zoom In Zoom Out Reset image size

The shock–cloud interaction problem has been extensively studied and used as a standard benchmark for the validation of MHD schemes and intercode comparisons (see Dai & Woodward 1994; Balsara 2001a; Lee & Deane 2009, and references therein). In an astrophysical context, it addresses the fundamental issue of the complex morphology of supernova remnants as well as their interaction with the interstellar medium. Energy, mass, and momentum exchange leading to cloud crushing strongly depends on the orientation of the magnetic field and the resulting anisotropy of thermal conduction.

Following Orlando et al. (2008), we consider the 2D Cartesian domain x ∈ [ − 4, 4], y ∈ [ − 1.4 × 6.6] (in parsecs) with the shock front propagating in the positive y-direction and initially located at y = −1. Ahead of the shock, for y > −1, the hydrogen number density n₁ has a radial distribution of the form:

$$\begin{equation} n_1 = n_a - \frac{n_a - n_c}{\cosh [\sigma (r/r_{\rm cl})^\sigma]}{,}\quad r = \sqrt{x^2+y^2}, \end{equation} \tag{ 64 }$$

where n_c = 1 cm⁻³ is the hydrogen number density at the cloud's center, n_a = 0.1n_c is the ambient number density, r_cl = 1 pc is the cloud's radius (σ = 10), and r is the radial distance from the center of the cloud. The ambient medium has a uniform temperature T_a = 10⁴ K in pressure equilibrium with the cloud, with a thermal pressure of p₁ = 2k_Bn_aT_a where k_B is the Boltzmann constant (we assume a fully ionized gas). Downstream of the shock, density, and transverse components of the magnetic field are compressed by a factor of n₂/n₁ = (Γ + 1)/(Γ − 1) while pressure and normal velocity are given by

$$\begin{equation} p_2 = 2n_2k_BT_s{,}\quad v_{y,2} = \sqrt{\frac{2}{\Gamma -1}\frac{2k_BT_s X}{m_{\rm H}}}, \end{equation} \tag{ 65 }$$

where X = 1/1.26 is the hydrogen mass fraction and T_s = 4.7 × 10⁶ K is the post-shock temperature. The normal component of the magnetic field (B_y) remains continuous through the shock front.

We perform two sets of simulations with a magnetic field strength of |B|  =  1.31 μG, initially parallel (case Bx4) or perpendicular (case By4) to the shock front. Adopting the same notations as Orlando et al. (2008), we solve in each case the MHD equations with (TN) and without (NN) thermal conduction effects for a total of four cases. Thermal conductivity coefficients along and across the magnetic field lines (see Equation (6)) are given, in cgs units, by

$$\begin{equation} k_\parallel = \displaystyle 9.2\times 10^{-7}T^{5/2} {,} \quad k_\perp = \displaystyle 5.4\times 10^{-16} \frac{n_{\rm H}^2}{\sqrt{T} {\bf B}^2} \,. \end{equation} \tag{ 66 }$$

The resolution of the base grid is set to 32² points and five levels of refinement are employed, yielding an equivalent resolution of 1024 × 1024. Upwind characteristic tracing (Section 3.2) together with PPM reconstruction (Equation (33)) and the HLLD Riemann solver are chosen to advance the equations in time. Open boundary conditions are applied on all sides, except at the lower y-boundary where we keep constant inflow values. The CFL number is C_a = 0.8. Refinement is triggered upon density with a threshold of χ_ref = 0.15.

In Figures 19 and 20, we show the cloud evolution for the four different cases at three different times t = 1, 2, 3 (in units of the cloud crushing time τ_cc = 5.4 × 10³ yr). Our results are in excellent agreement with those of Orlando et al. (2008), confirming the efficiency of thermal conduction in suppressing the development of hydrodynamic instabilities at the cloud borders. The topology of the magnetic field can be quite effective in reducing the efficiency of thermal conduction. For a field parallel to the shock front (TN-Bx4), the cloud's expansion and evaporation are strongly limited by the confining effect of the enveloping magnetic field. In the case of a perpendicular field (TN-By4), on the contrary, thermal exchange between the cloud and its surroundings becomes more efficient in the upwind direction of the cloud and promotes the gradual heating of the core and the consequent evaporation in a few dynamical timescales. In this case, heat conduction is strongly suppressed laterally by the presence of a predominantly vertical magnetic field.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Radiative jet models with a variable ejection velocity are commonly adopted to reproduce the knotty structure observed in collimated outflows from young stellar objects. The time variability leads to the formation of a chain of perturbations that eventually steepen into pairs of forward/reverse shocks (the "internal working surfaces"; see, for instance, de Colle & Raga 2006; Raga et al. 2007; Teşileanu et al. 2009) which are held responsible for the typical spectral signature of these objects. While these internal working surfaces travel down the jet, the emission from post-shocked regions occurs on a much smaller spatial and temporal scale, thus posing a considerable challenge even for AMR codes.

As an example application, we solve the MHD equations coupled to the chemical network described in Teşileanu et al. (2008), evolving the time-dependent ionization and collisionally excited radiative losses of the most important atomic/ionic species: H, He, C, N, O, Ne, and S. While hydrogen and helium can be at most singly ionized, we only consider the first three (instead of five) ionization stages of the other elements, which suffice for the range of temperatures and densities considered here. This amounts to a total of 19 additional non-homogeneous continuity/rate equations, such as Equation (3).

The initial condition consists of an axisymmetric supersonic collimated beam in cylindrical coordinates (r, z) in equilibrium with a stationary ambient medium. Density and axial velocity can be freely prescribed as

$$\begin{equation} \rho (r) = \rho _a + \frac{\rho _j-\rho _a}{\cosh \big(r^4\big/r_j^4\big)}\,,\quad v_z(r) = \frac{v_j}{\cosh \big(r^4\big/r_j^4\big)}{,}\quad \end{equation} \tag{ 67 }$$

where ρ_j and ρ_a = ρ_j/5 are the jet and ambient mass density, v_j = 200 km s⁻¹ is the jet velocity, and r_j = 2 × 10¹⁵ cm is the radius. The jet density is given by ρ_j = n_jμ_jm_a, where n_j = 10³ is the total particle number density, μ_j is the mean molecular weight, and m_a is the atomic mass unit. Both the initial jet cross section and the environment are neutral with the exception of C and S, which are singly ionized. The steady state results from a radial balance between the Lorentz and pressure forces and has to satisfy the time-independent momentum equation:

$$\begin{equation} \frac{dp}{dr} = - \frac{1}{2}\frac{1}{r^2}\frac{d(rB_\phi)^2}{dr}\,, \end{equation} \tag{ 68 }$$

where, for simplicity, we neglect rotations and assume a purely azimuthal magnetic field of the form

$$\begin{equation} B_\phi (r) = -\frac{B_m}{r}\sqrt{1 - \exp [-(r/a)^4]}{.} \end{equation} \tag{ 69 }$$

Here, a = 0.9r_j is the magnetization radius and B_m is a constant that depends on the jet and ambient temperatures. Direct integration of Equation (68) yields the equilibrium profile

$$\begin{equation} p(r) = p_j - \frac{1}{2}\frac{B_m^2\sqrt{\pi } \,\hbox{erf}(r^2/a^2)}{a^2}\,, \end{equation} \tag{ 70 }$$

where p_j = ρ_jk_BT_j/(μ_jm_a) is the jet pressure on the axis, μ_j is the mean molecular weight, m_a is the atomic mass unit, k_B is the Boltzmann constant, and T_j is the jet temperature. The value of B_m is recovered by solving Equation (70) in the r → ∞ limit given the jet and ambient temperatures T_j = 5 × 10³ K and T_a = 10³ K, respectively. Our choice of parameters is similar to the one adopted by Raga et al. (2007).

The computational domain is defined by r/r_j ∈ [0, 12.5], z/r_j ∈ [0, 25] with axisymmetric boundary conditions holding at r = 0. At z = 0 we keep the flow variables constant and equal to the equilibrium solution and add a sinusoidal time variability of the jet velocity with a period of 20 yr and an amplitude of 40 km s⁻¹ around the mean value. Free outflow is assumed on the remaining sides. We integrate the equations using linear reconstruction with the MC limiter, Equation (34), and the HLLC Riemann solver. The time step is computed from Equations (23) and (24) using a Courant number C_a = 0.6, a relative tolerance _c = 0.02 and following the considerations given in Section 4.2. Shock-adaptive hybrid integration is used by locally switching to the more diffusive MinMod limiter (Equation (35)) and the HLL Riemann solver, according to the mechanism outlined in Appendix B.

The base grid consists of 128 × 256 zones, and five additional levels with grid refinement jumps of 2:2:2:4:4 are used. At the effective resolution of 16,384 × 32,768 zones the minimum length scale that can be resolved corresponds to 1.526 × 10¹² cm (≈0.1 AU), the same as model M3 of Raga et al. (2007). Zones are tagged for refinement when the second derivative error norm of density exceeds the threshold value χ_ref = 0.15. In order to enhance the resolution of strongly emitting shocked regions, we prevent levels higher than two from being created if the temperature does not exceed 5 × 10³ × ℓ, where ℓ is the level number.

The results are shown, after ≈63.4 yr, in Figure 21 where the steepening of the perturbations leads to the formation of forward/reverse shock pairs. Radiative losses become strongly enhanced behind the shock fronts where temperature attains larger values and the compression is large. This is best illustrated in the close-up views of Figure 21 where we display density, temperature, hydrogen fraction, and magnetization for the second (upper panel) and third shocks (lower panel), respectively, located at z ≈ 14.2 and z ≈ 7.34 (in units of the jet radius). Owing to a much shorter cooling length, t_c ∼ p/Λ, a thin radiative layer forms the size of which, depending on local temperature and density values, becomes much smaller than the typical advection scale. An extremely narrow 1D cut shows, in Figure 22, the profiles of temperature, density, and ionization fractions across the central radiative shock, resolved at the largest refinement level. Immediately behind the shock wave, a flat transition region with constant density and temperature is captured on ∼6 grid points and has a very small thickness ∼0.005r_j. A radiatively cooled layer follows behind where temperature drops and the gas reaches the largest ionization degree. Once the gas cools below 10⁴ K, radiative losses become negligible and the gas is accumulated in a cold dense adiabatic layer. A proper resolution of the thin emission layers is thus crucial for correct and accurate predictions of emission lines and intensity ratios in stellar jet models (Teşileanu et al. 2011). Such a challenging computational problem can be tackled only by means of adaptive grid techniques.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our first test consists of an initial discontinuity located at x = 0.5 separating two regions of fluids characterized by

$$\begin{equation} (B_y,\, B_z, \, p)=\left\lbrace \begin{array}{@{}ll}(7,\, 7,\, 10^3) & \quad {\rm for}\quad x < 0.5 \\ (0.7,\, 0.7,\, 0.1) & \quad {\rm for}\quad x > 0.5 \, \end{array}\right. \end{equation} \tag{ 71 }$$

(see Balsara 2001b; Mignone & Bodo 2006; van der Holst et al. 2008, and references therein). The fluid is initially at rest with uniform density ρ = 1 and the longitudinal magnetic field is B_x = 10.

We solve the equations of relativistic MHD with the ideal gas law (Γ = 5/3) using the five-wave HLLD Riemann solver of Mignone et al. (2009) and C_a = 0.6. The computational domain [0, 1] is discretized using five levels of refinement with consecutive grid jump ratios of 4:2:2:2:2 starting from a base grid (level zero) of 400 grid zones (yielding an effective resolution of 25, 600 zones). Refinement is triggered upon the variable σ = (B²_y + B_z²)/(ργ) using Equation (50) with a threshold χ_r = 0.1. At t = 0.4 the resulting wave pattern (Figure 23) is comprised of two left-going rarefaction fans (fast and slow) and two right-going slow and fast shocks. The presence of magnetic fields makes the problem particularly challenging since the contact wave, slow, and fast shocks propagate very close to each other, resulting in a thin relativistic blast shell between x ≈ 0.88 and x ≈ 0.9.

Zoom In Zoom Out Reset image size

In Table 3, we compare, for different resolutions, the CPU timing obtained with the static grid version of the code versus the AMR implementation at the same effective resolution, starting from a base grid of 400 zones. In the unigrid computations, halving the mesh size implies approximately a factor of four in the total running time whereas it implies only a factor of two in the AMR approach. At the largest resolution employed here, the overall gain is approximately 50. This example confirms that the resolution of complex wave patterns in RMHD can largely benefit from the use of adaptively refined grids.

The generic Alfvén test (Giacomazzo & Rezzolla 2006; Mignone et al. 2009) consists of the following non-planar initial discontinuity:

$$\begin{equation} \left\lbrace \begin{array}{@{}lr}{\bf V}_L = \displaystyle \left(1, 0,0.3,0.4, 1, 6, 2, 5\right)^{\rm T} & \mathrm{for}\quad x_1 < 1/2 \,, \\ {\bf V}_R = \displaystyle \left(0.9, 0,0,0, 1, 5, 2, 5.3\right)^{\rm T} & \mathrm{for}\quad x_1 > 1/2 \,, \end{array}\right. \end{equation} \tag{ 72 }$$

where, as in Section 5.1.1, V = (ρ, v₁, v₂, v₃, B₁, B₂, B₃, p) is given in the frame of reference aligned with the direction of motion x₁. The ideal EoS (Equation (10)) with Γ = 5/3 is adopted.

Here, we consider a 2D version by rotating the discontinuity front by π/4 with respect to the mesh and following the evolution until t = 0.4cos π/4. The test is run on a coarse grid of 64 × 4 zones covering the domain x ∈ [0, 1], y ∈ [ − 1/32, 1/32] with six additional levels of refinement resulting in an effective resolution of 4096 × 256 zones (a factor of two in resolution is used between levels). In order to trigger refinement across jumps of a different nature, we define the quantity σ  =  B²/(ργ) together with Equation (50) and a threshold value χ_r = 0.03. The Courant number is C_a = 0.6 and the HLLD Riemann solver is used throughout the computation. Boundary conditions assume zero gradients at the x-boundaries while translational invariance is imposed at the top and bottom boundaries where, for any flow quantity q, we set q(i, j) = q(i ± 1, j∓1).

The breaking of the discontinuity, shown in Figure 24, leads to the formation of seven waves including a left-going fast rarefaction, a left-going Alfvén discontinuity, a left-going slow shock, a tangential discontinuity, a right-going slow shock, a right-going Alfvén discontinuity, and a right-going fast shock. Our results indicate that refined regions are created across discontinuous fronts which are correctly followed and resolved, in excellent agreement with the reference solution (obtained on a fine 1D grid). Note also that the longitudinal component of the magnetic field, B₁, does not present spurious jumps (Figure 24) and shows very small departures from the expected constant value. In the central region, for 0.4 ≲ x ≲ 0.6, rotational discontinuities and slow shocks are moving slowly and very adjacent to each other, thus demanding relatively high resolution to be correctly captured. With the given prescription, grid adaptation efficiently provides adequate resolution where needed. This is best shown in the close-up of the central region, Figure 25, showing σ together with the AMR block structure.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A comparison of the performance between static and adaptive grid computations is reported in Table 4 using different resolutions and mesh refinements. At the resolution employed here, the gain is approximately a factor of eight.

The 2D rotor problem (Del Zanna et al. 2003; van der Holst et al. 2008; Dumbser et al. 2008) consists of a rapidly spinning disk embedded in a uniform background medium (ρ = p = 1) threaded by a constant magnetic field B = (1, 0, 0). Inside the disk, centered at the origin with radius R = 0.1, the density is ρ = 10 and the velocity is prescribed as v = ω(− y, x, 0), where ω = 9.95 is the angular frequency of rotation. The computational domain is the square x, y ∈ [ − 1/2, 1/2] with outflow (i.e., zero gradient) boundary conditions. The ideal EoS with Γ = 5/3 is used. Computations are performed on a base grid with 64² grid points and six levels of refinement using the HLL Riemann solver and a CFL number C_a = 0.5. Zones are tagged for refinement whenever the second derivative error norm of ${\cal E}/(\rho \gamma)$, computed with Equation (50), exceeds χ_r = 0.15. The effective resolution amounts therefore to 4096² zones, the largest employed so far (to the extent of our knowledge) for this particular problem.

Results are shown at t = 0.4 in Figure 26 where one can see an emerging flow structure enclosed by a circular fast forward shock (r  ≈  0.425) traveling into the surrounding medium. An inward fast shock bounds the innermost oval region where density has been depleted to lower values. The presence of the magnetic field slows down the rotor, and the maximum Lorentz factor decreases from 10 to ≈2.2. The numerical solution preserves the initial point symmetry and density corrugations, present in lower resolution runs, seem drastically reduced at this resolution, in accordance with the findings of van der Holst et al. (2008). Divergence errors are mostly generated in proximity of the outer fast shock and are damped while being transported out of the domain at the speed of light in the GLM formulation. We monitor the growth of such errors by plotting the volume average of |∇ · B| as a function of time, showing (bottom panel in Figure 26) that the growth of monopole errors is limited to the same values of van der Holst et al. (2008; the factor of four comes from the fact that our computational domain is two times smaller).

Zoom In Zoom Out Reset image size

Parallel performance for this particular test is shown in Figure 27 plotting the speed-up, defined as $S=8T^{\rm u}_{8}/T_{N_{\rm CPU}}$, where T^u₈ is the execution time of the uniform grid run performed at the same effective resolution (4096²) on eight processors while $T_{N_{\rm CPU}}$ is the running time measured with N_CPU processors. AMR scaling has been quantified for computations using maximum patch sizes of 16 and 32 grid points resulting, respectively, in 11, 921 and 6735 total number of blocks at the end of integration. Approximately half of it belongs to the finest level, as reported in Figure 27. In general, the AMR approach offers a 5–10 speed-up gain in terms of execution time over the uniform, static grid approach. For N_CPU < 256 the parallel efficiency, measured as S/N_CPU, is larger than 0.7 while it tends to be reduced for more processors. Computations carried with fewer blocks (i.e., larger block sizes) tend to be more efficient when fewer CPUs are in use since inter-processor communications are reduced. However, this cost obviously becomes larger at 512 (or more) processors resulting in a loss of efficiency.

Zoom In Zoom Out Reset image size

Strong symmetric explosions in highly magnetized environments can become rather arduous benchmarks probing the code's ability to evolve strongly magnetized shocks (see, for instance, Komissarov 1999; Leismann et al. 2005; Mignone & Bodo 2006; Del Zanna et al. 2003; Beckwith & Stone 2011). Failures may lead to unphysical densities or pressures whenever the scheme does not introduce adequate dissipation across oblique discontinuities and/or if the divergence-free condition is not properly controlled.

The 2D setup considered here consists of the Cartesian square [ − 6, 6] × [ − 6, 6] initially filled with constant density and pressure values ρ = 10⁻⁴ and p = 5 × 10⁻³ and threaded by a constant horizontal magnetic field with strength B₀ = 1. A cylinder with radius r = 0.8 and centered at the origin delimits a higher density and pressure region where ρ = 10⁻², p = 1. The ideal EoS with Γ = 4/3 is used. Our choice of parameters results in a low β plasma (2p/B² = 10⁻²) and corresponds to the strongly magnetized cylindrical explosion discussed by Beckwith & Stone (2011). We point out that this configuration differs from the one discussed in Komissarov (1999) and Mignone & Bodo (2006) by having the ambient pressure be 10 times larger. This allows to run the computation without having to introduce any specific change or modification in the algorithm, such as shock flattening or redefinition of the total energy. We adopt a base grid of 48 zones in each direction with outflow boundary conditions holding on all sides. Integration proceeds until t = 4 using the HLLC Riemann solver (Mignone & Bodo 2006) using five levels of refinement, reaching an effective resolution of 1536² zones. Total energy triggers refinement with a threshold of χ_ref = 0.025.

The overpressurized region, Figure 28, sets an anisotropic blast wave delimited by a weak fast forward shock propagating almost radially. The explosion is strongly confined along the x-direction by the magnetic field where plasma is accelerated to γ_max ≈ 1.76. The inner structure is delimited by an elongated structure enclosed by a slow shock adjacent to a contact discontinuity. The two fronts blend together as the propagation becomes perpendicular to the field line. Since the problem is purely 2D, rotational discontinuties are absent. The block distribution, shown in the lower panel of Figure 28, shows that refinement clusters around the outer discontinuous wave as well as the multiple fronts delimiting the extended inner region. The total number of blocks is 2258 with the finest level giving a relative contribution of ∼0.55 and a corresponding volume filling factor of ∼0.12. The numerical solution retains the expected degree of symmetry and the agreement with earlier results demonstrate that the code can robustly handle strongly magnetized relativistic environments within the GLM-MHD formalism. Results obtained with four processors show that calculation done with AMR is faster than the static grid runs by a factor of ∼4.5.

Zoom In Zoom Out Reset image size

In the next example, we consider a 3D extension of the cylindrical blast wave problem presented in Del Zanna et al. (2007). The initial condition consists of a sphere with radius r = 0.8 centered at the origin filled with hot gas having ρ = 10⁻², p = 1 (the ideal EoS with Γ = 4/3 is used). Outside this region, density and pressure decrease linearly reaching the constant ambient values ρ = 10⁻⁴, p = 5 × 10⁻³ for r ⩾ 1. The plasma is at rest everywhere and a constant magnetic field is set oblique to the mesh, i.e., ${\bf B}=B_0(\hat{{\bf e}}_x + \hat{{\bf e}}_y)/\sqrt{2}$ with B₀ = 0.1. The computational domain is the Cartesian box [ − 6, 6]³ covered by a base grid of 40 zones in each direction with outflow boundary conditions holding on all sides. Using the HLLD Riemann solver, we follow integration until t = 4 using four levels of refinement (effective resolution 640³) and a CFL number C_a = 0.4. Zones are tagged for refinement upon density fluctuations with a threshold value of χ_ref = 0.15.

Results shown in Figure 29 reveal an oval-shaped explosion enclosed by an outer fast shock propagating almost radially. A strong reverse shock confines a prolate spheroidal region where the magnetic field has been drained and inside which expansion takes place radially. The largest Lorentz factor ∼6.3 is attained in the direction of magnetic field lines close to the ellipsoid poles.

Zoom In Zoom Out Reset image size

In order to measure parallel performance on a fixed number of blocks, we have integrated the solution array starting from t = 3.5 for a fixed number of steps by enforcing zero interface flux. Results are shown in Figure 30 where we compare two sets of "frozen-flux" computations with maximum block sizes of 20 and 40 mesh points corresponding, respectively, to 27,661 (22,598 on the finest level) and 11,229 (9462 on the finest level) total number of blocks. In the former case, the larger number of patches allows us to reach a better efficiency (more than 0.75) whereas the latter case performs worse for increasing number of processors. Even if the finest level has a number of blocks still larger than the maximum number of CPUs employed, many boxes are much smaller than the maximum block size allowed. The ideal workload per processor (i.e., the number of cells of the finest level divided by the number of processors) is comparable to the maximum block size in the 1024 CPUs run and becomes smaller in the 2048 CPUs run. In this latter case, the workload of the processors integrating the larger blocks is therefore larger than the ideal one and most of the CPUs must wait for these processors to complete the integration. This could explain the poor parallel performance of this test problem.

Zoom In Zoom Out Reset image size

In the next example (Bucciantini & Del Zanna 2006; Mignone et al. 2009), we consider a 2D planar domain with x ∈ [0, 1], y ∈ [ − 1, 1] initially filled with a (relativistically) hot gas with constant density ρ = 1 and pressure p = 20. An initially perturbed shear velocity profile of the form

$$\begin{equation} (v_x, v_y, v_z) = \left[V_0\tanh \frac{y}{\alpha }, \epsilon \sin (2\pi x)\exp \left(-\frac{y^2}{\beta ^2}\right), 0\right] \end{equation} \tag{ 73 }$$

is set across the domain, where V₀ = 1/4 is the nominal flow velocity, = V₀/100 is the amplitude of the perturbation while α = 0.01, β = 0.1. The TM EoS, Equation (11), is used during the evolution to recover the gas enthalpy from density and pressure. The magnetic field is initially uniform and prescribed in terms of the components parallel and perpendicular to the plane of integration:

$$\begin{equation} (B_x, B_y, B_z) = (\sqrt{2\sigma _{\rm pol}p},0, \sqrt{2\sigma _{\rm tor}p}), \end{equation} \tag{ 74 }$$

where σ_pol = 0.01, σ_tor = 1. We perform the integration on a base grid of 32 × 64 cells with six additional levels of refinement, reaching an effective resolution of 2048 × 4096 zones. The Courant number is set to C_a = 0.8 and refinement is triggered whenever Equation (50) computed with the total energy density exceeds χ_r = 0.14. Open boundary conditions hold at the lower and upper y-boundaries while periodicity is imposed in the x-direction.

Figure 31 shows density and magnetic field distributions at t = 5, which approximately marks the transition from the linear phase to the nonlinear evolution (see the discussion in Mignone et al. 2009). Both density and magnetic field are organized into filaments that wrap around, forming a number of vortices arranged symmetrically with respect to the central point. Refined levels concentrate mainly around the location of the fluid interface and accurately follow the formation of the central vortex and the more elongated adjacent ones. The computation preserves the initial point symmetry even if the grid generation algorithm does not necessarily do so. The lower resolution outside this region contributes to damp acoustic waves as they travel toward the outer boundaries and eventually reduces the amount of spurious reflections.

Zoom In Zoom Out Reset image size

In Figure 32, we compare the parallel speed-up between a number of adaptive grid calculations using either six levels (grid jump of 2) or three levels (grid jump of 4) and the equivalent uniform grid run with 2048 × 4096 zones. Generally, the AMR approach is ≈40 faster than the fixed grid run. In addition, the three-level computation seems to perform somewhat better than the six-level run although the number of blocks on the finest levels is essentially the same at the end of computation (1085 and 1021, respectively). Parallel scaling sensibly deteriorates when the number of blocks per processor becomes less than three or four, in accordance with the results of previous tests.

Zoom In Zoom Out Reset image size

In the next example, we consider a 3D extension of the planar relativistic shock–cloud interaction originally presented in Mignone & Bodo (2006). The initial condition consists of a high-density (ρ = 10) spherical clump (radius 0.15) centered at (0.8, 0, 0) adjacent to a shock wave located at x = 0 with downstream and upstream values given by (ρ, v_x, B_z, p) = (39.5052, 0, 1.9753, 129.72386) for x < 0 and $(\rho ,v_x,B_z,p)=(1,\, -\sqrt{0.99},\, 0.5,\, 10^{-3})$ for x ⩾ 0. The remaining quantities are set to zero. The computational domain is the box x ∈ [0, 2], y, z ∈ [ − 1/2, 1/2]. At the base level we fix the resolution to 64  ×  32  ×  32 zones and employ four levels of refinement equivalent to an effective resolution of 1024  ×  512  ×  512. The refinement criterion (Equation (50)) uses the conserved density for zone tagging with a threshold χ_r = 0.2. The equations of relativistic MHD are solved using the TM EoS (Mignone & McKinney 2007), the HLLD Riemann solver, and linear reconstruction with the harmonic limiter (Equation (41)). Shock-adaptive hybrid integration provides additional numerical dissipation in proximity of strong shock waves by switching the integration scheme to the HLL Riemann solver and the MinMod limiter (Equation (35)), according to the strategy described in Appendix B. For efficiency purposes, we take advantage of the symmetry across the xy- and xz-planes to reduce the computation to one quadrant only.

Figure 33 shows a 3D rendering of density and magnetic pressure through a cross-sectional view of the xy- and xz-planes at t = 1. The collision generates a fast, forward bow shock propagating ahead and a backward reverse shock transmitted back into the cloud. The cloud becomes gradually wrapped by the incident shock into a mushroom-shaped shell reaching a large compression (ρ_max ≈ 121, B²/2 ≈ 157). Grid refinement concentrates mainly around the incident and forward shocks and the tangential discontinuities bordering the edge of the cloud.

Zoom In Zoom Out Reset image size

As a final application, we discuss the 3D propagation of a relativistic magnetized jet carrying an initially purely azimuthal magnetic field. Indeed, the presence of a substantial toroidal component of the field is commonly invoked and held responsible for the acceleration and collimation of jets from active galactic nuclei. From an astrophysical point of view, relativistic jets appear to be quite stable although cylindrical MHD configurations are expected to be unstable to reflection, Kelvin–Helmholtz and current driven modes thus posing an unsolved issue, see Mignone et al. (2010), and reference therein.

Following Mignone et al. (2009), we set the box x, y ∈ [ − 12.5, 12.5], z ∈ [0, 50] as our computational domain, initially filled with constant uniform density and pressure ρ_a and p_a. The beam is injected at the lower boundary for z ⩽ 0, $r=\sqrt{x^2+y^2} \le 1$ with a density ρ_j, a longitudinal velocity corresponding to a Lorentz factor γ_j, and an azimuthal magnetic field given by

$$\begin{equation} B_\phi = \left\lbrace \begin{array}{@{}ll}\gamma _jb_m r/a &\quad {\rm for}\quad r < a {,}\\ \gamma _jb_m a/r &\quad {\rm for}\quad a < r < 1 {,}\\ 0 &\quad {\rm otherwise}{,} \end{array}\right. \end{equation} \tag{ 75 }$$

where b_m is the magnetic field strength in the fluid's frame and a = 0.5 is the magnetization radius. The pressure profile is found by imposing radial momentum balance across the beam giving

$$\begin{equation} p(r) = p_a + b_m^2\left[1 - \min \left(\frac{r^2}{a^2},1\right)\right]\,, \end{equation} \tag{ 76 }$$

where p_a is the ambient pressure defined in terms of the sonic Mach number M:

$$\begin{equation} p_a = \frac{\rho _jv_j^2(\Gamma -1)}{\Gamma (\Gamma -1)M^2 - \Gamma v_j^2}\,, \end{equation} \tag{ 77 }$$

with Γ = 5/3 although the TM EoS (Equation (11)) is used during the evolution. The actual value of b_m may be found by prescribing the magnetization parameter σ_ϕ defined as the ratio between beam-averaged magnetic energy density and thermal pressure. The final result yields (Mignone et al. 2009)

$$\begin{equation} b_m = \sqrt{\frac{-4p_a\sigma _\phi }{a^2(2\sigma _\phi -1+4\log a)}}\,. \end{equation} \tag{ 78 }$$

For the present simulation, we adopt γ_j = 7, M = 3, σ_ϕ = 1.3, ρ_j = 1, and ρ_a = 10³, thus corresponding to a magnetically dominated, underdense relativistic jet.

We follow the simulation up to t = 130 (in units of the speed of light crossing of the jet radius), starting with a base grid of 32 × 32 × 64 with four levels of refinement. The refinement ratio between consecutive levels is 2 so that, at the finest resolution, we have ≈20 points per beam radius. Zones belonging to levels zero, one, and two are tagged for refinement using the normalized second derivative Equation (50) of the total energy density with a threshold χ_r = 0.25. Zones at the finest level grid are instead created by using B²/(γρ) in place of the total energy density. The rationale for choosing this selective rule is to provide the largest resolution on the jet material only, still being able to track the sideway expansion of the cocoon at lower resolution. We use the HLLD Riemann solver with the harmonic limiter (Equation (41)) except at strong shocks where we employ the multidimensional shock dissipation switch outlined in Appendix B. The Courant number is C_a = 0.3.

The jet structure, Figure 34, shows that the flow maintains a highly relativistic central spine running through multiple recollimation shocks where jet pinching occurs. At the jet head the beam decelerates to sub-relativistic velocities favoring the formation of a strongly magnetized termination shock where B² ≈ 26γ²ρ. Here, the presence of a predominant toroidal magnetic field induces CD kink instabilities which are responsible for the observed jet wiggling and symmetry breaking (see Mignone et al. 2010). This test shows that PLUTO-CHOMBO can be effectively used in simulations involving relativistic velocities, highly supersonic flows, and strongly magnetized environments.

Zoom In Zoom Out Reset image size

At t = 130 the finest level grid has a volume filling factor of ≈11% with 10,761 blocks while level three occupies ≈41% of the total volume with 5311 blocks. For this particular application, the benefits offered by grid adaptivity are more evident at the beginning of the computation when most of the computational zones in the domain are unrefined. In this respect, we have found that the CPU time increases (for t ⩽ 130) with the number of steps n approximately as a + bn + cn² + dn³, where a ≈ 1.46, b ≈ 0.079, c ≈ 5.03 × 10⁻⁴, and d ≈ 2.38 × 10⁻⁶. This suggests that, as long as the filling factor remains well below 1, the AMR calculation with the proposed refinement criterion should be at least ∼2.5 times faster than an equivalent computation using static mesh refinement and considerably larger if compared to a uniform grid run at the effective resolution.