Initialising CHIMES

At the beginning of the simulation you will need to initialise the CHIMES chemistry solver. This will involve setting up the chemical network, based on which elements are to be included in the network as specified by the user, and then reading in the HDF5 data files and allocating memory to the various tables that store the reaction coefficients, photoionisation cross sections etc.

To initialise CHIMES, you need to call the init_chimes(*myGlobalVars) function. This takes as an argument a pointer to an instance of the globalVariables structure. Some of the variables in this structure will need to be set before the call to init_chimes(*myGlobalVars), as they define the paths where the data files are stored and details of the network such as which elements to include. Note that the speciesIndices[] and totalNumberOfSpecies variables will be set by init_chimes(*myGlobalVars) itself and do not need to be provided by the User. See the Data Structures section for details.

The data structures and arrays used to store the various reaction rate data are all defined internally within CHIMES. The call to init_chimes(*myGlobalVars) will automatically set up these structures, allocate memory to the arrays etc. as required.

Once you have initialised CHIMES, you will then need to set up the gasVariables structure for each gas particle/cell. When you first set up the abundance array for a given particle/cell, you will need to initiallise the abundances to a valid state, meaning that it needs to obey the various constraint equations, i.e. the sum of all ions/molecules of a given element needs to be equal to the total abundance of that element, and the total charge of all positive and negative ions needs to be neutral. For convenience, we provide a routine, initialise_gas_abundances(*myGasVars, *myGlobalVars), that will set all elements to a given ionisation state, as determined by the globalVariables.InitIonState parameter. For example, if globalVariables.InitIonState == 0, all elements will be set to fully neutral. If the maximum ionisation state for a given element is less than the InitIonState parameter, that element will be put in its maximum ionisation state. For example, if InitIonState == 2, all hydrogen will be singly ionised. The molecules are set to zero in this routine. It takes two arguments: a pointer to the instance of the gasVariables structure for the given particle/cell, and a pointer to the instance of the globalVariables structure.

The initialise_gas_abundances(*myGasVars, *myGlobalVars) function is very simple, as each element is put into just one ionisation state. The idea is to put the abundance array into a valid state that obeys the constraint equations. If you want your simulation to start from an equilibrium state, you can then evolve the abundances of each particle/cell for a long time (say, a Hubble time) at fixed temperature to compute the equilibrium state. Alternatively, you could compute the equilibrium state of each particle/cell beforehand (for example, using CHIMES Driver) and write the resulting equilibrium abundance arrays to the initial conditions (ICs) file for your simulation. These can then just be read in from the ICs at the beginning of the simulation. In this case you wouldn’t need to run initialise_gas_abundances(*myGasVars, *myGlobalVars) first as the equilibrium states would already obey the constraint equations.