ENZO: AN ADAPTIVE MESH REFINEMENT CODE FOR ASTROPHYSICS

The Eulerian equations of ideal magnetohydrodynamics (MHD) including gravity, in a coordinate systems comoving with the cosmological expansion, are given by

$$\begin{eqnarray} \frac{\partial \rho }{\partial t} + \frac{1}{a} \boldsymbol {\nabla } \cdot (\rho {\bf v}) & = & 0, \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} \frac{\partial \rho {\bf v}}{\partial t} + \frac{1}{a} \boldsymbol {\nabla } \cdot \left(\rho {\bf v}{\bf v}+ {\bf I}{p}^* - \frac{{\bf B}{\bf B}}{a} \right) & = & - \frac{\dot{a}}{a} \rho {\bf v}- \frac{1}{a} \rho \boldsymbol {\nabla }\phi, \quad \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} && \frac{\partial E}{\partial t} + \frac{1}{a} \boldsymbol {\nabla } \cdot \left[ (E + p^*) {\bf v}- \frac{1}{a} {\bf B}({\bf B}\cdot {\bf v}) \right] = - \frac{\dot{a}}{a} \left(2E - \frac{B^2}{2a} \right) \nonumber\\ && \quad\quad - \frac{\rho }{a} {\bf v}\cdot {\boldsymbol \nabla }{\phi} - { \Lambda} + {\Gamma} + \frac{1}{{a}^2} \boldsymbol {\nabla } \cdot {\bf F}_{\rm cond}, \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} \frac{\partial {\bf B}}{\partial t} - \frac{1}{a} \boldsymbol {\nabla }\times ({\bf v}\times {\bf B}) & = & 0. \end{eqnarray} \tag{ 4 }$$

In these equations, E, ρ, v, and B are the comoving total fluid energy density, comoving gas density, peculiar velocity, and comoving magnetic field strength, respectively. The matrix I is the identity matrix, and a is the cosmological expansion parameter (discussed in more detail below). The first equation represents conservation of mass, the second conservation of momentum, and the third conservation of total (kinetic plus thermal plus magnetic) fluid energy. They are respectively, the first, second, and third moments of the Boltzmann equation. The fourth equation is the magnetic induction equation. Terms representing radiative cooling (Λ) and heating (Γ) enter on the right-hand side of the energy equation (3), as does the flux due to thermal heat conduction (F_cond).

The comoving total fluid energy density E is given by

$$\begin{equation} E = e + \frac{\rho v^2}{2} + \frac{B^2}{2a}, \end{equation} \tag{ 5 }$$

where e is the comoving thermal energy density. The total comoving isotropic pressure p* is given by

$$\begin{equation} p^* = p + \frac{B^2}{2a}, \end{equation} \tag{ 6 }$$

and the quantity p is the thermal pressure. We use units such that the magnetic permeability is unity (μ₀ = 1). The equations are closed by an equation of state and Poisson's equation for the gravitational potential ϕ:

$$\begin{eqnarray} e &=& \frac{p}{(\gamma - 1)}, \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} \nabla ^2 \phi &=& \frac{4 \pi G}{a} \left(\rho _{\rm total} - \rho _0 \right). \end{eqnarray} \tag{ 8 }$$

The equation of state is shown here for an ideal gas with a ratio of specific heats γ. The gravitational potential ϕ is sourced by the total mass density contrast, where ρ_total = ρ_gas + ρ_dm + ρ_stars and ρ₀ is the mean density. Although we write the equations including both magnetic field terms and comoving coordinates, the code is frequently used both in the purely hydrodynamic limit (B = 0, referred to as HD below) and without cosmological expansion (a = 1, $\dot{a} = 0$).

For completeness, we note here that we have defined several key quantities in the comoving frame to make the previous equations more readable. Specifically, we define:

$$\begin{eqnarray} {\bf x}& = & {\bf x}^{\prime } / {a}, \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray} {\bf v}& = & a \dot{{\bf x}} = {\bf v}^\prime - \dot{a}{\bf x}, \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} \rho & = & a^3 \rho ^\prime, \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} p & = & a^3 p^\prime, \end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray} E & = & a^3 \left(E^\prime - \dot{a} \, {\bf x}\cdot {\bf v}^\prime - \frac{1}{2} \dot{a}^2 {\bf x}^2\right), \end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray} \phi & = & \phi ^\prime + \frac{1}{2} \, a \ddot{a} \, {\bf x}^2, \end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray} {\bf B}& = & a^2 {\bf B}^\prime, \end{eqnarray} \tag{ 15 }$$

where primes indicate quantities in proper coordinates (fixed frame). We note that the definition of comoving B used here is such that a uniform field is constant in a homogeneous expanding universe and is also the quantity used in the MHD-CT solver described below; however, it is not universal, and is slightly different from that used both in Li et al. (2008) and in the MHD-Dedner solver (see Section 4.2 for more details). The expansion parameter a ≡ 1/(1 + z) follows the expansion of a smooth, homogeneous background, where z, the redshift, is a function only of t. All spatial derivatives are determined with respect to the comoving position x, which removes the universal expansion from the coordinate system. The evolution of a(t) is governed by the second Friedmann equation for the expansion of a spatially homogeneous and isotropic universe

$$\begin{eqnarray} \frac{\ddot{a}}{a} & = & - \frac{4 \pi G }{3 a^3 } (\rho _0 + 3p_0/c^2) + \Lambda _c c^2 / 3. \end{eqnarray} \tag{ 16 }$$

Here ρ₀ is the mean comoving mass density (including both baryonic and dark matter), p₀ is the comoving background pressure contribution, and Λ_c is the cosmological constant. This system of equations is limited to the non-relativistic regime and assumes that curvature effects are not important—both assumptions are reasonable as long as the size of the simulated region is small compared to the radius of curvature and the Hubble length c/H, where c is the speed of light and $H = \dot{a}/a$ is the Hubble constant.

The comoving evolution equations are equivalent to the fixed coordinate version in the non-cosmological limit (a = 1, $\dot{a} = 0$). We include the cosmological terms with the understanding that Enzo is not restricted to cosmological applications—while historically the code was written with cosmological situations in mind (e.g., galaxy clusters), in more recent years it has been used to simulate a much broader range of astrophysical environments.

Any collisionless components (such as dark matter and stars) are modeled by N-body particles, whose dynamics are governed by Newton's equations in comoving coordinates:

$$\begin{eqnarray} \frac{d {\bf x}}{dt} & = & \frac{1}{a} {\bf v}, \end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray} \frac{d {\bf v}}{dt} & = & - \frac{\dot{a}}{a} {\bf v}- \frac{1}{a} \boldsymbol {\nabla }\phi. \end{eqnarray} \tag{ 18 }$$

The particles contribute to the gravitational potential through Poisson's equation (Equation (8)).

In addition, Enzo can solve the mass conservation equations for a set of chemical species and their reactions. For any species i with comoving number density n_i, these equations have the form:

$$\begin{equation} \frac{\partial n_i}{\partial t} + \frac{1}{a} \boldsymbol {\nabla } \cdot (n_i {\bf v}) = { \sum _j k_{ij}(T) \, n_i n_j + \sum _j \Gamma ^{\rm ph}_j n_j } \end{equation} \tag{ 19 }$$

where k_ij are the rate coefficients for the two-body reactions and are usually functions of only temperature (we will specifically note the cases where we either include three-body reactions or have density-dependent rates). The $\Gamma ^{\rm ph}_j$ are destruction/creation rates due to photoionizations and/or photodissociations. Currently the species Enzo can follow include H, H⁺, He, He⁺, He⁺⁺, and e⁻, and through additional options also H⁻, H₂, H$_2^+$, and HD, D, and D⁺. Lastly, Enzo can also track the advection of one or more comoving metal density fields, which can contribute to the radiative cooling and star formation processes.

The code can include either a homogeneous radiation background or evolve an inhomogeneous radiation field either by directly solving the radiative transfer equation along rays or by solving a set of moment equations derived from the radiative transfer equation. In comoving coordinates (e.g., Gnedin & Ostriker 1997) the radiative transfer equation reads

$$\begin{equation} \frac{1}{c} \; \frac{\partial I_\nu }{\partial t} + \frac{a_{\rm em}}{a} \, \hat{n} \cdot \boldsymbol {\nabla }I_\nu - \frac{H}{c} \; \left(\nu \frac{\partial I_\nu }{\partial \nu } - 3 I_\nu \right) = -\kappa _\nu I_\nu + j_\nu. \end{equation} \tag{ 20 }$$

Here I_ν ≡ I(ν, x, Ω, t) is specific intensity of the radiation, with dimensions of energy per time t per solid angle Ω per unit area per frequency ν. The second term represents the propagation of radiation, where the factor a_em/a accounts for cosmic expansion since the time of emission. The third term describes both the cosmological redshift and dilution of radiation. On the right hand side, the first term captures absorption (κ_ν ≡ κ_ν(x, ν, t)), and the second term emission (j_ν ≡ j_ν(x, ν, t)) from any sources of radiation (whether point or diffuse).

As an alternative to the ray-casting strategy, Enzo can also solve the angle-averaged radiative transfer equation through a flux-limited diffusion approximation, with couplings to both the gas energy and chemical number densities. The equations solved are

$$\begin{eqnarray} \frac{\partial E_r}{\partial t} + \frac{1}{a} \boldsymbol {\nabla } \cdot \left(E_r {\bf v}\right) &=&\frac{1}{a^2} {\boldsymbol \nabla } \cdot \left(D\boldsymbol {\nabla }E_r\right) - \frac{\dot{a}}{a} E_r - c\kappa E_r + \eta, \nonumber\\ \end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray} \frac{\partial e_c}{\partial t} &=& -\frac{2\dot{a}}{a} e_c + \Gamma - \Lambda, \end{eqnarray} \tag{ 22 }$$

where E_r = E_r(x, t) is a gray radiation energy density and e_c is the internal energy correction due to photo-heating and chemical cooling. Here, we define E_r through first assuming a fixed frequency spectrum, i.e., $E_{\nu }(\nu,{\bf x},{t}) = \tilde{E}_r({\bf x},{t}) \chi (\nu)$, and then defining the integrated quantity

$$\begin{equation} E_r({\bf x},{t}) \equiv \int _{\nu _0}^{\infty } {E}_{\nu }(\nu,{\bf x},{t})\,{d}\nu \ = \ \tilde{E}_{r}({\bf x},{t}) \int _{\nu _0}^{\infty } \chi (\nu)\,{d\nu}. \end{equation} \tag{ 23 }$$

The quantity D in Equation (21) is the Larsen square-root flux-limiter (see Morel 2000), κ is the total opacity, η is the field of radiation sources, Γ is the radiation-induced photo-heating rate, and Λ is the chemical cooling rate.

Finally, Enzo implements the equations of isotropic heat conduction in a manner similar to that of Parrish & Stone (2007), where the isotropic heat flux is given by

$$\begin{equation} {\bf F}_{\rm cond}= -\kappa _{\rm cond} \boldsymbol {\nabla }{T}. \end{equation} \tag{ 24 }$$

Here κ_cond = f_sp κ_sp is the heat conduction coefficient, given as a fraction f_sp of the Spitzer conductivity κ_sp (Spitzer 1962), and T is the gas temperature (with fluids explicitly assumed to be single-temperature). Saturation of the heat flux in high temperature, low density regimes (such as the intracluster medium in galaxy clusters) is taken into account.

Thermal conduction in a plasma can be strongly affected by the presence of magnetic field lines, which may strongly suppress heat flow perpendicular to the magnetic field. To include this effect, it is possible to allow heat transport only along (not across) magnetic field lines, as follows:

$$\begin{equation} {\bf F}_{\rm cond}= -\kappa _{\rm cond} \,{\bf b}({\bf b} \cdot \boldsymbol {\nabla }{T}), \end{equation} \tag{ 25 }$$

where b is the unit vector denoting the direction of the magnetic field.

In this section we briefly describe the numerical methods that are used to solve the equations outlined in Section 2.1. We proceed through the numerical methods in the same order as will be used in Sections 3 through 10 so that there is a one-to-one correspondence between each of the following overviews and the complete description provided in the later sections. The goal in this section is to introduce the reader to the basic principles of the methods without drowning them in detail.

2.2.1. Structured Adaptive Mesh Refinement

2.2.2. (Magneto)-hydrodynamic Solver Methods

2.2.3. Gravity

2.2.4. N-body Dynamics

2.2.5. Chemistry

2.2.6. Radiative Cooling and Heating

2.2.7. Homogeneous Radiation Backgrounds

2.2.8. Radiation Transport: Ray Tracing

2.2.9. Radiation Transport: Flux-limited Diffusion

2.2.10. Heat Conduction

2.2.11. Star Formation and Feedback

2.2.12. Timestep Constraints

In order to solve the hydrodynamic equations on a patch, it is necessary to specify the boundary conditions on that patch. There are two cases. The first is that the boundary is external to the computational domain. This occurs for all the boundaries of the root grid. The second possibility is that the boundary is within the computational domain, which is true for most, but not all, subgrids.

In the first case, we allow for four types of boundary conditions (for simplicity, represented here as being at x = 0, with positive x interior and negative x exterior to the domain):

• 1.  
  Reflecting. In this case, the boundary behaves like a mirror, with the solution on the boundary reflecting the solution in the computational domain (q(− x) = q(+ x)) but with the velocity normal to the boundary reversed: v_x(− x) = −v_x(+ x).
• 2.  
  Outflow. We approximate outflow boundary conditions by duplicating the solution at the edge of the computational domain: q(− x) = q(0).
• 3.  
  Inflow. In this case, the boundary values are fixed by a pre-determined function (q(− x) = q0(− x, t)). The code provides a simple way to set inflow values that are constant in both time and space over an entire face, or hooks are provided for special cases where the boundary must be time-varying.
• 4.  
  Periodic. The boundary solution is obtained from the other side of the grid: q(− x) = q(x_max − x).

In the second case, the value must be either (1) interpolated from a parent or (2) copied from a sibling. In practice, we solve this problem by first interpolating all boundaries from the grid's parent and then performing the sibling copy wherever possible, since these values will always be as good as, or superior to, interpolation.

We focus first on conserved, cell-centered quantities. The interpolation function should be conservative in the sense that the mean over the interpolated quantity, in a given cell, must return the cell-centered value. It should also be monotonic, so that new extrema are not introduced into a smoothly varying field.

In one dimension, such functions are well known; examples include donor cell, Van Leer, and piecewise parabolic interpolation, which are first-, second-, and third-order accurate respectively (e.g., Stone & Norman 1992a). However, we require fully multi-dimensional interpolation functions, and so focus on second-order schemes since they are tractable and relatively accurate.²⁵ In the Appendix, we explicitly write down functions in one, two and three dimensions, but here we discuss the general design considerations.

In two dimensions, it is relatively easy to build a linear, monotonic, conservative function by requiring, first, that the function at the center of the cell return the cell-centered quantity (and since it is linear, this implies that an average over the cell will also be equal to this value), and second, by treating the cell-diagonals as one-dimensional Van Leer-type functions. The reason to focus on the cell diagonals is that for a linear function, the extrema will occur at the cell corners, which are in turn linked by the diagonals. This process can be done in two dimensions since there are three constraints (one at the center and one along each diagonal), and three parameters in a two-dimensional linear function.

In three dimensions this process breaks down since the number of constraints (center + four diagonals = five) exceeds the number of parameters (four). One way to solve this problem would be to perform a least-squares fit under the given constraints; however, this is too computationally expensive. Therefore, in the Appendix we describe a heuristic approach that satisfies four constraints analytically and then adjusts the interpolation parameters in such a way as to satisfy the fifth. While the result is not necessarily optimal, it does appear to perform adequately in our tests.

Given the grid and its boundary conditions, the process of solving Equations (1)–(3) can be detached from the rest of the SAMR technique, so that it is possible to experiment with a number of different methods. The restrictions are that the method should be flux-conservative, or at least produce an accurate estimate of the flux of conserved quantities at cell interfaces. This is required in the flux-correction step (see Section 3.4). There is also the issue of centering, which refers to the position of physical quantities within a cell (the primary types are cell-centered, face-centered, and node-centered). The interpolation step, described in Section 3.1, must correctly interpolate with the selected centering.

Since a refined region is simulated with (at least) two different resolutions (one fine and one coarse), it has two or more "solutions" (i.e., values of density, velocity) that represent the same volume of space. This is one important aspect of the SAMR method—some physical regions are represented by grids at multiple resolutions. To maintain consistency, the coarse grid quantities are updated with the finer values, using a volume-weighted average of the conserved quantities:

$$\begin{equation} q^{\rm coarse} = r^{-d} \sum q^{\rm fine}_{i, j,k} \end{equation} \tag{ 26 }$$

where r is the refinement factor, and the sum is over all fine cells that "cover" the coarse cell. The exponent, d, refers to the weight of the average. For cell-centered quantities, this is the dimensionality of the problem. For face-centered quantities, this is one less than the dimensionality, and for edge-centered quantities this is one.

One of the significant advantages of flux-conservative methods is, of course, their ability to maintain conserved quantities to machine precision. Unfortunately, the projection step upsets this flux conservation since, at the boundary of a refined region, the coarse and fine cells on either side of the border are updated with different estimates of the flux across the boundary. At one refinement level this comes from the coarse grid computation, and in the other from the fine grid solution. Conservation can be restored by correcting all coarse cells that lie outside the boundary of a fine region with the difference between the coarse and fine estimates of the flux across the boundary:

$$\begin{equation} q^{\rm coarse} = \tilde{q}^{\rm coarse} - \Delta t \left(F^{\rm coarse} - \sum _{j,k} F^{\rm fine}_{j,k} \right) \end{equation} \tag{ 27 }$$

where F is the flux of quantity q_i in the i-direction through the cell interface. We use $\tilde{q}$ to indicate the uncorrected quantity on the coarse grid and the sum is over the r^{d − 1} fine cell interfaces in the perpendicular dimensions that share a face with the coarse cell being corrected. This correction is carried out for all coarse cells that share a face with a refined cell (and are not themselves covered with more refined cells).

Since the addition of more highly refined grids is adaptive, the conditions for refinement must be specified. Many different criteria can be simultaneously specified for a given simulation, and these refinement criteria are discussed in more detail in the following section (Section 3.6). Once all cells in a given grid that require additional refinement have been flagged, rectangular boundaries are determined that minimally encompass the flagged cells. This is done using an iterative machine-vision technique, as described in Berger & Rigoutsos (1992).

We first find the smallest region (or "proto-subgrid") that covers all of the refined cells on a grid. We compute an efficiency for this grid, defined as the ratio of flagged cells to total cells. If this efficiency is sufficiently large (typically over 30%), or if the grid is small enough, then we adopt the region as a new grid and move on to the next region (if any). If the region is not acceptable, then we split it into two using the following algorithm. First, we determine the longest dimension and compute the signature along that dimension. The signature is found by summing the number of refined cells with the same index along the longest dimension. Edges are typically indicated by inflection points in this signature, so we calculate the second derivative using a simple three-point finite difference estimator and look for the maximum in the second derivative, with ties going to the index closest to the center. If no inflection points are found along the longest axis, the next-longest axis is searched. If no inflection points are found, the longest axis is cut in half. The original region is then split into two new regions, or proto-subgrids. This process is then re-applied to each of the new regions.

This algorithm is not necessarily optimal, but experience shows it does a reasonable job of generating an efficient covering of the refined cells. Note that some cells are unnecessarily refined with the SAMR approach, resulting in wasted resources (unlike in cell-splitting AMR schemes). However, it has a number of advantages. First, traditional grid-based hydro methods can be simply added with only minor modifications. Second, it is easier to accommodate larger stencils (as in PPM, which requires three boundary zones), and third, it can be more efficient for some machines since there are more predictable and linear memory accesses, which more efficiently use a CPU's cache.

In this section, we turn to the refinement criteria themselves. Although generalized refinement criteria are available based on estimators for the error, such as the Richardson method (Atkinson & Han 2004), we have focused generally on methods that use physical estimators based on the assumption that the simulator has a specific quantity or quantities that indicate the need for refinement. The flagging methods currently built into Enzo are:

• 1.  
  Slope. If the normalized slope (q_{i + 1} − q_{i − 1})/(2q_i) of baryon quantities between adjacent cells is larger than a specified value (typically 0.3), then the cell is flagged.
• 2.  
  Baryon mass. This refinement criterion is designed to mimic a Lagrangian method in that it tries to keep a fixed mass resolution. In this method, a cell is refined if the baryonic mass in the cell (M_g = ρ (Δx)^d) is larger than a specified value:

  $$\begin{equation} M_g > \rho _{\rm flag}\, (\Delta x_{\rm root})^d r^{\epsilon _l l}, \end{equation} \tag{ 28 }$$

  where ρ_flag is the equivalent density on the root grid required for refinement, Δx_root is the root grid cell spacing, r is the refinement factor, l is the level and _l is a parameter that can be used to make the refinement more aggressive (_l < 0) or less aggressive (_l > 0) than pure Lagrangian refinement—in other words, super-Lagrangian or sub-Lagrangian refinement.
• 3.  
  Shocks. In this method, we identify and refine cells that contain shocks using the following criteria:

  $$\begin{eqnarray*} && (p_{i+1} - p_{i-1})/\min (p_{i+1}, p_{i-1}) > 0.33, \!\!\quad u_{i-1} - u_{i+1} > 0, \nonumber\\ && \quad {\rm and} \quad e_i / E_i > 0.1. \nonumber \end{eqnarray*}$$

  We note that the parameters 0.33 and 0.1 in the equations above were determined empirically, and can be modified by the user to refine more or less aggressively, as the situation warrants.
• 4.  
  Particle mass. This refinement criterion is precisely the same as the baryon mass criterion except that it uses the gridded (cloud-in-cell, or CIC) particle density. In general, this criteria is used in cosmological simulations, and thus refines on dark matter and stellar density.
• 5.  
  Jeans length. Following Truelove et al. (1998), the code can be forced to refine the Jeans length by a fixed number of cells, or, more specifically, whenever the following criterion is not met:

  $$\begin{equation} \Delta x < \left(\frac{\gamma \pi k_B T}{N_J^2 G \rho \mu m_{\rm H}} \right)^{1/2}, \end{equation} \tag{ 29 }$$

  where N_J is the required number of cells per Jeans length (four by default).
• 6.  
  Cooling time. Often in cosmological simulations, cooling is a rapid process and we do not necessarily want to be limited by the cooling time. However, in some cases it is desirable to resolve the cooling process, in which case we flag cells that have a cooling time shorter than their sound crossing time. The cooling time is given by t_cool = 1.5k_BT/(n Λ(T)) and the sound crossing time is t_cross = Δx/c_s, where c_s is the sound speed.
• 7.  
  Must-refine particles. In order to allow for refinement in moving regions, it is possible to generate "must-refine" particles that ensure local refinement of the 8 nearest cells to a specific level. The particles feel gravitational forces and their trajectories are followed as for any other collisionless particle. This feature can be used to ensure that a specific Lagrangian region stays refined to a specified level of resolution.
• 8.  
  Shear. Originally used in Kritsuk et al. (2006) and later modified into a more general form, turbulence calculations and Kelvin–Helmholtz instabilities can be successfully modeled with AMR if the shear is used as a refinement parameter; in particular, refinement occurs when the finite difference version of the following inequality holds:

  $$\begin{eqnarray} &&\left(\frac{\partial v_x}{\partial y} + \frac{\partial v_y}{\partial x} \right)^2 + \left(\frac{\partial v_z}{\partial y} + \frac{\partial v_y}{\partial z} \right)^2 \nonumber\\ &&\qquad+ \left(\frac{\partial v_x}{\partial z} + \frac{\partial v_z}{\partial x} \right)^2 > \left(\epsilon _s \frac{c_s}{\Delta x}\right)^2 \end{eqnarray} \tag{ 30 }$$

  where _s is a dimensionless user-defined parameter and Δx is the local cell size at the level where refinement is occurring.
• 9.  
  Optical depth. This refinement criterion, typically used for simulations using radiation transport, flags cells for refinement if they have an H i optical depth lower than unity. The optical depth of a cell is estimated to be:

  $$\begin{equation} \tau = \sigma _{{\rm H\,\scriptsize{I}}} n_{{\rm H\,\scriptsize{I}}} \Delta x \end{equation} \tag{ 31 }$$

  where $\sigma _{{\rm H\,\scriptsize{I}}}$ is the neutral atomic hydrogen absorption cross section, $n_{{\rm H\,\scriptsize{I}}}$ is the proper H i number density, and Δx is the cell width (which is assumed to be the same along each axis). If τ > 1, the cell is refined further, up to the maximum level of refinement.
• 10.  
  Resistive length. In a way similar to the Jeans length criterion described above, the code can be forced to locally resolve the resistive length scale by a minimum of a fixed number of cells. This is done by ensuring that the resistive length scale, defined as

  $$\begin{equation} L_R = |{B}| / |{\boldsymbol \nabla } \times {B}| \end{equation} \tag{ 32 }$$

  is always refined by at least N_R cells locally (i.e., we refine a given cell further if L_R ⩽ N_RΔx). Note that Enzo does not actually solve the resistive MHD equations—this is the "characteristic" local resistive length, rather than one calculated within the MHD solver.
• 11.  
  Must refine region. This refinement criterion is used to ensure that all cells in a given subvolume of the simulation are refined to at least some minimum level of refinement. The subvolume is user-specified, and may evolve over time (e.g., to follow a structure of interest that is traveling through the simulation volume).
• 12.  
  Metallicity. This refinement criterion ensures that all cells with a metallicity (defined as Z ≡ ρ_Z/ρ_b/0.022, where ρ_Z is the density of metals in the cell, ρ_b is the total baryon density (including metals), and 0.022 is the solar mass fraction of metals) above a user-defined value is refined to at least some minimum level of refinement.
• 13.  
  Second derivative. This is similar to the refinement criteria used in Paramesh and was originally developed by Loehner et al. (1987). This method calculates the second derivative of a user-defined field normalized by the first derivative, and sets a threshold for refinement. More background information can be found in the FLASH4 user guide,²⁶ but we reproduce the relevant equation here for completeness. The multi-dimensional version of the criteria takes in to account all cross-derivatives, and is expressed by,

  $$\begin{equation} S = \left\lbrace \displaystyle\frac{\sum _{ij}\left(\displaystyle\frac{\partial ^2f}{\partial x_i \partial x_j} \right)^2}{\sum _{ij}\left[ \displaystyle\frac{1}{2\Delta x_j}\left(\left|\displaystyle\frac{\partial f}{\partial x_i}\right|_{x_i+\Delta x_i} + \left|\displaystyle\frac{\partial f}{\partial x_i}\right|_{x_i-\Delta x_i} \right) + \epsilon \displaystyle\frac{\left|\bar{f_{ij}}\right|}{\Delta x_i \Delta x_j} \right]^2} \right\rbrace ^{1/2}, \end{equation} \tag{ 33 }$$

  where f is the user-defined field, $\left|\bar{f_{ij}}\right|$ is an average of f over cells in the i and j dimensions, and is a user-controlled value that provides a method of ignoring small fluctuations, set by default to 10⁻². If S is larger than a specified threshold (bounded by [0,1]) a cell is marked for refinement. This threshold is chosen at runtime and can vary for each of the fields under consideration.

Finally, it is also possible to specify a rectangular region that is statically refined to a certain level. These static refined regions are in place throughout the simulation regardless of the grid quantities and can be used for testing purposes, or for multiply-zoomed initial conditions (in cosmology simulations, for example) that would not otherwise be refined.

One (purely) hydrodynamic method used in Enzo is closely based on the PPM of Colella & Woodward (1984)—henceforth referred to as CW84—which has been modified for the study of cosmological fluid flows. CW84 describe two variants on this method: Lagrangian Remap and Direct Eulerian. We use the Direct Eulerian version in Enzo, which is better suited for AMR simulations. The Lagrangian Remap version has previously been adapted for cosmological use as described in Bryan et al. (1995), and the implementation we use here is closely based on that work.

The first term on the right hand side of Equations (2) and (3) comes from our choice of comoving coordinates (a similar term also appears in Equation (18), the velocity relation for the dark matter particles). A similar term does not appear in the mass conservation Equation (1) because of the comoving density definition. We note that these expansion terms could be eliminated entirely by the proper choice of variables (including time), although we have not done so here as they do not constitute a major source of error. We solve these terms using the same technique—they are split from the rest of the terms and solved using an (implicit) time-centered method, which is straightforward as there are no spatial gradients. Note that we use this method for all four of the hydro solvers.

The remainder of the terms in the fluid equations involve derivatives. Because we are interested in phenomena with no special geometry, we will restrict discussion to Cartesian coordinates. We can dimensionally split the equations and rewrite the one-dimensional (Eulerian) versions of Equations (1)–(3) without expansion terms in conservative form as,

$$\begin{eqnarray} \frac{\partial \rho }{\partial t} + \frac{1}{a} \frac{\partial \rho v }{\partial x} & = & 0 \end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray} \frac{\partial \rho v}{\partial t} + \frac{1}{a} \frac{\partial \rho v^2}{\partial x} + \frac{1}{a} \frac{\partial p}{ \partial x} & = & \rho \frac{g}{a} \end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray} \frac{\partial \rho E}{\partial t} + \frac{1}{a} \frac{\partial \rho v E}{\partial x} & = & \rho v \frac{g}{a}. \end{eqnarray} \tag{ 36 }$$

Here x and v refer to the one-dimensional comoving position and peculiar velocity of the baryonic gas, and g is the acceleration at the cell center. These equations are now in a form that can be solved by the split PPM scheme.

We now restrict ourselves to the solution of Equations (34)–(36) in one dimension. In the Direct Eulerian version of PPM, this is accomplished by a three-step process. First, we compute effective left and right states at each grid boundary by constructing a piecewise parabolic description of the primitive variables (ρ, u and E) and then averaging over the regions corresponding to the distance each characteristic wave can travel (u, u − c_s and u + c_s, where c_s is the sound speed). Second, a Riemann problem is solved using these effective left and right states, and finally the fluxes are computed based on the solution to this Riemann problem and the conserved quantities are updated. This is described in detail in CW84, but we will briefly outline the procedure here both for completeness and to put the changes we will make in context.

The Eulerian difference equations are:

$$\begin{equation} \rho _j^{n+1} = \rho _j^{n} + \Delta t \left(\frac{ \overline{\rho }_{j+1/2}\overline{v}_{j+1/2} - \overline{\rho }_{j-1/2}\overline{v}_{j-1/2} }{\Delta x_j} \right) \end{equation} \tag{ 37 }$$

$$\begin{eqnarray} \rho _j^{n+1} v_{j}^{n+1} &\,{=}\,& \rho _j^{n+1} v_j^n + \Delta t \nonumber\\ && \!\!\!\times\! \left(\!\frac{ \overline{\rho }_{j+1/2} \overline{v}_{j+1/2}^2 - \overline{\rho }_{j-1/2} \overline{v}_{j-1/2}^2 + \overline{p}_{j+1/2} - \overline{p}_{j-1/2}}{\Delta x_j} \!\right) \nonumber\\ && \!\!\!+ \frac{ \Delta t }{2} g_j^{n+1/2} \left(\rho _j^n + \rho _j^{n+1} \right), \end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray} &&\rho _j^{n+1} E_j^{n+1} = \rho _j^n E_j^n + \Delta t \nonumber\\ &&\quad \times \left(\displaystyle\frac{ \overline{\rho }_{j+1/2} \overline{v}_{j+1/2} \overline{E}_{j+1/2} - \overline{\rho }_{j-1/2} \overline{v}_{j-1/2} \overline{E}_{j-1/2} + \overline{v}_{j+1/2} \overline{p}_{j+1/2} - \overline{v}_{j-1/2} \overline{p}_{j-1/2} }{\Delta x_j} \right) \nonumber \\ && \quad + \displaystyle\frac{ \Delta t }{2} g_j^{n+1/2} \left(\rho _j^n v_j^n + \rho _j^{n+1} v_j^{n+1} \right). \end{eqnarray} \tag{ 39 }$$

We have used subscripts to indicate zone-centered (j) and face-centered (j + 1/2) quantities, while superscripts refer to position in time. The cell width is Δx_j. Although they have been discretized in space, the accuracy of the updates depend on how well we can compute the fluxes into and out of the cell during Δt. This in turn depends on our ability to compute the time-averaged values of p, ρ and v at the cell interfaces, denoted here by $\overline{p}_{j\pm 1/2}$, $\overline{\rho }_{j\pm 1/2}$, and $\overline{v}_{j\pm 1/2}$. We now describe the steps required to compute these quantities.

We first construct monotonic piecewise parabolic (third-order) interpolations in one dimension for each of p, ρ, and v. The pressure is determined from Equation (7), the equation of state. The interpolation formula for some quantity q is given by:

$$\begin{eqnarray} q_j(x) & = & q_{L,j} + \tilde{x}(\Delta q_j + q_{6,j}(1-\tilde{x})), \nonumber\\ \tilde{x} & \equiv & {x - x_{j-1/2} \over \Delta x_j}, \quad x_{j-1/2} \le x \le x_{j+1/2}. \end{eqnarray} \tag{ 40 }$$

This is Equation (1.4) of CW84 (in the spatial rather than mass coordinate, as is appropriate for the direct Eulerian approach). The quantities q_{L, j}, Δq_j, and q_{6, j} can be viewed simply as interpolation constants; however, they also have more intuitive meanings. For example, q_{L, j} is the value of q at the left edge of zone j, while Δq_j and q_{6, j} are analogous to the slope and first-order correction to the slope of q (see CW84 for a complete discussion):

$$\begin{equation} \Delta q_j \equiv q_{R,j} - q_{L,j} \quad q_{6,j} \equiv 6[q_j - 1/2(q_{L,j} + q_{R,j})]. \end{equation} \tag{ 41 }$$

We have reduced the problem to finding q_{L, j} and q_{R, j}. While this is simple in principle, it is complicated somewhat by the requirement that these values be of sufficient accuracy and that the resulting distribution be monotonic. That is, no new maxima or minima are introduced. The resulting formulae are straightforward but complicated and are not reproduced here, but see Equations (1.7)–(1.10) of CW84. We also optionally allow steepening as described in that reference.

Once we have the reconstruction, the primary quantities (p, ρ, v and E) are averaged over the domains corresponding to the three characteristics u−c, u, or u+c (where c is the sound speed in a cell). The linearized gas dynamics equations are then used to compute second-order accurate left and right states that take into account the multiple wave families. This process is described in CW84 and we use their Equations (3.6) and (3.7).

With these effective states, an approximation to the Riemann problem is found (see below for more detail about the Riemann solvers used), producing estimates for $\overline{p}_{j\pm 1/2}$, $\overline{\rho }_{j\pm 1/2}$, and $\overline{v}_{j\pm 1/2}$ that are third-order accurate in space and second-order accurate in time. These are then used to solve the difference Equations (37)–(39) for ρ^{n + 1}, v^{n + 1}, and E^{n + 1}.

We include an optional diffusive flux (and flattening for the parabolic curves) that can improve the solution in some cases. Our implementation follows that in the appendix of CW84.

In addition, as discussed earlier, the three-dimensional scheme is achieved by operator-splitting and repeating the above procedure in the other two orthogonal directions. The transverse velocities and any additional passive quantities are naturally and easily added to this system (see Equation (3.6) of CW84).

We note that the acceleration required in Equation (38) is actually the acceleration felt by the entire zone and not just at the zone center. Therefore, it is possible to find the mass-weighted average acceleration over the zone by expanding the density and acceleration distributions and retaining all terms up to second-order in Δx (Bryan et al. 1995), although we find that the potential is so slowly varying that this is unnecessary.

4.1.1. The Dual Energy Formulation for Very High Mach Flows

The system described thus far works well for gravitating systems with reasonable Mach numbers (<100) as long as the structures are well resolved. This section and the next detail changes that are required to correctly account for situations where one or both of these requirements are not met.

Large, hypersonic bulk flows appear to be very common in cosmological simulations and they present a problem because of the high ratio of kinetic energy E_k to gas internal energy e, which can reach as high as 10⁸. Inverted, we see that the internal energy consists of an extremely small portion of the total energy. In such a situation, the pressure, proportional to E − E_k, is the small difference between two large numbers: a disastrous numerical situation. This is not as large a problem as it may at first appear because it only occurs when the pressure is negligibly small. Therefore, even if we suffer large errors in the pressure distribution in these regions, the dynamics and total energy budget of the flow will remain unaffected. Nevertheless, if the temperature distribution is required for other reasons (e.g., for calculating radiative processes), a remedy is required.

To overcome this, we also solve the internal energy equation:

$$\begin{equation} \frac{\partial e}{\partial t} + \frac{1}{a} {\bf v}\cdot \boldsymbol {\nabla }{\rm e} = - 3 \frac{\dot{a}}{a} \frac{p}{\rho } - \frac{p}{a\rho } \boldsymbol {\nabla } \cdot {\bf v} \end{equation} \tag{ 42 }$$

in comoving coordinates. The structure is similar to the total energy equation; the second term on the left hand side represents transport, while the first term on the right is due to expansion of the coordinate system. It is differenced (again, in Eulerian form without the expansion term) as,

$$\begin{eqnarray} \rho _j^{n+1} e_j^{n+1} &\,{=}\,& \rho _j^n e_j^n \nonumber\\ && + \Delta t \left(\!\frac{ \overline{\rho }_{j+1/2} \overline{v}_{j+1/2} \overline{e}_{j+1/2} - \overline{\rho }_{j-1/2} \overline{v}_{j-1/2} \overline{e}_{j-1/2}}{\Delta x_j} \!\right) \nonumber\\ && - \Delta t p_j^{n} \left(\frac{ \overline{v}_{j+1/2} - \overline{v}_{j-1/2}}{\Delta x_j} \right) \end{eqnarray} \tag{ 43 }$$

Note that because of the structure of this equation, it is not in flux-conservative form. In particular, the pressure is evaluated at the cell center. Unfortunately, time-centering of this pressure has proved difficult to do without generating large errors in the internal energy and so we leave the pressure at the old time in this difference equation. This leads to some spreading of shocks; however, we note that this equation is only used in hypersonic flows.

It is necessary, however, to conserve the total energy so that the conversion of kinetic to thermal energy is performed properly. We must therefore combine the two formulations without allowing the separately-advected internal energy e to play a role in the gas dynamics. This is done by carrying both terms through the simulation and using the total energy E for hydrodynamic routines and the internal energy e when the temperature profile is required. One way to view this procedure is to treat e as enhanced precision (extra digits) for E that automatically "floats" to where it is needed. We only require that they be kept synchronized when the two levels of precision overlap.

When the pressure is required solely for dynamic purposes, the selection criterion operates on a cell by cell basis using,

$$\begin{equation} p = \left\lbrace \begin{array}{@{}ll}\rho (\gamma - 1)(E - {\bf v}^2/2),& (E - {\bf v}^2/2)/E > {\eta }_1, \\ \rho (\gamma -1)e,& (E - {\bf v}^2/2)/E < {\eta }_1 \end{array} \right.. \end{equation} \tag{ 44 }$$

It should be stressed that as long as the parameter η₁ is small enough the dual energy method will have no dynamical effect. We use η₁ = 10⁻³, which is consistent with the truncation error of the scheme for grid sizes that are typically used in our simulations. We are now free to select the method by which the internal energy field variable e is updated so that it will not become contaminated with errors advected by the total energy formulation but still give the correct distribution in shocked regions. Since we are concerned with the advection of errors, the selection criterion must look at each cell's local neighborhood. In one dimension, this is done with,

$$\begin{equation} e = \left\lbrace \begin{array}{@{}ll}(E - {\bf v}^2/2),& \rho (E - {\bf v}^2/2)/\max ({\rho }_{j-1}E_{j-1},{\rho }_j E_j,{\rho }_{j+1}E_{j+1}) > {\eta }_2,\\ e,& \rho (E - {\bf v}^2/2)/\max ({\rho }_{j-1}E_{j-1},{\rho }_j E_j,{\rho }_{j+1}E_{j+1}) < {\eta }_2 \end{array} \right.. \end{equation} \tag{ 45 }$$

Thus, η₂ determines when the synchronization (of e with E) occurs. Too high a value may mask relatively weak shocks, while spurious heating (via contamination) may occur if it is set too low. After some experimentation, we have chosen η₂ = 0.1, a somewhat conservative value. This scheme is optional and is generally only required in large-scale cosmological simulations where the gas cools due to the expansion of the universe but large bulk flows develop due to the formation of structure.

We note that others have independently developed a similar but distinct scheme for dealing with this problem, which is endemic to methods adopting the total energy equation. In Ryu et al. (1993), the two variables adopted are total energy and entropy (rather than total energy and thermal energy), with an analogous scheme for choosing which variable to employ.

4.1.2. Riemann Solvers and Fallback

The second method we describe is a MUSCL-based solver than can be used in both HD and MHD modes. The description here will be very brief both because the ideas are similar to those described in the previous section, and because this implementation has previously been described in more detail elsewhere (Wang et al. 2008; Wang & Abel 2009).

Much like the PPM solver, we have three basic steps: the first is reconstruction of the variables, the second is a solution of the Riemann problem, and the third is updating the conserved quantities with the fluxes as written above. For the reconstruction scheme we have implemented only the simple piecewise linear reconstruction (van Leer 1979; Colella & Glaz 1985), with options for both primitive and conservative variable reconstruction. The available Riemann solvers are HLL, HLLC, and HLLD, as described earlier.

To more clearly describe the Dedner divergence cleaning modifications, we write the equations of compressible inviscid hydrodynamics in the form of conservation laws as,

$$\begin{equation} \frac{\partial {U}}{\partial {t}} + \frac{\partial {F^x}}{\partial {x}} + \frac{\partial {F^y}}{\partial {y}} + \frac{\partial {F^z}}{\partial {z}}= 0, \end{equation} \tag{ 46 }$$

The conserved variable U is given by

$$\begin{equation} U = (\rho, \rho v_x, \rho v_y, \rho v_z, \rho E)^{T}, \end{equation} \tag{ 47 }$$

where ρ is density, v_i are the three components of velocity for i = x, y, z, E = v²/2 + e denotes the specific total energy and e the specific internal energy (note that in this section only, we use the specific energy).

For the generalized Lagrange multiplier (GLM) formulation of the MHD equations (Dedner et al. 2002), we consider these conserved variables

$$\begin{equation} U = (\rho, \rho v_x, \rho v_y, \rho v_z, \rho E+B^2/2, B_x, B_y, B_z, \psi)^{T}, \end{equation} \tag{ 48 }$$

where B_i with i = x, y, z are the three components of magnetic fields and ψ is the additional scalar field introduced in the GLM formulation for the divergence cleaning. The fluxes then are

$$\begin{eqnarray} F^x &=& \left(\rho v_x, \rho v_x^2+p+B^2/2-B_x^2, \rho v_yv_x-B_yB_x, \right. \nonumber \\ && \rho v_zv_x-B_zB_x, \rho \left({v^2\over 2} + h \right)v_x+B^2v_x-B_xB\cdot v, \nonumber \\ && \left. \psi, v_xB_y-v_yB_x, -v_zB_x+v_xB_z, c_h^2B_x \right)^T, \end{eqnarray} \tag{ 49 }$$

$$\begin{eqnarray} F^y &=& \left(\rho v_y, \rho v_xv_y-B_xB_y, \rho v_y^2+p+B^2/2-B_y^2,\right. \nonumber \\ && \rho v_zv_y-B_zB_y, \rho \left({v^2\over 2} + h \right)v_y+B^2v_y-B_yB\cdot v, \nonumber \\ && \left. v_yB_z-v_zB_y,\psi,-v_xB_y+v_yB_x, c_h^2B_y \right)^{T}, \end{eqnarray} \tag{ 50 }$$

$$\begin{eqnarray} F^z &=& \left(\rho v_z, \rho v_xv_z-B_xB_z, \rho v_yv_z-B_yB_z, \right. \nonumber \\ && \rho v_z^2+p+B^2/2-B_z^2, \rho \left({v^2\over 2} + h \right)v_z+B^2v_z-B_zB\cdot v, \nonumber \\ && \left. -v_yB_z+v_zB_y, v_zB_x-v_xB_z,\psi, c_h^2B_z \right)^{T}, \end{eqnarray} \tag{ 51 }$$

where c_h is a constant controlling the propagation speed and damping rate of ${\boldsymbol \nabla } \cdot {\boldsymbol B}$, and h = e + p/ρ denotes the specific enthalpy. All quantities are cell-centered.

The method is dimensionally un-split in that the fluxes are computed for all dimensions first and the conserved quantities are updated in one step, in contrast to the Strang splitting employed in the other fluid methods described in this paper. Also unlike the other schemes, time-integration is done with a second-order Runge–Kutta scheme (Shu & Osher 1988).

Finally, we note that for cosmological simulations, this solver uses a slightly different definition of the magnetic field than used in the rest of the paper. In particular, the field is defined as B = a^3/2B' (where B' is the proper field strength). This adds a source term of $-\dot{a}/(2a) {\bf B}$ on the right-hand side of Equation (4) and removes the factor of a in the B term in the energy equation (Equation (5)) and in the definition of the isotropic pressure p*.

The third solver we describe is an MHD method developed by Collins et al. (2010). Since a full description and suite of test problems can be found in that reference, we only describe the method briefly here.

The divergence of the magnetic field, $\boldsymbol {\nabla } \cdot {\bf B}$, is identically zero in reality due to the fact that the evolution of the magnetic field is the curl of a vector, and the divergence of the curl of a vector is identically zero. The CT method (Evans & Hawley 1988; Balsara & Spicer 1999) for MHD evolution employs this same vector property to evolve the magnetic field in a manner that preserves $\boldsymbol {\nabla } \cdot {\bf B}=0$. The electric field is computed using the fluxes from the Riemann solver. The curl of that electric field is then used to update the magnetic field. The advantage of this method is that it preserves $\boldsymbol {\nabla } \cdot {\bf B}$ to machine precision. The primary drawback is increased algorithm complexity. Note also that since only the update of the magnetic field is divergence free, any monopoles created by other numerical sources (such as ill-chosen initial conditions) persist.

The base Godunov method, described in Li et al. (2008), uses spatially and temporally second-order reconstruction (both MUSCL-Hancock and Piecewise Linear Method), and a selection of Riemann solvers including HLLC and HLLD (Mignone 2007), as described earlier. The CT methods are the first-order method described by Balsara & Spicer (1999) and the second-order methods described in Gardiner & Stone (2005). The AMR machinery is described by Balsara (2001) and Collins et al. (2010).

The increased complexity of the CT scheme comes in the form of area-averaged face and length-averaged edge-centered variables, while the rest of Enzo employs predominantly volume-averaged cell-centered variables. The magnetic field is represented by both a face-centered field, B_f, and a cell-centered field, B_c. The electric field is edge-centered. The magnetic field is updated in four steps: first, the Riemann problem is solved in the traditional manner, using the cell-centered field, second, an edge-centered electric field is computed using the fluxes from the Riemann solver, third, the curl of that electric field is used to update the face-centered field, finally, the cell-centered magnetic field is updated with an average of the face-centered field.

Divergence-free AMR is somewhat more complex than the AMR employed elsewhere in Enzo. First, the interpolation must be constrained to be divergence free. Thus, all three face-centered field components are interpolated in concert. Second, any magnetic field information from the previous timestep must be included in the interpolation, making the interpolation more complex than the simple parent-child relation used for other fields. Third, the flux correction involves more possible grid relations than traditional AMR. In order to circumvent this last complexity, the electric field is projected from fine grids to parent grids (rather than the magnetic field), and is then used to re-update the parent magnetic field. This is described in detail in Balsara & Spicer (1999) and Collins et al. (2010).

The dual energy formalism has also been incorporated in two possible ways—one that uses internal energy, as described in Section 4.1, and one that uses entropy (Ryu et al. 1993; Collins et al. 2010).

As an alternative to the previous Godunov methods, Enzo also includes an implementation of the finite difference hydrodynamic algorithm employed in the compressible magnetohydrodynamics code ZEUS (Stone & Norman 1992a, 1992b). Fluid transport is solved on a Cartesian grid using the upwind, monotonic advection scheme of van Leer (1977b) within a multistep (operator-split) solution procedure that is fully explicit in time. This method is formally second-order accurate in space but first-order accurate in time.

As discussed in the section describing the PPM (Section 4.1), operator-split methods break the solution of the hydrodynamic equations into parts, with each part representing a single term in the equations. Each part is evaluated successively using the results preceding it. In this method, in addition to operator-splitting the expansion terms (i.e., those terms in Equations (1)–(3) that depend on $\dot{a}$), we divide the remaining terms into source and transport steps. The terms to be solved in the source step are those on the right-hand side of Equations (1)–(3), while the transport terms are on the left-hand side of these equations and are responsible for the advection of mass, momentum and energy across the grid.

The ZEUS method uses a von Neumann–Richtmyer artificial viscosity to smooth shock discontinuities that may appear in fluid flows and can cause a breakdown of the finite difference equations. The artificial viscosity term is added in the source terms as:

$$\begin{eqnarray} \rho \frac{\partial {\bf v}}{\partial t} &=& {-} \boldsymbol {\nabla }p - \rho \boldsymbol {\nabla }\phi - \boldsymbol {\nabla } \cdot {\bf Q} \end{eqnarray} \tag{ 52 }$$

$$\begin{eqnarray} \frac{\partial e}{\partial t} &=& {-}p \boldsymbol {\nabla } \cdot {\bf v}- {\bf Q} : \boldsymbol {\nabla }{\bf v}. \end{eqnarray} \tag{ 53 }$$

Here Q is the artificial viscosity stress tensor, which we take to be diagonal with on-axis terms given by l²ρ(∂v/∂x)² as proposed by von Neumann & Richtmyer. The length scale l determines the width of shocks and is typically a few times the cell spacing.

In our Cartesian coordinate system, the finite difference version of these equations is particularly simple, although there is one important complication. In the ZEUS formalism, the velocity is a face-centered quantity—that is, the velocity is recorded on a grid that is staggered as compared to the density, pressure and energy, which are at the cell center. Therefore we must remember that v_j is at position x_{j − 1/2} (we use this notation, rather than writing v_{j + 1/2} both to match the original ZEUS paper and also to make it easier to compare these equations to what is actually in the code).

As in the original ZEUS paper, the source terms are added in three steps. First we add the pressure and gravity forces:

$$\begin{equation} v_j^{n+a} = v_j^n - \frac{\Delta t}{\Delta x_j} \frac{p^n_j - p^n_{j-1}}{\left(\rho ^n_j + \rho ^n_{j-1} \right)/2}. \end{equation} \tag{ 54 }$$

Partial updates are denoted by the n+a, n+b notation. We show updates in one dimension as the extension to the multi-dimensional case is straightforward (note that all dimensions are carried out for each substep before progressing to the next substep). We then add the artificial viscosity:

$$\begin{eqnarray} v_j^{n+b} & = & v_j^{n+a} - \frac{\Delta t}{\Delta x_j} \frac{q^{n+a}_j - q^{n+a}_{j-1}}{\left(\rho ^n_j + \rho ^n_{j-1} \right)/2} \end{eqnarray} \tag{ 55 }$$

$$\begin{eqnarray} e_j^{n+b} & = & e_j^n - \frac{\Delta t}{\Delta x_j} q^{n+a}_j \left(v^{n+a}_{j+1} - v^{n+a}_{j} \right). \end{eqnarray} \tag{ 56 }$$

The artificial viscosity coefficient q_j is given by:

$$\begin{equation} q_j = \left\lbrace \begin{array}{@{}ll}Q_{\rm AV} \rho _j (v_{j+1} - v_j)^2 & {\rm if}\ (v_{j+1} - v_j) < 0 \\ 0 & {\rm otherwise} \end{array} \right. \end{equation} \tag{ 57 }$$

where Q_AV is a constant with a typical value of 2. We refer the interested reader to Stone & Norman (1992a) and Anninos & Norman (1994) for more details. We also include the option (turned off by default) of adding a linear artificial viscosity as suggested in the ZEUS paper for stagnant flow regions. This is given by

$$\begin{equation} q_{{\rm lin},j} = Q_{\rm LIN} \rho c_j (v_{j+1} - v_{j}) \end{equation} \tag{ 58 }$$

where $c_j^2 = \gamma p/\rho$ is the adiabatic sounds speed.

Finally, the third source step is the compression term and is given by

$$\begin{equation} e^{n+c}_j = e^{n+b}_j \left(\frac{1 - (\Delta t/2) (\gamma - 1) (\boldsymbol {\nabla } \cdot {\bf v})_{j} }{1 + (\Delta t/2) (\gamma - 1) (\boldsymbol {\nabla } \cdot {\bf v})_{j} } \right). \end{equation} \tag{ 59 }$$

We have used the notation $(\boldsymbol {\nabla } \cdot {\bf v})_{\rm j}$ to indicate the (potentially) multi-dimensional velocity divergence evaluated at the cell center position x_j. This equation differs from the previous ones in that in the multi-dimensional case, still only one finite difference equation is evaluated, but the divergence becomes multi-dimensional.

We next examine the transport step, which is conservative. Once again, we dimensionally split the equations and present only the one-dimensional version. The finite difference equations actually solved are:

$$\begin{equation} \rho _j^{n+d} = \rho _j^{n} - \frac{\Delta t}{\Delta x} \left(v^{n+c}_{j+1/2} \rho ^{*}_{j+1/2} - v^{n+c}_{j-1/2} \rho ^{*}_{j-1/2} \right). \end{equation} \tag{ 60 }$$

Here $\rho ^*_j$ is the correctly upwinded value of ρ evaluated at the cell-face corresponding to v_j, making $\rho ^*_j v_j$ the mass flux at the cell boundary and guaranteeing mass conservation. This requires interpolating each cell-centered quantity to the cell edge. As recommended in Stone & Norman (1992a), we use the second-order van Leer scheme, which uses piecewise linear functions. These are given by Equations (48) and (49) of Stone & Norman (1992a). The transport steps for the other variables are similar. Note that we advect the specific energy and specific momenta using the mass flux, as dictated by the principle of consistent transport. This requires appropriate averaging for the momenta in the perpendicular directions as outlined in Equations (57)–(72) of Stone & Norman (1992a).

A limitation of a technique that uses an artificial viscosity is that, while the correct Rankine–Hugoniot jump conditions are achieved, shocks are broadened over three to four mesh cells. This may cause unphysical pre-heating of gas upstream of the shock wave, as discussed in Anninos & Norman (1994). On the other hand, it is much more robust than PPM and is easy to add additional physics. We also note that this method solves only the internal energy equation rather than total energy, so the dual energy formulation discussed in Section 4.1 is unnecessary.

Solving for the accelerations of the cells and particles on the grid due to self-gravity involves three steps: (1) computing the total gravitating mass, (2) solving for the gravitational potential field with the appropriate boundary conditions, and (3) differencing the potential to get the acceleration, and, if necessary, interpolating the acceleration back to the particles. These steps are described in detail below.

First, the massive (dark matter and star) particles are distributed onto the grids using the second-order cloud-in-cell (CIC) interpolation technique (Hockney & Eastwood 1988) to form a spatially-discretized density field ρ_DM. During the CIC interpolation, particle positions are (temporarily) advanced by 0.5vⁿΔt so that we generate an estimate of the time-centered density field. Particles on subgrids within the grid's volume are also added to its gravitating field using the same method. In addition, since the gravitating field for a grid is defined beyond the grid edges (see below), massive particles from sibling grids and sub grids that lie within the entire gravitating field are used. This step can involve communication.

Next, we add the baryonic grid densities in a similar fashion. In particular, we treat baryonic cells as virtual CIC particles that are placed at the grid center but are advanced by 0.5vⁿΔt in order to approximately time-center the gravitating mass field. Cells that are covered by further-refined grids are treated in a similar way (i.e., we also use the subgrid cells as virtual CIC particles). This procedure results in a total gravitating mass field $\rho _{\rm total}^{n+1/2}$.

To compute the potential field from this gravitating mass field on the root grid, we use a fast Fourier transform (FFT). For periodic boundary conditions, we can use either a simple Greens function kernel of −k⁻², or the finite difference equivalent (Hockney & Eastwood 1988):

$$\begin{equation} G({\bf k}) = - \frac{\Delta x}{2 (\sin (k_x \Delta x/2)^2 + \sin (k_y \Delta y/2)^2 + \sin (k_z \Delta z/2)^2 ) } \end{equation} \tag{ 61 }$$

where $k^2 = k_x^2 + k_y^2 + k_z^2$ is the wavenumber in Fourier space and the potential is calculated in k-space as usual with $\tilde{\phi }(k) = G(k) \tilde{\rho }(k)$.

For isolated boundary conditions, we use the James method (James 1977). In this case, the Greens function is generated in real-space so as to have the correct zero-padding properties and then transformed into the Fourier domain. In both cases, the potential is then transformed back into the real domain to get potential values at the cell centers. These are differenced with a two-point centered difference scheme to obtain accelerations at the cell centers (except if we are using the staggered ZEUS-like solver, in which case the accelerations are computed at the cell faces to match the velocities). Particle accelerations are obtained using a (linear) CIC interpolation from the grid.

In order to calculate more accurate potentials on the subgrids, Enzo uses a similar but slightly different technique from the root grid. The generation of the total gravitating mass field is essentially identical, using CIC interpolation for both the particles and baryons, including subgrids as before. To compute the potential on subgrids, however, we use the standard seven-point (in three dimensions) second-order finite difference approximation to Poisson's equation. Boundary conditions are then interpolated from the potential values on the parent grid. We use either tri-linear interpolation or a natural second-order spline for this: both methods give similar results, but the default is the tri-linear interpolation, which empirically provides a reasonable compromise between speed and accuracy. The potential equation on each subgrid is then solved with the given Dirichlet boundary conditions with a multigrid relaxation technique. This is applied to each subgrid separately.

The region immediately next to the boundary can contain unwanted oscillations (e.g., Anninos et al. 1994), and so we use an expanded buffer zone around the grid, of size three parent grid boundary zones (so typically six refined zones for a refinement factor of two). The density is computed in this region and the potential solved, but only the region that overlaps with the active region of the grid itself is used to calculate accelerations.

Simply interpolating the potential without feeding it back to higher levels leads to errors in the potential at more refined levels, due to the build-up of errors during the interpolation of coarse boundary values. In addition, neighboring subgrids are not guaranteed to generate the same potential values because of the lack of a coherent potential solve including the whole grid hierarchy. In an attempt to partially alleviate this problem, we allow for an iterative procedure across sibling grids, in which the potential values on the boundary of grids can be updated with the potential in "active" regions of neighboring subgrids. To prevent overshoot, we average the potential on the boundary and allow for (by default) four iterations, with the number of iterations determined by a parameter specified at runtime. This procedure can help in many cases, but does not necessarily produce a coherent solution across all grids and so does not completely solve the problem; we are working on a slower but more accurate method that does a multigrid solve across the whole grid (D. R. Reynolds et al., in preparation).

At this point it is useful to emphasize that the effective force resolution of an adaptive particle-mesh calculation is approximately twice as coarse as the grid spacing at a given level of resolution.

Enzo uses a particle-mesh N-body method to calculate the dynamics of collisionless systems (Hockney & Eastwood 1988). This method follows trajectories of a representative sample of individual particles that sample the phase space of the dark matter distribution, and is much more efficient than a direct solution of the Boltzmann equation in essentially all astrophysical situations for the levels of accuracy that are required for simulations of structure formation. As described earlier, the gravitational potential is computed by solving the elliptic Poisson's equation (Equation (8)) and differencing the potential to find accelerations, which are then interpolated back to particles. This acceleration is time-centered (because the underlying gravitating mass field is approximately time-centered), and so we have accelerations g^{n + 1/2} for each particle. These are used to update the particle positions and velocities starting from xⁿ and vⁿ using a standard drift-kick-drift technique (Hockney & Eastwood 1988):

$$\begin{eqnarray} {\bf x}^{n+1/2} & = & {\bf x}^{n} + \frac{\Delta {t}}{2 a^n} {\bf v}^{n} \nonumber \\ {\bf v}^{n+1} & = & {\bf v}^{n} \left(1 - \frac{\dot{a}^{n+1/2}}{{a}^{n+1/2}}\right) + \frac{\Delta {t}}{{a}^{ n+1/2}} {\bf g}^{ n+1/2} \nonumber\\ {\bf x}^{n+1} & = & {\bf x}^{n+1/2} + \frac{\Delta {t}}{2 a^{n+1}} {\bf v}^{n+1}. \end{eqnarray} \tag{ 62 }$$

Particles are stored in the most highly refined grid patch at the point in space where they exist, and particles that move out of a subgrid patch are sent to the grid patch covering the adjacent volume with the finest spatial resolution, which may be of the same spatial resolution, coarser, or finer than the grid patch from which the particles moved. This takes place in a communication process at the end of each timestep on a level.

To avoid unphysical point-mass effects, Enzo provides a parameter that governs the maximum level at which particles will be regarded as point masses. At higher levels, contributions from particles to the gravitating mass field will be smoothed over a spherical region centered at each particle's position.

While it is often safe to assume that species (both chemical and ionization) within a gas can be treated as being in equilibrium, in some regimes that are found in astrophysics this assumption leads to considerable error. For example, the cooling and collapse of primordial gas in Population III star formation is dominated by molecular hydrogen, which in the absence of dust forms via an inefficient pair of collisional processes that depend heavily on the local, highly non-equilibrium population of free electrons. As a result, when modeling primordial star formation it is critical to follow the non-equilibrium evolution of the chemical species of hydrogen, including molecular hydrogen and deuterium.

The primordial non-equilibrium chemistry routines used in Enzo were first described by Abel et al. (1997) and Anninos et al. (1997), but have since been extended with updated reaction rates and the inclusion of deuterium species (McGreer & Bryan 2008; Turk 2009). These routines follow the non-equilibrium chemistry of a gas of primordial composition with 12 total species: H, H⁺, He, He⁺, He⁺⁺, H⁻, H$_2^+$, H₂, e⁻, D, D⁺, and HD. Enzo also computes the radiative heating and cooling of the gas from atomic and molecular line excitation, recombination, collisional excitation, free–free transitions, Compton scattering of the cosmic microwave background, as well as several models for a metagalactic UV background that heat the gas via photoionization and photodissociation (see Section 6.2 for more details). The chemical and thermal states of the gas can be updated either at the same hydrodynamical timestep (i.e., decoupled and operator-split) or through the same subcycling system (i.e., a coupled chemical and thermal system). The default behavior of Enzo is to couple these two systems at subcycles of the hydrodynamic timestep; this results in updates to both the chemical and thermal states of the gas (which also inform the temperature, the reaction rate coefficients and the cooling functions of the gas) on timescales that are faster than those of the gas dynamics. The gamma used by Enzo to compute the temperature of the gas from the energy and density characteristics utilizes a variable gamma that includes effects of the rotational state of molecular hydrogen, enabling it to vary from 5/3 (fully-atomic) to 7/5 (fully molecular). This further coupling of the chemical and thermal states of the gas underscores the need for coupled chemistry and radiative cooling solutions.

Input parameters to Enzo govern the chemical species that are updated during the course of the simulation. This can include only the atomic species (H, H⁺, He, He⁺, He⁺⁺, and e⁻), those species relevant for molecular hydrogen formation (H₂, H$_2^+$, and H⁻), and can further include deuterium and its species (D, D⁺, and HD). A total of 9 kinetic equations are solved for the 12 species mentioned above, including 29 kinetic and radiative processes. See Table 2 for the collisional processes and Table 3 for the radiative processes solved.

Table 2. Collisional Processes Solved in the Enzo Nonequilibrium Primordial Chemistry Routines

Collisional Processes
(1)	H	+	e−	→	H+	+	2e−
(2)	H+	+	e−	→	H	+	γ
(3)	He	+	e−	→	He+	+	2e−
(4)	He+	+	e−	→	He	+	γ
(5)	He+	+	e−	→	He++	+	2e−
(6)	He++	+	e−	→	He+	+	γ
(7)	H	+	e−	→	H−	+	γ
(8)	H−	+	H	→	H2	+	e−
(9)	H	+	H+	→	H$_2^+$	+	γ
(10)	H$_2^+$	+	H	→	H2	+	H+
(11)	H2	+	H+	→	H$_2^+$	+	H
(12)	H2	+	e−	→	2H	+	e−
(13)	H2	+	H	→	3H
(14)	H−	+	e−	→	H	+	2e−
(15)	H−	+	H	→	2H	+	e−
(16)	H−	+	H+	→	2H
(17)	H−	+	H+	→	H$_2^+$	+	e−
(18)	H$_2^+$	+	e−	→	2H
(19)	H$_2^+$	+	H−	→	H2	+	H
(20)	2H	+	H2	→	2H2
(21)	2H	+	H	→	H2	+	H
(22)	H2	+	H2	→	H2	+	2H
(23)	3H			→	H2	+	H
(24)	D	+	e−	→	D+	+	2e−
(25)	D+	+	e−	→	D	+	γ
(26)	H+	+	D	→	H	+	D+
(27)	H	+	D+	→	H+	+	D
(28)	H2	+	D+	→	HD	+	H+
(29)	HD	+	H+	→	H2	+	D+
(30)	H2	+	D	→	HD	+	H
(31)	HD	+	H	→	H2	+	D
(32)	H−	+	D	→	HD	+	e−

Download table as:  ASCIITypeset image

Table 3. Radiative Processes Solved in the Enzo Nonequilibrium Primordial Chemistry Routines

Radiative Processes
(33)	H	+	γ	→	H+	+	e−
(34)	He	+	γ	→	He+	+	e−
(35)	He+	+	γ	→	He++	+	e−
(36)	H−	+	γ	→	H	+	e−
(37)	H2	+	γ	→	H$_2^+$	+	e−
(38)	H$_2^+$	+	γ	→	H	+	H+
(39)	H$_2^+$	+	γ	→	2H+	+	e−
(40)	H2	+	γ	→	H$_2^*$	→	2H
(41)	H2	+	γ	→	2H
(42)	D	+	γ	→	D+	+	e−

Download table as:  ASCIITypeset image

The chemical reaction equation network is technically challenging to solve due to the huge range of reaction timescales involved. The characteristic times for creation and destruction of the various species and reactions can differ by many orders of magnitude and are often very sensitive to the chemical and thermal state of the gas. This makes a fully-implicit scheme, with convergence criteria and error tolerance, strongly preferable for such a set of equations. However, most implicit schemes require an iterative procedure to converge, and for large networks (such as this one) an iterative, fully-implicit method can be very time-consuming and computationally costly for a relatively small increase in accuracy. At the present time, this makes fully-implicit methods somewhat undesirable for large, three-dimensional simulations.

Enzo solves the rate equations using a method based on a semi-implicit formulation in order to provide a stable, positive definite and first-order accurate solution. The update discretization splits chemical changes into formation components and destruction components and updates with a mixed set of time states, as described in Anninos et al. (1997). The formation components of species S_i are computed at the current subcycle time, where the contribution of species S_i to its own destruction components are computed at the updated time; all other contributions to the destruction component are computed at the current time. This mixed state improves accuracy, ensures species values are positive definite, and is equivalent to one Jacobi iteration of an implicit Euler solve. This technique is optimized by taking the chemical intermediaries H⁻and H$_2^+$, which have large rate coefficients and low concentrations, and grouping them into a separate category of equations. Due to their fast reactions, these species are very sensitive to small changes in the more abundant species and are (at almost all times in astrophysical calculations) close to their equilibrium values. Attempting to resolve their formation and destruction times would necessitate extremely small timesteps. Therefore, reactions governing these two species can be decoupled from the rest of the network and treated independently through analytic solutions for equilibrium values. This allows a significant speedup in the solution speed, as the timestepping scheme is applied only to the slower 7- or 10-species network (depending on whether deuterium is included or not), which will be much closer to the overall hydrodynamic timestep of the simulation.

Even so, the accuracy and stability of the scheme is maintained by subcycling the rate solver within a single hydrodynamic timestep. These subcycle timesteps are determined so that the estimated fractional change in the electron concentration is limited to no more than 10% per timestep; additional criteria may be applied based on the expected change in internal energy from radiative cooling and from chemical heating due to the formation of molecular hydrogen.

It is important to note the regime in which this model is valid. According to Abel et al. (1997) and Anninos et al. (1997), the reaction network is valid for temperatures between 10⁰ and 10⁸ K. The original model discussed in these two references was only claimed to be valid up to n_H ∼ 10⁴ cm⁻³. However, addition of the three-body H₂ formation process (Equation (20) in Table 2) allows correct modeling of the chemistry of the gas up until the point where collisionally induced emission from molecular hydrogen becomes an important cooling process, which occurs at n_H ∼ 10¹⁴ cm⁻³. A further concern is that the optically thin approximation for radiative cooling eventually breaks down, which occurs before n_H ∼ 10¹⁶–10¹⁷ cm⁻³ in gas of primordial composition. Beyond this point, modifications to the cooling function that take into account the non-negligible opacity in the gas must be made, as discussed by Ripamonti & Abel (2004), and was put into Enzo for the work published in (Turk et al. 2009; Turk 2009). The formation of molecular hydrogen as catalyzed by dust was recently added to Enzo to enable studies of low-metallicity gas, as well as the inclusion of appropriate timestepping criteria to account for the input of ionizing radiation. Even with these modifications, a completely correct description of the cooling of primordial gas at very high densities requires some form of radiation transport, which will greatly increase the cost of the simulations. Furthermore, at very high densities, the stiffness of the molecular hydrogen reaction rates may require better than a first-order accurate solution; as such, the transition to this regime will likely necessitate a fully-implicit, iterative solver.

Enzo has multiple methods for computing the energy change from radiative cooling and heating. All of them assume that the gas can be modeled either as completely optically thin or with a simple, local approximation to optical thickness. In this section, we describe the methods for computing the cooling rates from metal-free and metal-enriched gas. Sample cooling curves for each of Enzo's primary cooling methods are shown in Figure 3.

Zoom In Zoom Out Reset image size

6.2.1. Primordial Cooling

6.2.2. Metal Cooling

Enzo supports the use of a set of spatially uniform (but possibly time-varying) radiation fields that can interface with the chemistry and cooling/heating routines described in Section 6. Many of these use fits to the H, He and He⁺ ionizing and photo-ionization heating rates that are of the form

$$\begin{equation} {\rm rate} = k_0 (1+z)^{\alpha } \exp {\left(\frac{\beta (z-z_0)^2 }{1 + \gamma (z+z_1)^2} \right)} \end{equation} \tag{ 63 }$$

where the constant coefficients α, β, γ, z₀ and z₁ are fits from the literature. For radiation field types 1–3 (numbered as they appear in the Enzo source code and simulation parameter values), we give the coefficients for the photoionization (and photo-heating) rates of H, He, and He⁺ in Table 4. Radiation background types 1 and 2 are based on Haardt & Madau (1996) with two different intrinsic quasar spectra slopes, while Type 3 is from Haardt & Madau (2012), modified to match the normalized field found in Kirkman et al. (2005).

The remainder of the radiation field types either build on the fits in Types 1–3 or use a completely different form and therefore we describe them in the text rather than in Table 4. Type 4 is the same as Type 3 but also adds X-ray Compton heating from Madau & Efstathiou (1999), using Equations (4) and (11) of that paper.

Homogeneous radiation field types 5 and 6 start with a spectral shape that is then integrated against the appropriate H, He, and He⁺, cross-sections to compute the ionization and photo-ionization heating rates (we use 400 bins logarithmically spaced from 0.74 eV to 7.24 × 10⁹ eV). In particular, Type 5 has a hard, featureless quasar-like power law spectrum $f_{\nu } = f_{{\rm H\,\scriptsize{I}}} \nu ^{\alpha _0}$, where α₀ = 1.5 by default and the spectrum is normalized at the H i ionization edge. Type 6 has the same spectrum, but attenuated by a column density of 10²² cm⁻³ neutral hydrogen.

Types 8, 9, and 14 have only a photo-dissociating Lyman–Werner flux, with Type 8 being constant and Type 9 using the redshift-dependent results of Trenti & Stiavelli (2009). Type 14 uses a fit from Wise & Abel (2005) for the range 6 < z < 50, constant for z < 6, and zero for z > 50. Other types are either undefined or currently unused.

Stars and black holes strongly affect their surroundings through radiation. Radiation transport is a well-studied problem; however, its treatment in multidimensional calculations is difficult because of the dependence on seven variables—three spatial, two angular, frequency, and time. The non-local nature of the thermal and hydrodynamical response to radiation sources further adds to the difficulty. Here we briefly describe Enzo's ray tracing implementation Enzo+Moray, which is presented in full detail in Wise & Abel (2011) with seven code tests and six applications.

We solve the radiative transfer equation in comoving coordinates (given by Equation (20)). We can make some appropriate approximations to reduce the complexity of this equation in order to include radiation transport in numerical calculations. Typically timesteps in dynamic calculations are small enough so that Δa/a ≪ 1, therefore a_em/a ≈ 1 in any given timestep, reducing the second term to $\hat{n} \partial I_\nu /\partial {\bf x}$. To determine the importance of the third term, we evaluate the ratio of the third term to the second term. This is HL/c, where L is the simulation box length. If this ratio is ≪1, we can ignore the third term. For example at z = 5, this ratio is 0.1 when L = c/H(z = 5) = 53 proper Mpc. In large boxes where the light crossing time is comparable to the Hubble time, then it becomes important to consider cosmological redshifting and dilution of the radiation. Thus Equation (20) reduces to the non-cosmological form in this local approximation,

$$\begin{equation} \frac{1}{c} \frac{\partial I_\nu }{\partial t} + \hat{\bf n} \frac{\partial I_\nu }{\partial {\bf x}} = -\kappa _\nu I_\nu + j_\nu. \end{equation} \tag{ 64 }$$

We choose to represent the source term j_ν as point sources of radiation (e.g., stars, quasars) that emit radial rays that are propagated along the direction $\hat{n}$.

Ray tracing is an accurate method to propagate radiation from point sources on a computational grid as long as there are a sufficient number of rays passing through each cell. Along a ray, the radiation transfer equation reduces to

$$\begin{equation} \frac{1}{c} \frac{\partial P}{\partial t} + \frac{\partial P}{\partial r} = -\kappa P, \end{equation} \tag{ 65 }$$

where P is the photon number flux along the ray. To sample the radiation field at large radii, ray tracing requires at least N_ray = 4πR²/(Δx)² rays to sample each cell with one ray, where R is the radius from the source to the cell and Δx is the cell width. If one were to trace N_ray rays out to R, the cells at a smaller radius r would be sampled by, on average, (r/R)² rays, which is computationally wasteful because only a few rays per cell are required to provide an accurate calculation of the radiation field (as we will show later).

We avoid this inefficiency by utilizing adaptive ray tracing (Abel & Wandelt 2002), which is based on Hierarchical Equal Area isoLatitude Pixelation (HEALPix; Górski et al. 2005) and progressively splits rays when the sampling becomes too coarse. In this approach, the rays are traced along normal directions of the centers of the HEALPix pixels that evenly divide a sphere into equal areas. The rays are initialized at each point source with the photon luminosity (photon s⁻¹) equally spread across N_pix = 12 × 4^l rays, where l is the initial HEALPix level. We usually find l = 0 or 1 is sufficient because these coarse rays will usually be split before traversing the first cell.

The rays are traced through the grid in a typical fashion (e.g., Abel et al. 1999), in which we calculate the next cell boundary crossing. The ray segment length crossing the cell is

$$\begin{equation} dr = R_0 - \min _{i=1 \rightarrow 3} [({\bf x}_{ {{\rm cell}},i} - {\bf x}_{ {\rm src},i}) / \hat{\bf n}_{{\rm ray},i} ], \end{equation} \tag{ 66 }$$

where R₀, $\hat{\bf n}_{\rm ray}$, x_{cell, i}, and x_{src, i} are the initial distance traveled by the ray, normal direction of the ray, the next cell boundary crossing in the ith dimension, and the position of the point source that emitted the ray, respectively. However before the ray travels across the cell, we evaluate the ratio of the face area A_cell of the current cell and the solid angle Ω_ray of the ray,

$$\begin{equation} \Phi _c = \frac{A_{\rm cell}}{\Omega _{\rm ray}} = \frac{N_{\rm pix} (\Delta x)^2}{4\pi R_0^2}. \end{equation} \tag{ 67 }$$

If Φ_c is less than a pre-determined value (usually >3), the ray is split into 4 child rays. The pixel numbers of the child rays p' are given by the "nested" scheme of HEALPix at the next level, i.e., p' = 4 × p + [0, 1, 2, 3], where p is the original pixel number. The child rays (1) acquire the new normal vectors of the pixels, (2) retain the same radius of the parent ray, and (3) get a quarter of the photon flux of the parent ray. Afterward, the parent ray is discontinued.

A ray propagates and splits until at least one of the following conditions is met: (1) the photon has traveled c × dt_P, where dt_P is the radiative transfer timestep, (2) its photon flux is almost fully absorbed (>99.9%) in a single cell, which significantly reduces the computational time if the radiation volume filling fraction is small, (3) the photon leaves the computational domain with isolated boundary conditions, or (4) the photon travels $\sqrt{3}$ of the simulation box length with periodic boundary conditions. In the first case, the photon is halted at that position and saved, where it will be considered in the solution of I_ν at the next timestep. In the next timestep, the photon will encounter a different hydrodynamical and ionization state, hence κ, in its path. Furthermore any time variations of the luminosities will be retained in the radiation field. This is how this method retains the time derivative of the radiative transfer equation. The last restriction prevents our method from considering sources external to the computational domain. However, a uniform radiation background can be used in conjunction with ray tracing that adds the background intensity to the local radiation field.

The radiation field is calculated by integrating Equation (65) along each ray, which is done by considering the discretization of the ray into segments. In the following description, we assume the rays are monochromatic for simplicity. For convenience, we express the integration in terms of optical depth τ = ∫κ(r, t) dr, and for a ray segment

$$\begin{equation} d\tau = \sigma _{\rm abs}(\nu) n_{\rm abs} dr. \end{equation} \tag{ 68 }$$

Here σ_abs and n_abs are the cross section and number density of the absorbing medium, respectively. In the static case, Equation (65) has a simple exponential analytic solution, and the photon flux of a ray is reduced by

$$\begin{equation} dP = P (1 - e^{-\tau }) \end{equation} \tag{ 69 }$$

as it crosses a cell. We equate the photo-ionization rate to the absorption rate, resulting in photon conservation (Abel et al. 1999; Mellema et al. 2006). Thus the photo-ionization and photo-heating rates associated with a single ray (k_ph and Γ_ph, respectively) are

$$\begin{equation} k_{\rm ph} = \frac{P (1 - e^{-\tau })}{n_{\rm abs} \; V_{\rm cell} \; dt_P}, \end{equation} \tag{ 70 }$$

$$\begin{equation} \Gamma _{\rm ph} = k_{\rm ph} \; (E_{\rm ph} - E_i), \end{equation} \tag{ 71 }$$

where V_cell is the cell volume, E_ph is the photon energy, and E_i is the ionization energy of the absorbing material. In each cell, the photo-ionization and photo-heating rates from each ray in the calculation are summed. After the ray tracing is complete, these rates are used as inputs to the solver described in Section 2.2.5 to update the ionization, chemical, and energy states of the gas in each cell.

In addition to the ray-tracing approach for radiation transport described in Section 7.2, Enzo currently includes a field-based radiation transport solver for problems posed on uniform (i.e., non-AMR, non-static mesh refinement) grids, which has been tuned for large-scale simulations involving many ionizing sources. Detailed explanations of the model and solution approach may be found in Norman et al. (2008), Reynolds et al. (2009), Norman et al. (2009), and Norman et al. (2013), the salient features of which are reproduced here. In addition, comparisons of this solver with other astrophysical radiation transport solvers may be found in Iliev et al. (2009). Enzo's field-based radiation solver focuses on a flux-limited diffusion approximation for cosmological radiative transfer, with couplings to both the gas energy and chemical number densities.

The system of equations (Equations (21) and (22)) along with the chemical network (Equation (19)) is solved independently of Enzo's hydrodynamics, gravity, and dark-matter solvers (Sections 4.1–5.2), thereby allowing the advective portions of Equations (21) and (19) to be taken care of by the fluid solvers. Due to the disparate time scales between radiation transport and chemical ionization and heating, the remainder of these equations is solved using an operator-split algorithm. Within a given timestep to evolve $(E_r^n, e_c^n, {\mathrm{n}}_i^n) \rightarrow (E_r^{n+1}, e_c^{n+1}, {\mathrm{n}}_i^{n+1})$, we first evolve Equation (21): $(E_r^n, e_c^n, {\mathrm{n}}_i^n) \rightarrow (E_r^{n+1}, e_c^{n}, {\mathrm{n}}_i^{n})$. This uses an implicit Euler time discretization, and a second-order centered finite difference spatial discretization, resulting in a large linear system of equations. These are solved using a preconditioned conjugate gradient iteration, where the preconditioner consists of a geometric multigrid solver. Both of these linear solvers are provided by the HYPRE linear solver library (see Falgout & Yang 2002, 2012).

We then evolve the heating and chemistry system, Equations (22) and (19): $(E_r^{n+1}, e_c^n, {\mathrm{n}}_i^n) \rightarrow (E_r^{n+1}, e_c^{n+1}, {\mathrm{n}}_i^{n+1})$. Due to the lack of spatial derivatives (since advection is handled elsewhere), this system is a coupled system of nonlinear ordinary differential equations. This utilizes an implicit quasi-steady-state approximation, formulated as follows. We consider the modified equations,

$$\begin{eqnarray} \frac{\partial e_c}{\partial t} &=& {-}\frac{2\dot{a}}{a} e_c + \Gamma (\bar{E}_r,\bar{\mathrm{n}}_i) - \Lambda (\bar{E}_r,\bar{\mathrm{n}}_i), \end{eqnarray} \tag{ 72 }$$

$$\begin{eqnarray} \frac{\partial {\mathrm{n}}_i}{\partial t} &=& k_{i,j}(\bar{e}_c) {\mathrm{n}}_e \bar{\mathrm{n}}_j - {\mathrm{n}}_i \Gamma _i^{ph}(\bar{E}_r), \quad i=1,\ldots,N_s, \end{eqnarray} \tag{ 73 }$$

where we have defined the time-centered "background" states $\bar{E}_r = (E_r^{n}+E_r^{n+1})/2$, $\bar{\mathrm{n}}_i = ({\mathrm{n}}_i^{n}+{\mathrm{n}}_i^{n+1})/2$ and $\bar{e}_c = (e_c^{n}+e_c^{n+1})/2$. These equations may each be solved analytically for their solution at the time t, which we denote by

$$\begin{eqnarray} e_c(t) &=& {\rm sol}_{e}\left(t,\bar{E}_r,\bar{\mathrm{n}}_i,e_c^n\right), \end{eqnarray} \tag{ 74 }$$

$$\begin{eqnarray} {\mathrm{n}}_i(t) &=& {\rm sol}_{\mathrm{n}_i} \left(t,\bar{E}_r,\bar{e_c},\mathrm{n}_i^n\right), \quad i=1,\ldots,N_s. \end{eqnarray} \tag{ 75 }$$

We then define a nonlinear system of equations to compute the time-evolved solutions $(E_r^{n+1}, e_c^{n+1}, {\mathrm{n}}_i^{n+1})$ as

$$\begin{eqnarray} f_e\left(e_c^{n+1},{\mathrm{n}}_i^{n+1}\right) &\equiv & e_c^{n+1} - {\rm sol}_{e}\left(t^{n+1},\bar{E}_r,\bar{\mathrm{n}}_i,e_c^n\right) = 0, \end{eqnarray} \tag{ 76 }$$

$$\begin{eqnarray} f_{\mathrm{n}_i}\left(e_c^{n+1},{\mathrm{n}}_i^{n+1}\right) &\equiv & {\mathrm{n}}_i(t) - {\rm sol}_{\mathrm{n}_i} \left(t,\bar{E}_r,\bar{e_c},\mathrm{n}_i^n\right)=0, \nonumber\\ \quad i &=& 1,\ldots,N_s. \end{eqnarray} \tag{ 77 }$$

This system of N_s + 1 nonlinear equations is solved using a damped fixed point iteration,

$$U_i = U_i - \lambda f_i(U), \quad i=1,\ldots,N_s+1,$$

where U is a vector containing the solutions to Equations (76) and (77). In this iteration, for the first 50 sweeps we use λ = 1. For more challenging problems where this does not converge, we switch to a damping parameter of λ = 0.1.

Enzo implements the equations of isotropic heat conduction in a manner similar to that of Parrish & Stone (2007). The isotropic flux of heat is given by Equation (24) and we use a value for the Spitzer conduction coefficient, κ_sp = 4.6 × 10⁻⁷ T^5/2 erg s⁻¹ cm⁻¹ K⁻¹ (Spitzer 1962). In this situation we are using a value for the Coulomb logarithm, log Λ = 37.8, that is appropriate for the intracluster medium (Sarazin 1988)—in astrophysically relevant, fully ionized plasmas this value varies by no more than 50% (see, e.g., Smith et al. 2013). It is quite possible that the local heat flux computed in this way can become unphysically large in the high-temperature, low-density cluster regime when using this formulation; therefore, we take into account the saturation of the heat flux (Cowie & McKee 1977) at a maximum level of

$$\begin{equation} F_{{\rm sat}} \simeq 0.4 n_e k_b T \left(\frac{2 k_b T}{\pi m_e} \right)^{1/2}. \end{equation} \tag{ 78 }$$

To ensure a smooth transition between the Spitzer and saturated regimes, we define an effective conductivity using the formalism of Sarazin (1988)

$$\begin{equation} \kappa _{{\rm eff}} = \frac{\kappa _{\rm cond}}{1 + 4.2 \lambda _e / \ell _T} \;, \end{equation} \tag{ 79 }$$

where λ_e is the electron mean free path and $\ell _T \equiv T / |\boldsymbol {\nabla }T|$ is the characteristic length scale of the local temperature gradient. We also assume that the conductivity of the plasma can be described in terms of an effective conductivity, which can be expressed as a fraction f_sp of the Spitzer conductivity (where f_sp ⩽ 1.0 are considered physically realistic values). This takes into account physical processes below the resolution limit of the simulation, such as tangled magnetic fields, that can suppress heat transport.

Thermal conduction in a plasma can be strongly affected by the presence of magnetic field lines, which may suppress heat flow perpendicular to the magnetic field. In that case, we allow for heat transport only parallel to the magnetic field lines in magnetohydrodynamic simulations. Mathematically, this is given by Equation (25). As with the isotropic thermal conduction, we allow a multiplicative factor f_sp to take into account the possible suppression of magnetic fields below the resolution limit of the simulation.

Both isotropic and anisotropic thermal conduction in Enzo are treated in an operator-split manner. Furthermore, within the heat transport module, transport along the x, y, and z directions are computed in a directionally-split fashion, with heat flux along each direction calculated at the + and − faces of the cell using the arithmetic mean of the cell-centered temperature in cells n and n + 1 or n − 1 and n, respectively (empirically, this is more stable than taking the geometric mean of the cell-centered temperatures). The addition of transport along magnetic field lines requires the calculation of cross-terms in the temperature derivatives at cell faces, which can result in spurious oscillations in the temperature field in regions where the temperature gradient is strong in more than one spatial direction. Controlling these oscillations requires the addition of a flux limiter for calculations of the temperature field. In this case, we choose the monotonized central difference flux limiter (van Leer 1977a), which serves to maintain numerical stability without sacrificing substantial speed or accuracy.

8.2.1. Overview

8.2.2. Cen and Ostriker

The Cen & Ostriker method is a heuristic model of star formation on galactic scales. This method, first described by Cen & Ostriker (1992), assumes that stars form in substantially overdense, converging, and gravitationally unstable gas. Algorithmically, each cell at the locally highest level of refinement is examined at each timestep to see if it meets the criteria for star formation. Star particles are allowed to form in a cell if the following criteria are met:

$$\begin{eqnarray} \rho _b/\bar{\rho }_b & \ge & \eta, \end{eqnarray} \tag{ 80 }$$

$$\begin{eqnarray} \boldsymbol {\nabla } \cdot {\bf v}_b & < & 0, \end{eqnarray} \tag{ 81 }$$

$$\begin{eqnarray} t_{\rm cool} & < & t_{\rm dyn} \equiv \sqrt{3 \pi / 32G \rho _{\rm tot}}, \end{eqnarray} \tag{ 82 }$$

$$\begin{eqnarray} m_{b} & \, {>}\, & m_{J} \equiv G^{-3/2} \rho _{b}^{-1/2}c^{3} \left[ 1 + \frac{\delta \rho _{d}}{\delta \rho _{b}} \right]^{-3/2} \end{eqnarray} \tag{ 83 }$$

where η is the user-defined overdensity threshold, and m_b and m_j are the baryonic mass in the cell and the Jeans mass of the cell, and c is the isothermal sound speed in the cell. If all of these criteria are met, the mass of a star particle is calculated as m_* = m_b(Δt/t_dyn)f_*eff, where f_*eff is the star formation efficiency parameter.

If m_* is greater than a minimum star mass m_*min, a particle is created and given several attributes: mass, a unique index number, the time of formation t_form, the local dynamical free-fall time t_dyn and the metallicity fraction of the baryon gas in the cell f_Zb. There is a user-defined minimum dynamical time T_{dyn, min}, which is observationally motivated and affects the feedback rates (see below). The particle is placed in the center of the cell and given the same peculiar velocity as the gas in the cell, and is then treated in the same manner as the dark matter particles. An amount of gas corresponding to the new particle's mass is removed from the cell.

The star formation algorithm creates each star particle instantaneously. However, feedback should take place over a significant timescale, as all of the stars contained within the "star particle" would in reality form (and massive stars would die) over a substantial period of time. Therefore, we assume that for the purposes of feedback that the mass of stars formed at a time t with timestep Δt is:

$$\begin{eqnarray} \Delta m_{\rm sf} &=& \int _t^{t+\Delta t} \frac{dM}{dt} dt = \int _{\tau _0}^{\tau _1} m_* \tau e^{-\tau } d\tau \nonumber\\ &=& m_* [(1+\tau _0) e^{-\tau _0} - (1 + \tau _1) e^{-\tau _1} ] \end{eqnarray} \tag{ 84 }$$

where τ₀ = (t − t_form)/t_dyn and τ₁ = (t + Δt − t_form)/t_dyn.

During this timestep, the star particle returns metal-enriched gas and thermal energy from supernovae and from stellar winds. Since massive stars have very short lifetimes, we assume that there is an immediate feedback of some fraction f_SN of the rest energy from the stars into the gas, such that E_add = f_SNΔm_sfc², where c is the speed of light. In addition, a fraction f_Z* of the stellar mass is fed back in the form of metals. Finally, a fraction of the mass f_m* is added back into the gas along with momentum in order to simulate the mass ejection from all stars (not just supernovae).

There are six user-defined parameters in this algorithm: three deal with star formation (η, m_*min and t_{dyn, min}), and three deal with feedback (f_SN, f_Z* and f_m*). Some of these parameters are completely free, while others can be guided by observation or theory. For example, the supernova feedback parameter, f_SN, can be constrained by assuming that, for every 200 M_☉ of stars created, one supernova occurs, and this event feeds back approximately 10⁵¹ erg of thermal energy, giving:

$$\begin{equation} f_{\rm SN} = \frac{10^{51}\, \mathrm{erg}}{200\, M_\odot \, c^2} \simeq 3 \times 10^{-6}. \end{equation} \tag{ 85 }$$

The metal yield f_Z*, defined as the mass in metals produced per unit mass of stars created, can be constrained by, e.g., the theoretical model of Woosley & Weaver (1995). This model suggests that f_Z* = 0.02 is an appropriate number. The minimum dynamical time is set to be t_{dyn, min} = 10⁷ yr to agree with the SN timescales seen in nearby OB associations.

The other parameters, such as the overdensity threshold η, minimum star mass m_*min, and mass ejection fraction f_m* are not well constrained either theoretically or observationally. Indeed, m_*min is a purely numerical parameter designed to keep the code from producing too many star particles, and thus has no observational or theoretical counterpart. The η parameter, on the other hand, nominally has physical meaning (i.e., the density above which star formation must occur on a very short timescale in a self-gravitating cloud); however, in the vast majority of simulations the densities that are reachable are nowhere near the densities of protostellar clouds, and thus this parameter becomes a rough proxy for finding the densest environments in a given simulation.

8.2.3. Schmidt-law Method

8.2.4. H₂-regulated Method

The methods described in the previous two sections are generally used in simulations that have relatively poor resolution, Δx ≳ 1 kpc. At this physical scale, individual star forming regions are not resolved, so all uncertainty about the behavior of molecular clouds is folded into a density threshold for star formation. In calculations with much higher resolution, however (on the order of a few parsecs), individual molecular clouds can be resolved, thus rendering these approximations invalid. To that end, Kuhlen et al. (2012) have implemented a star formation algorithm that is specifically geared to high-resolution cosmological simulations of galaxy formation, where the formation of molecular hydrogen is followed directly and stars are allowed to form at the highest level of refinement when the local H₂ fraction exceeds a pre-determined threshold. Cells are examined every root grid timestep, Δt₀, and cells at the highest level that exceed the H₂ limit form star particles with masses proportional to the inferred mass of molecular hydrogen in the star-forming region (Krumholz et al. 2008, 2009; McKee & Krumholz 2010). The mass of the particle is calculated as:

$$\begin{equation} m_p = \epsilon \rho _{\rm gas} (\Delta x_{m})^3 \frac{\Delta t_0}{t_*} \end{equation} \tag{ 86 }$$

where Δx_m is the resolution of the maximum level of refinement, is an efficiency parameter with a standard value of 0.01 (as motivated by Krumholz & Tan 2007), and t_* is the local free-fall timescale.

The feedback method used in this method is identical to that described in Section 8.2.2.

8.2.5. Population III Star Formation

Unlike in the previous sections, under some circumstances it is both possible and desirable to simulate stars individually, rather than treating particles as ensembles. One particular example of this is Population III star formation (Abel et al. 2002; O'Shea & Norman 2007; Wise & Abel 2008; Turk et al. 2009), where a given halo may only form one star of primordial composition. To accommodate this, Enzo contains a star formation algorithm that forms individual Population III stars directly (Abel et al. 2007; Wise & Abel 2008; Wise et al. 2012). Using criteria similar to Cen & Ostriker (1992), a star particle forms when a cell meets all of the following conditions:

• 1.  
  A baryon overdensity of 5 × 10⁵ (corresponding to a hydrogen number density of roughly 10³ cm⁻³ at z = 10),
• 2.  
  A converging gas flow ($\boldsymbol {\nabla } \cdot {\bf v} < 0$), and
• 3.  
  A molecular hydrogen mass fraction $f_{\rm H_2} > 5 \times 10^{-4}$,

These are comparable to the conditions typical of a collapsing metal-free molecular cloud roughly 10 million years before the birth of a Pop III main-sequence star. If multiple neighboring cells are flagged as being able to form stars, a single star is created instead. This star has a mass that is randomly sampled from a stellar IMF with a functional form of

$$\begin{equation} f(\log M) dM = M^{\alpha } \exp \left [ - \left (\frac{M_{\rm char}}{M} \right)^{\beta } \right ]. \end{equation} \tag{ 87 }$$

Feedback from Population III stars created using this algorithm comes in multiple forms. Radiative feedback using the Moray radiation transport algorithm (Wise & Abel 2011) is available, using the mass-dependent hydrogen ionizing and Lyman–Werner photon luminosities and lifetimes of the Population III stars from Schaerer (2002). At the end of their main-sequence lifetimes, explosion energies, ejected gas mass, and ejected metal are calculated using a variety of sources that depend on the mass of the star, and are described in detail in Section 3.2.1 of Wise et al. (2012).

Detailed analysis of simulation results requires both the tools to ask sophisticated questions of the data and the ability to process vast quantities of data at high time-cadence. As simulations grow in size and complexity, storing data for post-processing simply becomes intractable. To cope with this, we have instrumented Enzo with the ability to conduct analysis during the course of a simulation. This enables analysis with extremely high time cadence (as often as every subcycle of the finest refinement level), without attempting to write an entire checkpoint output to disk. The current mechanism for conducting analysis in Enzo during the course of the simulation utilizes the same computational resources as are used by the simulation itself by transferring their usage from Enzo to the analysis routines; this is often referred to as in situ analysis or visualization. Utilizing a dynamically-scheduled second set of computing resources, often referred to as co-scheduled analysis or visualization, provides greater flexibility and overall throughput at the expense of simplicity.

We expose Enzo's mesh geometry and fluid quantities to the analysis platform yt (Turk et al. 2011b; Turk & Smith 2011). At compile time, Enzo is (dynamically or statically) linked against the Python and NumPy libraries necessary to create proxy objects exposing the mesh geometry, fluid quantities and particle arrays as NumPy arrays. This information is then passed to a special handler inside yt. yt interprets the mesh and fluid information and, without saving data to disk, constructs a native representation of the in-memory state of the simulation that appears identical to an on-disk simulation output. Enzo then executes a user-provided analysis script, which is able to access the in-memory simulation object. Once the analysis script has returned control to Enzo, the simulation proceeds. This process can occur at either the top of the main "EvolveHierarchy" loop or at the end of a timestep at the finest level, and the frequency with which it is called is adjustable by a run-time parameter. During the course of conducting analysis, the simulation is halted until the conclusion of the analysis.

Most analysis operations that can be performed on data sets that reside on disk can be performed on in-memory data sets in yt. This includes projections (i.e., line integrals, both on- and off-axis), slices, one-, two-, and three-dimensional fluid phase distributions, calculation of derived quantities and arbitrary data selection. As of version 2.5 of yt, the Rockstar phase-space halo finder (Behroozi et al. 2013) can be executed through yt on in-memory Enzo data, and so can the Parallel HOP halo finder (Eisenstein & Hut 1998; Skory et al. 2010). Operations that currently cannot be conducted on in-memory datasets are those that require spatial decomposition of data. For instance, calculating marching cubes on a data object with yt is a fully local operation and can be conducted in situ. However, calculating topologically-connected sets requires a spatial decomposition of data and thus cannot be conducted in situ. This prohibition extends to halo finding operations other than Rockstar, most multi-level parallelism operations, and volume rendering.

Where microphysical solvers or other operator-split physics calculations can be done in Python, yt can serve as a driver for these calculations. A major feature set that is currently being developed is to pass structured (i.e., non-fluid) information back from yt into Enzo. For instance, this could be the result of semi-analytic models of the growth and evolution of star clusters, galaxy particle feedback parameters that have been influenced by merger-tree analysis, or even spectral energy distributions that are calculated within yt and provided to Enzo as input into radiative transfer calculations. Future versions will include this, as well as the ability to dynamically allocate computational resources to yt such that the simulation may proceed asynchronously with analysis (co-scheduled analysis). With this functionality will also come the ability to dynamically partition data, such that spatially-decomposed operations such as volume rendering become feasible during the course of a simulation.

One of the inherent drawbacks of a grid-based fluid method is the inability to follow the evolution of a single parcel of fluid as it travels through the simulation volume. To address this, Enzo has the capability to introduce Lagrangian "tracer particles" into a calculation either at the beginning of the simulation or when restarting the calculation. These tracer particles are put into the simulation in a rectangular solid volume with uniform, user-specified spacing. Each particle's position and velocity is updated over the course of a single timestep Δt as follows:

$$\begin{eqnarray} x^{n+1/2} & = & x^n + (\Delta t/2) v^{{\rm interp},n} \nonumber \\ v^{n+1} & = & v^{{\rm interp},n+1} \nonumber\\ x^{n+1} & = & x^{n+1/2} + (\Delta t/2) v^{n+1}. \end{eqnarray} \tag{ 98 }$$

This is essentially a drift-kick-drift particle update from time n to time n + 1—however, instead of computing an acceleration at the half-timestep tⁿ + Δt/2 (as is done for massive particles—see Equation (62)), the particle's velocity is updated both at the beginning of each timestep and at the half-timestep by linearly interpolating the cell-centered baryon velocity to the position of the particle (v^interp), and assigning it to that value.

Enzo saves tracer particle data at user-specified intervals, independent of the intervals at which regular data sets are written. The data written out typically includes the tracer particle's unique ID and position, as well as the velocity, density, and temperature of the gas at its location, but is easily extensible to output any grid-based quantity that a user requires. This capability has been used quite effectively in several papers, including Silvia et al. (2010, 2012). It is crucial to keep in mind that tracer particles model a fixed number of Lagrangian fluid trajectories. The fluid on the grid, however, models the motion of all the mass and represents the average quantity in a grid cell's volume. Consequently, after a period of evolution, tracer particles—even if they initially had the same density distribution as the gas—will not have the same density distribution as the fluid. For example, they tend to accumulate at stagnation points of the flow, and care has to be taken in using these particles appropriately. Tracer particles are very useful in studies such as the variety of histories of the hydrodynamic quantities in Lagrangian fluid elements and when evaluating complex chemical and cooling models in regard to the simpler ones used in the actual numerical evolution.

Identification of shocks and their pre- and post-shock conditions can be accomplished through a combination of either temperature or velocity jumps with dimensionally split or unsplit search methods. The primary method used in Enzois the dimensionally unsplit temperature jump method described in detail in Skillman et al. (2008). We briefly outline the method here.

For every cell, we first determine whether it satisfies the following conditions necessary to be flagged as a shock:

$$\begin{eqnarray} &&\boldsymbol {\nabla } \cdot {\bf v} < 0, \end{eqnarray} \tag{ 99 }$$

$$\begin{eqnarray} &&\boldsymbol {\nabla }T \cdot \boldsymbol {\nabla }S > 0, \end{eqnarray} \tag{ 100 }$$

$$\begin{eqnarray} &&T_2 > T_1, \end{eqnarray} \tag{ 101 }$$

$$\begin{eqnarray} &&\rho _2 > \rho _1, \end{eqnarray} \tag{ 102 }$$

where v is the velocity field, T is the temperature, ρ is the density, and S = T/ρ^{γ − 1} is the entropy. T₂ and T₁ are the post-shock (downstream) and pre-shock (upstream) temperatures, respectively. An optional temperature floor may be chosen, which is useful for situations such as cosmological simulations without a radiation background where underdense gas in the IGM can cool adiabatically to unphysically low temperatures.

Once a cell is flagged, the local temperature gradient is calculated, which is then used to traverse cells parallel to the gradient to search for the first pre- and post-shock cell that do not satisfy the above conditions. If during the search a cell satisfying the conditions is found to have a more convergent flow, that cell is marked as the center, and the search is started again. Using the temperature values from each of these cells, the Mach number is then solved using the Rankine–Hugoniot temperature jump conditions:

$$\begin{equation} \frac{T_2}{T_1} = \frac{(5\mathcal {M}^2 - 1)(\mathcal {M}^2 + 3)}{16\mathcal {M}^2}, \end{equation} \tag{ 103 }$$

where $\mathcal {M}$ is the upstream Mach number.

Shock finding in the context of AMR is applied grid-by-grid. If a search for pre/post-shock cells goes outside the bounds of the grid ghost zones, the search is stopped for that particular shocked cell. In most situations this is adequate since the hydrodynamic shock is captured in fewer than the number of ghost zones. Shock finding can be run either upon data output or at each step in the evolution of an AMR level if the Mach and pre/post-shock quantities are needed for additional physics modules.

11.1.1. The Enzo Test Suite

11.1.2. Code Comparisons

The Enzo test suite (described in Section 11.1.1) contains hundreds of test problems that probe the code's behavior in a wide range of physical circumstances, and explicitly tests each physics package in the Enzo code, both individually and in combination. It is impractical to include a substantial fraction of these problems in a method paper; as a result, we have chosen to publish the results of only a small subset of particularly crucial test problems here. If the interested reader desires, they can download Enzo and yt, run the test suite, and see the results of any other test problem and its comparison to an analytic solution (if available), or to a "gold standard" solution from a stable version of Enzo. We also note that there are many other Enzo test problems in the associated papers which describe individual methods, including: MHD-CT (Collins et al. 2010), MHD-Dedner (Wang et al. 2008; Wang & Abel 2009), Radiation-flux-limited diffusion (Norman et al. 2008; Reynolds et al. 2009; Norman et al. 2009), Radiation-ray-tracing (Wise & Abel 2011).

The general structure of each test problem description is as follows. We will describe the construction of the test problem (including its initial and boundary conditions), the analytical or expected solutions, and motivate why we have included it in the paper. After that, we will show and describe Enzo's solution to the test problem. We remind the user that, as discussed at the end of Section 11.1.1, they can download a Mercurial repository containing this paper, the Enzo test problem parameter files that produced the simulation data used to generate the figures in this section, and the yt scripts necessary to create the figures themselves.

11.2.1. Sod Shock Tube

We begin the paper test suite with the classic one-dimensional Sod shock tube problem (Sod 1978), which provides a good test of a hydrodynamical solver's ability to resolve a clean Riemann problem with clear separation between the three resultant waves. These waves consist of a rarefaction fan, a contact discontinuity, and a moderate shock. The initial state is ρ_L, P_L = 1.0, 1.0 on the left of the boundary at x = 0.5 and ρ_R, P_R = 0.125, 0.1 on the right. All velocities are initially zero. In Figure 4, we show the density solution at the final time (t = 0.25) for three of our hydro solvers—the spatially third-order PPM, as well as the two second-order ZEUS and MUSCL schemes. We use 100 cells across the domain, which is a relatively standard choice in code method papers, and show solutions both with and without AMR. Without AMR (top row in Figure 4), it is clear that the PPM scheme produces by far the cleanest solution with all wave families crisply reproduced (in particular, the contact discontinuity and shock). ZEUS and MUSCL produce similar results, with MUSCL doing a slightly better job on the rarefaction fan. In Figure 4, we also show the integrated absolute deviation from the exact solution, ||E₁|| = ∑_iΔx_i|F(x_i) − Exact(x_i)|. These numbers confirm the qualitative differences noted previously.

Zoom In Zoom Out Reset image size

We also run the same simulation but with two levels of AMR (using a refinement factor of two), triggered based on a normalized slope greater than 10% in the density. This refinement criterion results in only the refinement of strong gradients, and does not include the rarefaction fan at late times. The results are shown in the bottom row of Figure 4. Using AMR, the results are much better for all three methods, with much sharper shocks and contact discontinuities and even a better representation of the rarefaction wave, which is only refined beyond the root grid at early times. Although the results are improved for all methods, PPM still produces the best result, as is shown clearly by the computed error norms (displayed again in each of the individual panels).

11.2.2. Wave Pool

In this one-dimensional test we pass a short wavelength linear wave through a static singly-refined region. The full domain is from 0 to 1 and the refined region is set from 0.25 to 0.75 at all times. We modify the left boundary consistent with a single linear sound wave with density given by ρ(x, t) = ρ₀(1 + Asin (kx − ωt)), where the amplitude is A = 0.01 and the wavelength λ = k/2π = 0.1. Similar expressions exist for the pressure and velocity. The initial, unperturbed density and pressure are set to unity and we adopt γ = 1.4. The (unrefined) root grid is covered with 100 cells, resulting in a wavelength for the linear wave of only 10 cells. This is, therefore, a challenging problem for hydro methods. We are particularly interested in any reflection or artifacts introduced by the wave entering and exiting the refined regions.

Figure 5 shows the evolution of the wave at three different times (t = 0.2, 0.3 and 0.8) for three of our solvers (PPM, ZEUS, and MUSCL). In all panels, the shaded regions denote the statically refined region. The leftmost column shows the wave just before it enters the refined region so that we can gauge how the solver is operating in the absence of AMR. The center column shows the wave after it has fully entered the refined region. The rightmost column shows the wave well after it has exited the refined region.

Zoom In Zoom Out Reset image size

With the PPM solver, we notice that even before the wave reaches the refined region it has been slightly damped. This is not unexpected for such a short wavelength mode, even for higher-order solvers such as PPM. The remaining panels demonstrate that the wave cleanly enters and exits the refined region. No significant reflection is seen on either entry or exit, and the amount of damping is mild.

For the ZEUS solver, we see that even before it has entered the refined region there are small oscillations excited behind the wave, although it is also worth noting that the wave itself is beautifully propagated without significant damping or phase errors. In this test we use only our standard, low amount of quadratic artificial viscosity. The trailing oscillations could be damped by additional viscosity, but we do not add any in order to be sensitive to any artifacts at the grid boundary. The remaining panels show that the trailing oscillations continue, but do not generate any additional noise. We note that the end result is similar to the case without any refined region.

Finally, the MUSCL solver also shows a very clean entry and exit from the refined region without any oscillations, although with mild damping on exit. However, because we use a piecewise linear reconstruction, the wave is spread more than with the other methods.

11.2.3. Shock Pool

The next problem is similar to the last one but instead examines how a shock with a Mach number of 2 passes in and out of a static refined region (again defined from z = 0.25 to 0.75). The density and pressure in the domain are initially set equal to 1.0 with zero velocity. At t = 0 the left boundary is set with the density, pressure and velocity appropriate for a $\mathcal {M}=2$ shock wave.

In Figure 6, we show the evolution of the shock wave at three times, again corresponding to just before entering the refine region (left column), after entering the refined region (center), and after exiting the refined region (right). We examine the same set of three solvers as in the previous test problem.

Zoom In Zoom Out Reset image size

Beginning with PPM (top row of Figure 6), we see that this method captures the shock in a small number of zones with only a small amount of oscillation. Upon entering the refined region, the shock finds itself broader than the natural width of the scheme (since the cell spacing is decreased by a factor of two), and so the shock front contracts, causing a slight entropy perturbation in the post-shock gas.

For ZEUS (middle row of Figure 6), the shock is broadened because of artificial viscosity and there are slightly more post-shock oscillations, although again quite mild. The impact of entering and exiting the refined region is somewhat larger than in PPM; however, the most noticeable difference is the incorrect position of the shock front. This is due to the fact that the scheme is not energy-conserving (see also the Sedov problem in Section 11.2.5), and has very little to do with the passage through the refined region.

Finally, the MUSCL scheme produces a shock that is intermediate in width between the two previous cases. This method is energy-conserving, and thus correctly reproduces the shock speed. The oscillations are mild except for the cell immediately outside of the refined region upon exiting.

11.2.4. Double Mach Reflection

The double Mach reflection test is a classic two-dimensional test of hydrodynamic algorithms originally described in Woodward & Colella (1984a) (and more recently in Stone et al. 2008). In this problem, shown in Figure 7, a shock is injected at an angle to a reflecting surface (the −y-boundary), and a jet appears along the reflecting surface. The ideal solution is self-similar, and the appearance of this solution is highly sensitive to numerical diffusion. If numerical noise is present, a Kelvin–Helmholz instability develops along this jet and breaks the self-similarity.

Zoom In Zoom Out Reset image size

In the test problem shown in Figure 7, a two-dimensional simulation with 960 × 240 cells was created with a domain of x = [0, 4] and y = [0, 1]. We use an ideal gas equation of state of γ = 1.4, a pre-shock density of 1.4, and a pre-shock specific internal energy of 2.5/1.4 (all in arbitrary units). A Mach 10 shock is initialized with a shock normal that is 30° from the x-axis and an initial position on the lower boundary of x = 1/6. The lower y-boundary and right x-boundary are reflecting; the left x boundary is inflowing, and the upper y-boundary has a time-dependent boundary condition that allows the shock to propagate into the domain as if it extends to infinity. The simulation starts at t = 0 and runs until t = 0.205 (arbitrary units), at which point the rightmost extent of the shock should be at roughly x = 3. In this simulation we use the direct Eulerian implementation of the piecewise parabolic hydrodynamic method with the diffusion, flattening, and shock steepening all enabled.

It is instructive to compare our Figure 7 to Figure 9 in Woodward & Colella (1984a). By the end of the simulation, a dense jet is apparent at the leading edge of the shock, propagating along the x-axis. The shape of this jet is sensitive to numerical diffusion, and our figure compares favorably to those shown in Woodward & Colella.

11.2.5. Sedov Explosion

The Sedov Blast test (Sedov 1959) models an intense explosion, initiated by depositing thermal energy into a homogenous distribution of gas. The result is a strong spherical shock wave centered on the point of energy injection. This problem is a popular test of astrophysical simulation codes for three reasons. First, this problem is representative of the astrophysical phenomenon of a supernova explosion. Second, an analytical solution to the evolution of the Sedov blast wave exists, allowing for a direct test of the accuracy of the code. The radius of the shock front as a function of time is given by

$$\begin{equation} r(t) = \left(\frac{E_0}{\alpha \rho _0}\right)^{1/5}t^{2/5} \end{equation} \tag{ 104 }$$

where E₀ is the initial energy injection, ρ₀ is the background density and α = 1.0 for cylindrical symmetry and an ideal gas with γ = 1.4; for the full derivation see Sedov (1959). Third, as the spherical shock expands, its symmetry, or lack thereof, serves to highlight any directional preferences of the hydrodynamics solvers. The test presented here is the two-dimensional version that is included in the Enzo distribution. The three-dimensional results from this test, both for Enzo and three other leading astrophysics codes, can be found in Tasker et al. (2008).

In the initial state, the box contains a homogenous distribution of gas at a density of 1 (note that, in the absence of gravity or radiative cooling, this problem is scale-free and thus without units). Thermal energy is deposited into a single cell at the center of the box with E₀ = 10.0. The problem is run in two dimensions with reflecting boundary conditions and a box having a length of 1 in both directions. We select a top grid of 100 × 100 cells and a maximum of four levels of adaptive refinement, refining by factors of two on shocks and the slopes of the density and total energy fields. The exception to this scheme is in the initial conditions, where additional refinement is placed directly around the injection point to better resolve spherical geometry. The results are assessed at t = 0.07, which corresponds to a time just before the shock reaches the box edge (see Figure 8).

Zoom In Zoom Out Reset image size

Figure 9 shows the radial profiles for the simulation run with the PPM hydro-solver (blue dashed line) and the ZEUS hydro-solver (red dash-dot-dot line) together with the analytical solution (black solid line). Clockwise from the top-left are density, velocity, internal energy and pressure. PPM matches the analytical solution extremely well for all quantities. However, the shock front in the ZEUS simulation lags behind the analytical position. This can also be seen in the slices shown in Figure 8. The cause of this discrepancy is that ZEUS shows a substantial energy loss during the first few timesteps and produces a diamond-shaped, rather than spherical, shockfront during this time. After this, the code correctly conserves energy but this initial energy loss remains clearly visible in the position of the shock at t = 0.07. This problem was addressed directly by Clarke (2010), who attributed the source of the issue to this version of ZEUS solving the internal, rather than total, energy equation. In situations with strong energy gradients, this choice caused an energy loss and the artificial viscosity produces the direction-dependent shockfront shape. In their paper, Clarke (2010) present results from an alternative version of ZEUS that conserves total energy. This problem is less marked for smaller energy gradients and it should be noted that the ZEUS hydro algorithm's stability and speed make it a highly competitive choice, despite the disagreements in this test.

Zoom In Zoom Out Reset image size

11.2.6. Point Source Gravity Test

This is a simple test of the ability of the Poisson solver to correctly reproduce the gravitational acceleration around a point source without any dynamics. We place a single particle in the center of the unit domain with root grid dimensions 32³ and add a static, nested grid refined by a factor of two that covers the central 1/8³ of the domain (i.e., so the refined grid is only 8³). We then place 5000 particles throughout the domain distributed uniformly in angle and in the logarithm of the radius (from the central point). The measured accelerations are shown in Figure 10. As is typical for particle-mesh based methods, the errors decrease at large distances and peak on roughly the grid scale. The accelerations are relatively smooth across the top grid/fine grid boundary with the most noticeable impact being a drop in the amplitude of tangential errors on the fine grid. The force law rolls over at about 1.6 fine cell widths (equivalent to 0.8 root grid widths in this figure), consistent with the CIC-deposition and interpolation used.

Zoom In Zoom Out Reset image size

11.2.7. Orbit Test

The "orbit test" simulates the orbit of an effectively massless particle around a massive primary particle and is used to test the accuracy of both the gravity solver and particle integration. In this problem (shown in Figure 11), a three-dimensional box with unit size and a 32³ grid is created, and a particle with mass 1 is placed at the center of the box (both box length and particle mass are in arbitrary units). A second "test" particle, with mass of 10⁻⁶, is placed at a distance of 0.3 L_box from the primary particle in the x–y plane (at z = 0.5 L_box). Both particles are given velocities such that they will orbit their common center of mass in the x–y plane. The simulation is run for 200 orbits, writing out a dataset approximately once per orbit (with the particle positions being written out several hundred times per orbit).

Zoom In Zoom Out Reset image size

The left panel of Figure 11 shows the path of the test particle in the x–y plane (at z = 0.5), plotted over the duration of the simulation. The test particle maintains a path that has effectively unnoticeable deviation from circularity, and with no discernable drift in either net radius or orbit center over time. The right panel displays the total specific energy of the test particle (defined as (1/2)v² + ϕ, where ϕ is the gravitational potential at the position of the test particle) as a function of time, plotted once per orbit. The mean total specific energy is −1.76797 in arbitrary units, with a standard deviation of 0.004885 (0.276% of the mean value), and with a maximum deviation from the mean of 0.546% of the mean value. Given that the potential is calculated on 32³ cells, the two particles are separated by between 10 and 14 cells depending on their relative positioning in the grid, and the particles are integrated using a second-order leapfrog method, this level of accuracy is expected and acceptable.

11.2.8. Self-similar Infall Test

This test problem is based on the self-similar solutions to the cosmological spherical infall problem found by Bertschinger (1985). It features a small density perturbation in a homogeneous Ω = 1 universe and is a strong test of the cosmological evolution, gravity, and hydrodynamic portions of the code. It is a close analog to cosmological halo formation.

For the initial conditions, we adopt a 32³ top grid with a single, initial subgrid (covering 1/8³ of the domain) with a refinement factor of two. An overdensity δ = 40 is placed in a single cell near the center of the domain. We begin at z = 199 and evolve to z = 0, which is a sufficient time for the evolution to largely forget its initial conditions and approach the self-similar result. An overdensity refinement criteria (δ_crit = 1.1 on the top grid) is used to add additional grids, going up to 5 additional levels beyond the root grid. We use the PPM solver without radiative cooling; only baryons are used for this problem. An ideal gas law with γ = 5/3 is adopted.

In Figure 12, we show the results, scaled according to the dimensionless variables as defined in Equations (2.9) and (3.2) of Bertschinger (1985). This demonstrates that we can quickly and easily obtain a good solution with only a fairly modest initial grid. The shock is sharply resolved and the asymptotic profiles at small λ are recovered. At very small values of λ, the initial conditions have not been fully forgotten and the self-similar solution is not recovered. This can be seen most clearly in the dimensionless mass. However, this is simply because of the limited amount of time for which we evolve the solution.

Zoom In Zoom Out Reset image size

11.2.9. Zel'dovich Pancake

The Zel'dovich pancake (Zel'dovich 1970) is particularly relevant to cosmological simulations, because it includes many features that are critical to structure formation: hydrodynamics, expansion, and self gravity. This problem represents the formation of an ideal, isolated caustic, and is thus a useful proxy for much more complicated structures in full three-dimensional cosmological simulations, such as the collapse of gas onto a cosmological halo or filament.

The initial conditions are simple, and we follow the prescription of Anninos et al. (1994). Assuming a geometrically flat cosmology, the density perturbation is given by

$$\begin{equation} \rho (x_l) = \rho _0 \left[ 1 - \frac{1+z_c}{1+z} \cos (k x_l) \right] \end{equation} \tag{ 105 }$$

with the internal energy of the gas set so that the entropy is constant throughout. The velocity perturbation is given by

$$\begin{equation} v(x_l) = -H_0 \frac{1 + z_c}{(1+z)^{1/2}} \frac{\sin (k x_l)}{k}. \end{equation} \tag{ 106 }$$

In the equations above, ρ₀ is the background density, z_c is a free parameter and is the redshift where the sheet forms a caustic (i.e., where it "pancakes"), z is the redshift of initialization, x_l is the Lagrangian mass coordinate, k = 2π/λ (where λ is the perturbation wavelength), and H₀ is the value of the Hubble constant at z = 0. Note that this solution is expressed in terms of Lagrangian positions, so one needs to convert this into the Eulerian coordinates, x_e, that are more useful to a grid-based calculation:

$$\begin{equation} x_e = x_l - \frac{1 + z_c}{1 + z} \frac{\sin (k x_l)}{k}. \end{equation} \tag{ 107 }$$

We note that the solution described above is exact up to the point of caustic formation. In Figure 13, we show the results of a test of the adaptive mesh version of Enzo's Zel'dovich Pancake test. A one-dimensional box of length 64 Mpc h⁻¹ is initialized at z = 20 in an Ω_M = 1 universe with h = 0.5 and a background temperature of 100 K, with a background density of ρ₀ = ρ_c. The simulation is initialized with 64 grid cells, refining by factors of four using criteria based on cell mass and the presence of shocks, for a maximum of two levels (i.e., an equivalent maximum resolution of 1024 grid cells). The simulation is evolved to z = 0 using the PPM hydro method.

Zoom In Zoom Out Reset image size

The final output of the calculation, with the key features of this test problem, is shown in Figure 13. The strong shocks and large density gradients are well-resolved, with density and velocity jumps being well-delineated and at the correct locations. The key features of this test problem can be resolved with far fewer cells—simulations including a mere 8 cells resolve the key features, as shown in Sections 3.3.4–3.3.5 of Bryan (1996)—but we choose a higher resolution here for illustrative purposes.

11.2.10. MHD: Brio & Wu Shock Tube

A standard test of MHD solvers is the MHD analog of the Sod shock tube, described by Brio & Wu (1988). This test, which uses γ = 2.0 and B_x = 0.75 throughout; ρ = 1.0, B_y = 1.0, and P = 1.0 in the left state; and ρ = 0.125, B_y = −1.0 and P = 0.1 in the right state. A total of 800 grid points were used, and run to a time of t = 0.08. Both tests were run using the HLLD Riemann solver of Miyoshi & Kusano (2005) and a Piecewise Linear Method (Colella & Glaz 1985) reconstruction. The results of the test are shown in Figure 14. In that figure, the black lines show the results of the CT solver while the gray lines show the reference solution, run at 10,000 points. A light gray line shows the Dedner results, but is completely indistinguishable from the points of the CT method. To highlight the differences between the CT and Dedner solvers, the lower pane of each plot shows the difference between the two solutions, which is at most a few percent. Most shocks are resolved by a few grid points, though the contact discontinuity (the large jump seen in density but not the other quantities) is slightly more smeared.

Zoom In Zoom Out Reset image size

11.2.11. MHD: Orszag–Tang Vortex

Figure 15 shows the Orszag–Tang vortex problem (Orszag & Tang 1979). The left panel shows the result using Enzo's CT MHD method, while the right panel shows the result using Enzo's implementation of Dedner MHD. The Orszag–Tang vortex test is a classic MHD test problem, and shows that significant small scale structure can be generated in MHD from large scale initial perturbations. It is often used to compare the effective resolution of different MHD schemes. The test begins with uniform density, ρ₀ = 25/(36π), and pressure, P₀ = 5/(12π) (as with other tests in this section, in the absence of gravity or chemistry/cooling we use dimensionless units). There is a single rotational mode in the velocity, and two in the magnetic field: ${\bf v}_0 = (-\sin (2\pi y) \hat{x},\, \sin (2\pi x) \hat{y})$, ${\bf B}_0 = (-\sin (2\pi y) \hat{x},\, \sin (4\pi x) \hat{y})$. The simulation is evolved to t = 0.48. One can see that the structures are accurately represented as compared to, for example, Tóth (2000), and that the resolution of shocks is comparable in both methods.

Zoom In Zoom Out Reset image size

11.2.12. One-zone Collapse Test

The one-zone collapse test simulates the collapse of a self-gravitating gas cloud using a semi-analytic model for the evolution of gas density and adiabatic heat input as a function of time. It is designed to test the chemistry and cooling modules over a wide range in densities and over physically motivated timescales. Because this test disables the hydrodynamic and gravity solvers and uses a simple model for the density evolution, it is far faster than running a true collapse simulation. The density evolution is based on the self-similar Larson–Penston solution for isothermal collapse (Larson 1969; Penston 1969) with a modification to account for the efficiency with which the heat introduced by compression can be radiated away (Yahil 1983). Our implementation, described briefly here, follows the work of Omukai et al. (2005), and we direct the interested reader to this paper for further details. The evolution of the gas density, ρ, is given by

$$\begin{equation} \frac{d\rho }{dt} = \frac{\rho }{t_{\rm col}}, \end{equation} \tag{ 108 }$$

where the collapse timescale, t_col, is

$$\begin{equation} t_{\rm col} = \frac{t_{\rm dyn}}{\sqrt{1 - f}}. \end{equation} \tag{ 109 }$$

In this equation, t_dyn is the dynamical time for the collapse of a spherical cloud, and is expressed as

$$\begin{equation} t_{\rm dyn} = \sqrt{\frac{3 \pi }{32 G \rho }}. \end{equation} \tag{ 110 }$$

The collapse timescale is altered from the dynamical time by a factor $1/\sqrt{1-f}$ in Equation (109), which is an approximation of the ratio of the gas pressure to the force of gravity. The value of f depends on the effective adiabatic index, γ_ef ≡ (∂ln p/∂ln ρ), which we linearly extrapolate from derivative values at the two previous timesteps. For the value of f in this test problem, we use the piecewise function of Omukai et al. (2005), given by

$$\begin{equation} f \,{=}\, \left\lbrace \begin{array}{@{}ll}0, & \gamma _{\rm ef} < 0.83,\\ 0.6 + 2.5 (\gamma _{\rm ef} - 1) - 6.0 (\gamma _{\rm ef} - 1)^{2}, & 0.83 < \gamma _{\rm ef} < 1,\\ 1.0 + 0.2 (\gamma _{\rm ef} - 4/3) - 2.9 (\gamma _{\rm ef} - 4/3)^{2}, & \gamma _{\rm ef} > 1. \end{array} \right. \end{equation} \tag{ 111 }$$

The specific energy evolves as

$$\begin{equation} \frac{de}{dt} = -p \frac{d}{dt} \frac{1}{\rho } - \Lambda, \end{equation} \tag{ 112 }$$

where Λ is the cooling rate in units of erg s⁻¹ g⁻¹ and energy, temperature, density, and pressure are related by the ideal gas law, including effects from molecular hydrogen as appropriate. Figure 16 shows an example of the one-zone collapse test performed with an initial number density of 1 hydrogen atom per cm⁻³ and an initial temperature of 100 K using the 12 species chemistry network with H, D, and He species and metal cooling rates calculated with the Cloudy code. The effects of metal cooling can be clearly seen; as the metallicity increases from zero to 10⁻² Z_☉, the gas rapidly and significantly deviates from the primordial result (black line). Our primordial results compare very well to those shown in Omukai et al. (2005); however, we use a different metal cooling method, so the lines describing the evolution of the metal-enriched gas are not directly comparable.

Zoom In Zoom Out Reset image size

11.2.13. Ion Balance Test

The one-zone collapse problem, described in Section 11.2.12, is a good test of the chemistry solver at low temperatures where it is difficult to calculate an analytical solution for the equilibrium fraction of each species. However, for a gas of primordial composition at high temperatures, a simple analytical solution exists for the ionization balance of H and He species. The equations balancing ionization and recombination for all ions of H and He result in six equations with six unknowns (e.g., Katz et al. 1996), are given by

$$\begin{equation} n_{{\rm H\,\scriptsize{I}}} = \frac{n_{\rm H} \alpha _{{\rm H}\,\scriptsize{II}}}{\alpha _{{\rm H}\,\scriptsize{II}} + \Gamma _{{\rm e}{\rm H\,\scriptsize{I}}} + \Gamma _{\gamma {\rm H\,\scriptsize{I}}} / n_{\rm e}}, \end{equation} \tag{ 113 }$$

$$\begin{equation} n_{{\rm H}\,\scriptsize{II}} = n_{\rm H} - n_{{\rm H\,\scriptsize{I}}}, \end{equation} \tag{ 114 }$$

$$\begin{equation} n_{{\rm He\,\scriptsize{II}}} = \frac{n_{\rm He}}{1 + \frac{\alpha _{{\rm He\,\scriptsize{II}}} + \alpha _{\rm d}}{\Gamma _{{\rm e{\rm He\,\scriptsize{I}}}} + \Gamma _{\gamma {\rm He\,\scriptsize{I}}} / n_{\rm e}} + \frac{\Gamma _{{\rm e}{\rm He\,\scriptsize{II}}} + \Gamma _{\gamma {\rm He\,\scriptsize{II}}} / n_{\rm e}}{\alpha _{{\rm He\,\scriptsize{III}}}}}, \end{equation} \tag{ 115 }$$

$$\begin{equation} n_{{\rm He\,\scriptsize{I}}} = n_{{\rm He\,\scriptsize{II}}} \frac{\alpha _{{\rm He\,\scriptsize{II}}} + \alpha _{\rm d}}{\Gamma _{{\rm e}{\rm He\,\scriptsize{I}}} + \Gamma _{\gamma {\rm He\,\scriptsize{I}}} / n_{\rm e}}, \end{equation} \tag{ 116 }$$

$$\begin{equation} n_{{\rm He\,\scriptsize{III}}} = n_{{\rm He\,\scriptsize{II}}} \frac{\Gamma _{{\rm e}{\rm He\,\scriptsize{II}}} + \Gamma _{ \gamma {\rm He\,\scriptsize{II}}} / n_{\rm e}}{\alpha _{{\rm He\,\scriptsize{III}}}}, \end{equation} \tag{ 117 }$$

$$\begin{equation} n_{\rm e} = n_{{\rm H}\,\scriptsize{II}} + n_{{\rm He\,\scriptsize{II}}} + 2 n_{{\rm He\,\scriptsize{III}}}, \end{equation} \tag{ 118 }$$

where n_H and n_He are the total number densities of all H and He ions, α_i is the recombination rate of species i, and Γ_ij is the ionization rate of species j by i, where i is either an electron for collisional ionization or a photon for photo-ionization.

In this test, we initialize a one-dimensional grid at a constant number density of 1 cm⁻³ with smoothly varying temperature from 10⁴ K to 10⁹ K. The gas is initially fully ionized and we then iterate the chemistry equations without allowing the gas to cool until convergence is reached. In Figure 17, we show the equilibrium ionization fractions for H and He ions and the resulting cooling rate. We assume no photo-ionization for this test. The results from this test are shown by the solid lines with the analytic solution described by Equations (113)–(118) using the rate equations in Enzo, referenced in Section 6.1, overplotted with the dashed lines and an analytical solution using the rates of Cen (1992), shown by the dotted lines. As can be seen, the Enzo chemistry solver and corresponding analytical solution are in perfect agreement. The results largely agree with those of Cen (1992), save for H i, He i, and He ii at high temperatures, where their ion fractions are quite low. The slightly higher abundance of He i in Cen (1992) results in a higher cooling rate by about a factor of two near T ∼ 10⁵ K due to the enhanced collisional excitation cooling of He i.

Zoom In Zoom Out Reset image size

11.2.14. Photo-evaporation of a Dense Clump

The photo-evaporation of dense clumps of gas is prevalent in radiation hydrodynamics simulations, and this test problem examines the ionization front propagation into a dense clump, shadowing effects behind the clump, and the hydrodynamic response on the clump from photo-heating, all using Enzo's Moray ray-tracing module. The problem setup is the same as Test 7 in the Cosmological Radiative Transfer Comparison Project (Iliev et al. 2009) and Wise & Abel (2011). The simulation domain is 6.6 kpc on a side with an ambient medium of pure neutral hydrogen of density n_H = 2 × 10⁻⁴ cm⁻³ and temperature T = 8000 K. We place a spherical overdensity in hydrostatic equilibrium with the ambient medium. It has a radius r = 0.8 kpc, hydrogen density n_H = 0.04 cm⁻³ (i.e., overdensity of 200), and temperature T = 40 K, and is centered at (x, y, z) = (5, 3.3, 3.3) kpc. In Iliev et al. (2009) all of the codes used a fixed 128³ grid to ease the comparison, but in this test to demonstrate a higher resolution AMR solution, we employ a 128³ grid with two additional levels of refinement by factors of two for cells with a baryon mass greater than 1.5 (method 2 in Section 3.6). This test is run for 15 Myr.

The cloud is subject to radiation from a point source at the center of the x = 0 boundary with an ionizing photon luminosity $\dot{N}_\gamma = 3 \times 10^{51}$ photons s⁻¹, corresponding to a flux F₀ = 10⁶ photons s⁻¹ cm⁻² at the clump surface closest to the radiation source. The radiation source has a spectrum of a T = 10⁵ K blackbody, and we use four energy groups with the following mean energies and relative luminosities: E_i = (17.98, 31.15, 49.09, 76.98) eV, L_i/L = (0.23, 0.36, 0.24, 0.06) that are optimized to reduce errors in the solutions with a full spectrum and energy discretization (Mirocha et al. 2012). (Note that this choice of energy groups is different from those used in Wise & Abel 2011.) We use a minimum angular resolution of 10 rays per cell and a constant radiative transfer timestep of 25 kyr.

The left panel in Figure 18 shows cuts of the density, temperature, neutral fraction, and pressure in a line connecting the source and the clump center at t = 1, 5, 10, 15 Myr. At 1 Myr, the ionization front has propagated through the leftmost 500 pc of the clump. This heated gas is now overpressurized, as seen in the pressure plot of Figure 18 that expands into the ambient medium. This photo-evaporative flow becomes apparent at 5 Myr with the increased density and pressure on the left side of the clump. The right panel in Figure 18 shows the volume-weighted distributions of the temperature and flow Mach number in the entire domain at t = 10 Myr. Overall, our results agree well with the codes presented in Iliev et al. (2009) and the uniform grid test of Enzo+Moray (Wise & Abel 2011).

Zoom In Zoom Out Reset image size

However, there are some small differences in the line cuts and distributions that arise from our choice of energy discretization that better samples the high-energy tail of a T = 10⁵ K blackbody spectrum and does not signify possible shortcomings of Enzo+Moray. First, in the temperature cut at t = 10 Myr, the high-energy radiation photo-heat the gas to T = 10⁴ K, but it remains primarily neutral. The temperature dip at x ≃ 0.88 is the gas has not come into thermal equilibrium with the incident radiation. A stronger temperature dip exists in the Iliev et al. results and varies from 1000 to 8000 K between codes. Another minor difference is the larger breadth in the temperature distribution feature centered at T = 4.5 × 10⁴ K that corresponds to the ambient medium. This difference also originates from the different sampling in the spectral energy distribution of the source. Figure 19 depicts the same quantities in a two-dimensional slice through the clump center at t = 10 Myr. At this time, it is apparent that the outer layers are expanding after being photo-heated. These slices also show the sharp shadowing effects of the dense clump in the neutral fraction plot that is representative of ray tracing techniques.

Zoom In Zoom Out Reset image size

11.2.15. Cosmological I-front Propagation

This test problem examines Enzo's flux-limited diffusion radiation transport and associated ionization chemistry solvers on a cosmological I-front test. This corresponds to test problem 4.6 from Reynolds et al. (2009), run using the cosmological deceleration parameter q₀ = 0.05 and initial redshift z_i = 10. In this test problem, the physics of interest is the expansion of an H ii region in a uniform medium around a single monochromatic ionizing source (with frequency hν = 13.6 eV). The I-front propagates rapidly at first, approaching 90% of the Strömgren radius by a scaled redshift of −log [(1 + z)/(1 + z_i)] ≈ 0.15, before cosmological expansion overtakes ionization, pushing the Strömgren radius outward faster than the I-front can propagate. The Strömgren radius is analytically given by

$$r_S(t) = \left[\frac{3\dot{N}_{\gamma }}{4\pi \alpha _{\rm B} n_{\rm H}(t)^2}\right]^{1/3},$$

where $\dot{N}_{\gamma } = 5 \times 10^{48}$ photon sec⁻¹ is the strength of the ionizing source, α_B = 2.52 × 10⁻¹³ s⁻¹ is the H ii recombination rate, and n_H(t) is the proper number density of hydrogen. We then define λ = α_B n_{H, i} / H₀ /(1 + z_i), where the subscript i indicates quantities at the initial redshift. The analytical solution for the I-front position as a function of time is then

$$\begin{eqnarray*} &&r_I(t) \;= r_{S,i} \left[\!\lambda e^{-\tau (a)} \displaystyle\int _1^{a(t)} e^{\tau (\tilde{a})} \left(\! 1 - 2q_0 + \displaystyle\frac{2q_0(1+z_i)}{\tilde{a}} \!\right)^{-1/2}d\tilde{a} \!\right]^{1/3}, \end{eqnarray*}$$

where

$$\begin{eqnarray*} \tau (a) &=& \lambda \displaystyle\frac{6q_0^2(1+z_i)^2}{F(a)-F(1)}, \\ F(a) &=& \left[2-4q_0 - \displaystyle\frac{2q_o(1+z_i)}{a}\right] \left[1-2q_0 + \displaystyle\frac{2q_0(1+z_i)}{a}\right]^{1/2}, \end{eqnarray*}$$

and a(t) is the cosmological expansion coefficient.

The simulation parameters are: box size L_i ≈ 27 kpc, Hubble constant H₀ = 100 km s⁻¹ Mpc⁻¹, and cosmological parameters Ω_m = 0.1, $\Omega _\Lambda =0$, and Ω_b = 0.1. The initial values for the simulation are a radiation energy density E_{r, i} = 10⁻³⁵ erg cm⁻³, temperature T_i = 10⁴ K, density ρ_{b, i} = 2.35 × 10⁻²⁸ g cm⁻³, and an ionized hydrogen fraction $n_{{\rm H}\,\scriptsize{II}}/n_{\rm H} = 0$.

In Figure 20, we plot profiles of the H i and H ii fractions at redshifts z = {6.24, 2.29, 1.02}, as well as the ratio of the I-front and Strömgren radii throughout the simulation. As seen in the left panel (species fractions as a function of radius), the I-front is initially quite narrow, but slowly becomes wider as the simulation proceeds due to the diffusion approximation used by the radiative transfer scheme. However, as seen in the I-front propagation history plot (right panel), the computed I-front location very accurately matches the exact value.

Zoom In Zoom Out Reset image size

11.2.16. Anisotropic Thermal Conduction

Figure 21 shows a test that demonstrates the correct behavior of anisotropic thermal conduction in Enzo. We initialize a two-dimensional, 256 × 256 cell simulation having a physical scale of 1 kpc on a side, a uniform density of 1 proton cm⁻³ and a background temperature of 10⁶ Kelvin. Magnetic fields with a strength of B₀ = 1 μG are initialized such that the field lines form circles around the center of the simulation volume, such that B_x = −B₀sin (θ) and B_y = B₀cos (θ), where θ is the angle measured from the +x direction in a counterclockwise manner. A Gaussian temperature pulse is injected at (0.75, 0.5) (in units of the box size), with a peak temperature of 10⁸ K and a FWHM of 1/64 of the box size. This initial setup is shown in the left panel of Figure 21. The simulation is then allowed to evolve with only anisotropic conduction turned on (e.g., no hydrodynamics, radiative cooling, or cosmological expansion), and with a Spitzer fraction of f_sp = 1.

Zoom In Zoom Out Reset image size

The right panel of Figure 21 shows the state of the simulation after 300 Myr. Heat has clearly been transported only along field lines—there has been no diffusion perpendicular to the magnetic field setup, which is critical for many studies involving anisotropic thermal conduction. No oscillations are seen in the temperature field in regions where the fields are not aligned with the grid, suggesting that the flux-limiter is operating as expected (see discussion in Section 8.1).

The current version of Enzo has been parallelized for distributed memory platforms using message passing interface (MPI). This is done using a single grid object as the basic unit of parallelization. Each grid object—including all cell and grid data—is fully contained on a single processor.²⁹ Parallelization is accomplished by distributing grids amongst processors. This is done on the root grid using a simple tiling system, where the root grid is split up into N_root tiles, with N_root typically equal to the number of processors, N_p.

Load balancing on levels other than the root level (i.e., grids for which the level l > 0) is different, as the refined patches are not generally uniformly distributed. Grids on refined levels are first placed on the same processor as their parent to minimize communication; however, this generally does not result in a well-balanced computational load. Therefore, the code has a number of options for load balancing the grids on a given level. Each grid is assigned an estimated computational load, generally equal to the number of cells in the grid (which, empirically, is a good estimate of computational cost). The first load-balancing option is to move a grid from the processor with the highest computational load to the processor with the lowest load, with the proviso that only grids with load factors less than half the difference between the highest and lowest loaded processors will be moved. This continues until the load ratio between the most-to-least loaded processor is below 1.05 or until no suitable grid can be found to transfer. A second load balancing option uses a space-filling Hilbert curve to order the grids by their approximate spatial position. Then, once the grids have a specific one-dimensional ordering, we can divide up the grids into N_p groups (with the division taking place as equally as possible). Load balancing is done separately for each level. Clearly load balancing is most successful if there are significantly more grids than processors; however, small grids are less efficient because of the large number of ghost zones compared to active zones, and so the code uses a simple heuristic in order to split up grids until there are of order 10 grids per processor. This generally results in good load-balancing while not producing grids that are wastefully small.

Communication between processors is done using a non-blocking communication strategy that allows overlap of communication and computation. This can be done efficiently because each processor retains a copy of the entire hierarchy of grids, except that grids that do not "live" on a given processor only contain meta data (essentially just the location and size of the grid; such grids are denoted as "ghost" grids). Replicating the hierarchy means that all communication of data from one grid to another can be identified by each processor independently. The metadata for "ghost" grids are quite small and so the extra memory required is generally not onerous unless very large numbers of grids are used (more than a few hundred thousand grids). A schematic of this distribution is shown in Figure 1.

Data is transferred through a three-step procedure that takes advantage of the capabilities of the MPI library: (1) as the code progresses and data is needed from another grid on another processor, the receiving processor posts an MPI non-blocking receive indicating that it is expecting data, this outstanding receive is recorded in a table, (2) the sending processor calls the MPI non-blocking send function, and then finally (3) the receiving processor, after it has carried out all the computation it can, waits for any MPI message to arrive. Each message is coded so that it can be matched with the appropriate receive posted in the first step, and based on that, the appropriate routine is called to processes the data. Step (3) is repeated until there are no outstanding receives.

12.2.1. Performance Measurement and Instrumentation

12.2.2. Unigrid Scaling

It is advantageous to use Enzo in its "unigrid" (i.e., non-adaptive mesh) mode for a variety of problems, including non-gravitating turbulence (e.g., Kritsuk & Norman 2002, 2004; Kritsuk et al. 2007, 2009b), the Lyα forest (Jena et al. 2005; Paschos et al. 2009), or feedback of metal-enriched gas from galaxies into the IGM (Norman et al. 2004; Smith et al. 2011). Achieving good scaling of the code in unigrid mode is relatively straightforward—upon initialization, unigrid simulations are decomposed such that each MPI process has a roughly equal subvolume (and thus an equal number of grid cells), meaning that work is evenly distributed among computational elements. Communication patterns for both the gravity solve (which uses an FFT) and the fluid solves (which transfer boundary information between subvolumes) are predictable and straightforward, and rebuilding of the grid hierarchy does not take place, removing a substantial global operation and a great deal of communication.

Figure 22 shows Enzo weak scaling results for a sequence of scaled unigrid Lyα forest calculations. These calculations include dark matter dynamics, hydrodynamics using the PPM solver, six-species non-equilibrium chemistry and radiative cooling, and a uniform metagalactic ultraviolet background. In this sequence of test calculations, we perform a weak scaling test on up to 13,824 MPI tasks on the NICS Kraken XT4 and ORNL Jaguar XT4 supercomputers.³⁰ In this test, each MPI task was given a 128³ root grid tile (i.e., 128³ grid cells containing baryon quantities) and, initially, approximately 128³ dark matter particles. The number of grid cells on each processor was constant throughout the calculation; the number of dark matter particles on each processor varies as they are moved from subvolume to subvolume as structure evolves. The grid resolution was kept at a constant comoving size of ≃ 40 kpc h⁻¹, and as the core count was increased, so was the simulation volume. On each machine, a compute node contained a single AMD Opteron quad-core chip (2.1 GHz on Jaguar; 2.3 GHz on Kraken) with 2 GB/memory per core (8 GB/total per node). Both machines used the SeaStar2 interconnect. In the scaling study, calculations were run with 1, 2, or 4 MPI tasks per node. The figure shows cell updates per second per MPI process; perfect scaling would be a horizontal line.

Zoom In Zoom Out Reset image size

As can be seen in Figure 22, the unigrid weak scaling performance of the code is extremely good for this problem, with only a 20% decrease in cell updates per second per MPI task as the code is scaled from 8 to 4096 MPI tasks, and a 40% decrease in performance overall going from 8 to 13,824 (or 24³) MPI tasks. We speculate that this decrease is likely to be partially due to global MPI communications used to, e.g., calculate the overall timestep of the simulation, and also likely due to load imbalances due to increasing cosmological power (and thus an increasingly uneven distribution of dark matter particles between MPI tasks at late times) as the simulation volume grows. We also observe that a systematic difference in speed can be seen between the two machines, which can be attributed primarily to the slightly faster CPUs on Kraken at the time (2.3 GHz, versus 2.1 GHz on Jaguar). The difference in speed when using different numbers of MPI tasks per node can be attributed primarily due to differences in competing usage of shared cache on the quad-core chips used on this machine.

Broadly, excellent scaling in Enzo's unigrid mode is seen for a variety of problems as long as each compute core is given an adequate amount of work to do. For cosmological simulations, this value has been empirically determined to be roughly 128³ cells per core. If fewer cells per core are used, the CPU is essentially data-starved and poor scaling is observed due to computing units being idle while waiting for information to be communicated from other processes (for, e.g., boundary information or gravity solves). Substantially larger cell counts per core would in principle help scaling by reducing the amount of inter-process communication needed.

As a final point, we observe that scaling at larger core counts has been measured, but only with an experimental hybrid-parallel (MPI + OpenMP) version of Enzo. Using this version, scaling comparable to that shown in Figure 22 was seen on up to 98,304 cores on the NICS Kraken XT5 (an upgraded version of the XT4 machine used for the scaling study shown in the figure), using 2–8 OpenMP threads per MPI process. Hybrid parallelism has the potential to substantially improve scaling by reducing the amount of communication per grid tile, as described in the previous paragraph.

12.2.3. AMR Scaling

Many astrophysical problems, such as cosmological galaxy formation (Hummels & Bryan 2012), high-resolution disk galaxy simulations (Kim et al. 2011), high-redshift (Turk et al. 2009), and present-day star formation (Collins et al. 2012), involve multi-scale physics that span several orders of magnitude in both space and time. In these situations, using Enzo in its AMR mode is beneficial. Because the adaptive grid hierarchy is dynamic, grid boundaries and, thus, communication patterns can be unpredictable, hindering strong scaling to high core counts.

Figure 23 shows strong scaling results from a single 50 Mpc h⁻¹ cosmology simulation run on N_core compute cores, ranging from 128 to 16,384 cores at power-of-two intervals. It was run on the NICS Kraken XT5 supercomputer, which has two AMD Opteron hexa-core processors with 16 GB of memory per compute node. The simulations that utilized 128, 256, and 512 cores were executed on 128 nodes because of the memory requirements. The higher core count simulations were run with 8 cores per node. This simulation would not run with 12 cores per node because of the memory overhead associated with the grid hierarchy being duplicated on each MPI process. However, this overhead is greatly diminished if a hybrid-parallel (MPI + OpenMP) approach is used.

Zoom In Zoom Out Reset image size

This simulation includes dark matter dynamics, hydrodynamics using the PPM, six-species non-equilibrium chemistry and radiative cooling, and a uniform metagalactic ultraviolet background. The simulation uses the space-filling curve method for load balancing the AMR grids. It has a 512³ root grid that is divided into 512 tiles, and 6 additional AMR levels are used. We perform these scaling tests when the simulation reaches z = 4, where the AMR grid hierarchy is well-developed and is thus a reasonable representation of AMR simulation behavior in general. The results shown in Figure 23 come from a single root-level timestep of Δt = 2.1 Myr. At this time, there are 3.04 × 10⁵ AMR grids, 9.03 × 10⁸ (∼967³) computational cells, and 1.34 × 10⁸ dark matter particles in total. The breakdown of the number of AMR grids, cells, timesteps, and number of cell updates on each level is shown in Table 5.

Table 5. Strong Scaling Test Computational Details

Level	Ngrid	Ncells	Nup	Ntimesteps
0	512	1.34 × 108	1.34 × 108	1
1	61,868	4.01 × 108	4.01 × 108	1
2	91,182	1.99 × 108	5.96 × 108	3
3	59,932	7.62 × 107	5.34 × 108	7
4	40,700	3.32 × 107	5.65 × 108	17
5	28,346	2.76 × 107	1.36 × 109	49
6	19,771	2.80 × 107	5.25 × 109	187
Total	302,311	9.03 × 108	8.83 × 109	⋅⋅⋅

Notes. Data shown at z = 4 for a root grid timestep of 2.1 Myr. Column 1: AMR Level. Column 2: Number of grids. Column 3: Number of computational cells. Column 4: Number of cell updates. Column 5: Number of timesteps.

Download table as:  ASCIITypeset image

The left panel in Figure 23 shows the computational and communication time spent on each level. In the AMR levels, there exists good strong scaling up to 2048 cores, and marginal speed-ups are found at 4096 cores. On the root-grid level, there exists good scaling up to 1024 cores, but the performance dramatically decreases at higher core counts. This occurs because the root grid is not re-partitioned into N_core tiles when the simulation is restarted with a different core count. This feature can be easily implemented and is planned in the next major release of Enzo, where scaling results would be similar to the weak scaling shown in Section 12.2.2. The right panel in Figure 23 shows the time spent in some representative major routines in Enzo. The local physics routines, for example SolveHydroEquations, exhibit perfect strong scaling because they involve no communication. By investigating the scaling behavior in each routine, it is clear that the communication in the SetBoundaryConditions, SolveForPotential (the multi-grid solver in AMR subgrids), and RebuildHierarchy are responsible for the lack of strong scaling at N_core ≳ 4096 in this particular simulation. The transpose of the root grid tiles are responsible for the performance decrease because it is not optimized for situations where the number of tiles is greater than the number of MPI processes. These results are directly applicable to simulations with similar computational demands. In simulations with fewer computational cells, strong scaling will cease at a smaller core count because the CPUs will become data-starved more quickly, and the opposite occurs with larger simulations.

12.2.4. An Approximate Time and Memory Scaling Model

In this section, we develop a simple, approximate model to estimate the computational time and memory required to complete a given calculation. Note that we are not trying to model parallel performance (which was discussed in the previous section), but have an even simpler goal: to determine how to scale computational time estimates as we increase the simulation resolution. This is a straightforward thing to do in codes that use static grids, as the computational effort per timestep is constant. However, for an AMR calculation, the number of grids that will be generated during the run is not known in advance, and so the CPU time per problem time can vary drastically throughout the calculation. For example, at the beginning of a cosmological simulation, when the densities are nearly uniform, only the static grid is required and the calculation progresses rapidly. However, as structure forms and dense clumps are generated, the number of grid points swells by orders of magnitude (an increase of 10³ is not uncommon) and most of the CPU time is consumed at late physical time in the simulation. Therefore, simply performing a few steps at the beginning of the calculation does not produce a good estimate of the required CPU time.

Nevertheless, we can try to determine how the compute time scales for a given run as we increase the resolution. First, we neglect particles and concentrate on the time taken by the grids, which can be justified both theoretically and empirically.³¹ Second, we break the problem down slightly, and examine the scaling over a short enough period of time that the grid structure does not change significantly (i.e., the number of grids at each level remains approximately constant). We then assume that the whole run scales in the same way, or in other words, that changes in the resolution affect the simulation in the same way at each timestep. This is usually a good approximation, as the most costly parts of the calculation typically don't change their characteristics substantially between timesteps.

To make progress, we assume that the computational cost to advance a single cell by one timestep is a constant. This can be incorrect if the chemistry solver requires many iterations, but is usually fairly accurate. For a unigrid calculation, the time would by proportional to $N_{\rm root}^4$ since the number of cells scales as $N_{\rm root}^3$ and, assuming the Courant condition is the factor that controls the timestep, the number of steps to advance the calculation over a given time interval is proportional to N_root. Therefore, to advance a hierarchy a given time interval, we find, accounting for all levels and using a refinement factor of two,

$$\begin{equation} t_{\rm SU} = C_1 \sum _l f_l N_{\rm root}^4 2^{4l} \end{equation} \tag{ 119 }$$

where C₁ is a constant, which can be thought of as the time taken to advance a single cell. The factor f_l is the fraction of the volume on a given level that is actually refined. By definition f₀ = 1, and f_l ⩽ 1. In writing down this equation, we neglect a number of costs that are not directly proportional to cell count, including communication between processors, the log N factor for the root grid FFT, cache misses, optimizations, and other costs associated with processing the hierarchy. The first item on this list, in particular, is clearly important for large processor counts; however, we neglect parallel considerations in this section.

Note that unless f_l/f_{l + 1} < 1/16 (i.e., if less than about 6% of a given level is further refined), the cost per level will increase with level. A key question, therefore, is the value of f_l for each level. Unfortunately, this depends strongly on the simulation being run. In Figure 24, we show the values of f_l for three simulations: two cosmological runs with varying box size, and a third simulation focusing on a single disk galaxy.

Zoom In Zoom Out Reset image size

As, the figure demonstrates, the f_l values all show a sharp decline with level, dropping as power-laws with the level l. Focusing first on the "refine-everywhere" cosmological runs, which both show similar behavior despite the different box sizes and redshifts at which we collect the data. In fact, we find these results are very typical for cosmological runs, with only slight variations depending on the exact problem that is being simulated. The top-right panel shows the ratio of f_l/f_{l + 1} which, remarkably, hovers around the critical value of 0.0625 determined earlier. This implies that the relative CPU usage of each level, shown in the lower-left panel, is nearly flat, and that going to additional levels of refinement only increases the amount of calculation by a factor of roughly 1/l_max. Even in the most extreme case, the L80 run, the last level adds only 25% to the computation time.

The memory used by each level, on the other hand, scales as

$$\begin{equation} {\rm memory} = C_2 \sum _l f_l N_{\rm root}^3 2^{3l} \end{equation} \tag{ 120 }$$

and the terms in this sum are shown in the bottom-right panel of Figure 24. This demonstrates that the memory usage is dominated by the top three levels, and adding additional levels adds minimally to the memory usage of the run. We have empirically tested this scaling for cosmological simulations and found it to be reasonably accurate (although the increase is typically slightly larger than found here, usually 30%–50%, probably because of less-than-ideal parallel scaling).

The physical reason for the f_{l + 1}/f_l ratios found in the cosmological simulations appears to be due to a combination of the density structure of individual clouds and to the distribution of the clouds themselves. In particular, we note that for a ρ∝r⁻² density profile, the Lagrangian refinement criteria typically used in such simulations produces an f_{l + 1}/f_l ratio of 1/8 for a single, resolved cloud. Of course, at some point we would resolve the flat density of individual clouds and the ratio would climb, making further refinement more costly; however, we are not yet in this regime in this example.

These conclusions are, however, highly dependent on the type of run being performed. We contrast this cosmological case to the simulation of a galaxy formation run, as shown by the dot-dashed curves in Figure 24. Note that grid levels 1 and 2 in this simulation are statically refined to ensure that the Lagrangian volume of the galaxy halo is resolved to at least level 2 at all times. More importantly, the f_{l + 1}/f_l ratio of the highest levels are systematically slightly larger than the critical 0.0625 value, indicating that the CPU time is dominated by the most refined levels, as shown in the bottom-left panel. The memory usage is dominated by levels 2 and 3, the static levels, as the bottom-right panel demonstrates. In this case, adding more levels of refinement would substantially increase the cost of the simulation.

We have focused thus far on scaling when increasing l_max, but without changing the refinement criteria. However, we can also keep the levels fixed and modify the refinement criteria. For cosmological runs, this is typically done by increasing the mass resolution (and simultaneously decreasing the dark matter particle mass), and adding additional linear small-scale modes to the initial conditions. The effect of this is to boost the f_l values at all levels (except for the root grid where f₀ is already 1) by a factor linearly proportional to the increase in the mass resolution. This is because the refinement criteria typically boosts the number of refined cells on the first and subsequent levels by this factor. Again, we have empirically tested that this scaling is approximately valid, provided that we keep the maximum level constant.

Putting this together, we find an approximate scaling for the computational time required for cosmological simulations:

$$\begin{equation} t_{\rm SU} \propto M_{\rm res}^{-1} l_{\rm max}. \end{equation} \tag{ 121 }$$

This is approximately accurate, and gives users some way to estimate compute times for Enzo calculations.

Many of today's leading supercomputers use a heterogeneous computing platform: on a single node of a distributed-memory platform, a multi-core CPU is often paired with one or more many-core accelerators. One programming model that has been shown to successfully take advantage of these hardware accelerators is to run the serial component of the algorithm on the CPU, and the vector-parallel part of the algorithm on the many-core hardware accelerator. Due to the much higher computational performance of the many-core accelerators, if a code can be ported to effectively utilize this heterogeneous architecture, massive speedup in simulation performance may be achieved. This heterogeneous trend is likely to continue in future supercomputers, due to the relatively low energy needs of hardware accelerators as well as other factors. Thus, to ensure Enzo's efficiency on future high-end computational platforms, some of the most time-consuming parts of Enzo have been ported to many-core architecture. Since GPUs are currently the most popular many-core accelerator, NVIDIA's CUDA C/C++ programming technology was chosen to port some Enzo modules to NVIDIA GPUs.

Thus far, the PPM and Dedner MHD solvers (described in Sections 4.1 and 4.2, respectively), have been ported to take advantage of GPUs using CUDA. Since the CPU and GPU on a node have access to separate pools of RAM, fields updated by other modules will be transferred to the GPU before calling the GPU fluid solver. Correspondingly, after calling the GPU solver, the updated fields will be transferred back to CPU. This ensures that the GPU-parallelized solvers can work correctly with other parts of the code that have not yet been ported to the GPU, as well as with the communication infrastructure within Enzo. In addition, the GPU solver supports SAMR, which is not a trivial task. As discussed in Section 3, one of the key steps in Enzo's AMR implementation is flux correction, which is required when each level of resolution is allowed to take its own time step. In the GPU version of the solvers, the fluxes are calculated on the GPU and only the fluxes required for flux correction are transferred back to the CPU. This reduces data transfer overhead, which can be substantial in a heterogeneous architecture of this sort.

The key step in porting to many-core architectures such as the GPU is exposing the massive parallelism inherent in the algorithm. Due to the explicit, directionally-split stencil pattern of both the PPM and Dedner MHD solvers, they are inherently massively parallel and thus should be a good fit for hardware acceleration. Both solvers essentially contain two parts—computation of fluxes, and a cell update. When computing fluxes, the basic procedure is to compute the flux at each cell-interface given the inputs from neighboring cells. In the cell update part, the fundamental computation is updating the cell-centered values using the previously-computed fluxes. In the CPU serial code, a loop over the grid is used, where the loop body contains these basic computations. In both parts, algorithmically-different loop iterations are completely independent of each other. Thus, the natural parallelization scheme is to map one GPU thread to one iteration of the loop. However, the original CPU code, which operates on each grid as a serial process, does not completely expose this parallelism as some small temporary arrays are reused among loop iterations. This re-use of arrays introduces data dependency among loop iterations, which is undesirable for GPU parallelization. Because of this re-use of arrays, the main change in porting the serial CPU solvers to allow massively parallel GPU computation was replacing those shared temporary arrays by larger temporary arrays that are not shared among loop iterations. This change exposed the massive parallelism in the algorithm, which could then be accelerated in a straightforward manner using CUDA.

To illustrate the speedup provided by porting two solvers to GPUs, we show the results of a weak scaling test of driven MHD turbulence on a uniform mesh in Figure 25. This particular problem type contains no physics other than the equations of ideal MHD, and thus is representative of the type of calculation than can currently benefit from the GPU-optimized solvers in Enzo. In this scaling test, we use the Dedner MHD solver, which runs on both CPUs and GPUs and has been tuned to maximize performance on both platforms. The benchmarking platform is a Cray XK6 system with a Gemini interconnect. Each node has a 16-core AMD Opteron 6272 CPU and a single NVIDIA Tesla K20 GPU, which has 2496 CUDA cores and a maximum theoretical speed of 1.17 Tflops for double-precision computations. For the CPU run, all the 16 cores in each node are used by launching 16 MPI processes per node. For the GPU run, 16 MPI processes are also launched per node, and all processes on this node share the single GPU on that node. This can work because NVIDIA's MPS (Multi-Process Service) technology allows multiple processes to concurrently use the same GPU. A series of weak scaling tests were run, varying from 1 to 8 nodes (16 to 128 MPI processes), with each node containing a cube of 256³ cells (so, a 4-node computation would have a domain that is 512 × 512 × 256 cells). Figure 25 displays the scaling results in terms of cell updates per second per node—as a result, ideal weak scaling is a horizontal line. In general, the simulations using the GPU-accelerated MHD solver perform roughly 5 times better overall on this system than the equivalent CPU-only calculation run on the same number of nodes. We note, however, that this somewhat mis-represents the speedup obtained by porting the MHD code to GPUs—in the GPU simulation case, the 16 CPU cores on each node are mostly idle, as they are used only for timestep calculation and boundary condition transmission (which are both computationally inexpensive compared to the fluid solve). This means that, effectively, the GPU simulation is 80 times faster than a single CPU core, so a user would likely see much greater effective speedup on a system where the ratio of CPU cores to GPUs is lower.

Zoom In Zoom Out Reset image size

This interpolation algorithm is generally second-order, but has monotonicity constraints as described below. In one dimension, we define the parent values as Q₋₁, Q₀, and Q₊₁, where the central parent cell has a left edge at x₀ and width Δx. We first linearly interpolate the parent values to the cell edges: Q_−1/2 = (Q₀ + Q₋₁)/2, and similarly for Q_1/2 = (Q₀ + Q₁)/2, and then compute a monotonic slope: ΔQ₀ = minmod(Q_1/2 − Q₀, Q₀ − Q_−1/2) where

$$\begin{equation} {\rm minmod} \left({a},{b} \right) = \left\lbrace \begin{array}{@{}ll}0 & {\rm if} \quad ab < 0 \\ \min { \left(| a |, | b | \right) } {\rm sign} (a) & {\rm otherwise} \end{array}\right. \end{equation} \tag{ A1 }$$

This slope is then used to compute the interpolated subgrid values. Defining q_i to be the subgrid values at cell centers for the r cells corresponding to parent cell Q₀ (for refinement factor r), we can write:

$$\begin{equation} q_i = Q_0 + \frac{i+(1-r)/2}{r} f_x \end{equation} \tag{ A2 }$$

where f_x = 2ΔQ₀ (this notation is used for consistency with the two- and three-dimensional cases).

In two dimensions, the procedure is very similar in that we linearly interpolate parent values to the four cell corners: Q_{−1/2, −1/2}, Q_{−1/2, 1/2}, Q_{1/2, −1/2}, Q_{1/2, 1/2} (by averaging parent cell-centered parent values). We then compute monotonic slopes across the two diagonals, since for a linear function extrema occur at corners:

$$\begin{eqnarray} \Delta Q_0 & = & {\rm minmod} ({Q_{1/2, 1/2} - Q_0},{Q_0 - Q_{-1/2, -1/2}} ) \quad \end{eqnarray} \tag{ A3 }$$

$$\begin{eqnarray} \Delta Q_1 & = & {\rm minmod} ({Q_{-1/2, 1/2} - Q_0},{Q_0 - Q_{1/2, -1/2}} ). \quad \end{eqnarray} \tag{ A4 }$$

We then translate these slopes into the grid axes with f_x = ΔQ₀ − ΔQ₁ and f_y = ΔQ₀ + ΔQ₁ so that the interpolation itself can be written simply as

$$\begin{equation} q_{i,j} = Q_0 + \frac{i+(1-r)/2}{r} f_x + \frac{j+(1-r)/2}{r} f_y. \end{equation} \tag{ A5 }$$

Finally, we write down the three-dimensional version—unfortunately, here the number of monotonicity constraints is four (the four diagonals across the eight opposing corners of the cube), while the number of slopes is three, so the problem is over-constrained. Somewhat arbitrarily, we adopt the following procedure. As before, we compute ΔQ₀, ΔQ₁, ΔQ₂, and ΔQ₃ with the minmod limiter across the 4 diagonals. We then define

$$\begin{equation} s =\frac{\Delta Q_1 + \Delta Q_2 + \Delta Q_3}{\Delta Q_0} \end{equation} \tag{ A6 }$$

which is the value of ΔQ₀ that the other slopes (ΔQ₁, ΔQ₂, and ΔQ₃) imply, normalized by the desired value of ΔQ₀ itself. If 0 < s < 1, then no adjustment needs to be made, as the monotonicity constraints are met. If these equalities are not met, then, if s < 0 we define χ_n = 1 if ΔQ_n/ΔQ₀ < 0, and χ_n = 0 otherwise (if s > 1, this is reversed, so χ_n = 1 if ΔQ_n/ΔQ₀ > 0, and 0 otherwise). These weights are used to determine which slopes to modify to match the ΔQ₀ constraint. We then compute the amount of adjustment required:

$$\begin{equation} f = - \frac{\Delta Q_0 s^\prime + \sum _n (1-\chi _n) \Delta Q_n}{\sum \chi _n \Delta Q_n + \epsilon } \end{equation} \tag{ A7 }$$

where s' = min (max (s, 0), 1) and is a small number to prevent numerical errors. We then use this adjustment fraction to compute the new, adjusted ΔQ_n (for n = 1, 2, 3) that match the ΔQ₀ constraint as closely as possible with

$$\begin{equation} \Delta Q_n = f \chi _n \Delta Q_n + (1-\chi _n) \Delta Q_n. \end{equation} \tag{ A8 }$$

Finally, these are converted to the grid axes with f_x = ΔQ₂ + ΔQ₃, f_y = ΔQ₁ + ΔQ₃, f_z = ΔQ₁ + ΔQ₂. We then do the interpolation itself with

$$\begin{eqnarray} q_{i,j,k} &=& Q_0 + \frac{i+(1-r)/2}{r} f_x + \frac{j+(1-r)/2}{r} f_y \nonumber\\ && + \frac{k+(1-r)/2}{r} f_z. \end{eqnarray} \tag{ A9 }$$

We also provide a variant on the above procedure, with two changes. The first is that the slopes across the diagonals (ΔQ₀, etc.) are computed directly rather than with the minmod limiter (e.g., ΔQ₀ = (Q_{1/2, 1/2, 1/2} − Q_{−1/2, −1/2, −1/2})/2). To ensure positivity in the resulting interpolation when applied to positive conserved quantities, the slopes are limited so that the smallest corner value is 0.2 of the cell center value. The procedure described above is then applied to turn these four slopes into a linear interpolation.

This is a completely different second-order interpolation scheme, which is based on Cloud-In-Cell (CIC) interpolation (Hockney & Eastwood 1988). In one dimension, we define the parent values as Q₀, and Q₊₁, where the left parent cell has a cell center at x₀ and width Δx. Then, the interpolated value for a subgrid cell q_i with a cell left edges at x_i = x₀ + iΔx^p/r, where i runs from 0 to r − 1 for refinement factor r, is given by:

$$\begin{equation} q_i = \frac{2r - 1 - 2i}{2r} Q_0 + \frac{1+2i}{2r} Q_{+1}. \end{equation} \tag{ A10 }$$

The extension to two and three dimensions is straightforward, with weights in the other dimensions computed in the same way and then multiplied to get a total of four and eight weights for the two- and three-dimensional cases, respectively. This scheme preserves monotonicity, but is not conservative.

This interpolation method provides third-order accuracy based on the Triangular Shaped Cloud methodology (Hockney & Eastwood 1988). As usual, in one dimension, we define the parent values as Q₋₁, Q₀, and Q₊₁, where the central parent cell has a left edge at x₀ and width Δx. Then, the interpolated value for a subgrid cell q_i with a cell left edges at x_i = x₀ + iΔx^p/r where i runs from 0 to r − 1 for refinement factor r, is given by:

$$\begin{equation} q_i = a_i Q_{-1} + b_i Q_0 + c_i Q_{+1} \end{equation} \tag{ A11 }$$

and the weights are given by:

$$\begin{equation} a(i) = \frac{(r-i)^3 - (r-i-1)^3}{6r^3}; \quad c(i) = \frac{3i^2 + 3i + 1}{6r^3} \end{equation} \tag{ A12 }$$

and b(i) = 1/r − a(i) − c(i). The extension to two and three dimensions is straightforward with the weights computed in the same way as for the one-dimensional case but then multiplied to determine the 9 and 27 factors necessary for the two- and three-dimensional cases, respectively.

One problem with this interpolation technique is that it is not conservative. In particular, the sum of the interpolated subgrid values:

$$\begin{equation} \tilde{Q_0} = \sum _{i=0}^{r-1} q_i \end{equation} \tag{ A13 }$$

is not, in general, equal to Q₀. We can retain conservation by adding the factor $(Q_0 - \tilde{Q_0})/r$ to the interpolated values (or by multiplying the interpolated values by the ratio $Q_0/\tilde{Q_0}$). Unfortunately, the result of this procedure does not preserve monotonicity—it can introduce local minima and maxima at parent cell boundaries. We can then attempt to correct that by taking a weighted average between the interpolated values on the subgrid and the parent value, with weights computed such that the interpolated values at the edge of the parent cell are not local maxima compared to the interpolated values in the neighboring parent cell.

Finally, for completeness, we include a first-order accurate piecewise constant interpolator for which, using the same definitions as in the previous case, we take q_i = Q₀, for i = 0 to r − 1.