Enzo-E Methods
[ This page is under development ]
This section decribes all Methods available in Enzo-E: what they do, what method-specific parameters they have, what fields / particles they access or update, what solvers if any used, and how different Methods are intended to be used with each other.
Current Enzo-E methods available include those listed below.
"ppm": hydrodynamics
This implements the modified piecewise parabolic method (PPM) in Enzo.
parameters
Parameter |
Type |
Default |
Description |
|---|---|---|---|
|
logical |
false |
PPM diffusion parameter |
|
integer |
3 |
PPM flattening parameter |
|
integer |
100 |
Enzo’s MinimumPressureSupportParameter |
|
logical |
false |
Pressure-free flag |
|
logical |
false |
PPM steepening parameter |
|
logical |
false |
Minimum pressure support |
fields
Field |
Type |
Read/Write |
Description |
|---|---|---|---|
|
|
[rw] |
|
|
|
[rw] |
|
|
|
[rw] |
if rank ≥ 2 |
|
|
[rw] |
if rank ≥ 3 |
|
|
[rw] |
|
|
|
[r] |
if gravity |
|
|
[r] |
if gravity and rank ≥ 2 |
|
|
[r] |
if gravity and rank ≥ 3 |
|
|
[w] |
computed from |
fluid_props compatability
This method is also compatible with the "bryan95" dual-energy formalism.
See dual-energy formalism for additional details.
This method currently ignores all of the floor parameters that are set in the physics:fluid_props:floors section of the parameter file.
"ppml": MHD
PPML ideal MHD solver
"mhd_vlct": hydrodynamics/MHD
This implements the VL + CT (van Leer + Constrained Transport) unsplit
Godunov method described by Stone & Gardiner (2009)
. This solver operates in 2 modes (designated by the required
Method:mhd_vlct:mhd_choice parameter):
a pure hydrodynamical mode (it cannot handle magnetic fields)
a MHD mode.
The algorithm is a predictor-corrector scheme, with attributes similar to the MUSCL-Hancock method. For both modes, the method includes two main steps. First the values are integrated to the half time-step. Then, the values at the half time-step are used to caluclate the change in the quantities over the full time-step.
In the MHD mode, the algorithm is combined with the constrained transport method. The primary representation of the magnetic field, \(\vec{B}\), components are stored in 3 face-centered Cello fields. In more detail:
\(B_x\) is stored on the x-faces
\(B_y\) is stored on the y-faces
\(B_z\) is stored on the z-faces
The method also tracks components of the magnetic fields at the cell-centers, which just store the average of the values from the cell’s faces.
Currently, this method only support 3-dimensional problems. In the future, alternative modes supporting MHD could easily be implemented that use divergence-cleaning schemes instead of constrained transport.
parameters
Note that the courant factor (specified by the "courant"
parameter) should be less than 0.5.
Parameter |
Type |
Default |
Description |
|---|---|---|---|
|
string |
none |
Specifies handling of magnetic fields (or lack thereof) |
|
string |
hlld |
name of the Riemann Solver to use |
|
string |
nn |
Reconstruction method for half timestep |
|
string |
plm |
Reconstruction method for full timestep |
|
float |
1.5 |
controls dissipation of the “plm”/”plm_enzo” reconstruction method. |
fields
Field |
Type |
Read/Write |
Description |
|---|---|---|---|
|
|
[rw] |
|
|
|
[rw] |
|
|
|
[rw] |
|
|
|
[rw] |
|
|
|
[rw] |
specific total energy |
|
|
[rw] |
Cell-centered bfield (average of the corresponding |
|
|
[rw] |
Cell-centered bfield (average of the corresponding |
|
|
[rw] |
Cell-centered bfield (average of the corresponding |
|
|
[rw] |
Primary representation of x-component of bfield (lies on x-faces). |
|
|
[rw] |
Primary representation of y-component of bfield (lies on y-faces). |
|
|
[rw] |
Primary representation of z-component of bfield (lies on z-faces). |
|
|
[w] |
computed from |
|
|
[rw] |
if dual-energy |
In hydro-mode, none of the 6 fields used to store the magnetic field should be defined.
At initialization the face-centered magnetic field should be
divergence free. Trivial configurations (e.g. a constant magnetic
field everywhere) can be provided with the "value"
initializer. For non-trivial configurations, we have provide the
"vlct_bfield" initializer which can initialize the magnetic fields
(face-centered and cell-centered) from expression(s) given in the
parameter file for component(s) of the magnetic vector potential.
fluid_props compatability
This method makes use of the density and pressure floor parameters that are set in the physics:fluid_props:floors section of the parameter file.
See Floors for more details about specifying these parameters.
This method requires that both parameters are specified and that they have positive values.
This method is also compatible with the "modern" dual-energy formalism.
See dual-energy formalism for additional details.
reconstruction
This subsection details the available interpolation methods for reconstructing the left and right states of the cell-centered interfaces. Presently, all available methods perform reconstruction on cell-centered primitive quantites, \({\bf w} = (\rho, {\bf v}, p, {\bf B})\)
To simplify the determination of the necessary number of ghost zones for a given combination of reconstruction algorithms on a unigrid mesh, we define the concepts of “stale depth” and “staling rate”. We define a “stale” value as a value that needs to be refreshed. “Stale depth” indicates the number of field entries, starting from the outermost field values on a block, that the region encompassing “stale” values extends over. Every time quantities are updated over a (partial/full) timestep, the stale depth increases. We define the amount by which it increases as the “staling rate” (which depends on the choice of interpolation method).
For a unigrid simulation, the number of required ghost zones is given by the sum of the staling rates for each selected reconstruction method.
We provide the names used to specify each available method in the input file, the associated staling depth, and a brief description.
Name |
Staling Depth |
Description |
|---|---|---|
|
1 |
Nearest Neighbor - (1st order) reconstruction of primitives |
|
2 |
Piecwise Linear Method - (2nd order) reconstruction of
primitives using the slope limiter from Enzo’s Runge–Kutta
integrator. This is tuned by the |
|
2 |
Piecwise Linear Method - (2nd order) reconstruction of primitives using the slope limiter from Athena (& Athena++). For some primitive variable, \({\bf w}_{i}\), the limited slope is defined in terms of the left- and right-differences: \(\delta{\bf w}_{L,i}={\bf w}_{i}-{\bf w}_{i-1}\) and \(\delta{\bf w}_{R,i}={\bf w}_{i-1}-{\bf w}_{i+1}\). If the signs of the differences don’t match (or at least 1 is 0), then the limited slope is 0. Otherwise the limited slope is the harmonic mean of the differences. |
We provide a few notes about the choice of interpolator for this algorithm:
The recommended choices of reconstruction algorithms are
"nn"for the half-timestep and then piecewise-linear reconstruction for the full-timestep (most test problems have been run usingplmwiththeta_limiter=2, matching the integrator description in Stone & Gardiner 2009 ). Using"nn"both times also works, however tests show that errors arise when piecewise linear reconstruction is used both times.It is supposed to be possible to reconstruct the characteristic quantities for this method or to use higher order reconstruction in place of
"plm"Reconstruction is always performed on the cell-centered magnetic fields. After reconstructing values along a given axis, the values of the reconstructed magnetic field component for that axis are replaced by the face-centered magnetic field values.
riemann solvers
This subsection details the available Riemann Solvers. Currently all available Riemann Solvers are defined to use magnetic fields, however, they all appropriately handle the cases where the magnetic fields are unformly 0. We provide a list of the names used to specify each Riemann Solver in the input file, and a brief description for each of them:
"hll"The HLL approximate Riemann solver with wavespeeds bounds estimated as \(S_L = \min(u_L - a_L, u_R - a_R)\) and \(S_R = \max(u_L + a_L, u_R + a_R)\). This is one of the proposed methods from Davis, 1988, SIAM J. Sci. and Stat. Comput., 9(3), 445–473. The same wavespeed estimator was used in MHD HLL solver implemented for Enzo’s Runge Kutta solver. Currently, this has only been implemented for MHD mode and it will raise an error as it isn’t tested.
"hlle"The HLLE approximate Riemann solver - the HLL solver with wavespeed bounds estimated according to Einfeldt, 1988, SJNA, 25(2), 294–318. This method allows the min/max eigenvalues of Roe’s matrix to be wavespeed estimates. For a description of the procedure for MHD quantities, see Stone et al. (2008) . If using an HLL Riemann Solver, this is the recommended choice. Currently, this has only been implemented for MHD mode.
"hllc"The HLLC approximate Riemann solver. For an overview see Toro, 2009, Riemann Solvers and Numerical Methods for Fluid Dynamics, Springer-Verlag. This is a solver for hydrodynamical problems that models contact and shear waves. The wavespeed bounds are estimated according to the Einfeldt approach. This can only be used in hydro mode.
"hlld"The HLLD approximate Riemann solver described in Miyoshi & Kusano, 2005. JCP, 315, 344. The wavespeed bounds are estimated according to eqn 67 from the paper. This reduces to an HLLC Riemann Solver when magnetic fields are zero (the wavespeed bounds will differ from"hllc). This can only be used in MHD mode.
Note
When the dual-energy formalism is in use, all of the solvers treat the internal energy as a passively advected scalar.
This is not completely self-consistent with the assumptions made by the HLLD solver. Unlike the other HLL-solvers which assume constant pressure in the intermediate regions of the Riemann Fan the HLLD solver assumes constant total pressure. It is unclear whether this causes any problems.
"pm_deposit": particle-mesh
Particle-mesh (“PM”) method component to deposit of field and particle mass into a “total density” field
parameters
Parameter |
Type |
Default |
Description |
|---|---|---|---|
|
float |
0.5 |
Deposit mass at time t + alpha * dt |
fields
particles
For a given particle type to be deposited to the total density field,
it must be part of the "is_gravitating" group, and must have either
an attribute called "mass", or a constant called "mass", but
not both.
If "mass" is an attribute, we loop through the mass attribute array
to get the mass of each particle; and if "mass" is a constant with a
value specified in the input parameter file, the mass of each particle is
equal to this value. In either case, the value of the divided by the cell
volume to get a density quantity, which is deposited on to the grid via
a CIC interpolation scheme.
"pm_update": particle-mesh
Particle-mesh (“PM”) method component to update particle positions
given acceleration fields. Only particle types in the "is_gravitating"
group are updated.
"heat": heat equation
A sample Method for implementing forward-euler to solve the heat equation.
"grackle": chemistry/cooling
Calls methods provided by the external Grackle 3.0 chemistry and cooling library.
Compatability with hydro/mhd solvers
The "grackle" method is compatible with both the "ppm" and the
"mhd_vlct" methods. The convention is to list the hydro method
before "grackle" in the Field:list parameter. This
configuration performs advection and radiative cooling in an
operator-split manner (Note: there isn’t currently support for
performing radiative cooling during the predictor step of the
VL+CT solver).
Integration with hydro-solvers is self-consistent when
Method:Grackle:primordial_chemistry has values of 0 or 1.
However, the integration is somewhat inconsistent when the parameter
exceeds 1. While users shouldn’t be too concerned about this
latter scenario unless they are simulating conditions where
\({\rm H}_2\) makes up a significant fraction of the gas density,
we describe the inconsistencies in greater detail below.
When Method:Grackle:primordial_chemistry > 1, the Grackle library
explicitly models chemistry involving \({\rm H}_2\) and how it
modifies the adiabtic index. Grackle’s routines treat
\(\gamma_0\), the “nominal adiabatic index” specified by
Physics:fluid_props:eos:gamma, as the adiabatic index for all
monatomic species (this should be 5.0/3.0). To that end, Grackle
supplies functions that can effectively be represented as
\(\gamma(e, n_{{\rm H}_2}, n_{\rm other})\) and \(p(\rho, e,
n_{{\rm H}_2}, n_{\rm other})\). In these formulas:
\(p\), \(\rho\) and \(e\) correspond to the quantities held by the
pressure,densityandinternal_energyfields. (Note: the \(\gamma\) function’s dependence on \(e\) accounts for the dependence of \(\gamma_{{\rm H}_2}\) on temperature)\(n_{{\rm H}_2}\) specifies the number density of \({\rm H}_2\). \(n_{\rm other}\) specifies a selection of the other primordial species (that roughly approximate the total number density). In practice, these are computed from passively advected species fields.
There are a handful of locations within the "ppm" and
"mhd_vlct" methods where this treatment is relevant:
Computing the timestep: each hydro/mhd method uses the \(p(\rho, e, n_{{\rm H}_2}, n_{\rm other})\) function for the pressure values. However, they both use \(\gamma_0\) in other places (such as the occurence of adiabatic index in the sound speed formula).
Pre-reconstruction pressure calculation: each hydro/mhd solver internally computes the pressure that is to be reconstructed with \(p=(\gamma_0 - 1)e\rho\).
Riemann Solver: in each hydro/mhd solver, the Riemann Solver completely ignore the grackle supplied functions.
VL+CT Energy floor and DE synchronization: the internal energy floor is computed from the pressure floor using: \(e_{\rm floor} = \frac{p_{\rm floor}}{(\gamma_0 - 1)\rho}\) (thus, \(p_{\rm floor}\) may exceed \(p(\rho, e_{\rm floor}, \ldots)\)). Additionally, synchronizing the internal energy with total energy relies on \(\gamma_0\).
PPM reconstruction: uses \(\gamma_0\).
"comoving_expansion": comoving expansion
Adds the comoving expansion terms to the physical variables.
"turbulence": driving
Turbulence driving.
"gravity": particle-mesh
Particle-mesh (“PM”) method component to compute gravitational potential given a total density field, and calculate associated acceleration fields.
"trace": tracer particles
Moves tracer particles given the velocity field.
"merge_sinks": merge sinks
Merges together sink particles which are separated by less than a given “merging radius”. This is done by copying all sink particles to / from all neighbouring blocks. A Friend-of-Friends algorithm is used to partition particles into groups, where all particles within a given group are separated by less than a merging radius. If a group has more than one particle, one of the particles has its properties changed: its position becomes that of the centre-of-mass of the group, and it takes the total mass, momentum and mass fraction of the whole group. In addition, its ‘lifetime’ attribute is set to be the maximum of the group, its ‘creation_time’ attribute is set to be the minimum of the group, and its ‘id’ attribute is set to the minimum of the group. Other particles in the group are marked for deletion. The final step is for each block to delete all the remaining sink particles which are ‘out-of-bounds’ of the block.
This method requires sink particles to have the following attributes: "mass", "x",
"y", "z", "vx", "vy", "vz", "is_copy", "id", "lifetime",
and "creation_time". All these attributes must be of type "default", except for
"is_copy" and "id" which must be of type "int64". Furthermore, "is_copy"
must be initialized to 0 for all particles.
This method also requires that the number of root blocks across all axes is greater than
2, i.e., that "Mesh:root_blocks" = [a,b,c], where a, b, and c are all
greater than 2.
This procedure cannot handle the case where particles originally from non-neighbouring blocks are put into the same FoF group. If this is found to occur, the program stops and prints an error message. This situation is unlikely to happen, unless the merging radius is too large relative to the block size.
Currently this will only run in unigrid mode. This is because this method
will only work correctly if all blocks containing sink particles are of the
same size, or equivalently, on the same refinement level.
For this reason, there is a check in the constructor of EnzoMethodMergeSinks
for whether "Adapt: max_level" is equal to zero. In future, we plan to
implement a refinement condition that
any block containing a sink particle needs to be on the highest level of
refinement. In this case, the assumption that
blocks containing sink particles are all on the same level of refinement
would be valid.
WARNING: there is currently a memory leak issue when running with this method
which can cause Enzo-E to crash in mysterious ways. If this problem is
encountered, it is advised to increase the batch size parameter
("Particle:batch_size") by a factor of a few
before attempting to run again. To be completely safe, the user can set a
batch size larger than the total number of sink particles in the whole
simulation, which should be feasible for small test problems.
parameters
Parameter |
Type |
Default |
Description |
|---|---|---|---|
|
float |
8.0 |
The merging radius in units of the minimum cell width (i.e., the minimum across all 3 dimensions), at the highest refinement level. |
"accretion": accretion
For cells within a spherical accretion zone around a sink particle, mass is removed
(i.e., the values of the density field are reduced) and added to the sink particle.
The momentum change of gas is in the accretion zone due to the mass loss is accounted
for by changing the momentum of the sink particle, so that total momentum is
conserved. The amount of mass removed is determined by which “flavor” of accretion is
chosen (specified by the "accretion:flavor" parameter), as well as the values
of the “density threshold” (specified by "accretion:physical_density_threshold_cgs") and the
“maximum mass fraction” (specified by "accretion:max_mass_fraction").
In "threshold" flavor accretion, the change in density of each cell is zero if the current
density is below the density threshold. If the current density is above the density threshold,
the change in density is the current density minus the density threshold, or the maximum mass
fraction times the current density, whichever is smaller.
In "bondi_hoyle" flavor accretion, the density change in each cell is calculated according
to the method described in Mark R. Krumholz et al 2004, ApJ, 611, 399. Furthermore, the
density change is limited in the same way as in "threshold" accretion.
In "flux" flavor accretion, the density change in each cell is calculated according to the
method described in Andreas Bleuler & Romain Teyssier 2004, MNRAS, 445, 4015-4036.
Furthermore, the density change is limited in the same way as in "threshold" accretion.
In "dummy" flavor accretion, no accretion is done (essentially, the accretion rate is zero).
This can be useful for testing purposes.
This method can only be used if "merge_sinks" is also used, with "merge_sinks" preceding
"accretion". In addition, this method requires the use of three spatial dimensions.
This method requires the following fields (in addition to the fields required by the hydro
method): "density_source", "density_source_accumulate", "mom_dens_x_source",
"mom_dens_x_source_accumulate", "mom_dens_y_source", "mom_dens_y_source_accumulate",
"mom_dens_z_source", and mom_dens_z_source_accumulate". In addition, if sink particles
have a "metal_fraction" attribute, there must be a "metal_density" field.
This method also requires sink particles to have the following attributes: "mass", "x",
"y", "z", "vx", "vy", "vz", and "accretion_rate", which must all be
of type "default".
parameters
Parameter |
Type |
Default |
Description |
|---|---|---|---|
|
float |
4.0 |
The accretion radius (i.e., the radius of the spherical accretion zone) in units of the minimum cell width (i.e., if the cell width along all the x, y, and z-axes are hx, hy, and hz, then the minimum cell width is the minimum of hx, hy, and hz), at the highest refinement level. Its value must be less than one fewer than the minimum ghost depth for “flux” accretion, and less than the minimum ghost depth for other flavors of accretion. The ghost depth is 4 (along all axes) by default. |
|
string |
|
The flavor of accretion used, which can be either “threshold”, “bondi_hoyle”, or “flux”. If this parameter is not set in the parameter file, or if some other string is provided, then the accretion method will be called but will do nothing. |
|
float |
1.0e-24 |
The value of the physical density threshold in cgs units. The density in each cell in the accretion zone cannot go below this value during the accretion process. The value of the density threshold in code units must be greater than or equal to the value of the density floor imposed by the hydro method. |
|
float |
0.25 |
This parameter specifies the maximum fraction of mass which can be accreted from a cell in one timestep. This value of this parameter must be between 0 and 1. |
"sink_maker": sink maker
This method runs on blocks at the highest level of refinement, and forms sink particles in cells which satisy certain criteria.
First, the gas density in the cell must be larger than the
density threshold, which is specified by the "sink_maker:physical_density_threshold_cgs"
parameter. If so, the mass of the potential sink particle is
\(V_{cell} \times \max(\rho - \rho_{thresh}, f_{max} \, \rho)\), where \(V_{cell}\) is the
cell volume, \(\rho\) is the cell gas density, \(\rho_{thresh}\) is the density
threshold, and \(f_{max}\) is the maximum fraction of the cell mass which can be turned
into a sink particle in one timestep, which is specified by the
"sink_maker:maximum_mass_fraction" parameter. This mass must be greater
than the minimum sink mass, which is specified by "sink_maker:min_sink_mass_solar" (in solar
mass units).
Next, the local Jeans length \(\lambda_J\) is calculated, where
\(\lambda_J = \frac{\pi c_s^2}{G \rho}\), where \(c_s\) is the sound speed of the gas
in the cell, and \(G\) is the gravitational constant. It is then checked whether
\(\lambda_J < N_J \times h_{max}\), where \(N_J\) is specified by
"sink_maker:jeans_length_resolution_cells", and \(h_{max} = max(h_x, h_y, h_z)\), where
\(h_x\), \(h_y\), and \(h_z\) are the cell widths along the x-, y- and z-axes,
respectively.
The next check is that the flow is converging. This is done by computing the strain tensor, given by \(A_{ij} = \frac{1}{2} \, \left( \frac{dv_i}{dx_j} + \frac{dv_j}{dx_i} \right)\). Since this tensor / matrix is real and symmetric, it has three real eigenvalues, and the check is equivalent to checking that all three eigenvalues are negative.
The final check is optional, i.e., it is only done if "sink_maker:check_density_maximum"
is “true”, and a cell will pass this check if it is a local density maximum, that is, its
density is larger than the density in all 26 neighboring cells.
If a cell passes all the checks that are performed, a sink particle is created. Its position
is the coordinates of the center of the cell, plus a small random offset. The maximum size
of the random offset is controlled by "sink_maker:max_offset_cell_fraction".
This method requires sink particles to have the following attributes: "mass", "x",
"y", "z", "vx", "vy", "vz", and "creation_time", which must all be
of type "default"; and "id" and "is_copy", which must be of type "int64".
If sink particles have a "metal_fraction" attribute, there must be a
"metal_density" field.
parameters
Parameter |
Type |
Default |
Description |
|---|---|---|---|
|
float |
4.0 |
If the local Jeans length in a cell is less than this quantity multiplied by the maximum cell width, then the cell is a candidate for forming a sink. The maximum cell width is maximum value out of hx, hy, and hz, where hx, hy, and hz are the cell widths across the x-, y- and z-axes, respectively. |
|
float |
1.0e-24 |
The value of the physical density threshold in cgs units. The density in a cell must be greater than the density threshold to be able to form a sink. The density in a cell after sink formation will be no less than the density threshold. The value of the density threshold in code units must be greater than or equal to the value of the density floor imposed by the hydro method. |
|
float |
0.25 |
The mass of a newly-formed sink is bounded above by this parameter multiplied by the cell density multiplied by the cell volume. The value of this parameter must be between 0 and 1. |
|
float |
0.0 |
The minimum mass of a newly-formed sink particle, in solar mass units. |
|
logical |
true |
Determines whether a cell is required to be a local density maximum in order to form a sink particle. |
|
float |
0.0 |
When a cell creates a sink particle, the x/y/z coordinate of its initial position will be the x/y/z coordinate of the center of the cell, plus a random value generated from a uniform distribution on the interval [-A,A], where A is equal to this parameter multiplied by the cell width along the x/y/z axis. |
|
integer |
0 |
When computing the random offset for the initial position of a sink particle, we compute an unsigned 64 bit integer value from the cycle number, the block index, and the cell index, and then add on this value to give the seed for the random number generator. |