Chemistry
If you use this chemistry module in your publication, please cite the accompanying code paper Gong et al. (2023). The paper also provides some more details beyond this documentation.
Because the chemical reaction rates for different species often differ by orders of magnitude, implicit solvers are better suited to solve the rate equations. We use the operator split method: at the last stage of the passive scalar flux integration task, before the conserved to primitive conversion, we advance the chemical abundances and the internal energy for the hydrodynamic time step dt, by solving a set of coupled ODEs of chemical rate equations and heating/cooling implicitly:
In the code, the reaction rate RHS() (right-hand-side) function and the Edot. Because the chemical reaction rate T, the user often need to have some expression of T as a function of the internal energy e in their RHS. Currently, chemistry works with a simple isothermal or adiabatic equation of state with a fixed adiabatic index of 5/3. General EOS dependence on chemistry has not been included yet.
Note that the passive scalar advection (and diffusion) is solved separately prior to the chemical rate equations. This operator split method is only first-order accurate. The conservation of elemental abundances (such as by rescaling the flux of different species) has not been implemented yet.
./configure.py --prob=[problem] --chemistry=[network] --ode_solver=[ode_solver]
-
[network]is the name of the chemical network insrc/chemistry/network/directory. -
[ode_solver]is the choice of the ODE solver.
Sample input files are in the input/chemistry/ directory. Several input parameters related to the CVODE solver in the <chemistry> block are with noting:
-
h_init: the initial timestep for the CVODE solver. Default uses the estimation from CVODE library, which usually works well. -
use_previous_h: choose whether to use the last CVODE solver internal subcycling timestep in the previous MHD timestep for the initial guess of CVODE solver in the next MHD timestep. Default is set tofalse, which works better for out-of-equilibrium chemistry. -
reltol, abstol: relative and absolute tolerance. More relaxed tolerance runs faster in general, andreltoloften has a larger impact on the solver speed thanabstol. Also note thatabstolshould not be too large, since the solver might be unstable and the abundance of trace species may not be correct otherwise. -
abstol_E: You can set the relative and absolute tolerance for each chemical species and internal energy separately. For example,abstol_Esets the absolute tolerance for internal energy, which correspond to about 0.01 K in thechem_six_rayregression test. -
maxorder: Maxiumum order of the backward differentiation formulas used by CVODE. The default value is 2, which gives a reasonable balance between speed and stability.
In addition, the chemistry module can be used for post-processing by setting the parameter active=fixed in the <hydro> block. In this case, the chemistry and radiation is solved iteratively at a timestep of dt=tlim/nlim, where tlim and nlim can be set in the <time> input block.
There is a set of regression tests available for chemistry. Exploring the regression tests is also a good way to get familiar with the chemistry module. You can run the regression tests in the tst/regression directory:
python run_test.py chemistry
for all chemistry regression tests, or
python run_test.py chemistry/[test_name]
for one specific test.
Most tests (except for chem_H2 and read_vtk) require installing the CVODE library. Below, the existing tests are described briefly.
The simplest test chem_H2 uses a 2-species HI-H2 reaction network without hydrodynamics. It compares the numerical results of chemical evolution with the existing analytic solution. The simple forward Euler ODE solver is used.
The test chem_H2_gaussian combines chemistry with hydrodynamics. An initial Gaussian profile of HI is advected and converted to H2 through the HI-H2 reaction network. The test is run at different numerical resolutions, checking the numerical error and convergence.
Several tests run chemistry with radiation. They all use the GOW17 network appropriate for the ISM chemistry and heating/cooling. These tests do not have an analytic solution and use a numerical result as a reference solution. The chem_gow17 and chem_kida_gow17 tests use a constant radiation field. The later uses the general KIDA network format to implement the GOW17 chemistry. The chem_six_ray and chem_pdr_static uses the six-ray radiation in 3D and 1D setup, respectively. The read_vtk actually does not include chemistry, and only is used as a demonstration for reading vtk files, which is useful for chemistry post-processing of numerical simulation outputs (see also the chem_six_ray test).
The chemical species uses the passive scalar structure in the code. By default, the number of species is the number of scalars. If the user wants additional passive scalars to be added (for example, tracer scalars that is not involved in chemical reactions), the number of scalars can be set separately in the --nscalars=[number of passive scalars] configure option, and should be larger than the number of chemical species. Currently, the program does not support the option of setting --nscalars=0 and only use the chemistry solver for heating and cooling.
This is a 2-species HI-H2 network. Although this network is very simplified and does not represent the realistic HI-H2 transition in the ISM, it has the advantage of having an analytic solution, and is used in regression tests.
The GOW17 network is from the work of Gong, Ostriker and Wolfire (2017). It is a widely used and well-tested network for carbon and oxygen chemistry in the atomic and molecular ISM. It has the advantage of being relatively simple with only 18 species and about 50 reactions, and is still accurate in capturing the most important chemical and thermal processes compared to more sophisticated networks.
This primordial chemistry network is originally implemented in Grassi et al. (2014), and used in the modified shock tube test.
The KIDA network allows flexible chemical networks to be constructed by the user simply by providing some text files in the KIDA format Wakelam et al. (2012). We provide an example of the GOW17 network written in the KIDA format. This general KIDA network has the advantage of being flexible and easy to use. The user can simply look up chemical reactions from the KIDA website and construct a network without modifying the code. It has the disadvantage of being slower than the hard-coded chemical networks, since the code structure is optimized toward flexibility instead of efficiency.
If you want to write your own chemistry network, you have two choices: using the KIDA network interface, or hard-code in your own network. If you choose the later, you need to provide the chemical reaction rates and heating/cooling rates. Examples are provided in the the directory athena/src/chemistry/network/ (see e.g. H2.hpp and H2.cpp). You must also add the configure option and number of species in the configure.py file.
The chemical reaction and energy (heating/cooling) source terms are solved as coupled ordinary differential equations (ODEs) in the operator-split manner. Two ODE solvers are implemented here, the CVODE solver and the forward Euler solver.
This solver uses the open-source CVODE library, which is part of the SUNDIALS suite. CVODE uses the backward differentiation formulas and variable-order, variable-step multistep methods to solve ODEs implicitly. It is especially suited for stiff ODE systems often encountered in chemistry. The user must install the CVODE library separately, and add the installation path to the configure option --cvode_path=[cvode_path]. Many chemistry regression tests requires the CVODE library to be installed and the environment variable $CVODE_PATH to be specified. Currently, we support the CVODE version 6.2.0. Here is an example of the configure options you can use to install CVODE, on a Slurm-based cluster using the classic Intel C++ compiler with MPI and OpenMP:
wget https://github.com/LLNL/sundials/releases/download/v6.2.0/cvode-6.2.0.tar.gz
tar -xvf cvode-6.2.0.tar.gz
mkdir cvode_build
mkdir cvode_install
export CVODE_PATH=$PWD/cvode_install
cd cvode_build
cmake -D CMAKE_INSTALL_LIBDIR=lib -D CMAKE_INSTALL_PREFIX=${CVODE_PATH} -D ENABLE_MPI=ON -D ENABLE_OPENMP=ON -D CMAKE_C_COMPILER=icc -D MPI_C_COMPILER=mpiicpc -D MPIEXEC_EXECUTABLE=srun -D CMAKE_BUILD_TYPE=RelWithDebInfo ../cvode-6.2.0
make -j16
make install
cd ..
For more details, consult the CVODE library documentation.
The simple forward Euler solver uses the explicit first-order differentiation formula to solve the ODEs by sub-cycling. Since the ODE systems are often stiff in chemistry, the forward Euler solver is often unsuitable for realistic problems. We implement the forward Euler solver as a demonstration for the user to write their own ODE solver that is not dependent on the CVODE library. It is used in the regression test chem_H2.
Currently, the overall chemical timestep is the same as the hydro timestep: the ODE solver advances the chemical abundances and energy for a single hydro timestep. Note that the ODE solver uses subcycling internally, and the subcycling timestep is adjusted to ensure the solution converges, but this process is hidden from the user. In some applications, the user may want to limit the hydro timestep so that some chemical abundance or gas temperature does not change too much in a single timestep. Because this depends on the application, we leave this implementation to the user. The user can provide timestep criteria in the problem generator.
Chemistry can depend on radiation, such as in the GOW17 network. The radiation transfer method can be chosen in the configure option --chem_radiation=[choice]. Two methods are currently implemented. The trivial constant radiation and the six-ray radiation transfer. We note that other radiation modules such as non-relativistic radiation transport is also implemented in Athena++, but it has not ben coupled with chemistry yet.
You can set the radiation field to be constant in time. This simplest choice gives a template for the user to write their own radiation transfer method, and is useful for testing.
In the six-ray approximation, the radiation field is calculated by ray tracing and averaged over the six directions along the Cartesian axes. The incident radiation field is assumed to come from the edge of the computational domain along each ray. We allow the incident radiation field to be separately specified in different directions by the user, and thus the six-ray method can also be used in calculating PDR structures.
Currently, the six-ray radiation is only compatible with the GOW17 chemical network, and not the general KIDA network. Since the six-ray assume the incident radiation field comes from the physical boundary, it is currently compatible with outflow, reflective, and user-specified boundary conditions, but not the periodic boundary condition.