Installation
The compilation process for grackle is very similar to that for
Enzo. For more details on the Enzo build
system, see the Enzo build documentation.
Dependencies
In addition to C/C++ and Fortran compilers, the following dependency must
also be installed:
- HDF5, the hierarchical data format.
HDF5 also may require the szip and zlib libraries, which can be
found at the HDF5 website. Compiling with HDF5 1.8 or greater
requires that the compiler directive
H5_USE_16_API be specified.
This can be done with -DH5_USE_16_API, which is in the machine
specific make files.
Downloading
Grackle is available in a git repository
here. Excellent guides
to git and GitHub are available at
guides.github.com. To clone the Grackle
repo, do the following:
~ $ git clone https://github.com/grackle-project/grackle
Building
- Initialize the build system.
~ $ cd grackle
~/grackle $ ./configure
- Proceed to the source directory.
- Configure the build system.
Note
As of version 2.1, Grackle uses libtool for building and installation.
As such, both shared and static libraries will be built automatically and
it is not necessary to add the -fPIC compiler flag.
Compile settings for different systems are stored in files starting with
“Make.mach” in the source directory. Grackle comes with three sample make
macros: Make.mach.darwin for Mac OSX, Make.mach.linux-gnu for
Linux systems, and an unformatted Make.mach.unknown. If you have a make
file prepared for an Enzo install, it cannot be used straight away, but is a
very good place to start.
Once you have chosen the make file to be used, a few variables should be set:
MACH_LIBTOOL - path to libtool executable. Note, on a Mac,
this should point to glibtool, which can be installed with macports
or homebrew.LOCAL_HDF5_INSTALL - path to your hdf5 installation.LOCAL_FC_INSTALL - path to Fortran compilers (not including the bin
subdirectory).MACH_INSTALL_PREFIX - path where grackle header and library files
will be installed.MACH_INSTALL_LIB_DIR - path where libgrackle will be installed (only
set if different from MACH_INSTALL_PREFIX/lib).MACH_INSTALL_INCLUDE_DIR - path where grackle header files will be
installed (only set if different from MACH_INSTALL_PREFIX/include).
Once the proper variables are set, they are loaded into the build system by
doing the following:
~/grackle/src/clib $ make machine-<system>
Where system refers to the make file you have chosen. For example, if you
chose Make.mach.darwin, type:
~/grackle/src/clib $ make machine-darwin
Custom make files can be saved and loaded from a .grackle directory in the
home directory.
Compiler Settings
There are three compile options available for setting the precision of
baryon fields, compiler optimization, and enabling OpenMP. To see them,
type:
~/grackle/src/clib $ make show-config
MACHINE: Darwin (OSX)
MACHINE-NAME: darwin
CONFIG_PRECISION [precision-{32,64}] : 64
CONFIG_OPT [opt-{warn,debug,high,aggressive}] : high
CONFIG_OMP [omp-{on,off}] : off
For example, to change the optimization to high, type:
~/grackle/src/clib $ make opt-high
Custom settings can be saved for later use by typing:
~/grackle/src/clib $ make save-config-<keyword>
They will be saved in the .grackle directory in your home directory. To
reload them, type:
~/grackle/src/clib $ make load-config-<keyword>
For a list of all available make commands, type:
~/grackle/src/clib $ make help
========================================================================
Grackle Makefile Help
========================================================================
make Compile and generate librackle
make install Copy the library somewhere
make help Display this help information
make clean Remove object files, executable, etc.
make dep Create make dependencies in DEPEND file
make show-version Display revision control system branch and revision
make show-diff Display local file modifications
make help-config Display detailed help on configuration make targets
make show-config Display the configuration settings
make show-flags Display specific compilation flags
make default Reset the configuration to the default values
- Compile and Install
To build the code, type:
~/grackle/src/clib $ make
Updating DEPEND
Compiling calc_rates.F
Compiling cool1d_multi.F
....
Linking
Success!
Then, to install:
~/grackle/src/clib $ make install
- Test your Installation
Once installed, you can test your installation with the provided example to
assure it is functioning correctly. If something goes wrong in this process,
check the out.compile file to see what went wrong during compilation,
or use ldd (otool -L on Mac) on your executable to determine what went
wrong during linking.
~/grackle/src/clib $ cd ../example
~/grackle/src/example $ make clean
~/grackle/src/example $ make
Compiling cxx_example.C
Linking
Success!
~/grackle/src/example $ ./cxx_example
The Grackle Version 2.2
Mercurial Branch default
Mercurial Revision b4650914153d
Initializing grackle data.
with_radiative_cooling: 1.
primordial_chemistry: 3.
metal_cooling: 1.
UVbackground: 1.
Initializing Cloudy cooling: Metals.
cloudy_table_file: ../../input/CloudyData_UVB=HM2012.h5.
Cloudy cooling grid rank: 3.
Cloudy cooling grid dimensions: 29 26 161.
Parameter1: -10 to 4 (29 steps).
Parameter2: 0 to 14.849 (26 steps).
Temperature: 1 to 9 (161 steps).
Reading Cloudy Cooling dataset.
Reading Cloudy Heating dataset.
Initializing UV background.
Reading UV background data from ../../input/CloudyData_UVB=HM2012.h5.
UV background information:
Haardt & Madau (2012, ApJ, 746, 125) [Galaxies & Quasars]
z_min = 0.000
z_max = 15.130
Setting UVbackground_redshift_on to 15.130000.
Setting UVbackground_redshift_off to 0.000000.
Cooling time = -1.434987e+13 s.
Temperature = 4.637034e+02 K.
Pressure = 3.345738e+34.
gamma = 1.666645e+00.
In order to verify that Grackle is fully functional, try running the
test suite.
Running the Test Suite
Grackle contains a number of unit and answer tests to verify that everything is
working properly. These will verify that:
- proper and comoving unit systems are consistent
- atomic, primordial collisional ionization equilibrium agrees with
the analytical solution
- all code examples build and run
- all python examples run and give correct results
- all Python code conforms to PEP 8
Once you have installed pygrackle, the tests can be run from the
src directory by typing make test:
~ $ cd grackle/src
~/grackle/src $ make test
or from the src/python directory by typing py.test:
~ $ cd grackle/src/python
~/grackle/src $ py.test
===================================== test session starts ======================================
platform darwin -- Python 2.7.11, pytest-2.8.1, py-1.4.30, pluggy-0.3.1
rootdir: /Users/britton/Documents/work/simulation/grackle/grackle/src/python, inifile:
collected 13 items
tests/test_chemistry.py ...
tests/test_code_examples.py ....
tests/test_examples.py ........
tests/test_flake8.py .
tests/test_primordial.py .
================================== 17 passed in 68.83 seconds ==================================
Now it’s time to integrate grackle into your simulation code.
Adding Grackle to Your Simulation Code
The majority of this document follows the implementation of Grackle in
a C++ simulation code. Full implementation examples for
C, C++, and Fortran are also available in the Grackle source. See
Example Executables for more information. For a list of all available
functions, see the API Reference.
Example Executables
The grackle source code contains examples for C, C++, and Fortran codes.
They are located in the src/example directory and provide examples
of calling all of grackle’s functions.
- c_example.c - C example
- cxx_example.C - C++ example
- cxx_omp_example.C - C++ example using OpenMP
- fortran_example.F - Fortran example
Once you have already installed the grackle library, you can build the examples
by typing make and the name of the file without extension. For example, to
build the C++ example, type:
To run the example, make sure to add the path to the directory containing
the installed libgrackle.so to your LD_LIBRARY_PATH (or
DYLD_LIBRARY_PATH on Mac).
Data Types
The grackle library provides a configurable variable type to control the
precision of the baryon fields passed to the grackle functions. For C and
C++ codes, this is gr_float. For Fortran codes, this is
R_PREC. The precision of these types can be configured with the
precision compile option. Compile with precision-32 to make
gr_float and R_PREC a 4 byte float (float for C/C++
and real*4 for Fortran). Compile with precision-64 to make
gr_float and R_PREC an 8 byte float (double for C/C++
and real*8 for Fortran).
-
gr_float
Floating point type used for the baryon fields. This is of type float
if compiled with precision-32 and type double if compiled with
precision-64.
-
R_PREC
The Fortran analog of gr_float. This is of type real*4 if
compiled with precision-32 and type real*8 if compiled with
precision-64.
Enabling Output
By default, grackle will not print anything but error messages. However,
a short summary of the running configuration can be printed by setting
grackle_verbose to 1. In a parallel code, it is recommended that
output only be enabled for the root process.
// Enable output
grackle_verbose = 1;
Code Units
It is strongly recommended to use comoving coordinates with any
cosmological simulation. The code_units structure contains
conversions from code units to CGS. If comoving_coordinates is set to
0, it is assumed that the fields passed into the solver are in the
proper frame. All of the units (density, length, time, velocity, and
expansion factor) must be set. When using the proper frame, a_units
(units for the expansion factor) must be set to 1.0.
-
code_units
This structure contains the following members.
-
int
comoving_coordinates
If set to 1, the incoming field data is assumed to be in the comoving
frame. If set to 0, the incoming field data is assumed to be in the
proper frame.
-
double
density_units
Conversion factor to be multiplied by density fields to return
densities in proper g/cm3.
-
double
length_units
Conversion factor to be multiplied by length variables to return
lengths in proper cm.
-
double
time_units
Conversion factor to be multiplied by time variables to return
times in s.
-
double
velocity_units
Conversion factor to be multiplied by velocities to return proper cm/s.
-
double
a_units
Conversion factor to be multiplied by the expansion factor such that
atrue = acode* a_units.
-
double
a_value
The current value of the expansion factor in units of a_units.
The conversion from redshift to expansion factor in code units is given
by a_value = 1 / (1 + z) / a_units. If the
simulation is not cosmological, a_value should be set to 1.
Note, if a_value is set to something other than 1 in a
non-cosmological simulation, all redshift dependent chemistry and
cooling terms will be set corresponding to the redshift given.
code_units my_units;
my_units.comoving_coordinates = 0; // 1 if cosmological sim, 0 if not
my_units.density_units = 1.67e-24; // 1 m_H/cc
my_units.length_units = 3.086e21; // 1 kpc
my_units.time_units = 3.15569e13; // 1 Myr
my_units.velocity_units = my_units.length_units / my_units.time_units;
my_units.a_units = 1.0; // units for the expansion factor
my_units.a_value = 1. / (1. + current_redshift) / my_units.a_units;
If comoving_coordinates is set to 1, it is assumed that the fields being
passed to the solver are in the comoving frame. Hence, the units must
convert from code units in the comoving frame to CGS in the proper
frame.
Note
With comoving_coordinate set to 1, velocity units need to be
defined in the following way.
my_units.velocity_units = my_units.a_units *
(my_units.length_units / a_value) / my_units.time_units; // since u = a * dx/dt
For an example of using comoving units, see the units system in the
Enzo code. For cosmological simulations, a
comoving unit system is preferred, though not required, since it allows the
densities to stay close to 1.0.
Chemistry Data
The main Grackle header file contains a structure of type
chemistry_data called grackle_data, which contains all of the
parameters that control the behavior of the solver. The routine,
set_default_chemistry_parameters() is responsible for the initial setup
of this structure and for setting of all the default parameter values. This
function must be handed a pointer to an instance of chemistry_data,
which will then be attached to grackle_data. The function will return an
integer indicating success (1) or failure (0). After this, parameters can then
be set to their desired values by accessing grackle_data. See
Parameters and Data Files for a full list of the available parameters.
-
chemistry_data
This structure holds all grackle run-time parameters, which are listed in
Parameters and Data Files.
-
chemistry_data_storage
This structure holds all chemistry and cooling rate arrays. All functions
described here make use of an internally stored instance of this type.
The user will not normally encounter this data type, except when using the
Internal Functions.
chemistry_data *my_grackle_data;
my_grackle_data = new chemistry_data;
if (set_default_chemistry_parameters(my_grackle_data) == 0) {
fprintf(stderr, "Error in set_default_chemistry_parameters.\n");
}
// Set parameter values for chemistry.
// Now access the global copy of the chemistry_data struct (grackle_data).
grackle_data->use_grackle = 1; // chemistry on
grackle_data->with_radiative_cooling = 1; // cooling on
grackle_data->primordial_chemistry = 3; // molecular network with H, He, D
grackle_data->metal_cooling = 1; // metal cooling on
grackle_data->UVbackground = 1; // UV background on
grackle_data->grackle_data_file = "CloudyData_UVB=HM2012.h5"; // data file
Once the desired parameters have been set, the chemistry and cooling rates
must be initialized by calling initialize_chemistry_data() with a
pointer to the code_units struct created earlier. This function
will return an integer indicating success (1) or failure (0).
// Finally, initialize the chemistry object.
if (initialize_chemistry_data(&my_units) == 0) {
fprintf(stderr, "Error in initialize_chemistry_data.\n");
return 0;
}
The Grackle is now ready to be used.
Running with OpenMP
As of version 2.2, Grackle can be run with OpenMP parallelism. To do this,
the library must first be compiled with OpenMP support enabled by issuing the
command, “make omp-on”, before compiling. See Compiler Settings for
more information on how to change settings.
For an example of how to compile your code with OpenMP, see the
cxx_omp_example.C code example (Example Executables). Once your code has
been compiled with OpenMP enabled, the number of threads used can be controlled
by setting the omp_nthreads parameter, stored in the grackle_data
struct.
// 8 threads per process
grackle_data->omp_nthreads = 8;
If not set, this parameter will be set to the maximum number of threads
possible, as determined by the system or as configured by setting the
OMP_NUM_THREADS environment variable.
Creating the Necessary Fields
As of version 3.0, the various density and energy fields are passed to
Grackle’s functions using a struct of type grackle_field_data.
The struct contains information about the size and shape of the field arrays
and pointers to all field arrays.
-
grackle_field_data
This structure is used to pass field data to Grackle’s functions. It
contains the following members:
-
int
grid_rank
The active dimensions (not including ignored boundary zones) of the field
arrays.
-
int*
grid_dimension
This should point to an array of size grid_rank. This stores
the size of the field arrays in each dimension.
-
int*
grid_start
This should point to an array of size grid_rank. This stores
the starting value in each dimension for the field data. This can be
used to ignore boundary cells in grid data.
-
int*
grid_end
This should point to an array of size grid_rank. This stores
the end value in each dimension for the field data. This can be used
to ignore boundary cells in grid data.
-
gr_float*
grid_dx
This is the grid cell width in length_units. This is currently
used only in computing approximate H2 self-shielding when H2 is tracked
(primordial_chemistry >= 2) and H2_self_shielding is
set to 1.
-
gr_float*
density
Pointer to the density field array.
-
gr_float*
HI_density
Pointer to the HI density field array. Used when
primordial_chemistry is set to 1, 2, or 3.
-
gr_float*
HII_density
Pointer to the HII density field array. Used when
primordial_chemistry is set to 1, 2, or 3.
-
gr_float*
HM_density
Pointer to the H- density field array. Used when
primordial_chemistry is set to 2 or 3.
-
gr_float*
HeI_density
Pointer to the HeI density field array. Used when
primordial_chemistry is set to 1, 2, or 3.
-
gr_float*
HeII_density
Pointer to the HeII density field array. Used when
primordial_chemistry is set to 1, 2, or 3.
-
gr_float*
HeIII_density
Pointer to the HeIII density field array. Used when
primordial_chemistry is set to 1, 2, or 3.
-
gr_float*
H2I_density
Pointer to the H2 density field array. Used when
primordial_chemistry is set to 2 or 3.
-
gr_float*
H2II_density
Pointer to the H2+ density field
array. Used when primordial_chemistry is set to
2 or 3.
-
gr_float*
DI_density
Pointer to the DI density field array. Used when
primordial_chemistry is set to 3.
-
gr_float*
DII_density
Pointer to the DII density field array. Used when
primordial_chemistry is set to 3.
-
gr_float*
HDI_density
Pointer to the HD density field array. Used when
primordial_chemistry is set to 3.
-
gr_float*
e_density
Pointer to the electron density field array. Used when
primordial_chemistry is set to 1, 2, or 3. Note,
the electron mass density should be scaled by the ratio of the
proton mass to the electron mass such that the electron density
in the code is the electron number density times the proton
mass.
-
gr_float*
metal_density
Pointer to the metal density field array. Used when
metal_cooling is set to 1.
-
gr_float*
internal_energy
Pointer to the internal energy field array.
-
gr_float*
x_velocity
Pointer to the x-velocity field array. Currently not used.
-
gr_float*
y_velocity
Pointer to the y-velocity field array. Currently not used.
-
gr_float*
z_velocity
Pointer to the z-velocity field array. Currently not used.
-
gr_float*
volumetric_heating_rate
Pointer to values containing volumetric heating rates. Rates
should be in units of erg/s/cm3. Used when
use_volumetric_heating_rate is set to 1.
-
gr_float*
specific_heating_rate
Pointer to values containing specific heating rates. Rates
should be in units of erg/s/g. Used when
use_specific_heating_rate is set to 1.
-
gr_float *
RT_heating_rate
Pointer to the radiation transfer heating rate field. Rates
should be in units of (erg/s/cm3) / nHI, where
nHI is the neutral hydrogen number density. Heating rates
for additional species are currently not yet supported.
Used when use_radiative_transfer is set to 1.
-
gr_float *
RT_HI_ionization_rate
Pointer to the HI photo-ionization rate field used with
radiative transfer. Rates should be in units of
1/time_units. Used when
use_radiative_transfer is set to 1.
-
gr_float *
RT_HeI_ionization_rate
Pointer to the HeI photo-ionization rate field used with
radiative transfer. Rates should be in units of
1/time_units. Used when
use_radiative_transfer is set to 1.
-
gr_float *
RT_HeII_ionization_rate
Pointer to the HeII photo-ionization rate field used with
radiative transfer. Rates should be in units of
1/time_units. Used when
use_radiative_transfer is set to 1.
-
gr_float *
RT_H2_dissociation_rate
Pointer to the H2 photo-dissociation rate field
used with radiative transfer. Rates should be in units of
1/time_units. Used when
use_radiative_transfer is set to 1 and
primordial_chemistry is either 2 or 3.
-
gr_float *
H2_self_shielding_length
Pointer to a field containing lengths to be used for
calculating molecular hydrogen column denisty for
H22self-shielding. Used when
H2_self_shielding is set to 2. Field data
should be in length_units.
It is not necessary to attach a pointer to any field that you do
not intend to use.
// Create struct for storing grackle field data
grackle_field_data my_fields;
// Set grid dimension and size.
// grid_start and grid_end are used to ignore ghost zones.
int field_size = 1;
my_fields.grid_rank = 3;
my_fields.grid_dimension = new int[3];
my_fields.grid_start = new int[3];
my_fields.grid_end = new int[3];
my_fields.grid_dx = 1.0; // only matters if H2 self-shielding is used
for (int i = 0;i < 3;i++) {
my_fields.grid_dimension[i] = 1;
my_fields.grid_start[i] = 0;
my_fields.grid_end[i] = 0;
}
my_fields.grid_dimension[0] = field_size;
my_fields.grid_end[0] = field_size - 1;
// Set field arrays.
my_fields.density = new gr_float[field_size];
my_fields.internal_energy = new gr_float[field_size];
my_fields.x_velocity = new gr_float[field_size];
my_fields.y_velocity = new gr_float[field_size];
my_fields.z_velocity = new gr_float[field_size];
// for primordial_chemistry >= 1
my_fields.HI_density = new gr_float[field_size];
my_fields.HII_density = new gr_float[field_size];
my_fields.HeI_density = new gr_float[field_size];
my_fields.HeII_density = new gr_float[field_size];
my_fields.HeIII_density = new gr_float[field_size];
my_fields.e_density = new gr_float[field_size];
// for primordial_chemistry >= 2
my_fields.HM_density = new gr_float[field_size];
my_fields.H2I_density = new gr_float[field_size];
my_fields.H2II_density = new gr_float[field_size];
// for primordial_chemistry >= 3
my_fields.DI_density = new gr_float[field_size];
my_fields.DII_density = new gr_float[field_size];
my_fields.HDI_density = new gr_float[field_size];
// for metal_cooling = 1
my_fields.metal_density = new gr_float[field_size];
// volumetric heating rate (provide in units [erg s^-1 cm^-3])
my_fields.volumetric_heating_rate = new gr_float[field_size];
// specific heating rate (provide in units [egs s^-1 g^-1]
my_fields.specific_heating_rate = new gr_float[field_size];
// heating rate from radiative transfer calculations (provide in units [erg s^-1 cm^-3]
my_fields.RT_heating_rate = new gr_float[field_size];
// HI ionization rate from radiative transfer calculations (provide in units of [ 1/time_units ]
my_fields.RT_HI_ionization_rate = new gr_float[field_size];
// HeI ionization rate from radiative transfer calculations (provide in units of [1/time_units]
my_fields.RT_HeI_ionization_rate = new gr_float[field_size];
// HeII ionization rate from radiative transfer calculations (provide in units of [1/time_units]
my_fields.RT_HeII_ionization_rate = new gr_float[field_size];
// H2 dissociation rate from radiative transfer calculations (provide in units of [1/time_units]
my_fields.RT_H2_dissociation_rate = new gr_float[field_size];
Note
The electron mass density should be scaled by the ratio of the
proton mass to the electron mass such that the electron density in the
code is the electron number density times the proton mass.
Calling the Available Functions
There are five functions available, one to solve the chemistry and cooling
and four others to calculate the cooling time, temperature, pressure, and the
ratio of the specific heats (gamma). The arguments required are the
code_units structure and the grackle_field_data struct.
For the chemistry solving routine, a timestep must also be given. For the
four field calculator routines, the array to be filled with the field values
must be created and passed as an argument as well.
The examples below make use of Grackle’s Primary Functions, where
the parameters and rate data are stored in instances of the
chemistry_data and chemistry_data_storage structs
declared in grackle.h. Alternatively, a set of Local Functions
require these structs to be provided as arguments, allowing for explicitly
thread-safe code.
Solve the Chemistry and Cooling
// some timestep (one million years)
double dt = 3.15e7 * 1e6 / my_units.time_units;
if (solve_chemistry(&my_units, &my_fields, dt) == 0) {
fprintf(stderr, "Error in solve_chemistry.\n");
return 0;
}
Calculating the Cooling Time
gr_float *cooling_time;
cooling_time = new gr_float[field_size];
if (calculate_cooling_time(&my_units, &my_fields,
cooling_time) == 0) {
fprintf(stderr, "Error in calculate_cooling_time.\n");
return 0;
}
Calculating the Temperature Field
gr_float *temperature;
temperature = new gr_float[field_size];
if (calculate_temperature(&my_units, &my_fields,
temperature) == 0) {
fprintf(stderr, "Error in calculate_temperature.\n");
return EXIT_FAILURE;
}
Calculating the Pressure Field
gr_float *pressure;
pressure = new gr_float[field_size];
if (calculate_pressure(&my_units, &my_fields,
pressure) == 0) {
fprintf(stderr, "Error in calculate_pressure.\n");
return EXIT_FAILURE;
}
Calculating the Gamma Field
gr_float *gamma;
gamma = new gr_float[field_size];
if (calculate_gamma(&my_units, &my_fields,
gamma) == 0) {
fprintf(stderr, "Error in calculate_gamma.\n");
return EXIT_FAILURE;
}
Cleaning the memory
_free_chemistry_data(my_grackle_data, &grackle_rates);
Grackle is using global structures and therefore the global structure grackle_rates needs also to be released.
Parameters and Data Files
Parameters
For all on/off integer flags, 0 is off and 1 is on.
-
int
use_grackle
Flag to activate the grackle machinery. Default: 0.
-
int
with_radiative_cooling
Flag to include radiative cooling and actually update the thermal
energy during the chemistry solver. If off, the chemistry species
will still be updated. The most common reason to set this to off
is to iterate the chemistry network to an equilibrium state.
Default: 1.
-
int
primordial_chemistry
Flag to control which primordial chemistry network is used.
Default: 0.
- 0: no chemistry network. Radiative cooling for primordial
species is solved by interpolating from lookup tables
calculated with Cloudy.
- 1: 6-species atomic H and He. Active species: H, H+,
He, He+, ++, e-.
- 2: 9-species network including atomic species above and species
for molecular hydrogen formation. This network includes
formation from the H- and H2+
channels, three-body formation (H+H+H and H+H+H2),
H2 rotational transitions, chemical heating, and
collision-induced emission (optional). Active species: above +
H-, H2, H2+.
- 3: 12-species network include all above plus HD rotation cooling.
Active species: above + D, D+, HD.
Note
In order to make use of the non-equilibrium chemistry
network (primordial_chemistry options 1-3), you must add
and advect baryon fields for each of the species used by that
particular option.
-
int
h2_on_dust
Flag to enable H2 formation on dust grains, dust cooling, and
dust-gas heat transfer follow Omukai (2000). This assumes
that the dust to gas ratio scales with the metallicity. Default: 0.
-
int
metal_cooling
Flag to enable metal cooling using the Cloudy tables. If enabled, the
cooling table to be used must be specified with the
grackle_data_file parameter. Default: 0.
Note
In order to use the metal cooling, you must add and advect a
metal density field.
-
int
cmb_temperature_floor
Flag to enable an effective CMB temperature floor. This is implemented
by subtracting the value of the cooling rate at TCMB from the
total cooling rate. Default: 1.
-
int
UVbackground
Flag to enable a UV background. If enabled, the cooling table to be
used must be specified with the grackle_data_file parameter.
Default: 0.
-
float
UVbackground_redshift_on
Used in combination with UVbackground_redshift_fullon,
UVbackground_redshift_drop, and
UVbackground_redshift_off to set an attenuation factor for the
photo-heating and photo-ionization rates of the UV background model.
See the figure below for an illustration its behavior. If not set,
this parameter will be set to the highest redshift of the UV background
data being used.
-
float
UVbackground_redshift_fullon
Used in combination with UVbackground_redshift_on,
UVbackground_redshift_drop, and
UVbackground_redshift_off to set an attenuation factor for the
photo-heating and photo-ionization rates of the UV background model.
See the figure below for an illustration its behavior. If not set,
this parameter will be set to the highest redshift of the UV background
data being used.
-
float
UVbackground_redshift_drop
Used in combination with UVbackground_redshift_on,
UVbackground_redshift_fullon, and
UVbackground_redshift_off to set an attenuation factor for the
photo-heating and photo-ionization rates of the UV background model.
See the figure below for an illustration its behavior. If not set,
this parameter will be set to the lowest redshift of the UV background
data being used.
-
float
UVbackground_redshift_off
Used in combination with UVbackground_redshift_on,
UVbackground_redshift_fullon, and
UVbackground_redshift_drop to set an attenuation factor for the
photo-heating and photo-ionization rates of the UV background model.
See the figure below for an illustration its behavior. If not set,
this parameter will be set to the lowest redshift of the UV background
data being used.
-
char*
grackle_data_file
Path to the data file containing the metal cooling and UV background
tables. Default: “”.
-
float
Gamma
The ratio of specific heats for an ideal gas. A direct calculation
for the molecular component is used if primordial_chemistry
> 1. Default: 5/3.
-
int
three_body_rate
Flag to control which three-body H2 formation rate is used.
The first five options are discussed in Turk et. al. (2011). Default: 0.
-
int
cie_cooling
Flag to enable H2 collision-induced emission cooling from
Ripamonti & Abel (2004). Default: 0.
-
int
h2_optical_depth_approximation
Flag to enable H2 cooling attenuation from Ripamonti &
Abel (2004).
Default: 0.
-
int
photoelectric_heating
Flag to enable a spatially uniform heating term approximating
photo-electric heating from dust from Tasker & Bryan (2008). Default: 0.
-
int
photoelectric_heating_rate
If photoelectric_heating is enabled, the heating rate in
units of (erg cm-3 s-1) n-1, where n is
the total hydrogen number density. In other words, this is the
volumetric heating rate at a hydrogen number density of
n = 1 cm-3. Default: 8.5e-26.
-
int
Compton_xray_heating
Flag to enable Compton heating from an X-ray background following
Madau & Efstathiou (1999). Default: 0.
-
float
LWbackground_intensity
Intensity of a constant Lyman-Werner H2 photo-dissociating
radiation field in units of 10-21 erg s-1 cm-2 Hz-1 sr-1. Default: 0.
-
int
LWbackground_sawtooth_suppression
Flag to enable suppression of Lyman-Werner flux due to Lyman-series
absorption (giving a sawtooth pattern), taken from Haiman & Abel,
& Rees (2000).
Default: 0.
-
float
HydrogenFractionByMass
The fraction by mass of Hydrogen in the metal-free portion of the
gas (i.e., just the H and He). In the non-equilibrium solver, this is
used to ensure consistency in the densities of the individual species.
In tabulated mode, this is used to calculate the H number density from
the total gas density, which is a parameter of the heating/cooling tables.
When using the non-equilibrium solver, a sensible default is 0.76.
However, the tables for tabulated mode were created assuming
nHe/nH = 0.1, which corresponds to an H mass fraction of
about 0.716. When running in tabulated mode, this parameter will automatically
be changed to this value. Default: 0.76.
-
float
DeuteriumToHydrogenRatio
The ratio by mass of Deuterium to Hydrogen. Default: 6.8e-5 (the value
from Burles & Tytler (1998)
multiplied by 2 for the mass of Deuterium).
-
float
SolarMetalFractionByMass
The fraction of total gas mass in metals for a solar composition.
Default: 0.01295 (consistent with the default abundances in the Cloudy code).
-
int
use_volumetric_heating_rate
Flag to signal that an array of volumetric heating rates is being
provided in the volumetric_heating_rate field of the
grackle_field_data struct. Default: 0.
-
int
use_specific_heating_rate
Flag to signal that an array of specific heating rates is being
provided in the specific_heating_rate field of the
grackle_field_data struct. Default: 0.
-
int
use_radiative_transfer
Flag to signal that arrays of ionization and heating rates from
radiative transfer solutions are being provided. Only
available if primordial_chemistry is greater than 0. HI, HeI,
and HeII ionization arrays are provided in RT_HI_ionization_rate,
RT_HeI_ionization_rate, and RT_HeII_ionization_rate
fields, respectively, of the grackle_field_data struct.
Associated heating rate is provided in the RT_heating_rate
field, and H2photodissociation rate can also be provided in the
RT_H2_dissociation_rate field when
primordial_chemistry is set to either 2 or 3. Default: 0.
-
int
radiative_transfer_coupled_rate_solver
When used with use_radiative_transfer set to 1, this flag
makes it possible to solve the chemistry and cooling of the
computational elements for which the radiation field is non-zero
separately from those with no incident radiation. This allows radiation
transfer calculations to be performed on a smaller timestep than the
global timestep. The parameter,
radiative_transfer_intermediate_step, is then used to toggle
between updating the cells/particles receiving radiative input and those
that are not. Default: 0.
-
int
radiative_transfer_intermediate_step
Used in conjunction with radiative_transfer_coupled_rate_solver
set to 1, setting this parameter to 1 tells the solver to only update
cells/particles where the radiation field is non-zero. Setting this
to 0 updates only those elements with no incident radiation. When
radiative_transfer_coupled_rate_solver is set to 0, changing
this parameter will have no effect. Default: 0.
-
int
radiative_transfer_hydrogen_only
Flag to only use hydrogen ionization and heating rates from the
radiative transfer solutions. Default: 0.
-
int
H2_self_shielding
Switch to enable approximate H2 self-shielding from both the UV
background dissociation rate and the H2 dissociation rate
given by RT_H2_dissociation_rate (if present). Three options
exist for the length scale used in calculating the H2 column
density. Default: 0.
-
int
self_shielding_method
Switch to enable approximate self-shielding from the UV background.
All three of the below methods incorporate Eq. 13 and 14 from
Rahmati et. al. 2013.
These equations involve using the spectrum averaged photoabsorption cross
for the given species (HI or HeI). These redshift dependent values are
pre-computed for the HM2012 and FG2011 UV backgrounds and included in
their respective cooling data tables. Default: 0
Care is advised in using any of these methods. The default behavior is to
apply no self-shielding, but this is not necessarily the proper assumption,
depending on the use case. If the user desires to turn on self-shielding,
we strongly advise using option 3. All options include HI self-shielding, and
vary only in treatment of HeI and HeII. In options 2 and 3, we approximately
account for HeI self-shielding by applying the Rahmati et. al. 2013 relations,
which are only strictly valid for HI, to HeI under the assumption that it behaves
similarly to HI. None of these options are completely correct in practice,
but option 3 has produced the most reasonable results
in test simulations. Repeating the analysis of Rahmati et. al. 2013 to
directly parameterize HeI and HeII self-shielding behavior would be a valuable
avenue of future research in developing a more complete self-shielding model.
Each self-shielding option is described below.
- 0: No self shielding. Elements are optically thin to the UV background.
- 1: Not Recommended. Approximate self-shielding in HI only.
- HeI and HeII are left as optically thin.
- 2: Approximate self-shielding in both HI and HeI. HeII remains
- optically thin.
- 3: Approximate self-shielding in both HI and HeI, but ignoring
- HeII ionization and heating from the UV background entirely
(HeII ionization and heating rates are set to zero).
These methods only work in conjunction with using updated Cloudy
cooling tables, denoted with “_shielding”. These tables properly account
for the decrease in metal line cooling rates in self-shielded regions,
which can be significant.
For consistency, when primordial_chemistry > 2, the self-shielding
attenutation factors calculated for HI and HeI are applied to the
H2ionization (15.4 eV) and H2+ dissociation
rates (30 eV) respectively. These reaction rates are distinct from the
H2self-shielding computed using the H2_self_shielding
flag.
-
int
omp_nthreads
Sets the number of OpenMP threads. If not set, this will be set to
the maximum number of threads possible, as determined by the system
or as configured by setting the OMP_NUM_THREADS environment
variable. Note, Grackle must be compiled with OpenMP support
enabled. See Running with OpenMP.
Data Files
All data files are located in the input directory in the source.
The first three files contain the heating and cooling rates for both
primordial and metal species as well as the UV background photo-heating
and photo-ionization rates. For all three files, the valid density and
temperature range is given below. Extrapolation is performed when
outside of the data range. The metal cooling rates are stored for
solar metallicity and scaled linearly with the metallicity of the gas.
Valid range:
- number density: -10 < log10 (nH / cm-3) < 4
- temperature: the temperature range is 1 < log10 (T / K) < 9.
Data files:
- CloudyData_noUVB.h5 - cooling rates for collisional ionization
equilibrium.
- CloudyData_UVB=FG2011.h5 - heating and cooling rates and UV
background rates from the work of Faucher-Giguere et. al. (2009), updated in 2011.
The maxmimum redshift is 10.6. Above that, collisional ionization
equilibrium is assumed.
- CloudyData_UVB=HM2012.h5 - heating and cooling rates and UV
background rates from the work of Haardt & Madau (2012). The maximum
redshift is 15.13. Above that, collisional ionization equilibrium is
assumed.
To use the self-shielding approximation (see self_shielding_method),
one must properly account for the change in metal line cooling rates in
self-shielded regions. Using the optically thin tables described above can
result in an order of magnitude overestimation in the net cooling rate at
certain densities. We have re-computed these tables by constructing
Jeans-length depth models in Cloudy at each density - temperature pair,
tabulating the cooling and heating rates from the core of each of these
clouds. These models enforce a maximum depth of 100 pc.
In addition, these tables contain the spectrum averaged absorption
cross sections needed for the Rahmati et. al. 2013 self-shielding
approximations. Currently only the HM2012 table has been recomputed.
- CloudyData_UVB=HM2012_shielded.h5 - updated heating and cooling
rates with the HM2012 UV background, accounting for self-shielding.
- CloudyData_UVB=FG2011_shielded.h5 - updated heating and cooling
rates with the FG2011 UV background, accounting for self-shielding.
The final file includes only metal cooling rates under collisional
ionization equilibrium, i.e., no incident radiation field. This table
extends to higher densities and also varies in metallicity rather than
scaling proportional to the solar value. This captures the
thermalization of metal coolants occuring at high densities, making this
table more appropriate for simulations of collapsing gas-clouds.
Valid range:
- number density: -6 < log10 (nH / cm-3) < 12
- metallicity: -6 < log10 (Z / Zsun) < 1
- temperature: the temperature range is 1 < log10 (T / K) < 8.
Data file:
- cloudy_metals_2008_3D.h5 - collisional ionization equilibrium,
metal cooling rates only.
Pygrackle: Running Grackle in Python
Grackle comes with a Python interface, called Pygrackle, which provides
access to all of Grackle’s functionality. Pygrackle requires the following
Python packages:
The easiest thing to do is follow the instructions for installing yt,
which will provide you with Cython, matplotlib, and NumPy. Flake8 and
py.test can then be installed via pip.
Installing Pygrackle
Once the Grackle library has been built and the above dependencies have been
installed, Pygrackle can be installed by moving into the src/python
directory and running python setup.py install.
~/grackle $ cd src/python
~/grackle/src/python $ python setup.py install
Note
Pygrackle can only be run when Grackle is compiled without OpenMP.
See Running with OpenMP.
Running the Example Scripts
A number of example scripts are available in the src/python/examples
directory. These scripts provide examples of ways that Grackle can be
used in simplified models, such as solving the temperature evolution of
a parcel of gas at constant density or in a free-fall model. Each example
will produce a figure as well as a dataset that can be loaded and analyzed
with yt.
Cooling Cell Example
This sets up a single grid cell with an initial density and temperature and solves
the chemistry and cooling for a given amount of time. The resulting dataset gives
the values of the densities, temperatures, and mean molecular weights for all times.
>>> import yt
>>> ds = yt.load("cooling_cell.h5")
>>> print ds.data["time"].to("Myr")
YTArray([ 0.00000000e+00, 6.74660169e-02, 1.34932034e-01, ...,
9.98497051e+01, 9.99171711e+01, 9.99846371e+01]) Myr
>>> print ds.data["temperature"]
YTArray([ 990014.56406726, 980007.32720091, 969992.99066987, ...,
9263.81515866, 9263.81515824, 9263.81515865]) K
Free-Fall Collapse Example
This sets up a single grid cell with an initial number density of 1 cm-3.
The density increases with time following a free-fall collapse model. As the density
increases, thermal energy is added to model heating via adiabatic compression.
This can be useful for testing chemistry networks over a large range in density.
The resulting dataset can be analyzed similarly as above.
>>> import yt
>>> ds = yt.load("freefall.h5")
>>> print ds.data["time"].to("Myr")
[ 0. 0.45900816 0.91572127 ..., 219.90360841 219.90360855
219.9036087 ] Myr
>>> print ds.data["density"]
[ 1.67373522e-25 1.69059895e-25 1.70763258e-25 ..., 1.65068531e-12
1.66121253e-12 1.67178981e-12] g/cm**3
>>> print ds.data["temperature"]
[ 99.94958248 100.61345564 101.28160228 ..., 1728.89321898
1729.32604568 1729.75744287] K
Simulation Dataset Example
This provides an example of using the grackle library for calculating chemistry and
cooling quantities for a pre-existing simulation dataset. To run this example, you
must also download the IsolatedGalaxy dataset from the yt sample data page.
How to Develop Grackle
Grackle is a community project!
We are very happy to accept patches, features, and bugfixes from any member of
the community! Grackle is developed using Git, primarily because of how well
it enables open-source, community contribution. We’re eager to hear from you.
Note
If you are already familiar with Git and GitHub,
the best way to contribute is to fork the main Grackle repository, make your changes, push them
to your fork, and issue a pull request. The rest of this document is just an
explanation of how to do that.
Keep in touch, and happy hacking!
Open Issues
If you’re interested in participating in Grackle development, take a look at the
issue tracker on GitHub.
If you are encountering a bug that is not already tracked there, please open a
new issue.
Contributing to Grackle with Git and Github
We provide a brief introduction to submitting changes here. We encourage
contributions from any user. If you are new to Git and/or GitHub, there are
excellent guides available at guides.github.com,
specifically the Git Handbook, and the GitHub
Hello World. We are also
happy to provide guidance on the mailing list or in our slack channel.
Licensing
Grackle is under the Enzo public license, a BSD-like license.
All contributed code must be BSD-compatible. If you’d rather not license in
this manner, but still want to contribute, please consider creating an external
package, which we’ll happily link to in the Grackle documentation.
How To Get The Source Code For Editing
Grackle is hosted on GitHub, and you can see all of the Grackle repositories at
https://github.com/grackle-project/. In order to modify the source code for Grackle,
we ask that you make a “fork” of the main Grackle repository on GitHub. A
fork is simply an exact copy of the main repository (along with its history)
that you will now own and can make modifications as you please. You can create
a personal fork by visiting the Grackle GitHub webpage at
https://github.com/grackle-project/grackle/. After logging in, you should see an
option near the top right labeled “fork”. You now have a forked copy of
the Grackle repository for your own personal modification.
This forked copy exists on GitHub under your username, so in order to access
it locally, follow the instructions at the top of that webpage for that
forked repository:
$ git clone http://bitbucket.org/<USER>/<REPOSITORY_NAME>
This downloads that new forked repository to your local machine, so that you can
access it, read it, make modifications, etc. It will put the repository in a
local directory of the same name as the repository in the current working
directory.
Verify that you are on the master branch of Grackle by running:
If you’re not on the master branch, you can get to it with:
You can see any past state of the code by using the git log command.
For example, the following command would show you the last 5 revisions
(modifications to the code) that were submitted to that repository.
Using the revision specifier (the number or hash identifier next to each
changeset), you can update the local repository to any past state of the
code (a previous changeset or version) by executing the command:
$ git checkout revision_specifier
Making and Sharing Changes
The simplest way to submit changes to Grackle is to do the following:
- Fork the main repository.
- Clone your fork.
- Make some changes and commit them.
- Push the changesets to your fork.
- Issue a pull request.
Here’s a more detailed flowchart of how to submit changes.
Fork Grackle on GitHub. (This step only has to be done once.) You can do
this by clicking on the fork button in the top-right corner of the main
repository.
Create a new branch in which to make your changes by doing git
checkout -b <new branch name>. This will make it easy to move back and
forth between the main branch of the code and your changes.
Edit the source file you are interested in and test your changes.
Use git add <files> to stage files to be committed.
Commit your changes with git commit. This will open a text editor so you
can write a commit message. To add your message inline, do
git commit -m "<commit message>". You can list specific file to be
committed.
Remember that this is a large development effort and to keep the code
accessible to everyone, good documentation is a must. Add in source code
comments for what you are doing. Add documentation to the appropriate
section of the online docs so that people other than yourself know how
to use your new code.
If your changes include new functionality or cover an untested area of the
code, add a test. Commit these changes as well.
Push your changes to your new fork using the command:
$ git push origin <branch name>
Note
Note that the above approach uses HTTPS as the transfer protocol
between your machine and GitHub. If you prefer to use SSH - or
perhaps you’re behind a proxy that doesn’t play well with SSL via
HTTPS - you may want to set up an SSH key
on GitHub. Then, you use
the syntax ssh://git@github.com/<USER>/grackle, or equivalent, in
place of https://github.com/<USER>/grackle in git commands.
For consistency, all commands we list in this document will use the HTTPS
protocol.
Issue a pull request by going to the main repository and clicking on the
green button that says Compare & pull request. This will open up a page
that will allow you to enter a description of the changes to be merged. Once
submitted, a series of automated tests will run and their status will be
reported on the pull request page.
During the course of your pull request you may be asked to make changes. These
changes may be related to style issues, correctness issues, or requesting
tests. The process for responding to pull request code review is relatively
straightforward.
Make requested changes, or leave a comment on the pull request page on
GitHub indicating why you don’t think they should be made.
Commit those changes to your local repository.
Push the changes to your fork:
$ git push origin <branch name>
Your pull request will be automatically updated.
Once your pull request has been accepted, you can safely delete your
branch:
$ git branch --delete <branch name>
Updating Your Branch
If your branch or pull request has been open for some time, it may be useful
to keep it up to date with the latest changes from the main repository. This
can be done by rebasing your changes.
Before doing this, you will need to be able to pull the latest changes from
the main repository.
Add the main repository as a remote:
$ git remote add grackle https://github.com/grackle-project/grackle
You can verify that it has been added by doing git remote -v. This
only needs to be done once.
Go back to the master branch and pull the changes:
$ git checkout master
$ git pull grackle master
Return to your branch and rebase your changes onto the head of the master
branch:
$ git checkout <branch name>
$ git rebase master
This should go smoothly unless changes have been made to the same lines in
the source, in which case you will need to fix conflicts. After rebasing,
you will get an error when trying to push your branch to your fork. This is
because you have changed the order of commits and git does not like that.
In this case, you will need to add “-f” to your push command to force
the changes to be accepted.:
$ git push -f origin <branch name>
Have fun!
