chemistry

Currently, the chemistry module is on the developing branch "chemistry_scalar". In order to use it, checkout the branch:

$git checkout chemistry_scalar

Note that the chemistry module is not fully tested yet. Use it at your own risk, and report any problems to Munan Gong.

Formulation

Chemistry uses the basic functionality of the passive scalar module. The reaction rates act as a source term for the passive scalars. It can also interact with other parts of the code, e.g. heating and cooling rates in the energy equation, the equation of state, opacity in radiation transfer, etc.

code implementation

Because the chemical reaction rates for different species often differ by orders of magnitude, implicit solvers are suited better 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 "dt", the simulation cycle timestep, by solving a set of coupled ODEs implicitly:

In the code, the reaction rate K_i is referred to as the RHS (right-hand-side) and the de/dt (heating and cooling) is refered to as Edot. Because the chemical reaction rate Ki usually depends on the temperature 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. 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. No conservation of elemental abundances (such as by rescaling the flux of different species) has been implemented yet.

Running the Code

Prerequisite: CVODE library

The chemistry module uses the CVODE library, which the user must install on their system: https://computing.llnl.gov/projects/sundials/cvode

regression tests

Currently, there are two chemistry regression tests: One of them tests the equilibrium solution of the chemical network in Gong, Ostriker and Wolfire (2017). The other compares the analytic solution to the numerical one in a simple H2 formation network (figure below).

The regression test can be run from:

~/athena/tst/regression$python run_test.py chemistry

In order to run the regression test successfully, one needs to configure the $CVODE_PATH environment variable. This can be done by adding the environmental variable in the ~/.bashrc file

export CVODE_PATH=/usr/local/sundials (or another path where CVODE is installed)

Another interesting regression test is the convergence test (Figure above, regression test chem_H2_gaussian). On the left panel, the relative tolerance for the chemistry ODE solver (reltol) is set to be very small (1e-15). Thus, the passive scalar advection calculations dominate the error. On the right panel, reltol is set to a larger value (1e-5), and the ODE solver error dominates at large N (high resolution). Interestingly, the L1 error increases slightly with resolution in this case. This is perhaps due to the fact that passive scalar does not guarantee the elemental abundance conservation (yellow line). Future implementation of elemental abundance conservation will be tested to see whether this issue can be improved.

configuruation

./configure.py --prob=[problem] --chemistry=[network] --cvode_path=[cvode_path]

[network] is the name of the network in athena/src/chemistry/network/ directory.

input file

Sample input files can be found in the athena/input/chemistry/ directory.

output

The output names are related to the species names set in the chemistry network. For example, for a species with the name 'H2', in the vtk and hdf5 output, the conservative scalar will be 'sH2', and the primitive scalar (relative abundance Ci) will be 'rH2'. In the history output, the total mass will be 'H2'. Currently, restart output is not yet implemented.

write your own chemistry network

See example of H2.hpp and H2.cpp in the directory athena/src/chemistry/network/. You must also add the configure option and number of species in the configure.py file.