EnzoE Methods¶
[ This page is under development ]
This section decribes all Methods available in EnzoE: what they do, what methodspecific 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 EnzoE 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 
Pressurefree 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"
dualenergy formalism.
See dualenergy 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 predictorcorrector scheme, with attributes similar to the MUSCLHancock method. For both modes, the method includes two main steps. First the values are integrated to the half timestep. Then, the values at the half timestep are used to caluclate the change in the quantities over the full timestep.
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 facecentered Cello fields. In more detail:
\(B_x\) is stored on the xfaces
\(B_y\) is stored on the yfaces
\(B_z\) is stored on the zfaces
The method also tracks components of the magnetic fields at the cellcenters, which just store the average of the values from the cell’s faces.
Currently, this method only support 3dimensional problems. In the future, alternative modes supporting MHD could easily be implemented that use divergencecleaning 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] 
Cellcentered bfield (average of the corresponding 


[rw] 
Cellcentered bfield (average of the corresponding 


[rw] 
Cellcentered bfield (average of the corresponding 


[rw] 
Primary representation of xcomponent of bfield (lies on xfaces). 


[rw] 
Primary representation of ycomponent of bfield (lies on yfaces). 


[rw] 
Primary representation of zcomponent of bfield (lies on zfaces). 


[w] 
computed from 


[rw] 
if dualenergy 
In hydromode, none of the 6 fields used to store the magnetic field should be defined.
At initialization the facecentered magnetic field should be
divergence free. Trivial configurations (e.g. a constant magnetic
field everywhere) can be provided with the "value"
initializer. For nontrivial configurations, we have provide the
"vlct_bfield"
initializer which can initialize the magnetic fields
(facecentered and cellcentered) 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"
dualenergy formalism.
See dualenergy formalism for additional details.
reconstruction¶
This subsection details the available interpolation methods for reconstructing the left and right states of the cellcentered interfaces. Presently, all available methods perform reconstruction on cellcentered 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 rightdifferences: \(\delta{\bf w}_{L,i}={\bf w}_{i}{\bf w}_{i1}\) and \(\delta{\bf w}_{R,i}={\bf w}_{i1}{\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 halftimestep and then piecewiselinear reconstruction for the fulltimestep (most test problems have been run usingplm
withtheta_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 cellcentered magnetic fields. After reconstructing values along a given axis, the values of the reconstructed magnetic field component for that axis are replaced by the facecentered 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, SpringerVerlag. 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 dualenergy formalism is in use, all of the solvers treat the internal energy as a passively advected scalar.
This is not completely selfconsistent with the assumptions made by the HLLD solver. Unlike the other HLLsolvers 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"
: particlemesh¶
Particlemesh (“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"
: particlemesh¶
Particlemesh (“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 forwardeuler 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
operatorsplit manner (Note: there isn’t currently support for
performing radiative cooling during the predictor step of the
VL+CT solver).
Integration with hydrosolvers is selfconsistent 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
,density
andinternal_energy
fields. (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).
Prereconstruction 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"
: particlemesh¶
Particlemesh (“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 FriendofFriends 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 centreofmass 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 ‘outofbounds’ 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 nonneighbouring 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 EnzoE 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, 40154036.
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 zaxes 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.0e24 
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 zaxes,
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 zaxes, respectively. 

float 
1.0e24 
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 newlyformed 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 newlyformed 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. 