Rate Functions

Grackle supports the calculation of individual rate coefficients at specific temperatures. When running Grackle’s chemistry solving schemes, these coefficients are computed internally, at the temperatures required by the solver. However, there may be cases where obtaining a specific rate coefficient at a given temperature is desirable – with Grackle’s initialisation routines ported from Fortran to C, this is now possible. Documented here is an overview of how to use these functions, alongside relevant examples.

All rate coefficients are designated an individual function for their calculation; their definitions are located within rate_coefficients.c and their prototypes in rate_coefficients.h, both found within the grackle/src/clib/ directory. The vast majority of the rate coefficient’s functions are of the same general form, however there are exceptions which will be discussed individually.

The General Rate Function

Structure

As mentioned previously, most rate functions have the same structure and can therefore be called identically as below, where RATE_NAME is the name of the coefficient you wish to calculate.

double RATE_NAME_rate(double T, double units, chemistry_data *my_chemistry);

This returns the value of the rate at a given temperature. The calculated rate will be divided through by the value given in the units argument. If set to 1, the rate will be returned in CGS units.

Parameters:

  • T (double) – gas temperature

  • units (double) – units

  • *my_chemistry (chemistry_data) – pointer to chemistry_data struct containing parameters.

Return type:

double

Returns:

value of rate coefficient

Examples

Example 1

Suppose we wish to print the value of the k1 rate coefficient at a temperature of 1000 K. Given we have already configured our chemistry parameters within a chemistry_data struct named my_chemistry, we can obtain the result by simply calling the function as follows:

#include "grackle_rate_functions.h"

double result = k1_rate(1e3, 1., &my_chemistry);
printf("k1 %e", result);

where our result is in cgs units.

Example 2

The rate coefficient reHII can be calculated via two different methods depending on the status of the CaseBRecombination parameter within the chemistry_data struct (named my_chemistry in this example), which can take values of either 0 or 1. Suppose we wish to compare the results of each calculation method over an arbitrary temperature range, we can achieve this by the following:

#include "grackle_rate_functions.h"

// Define temperature range to calculate coefficients over
double tempStart = 1e1;
double tempEnd = 1e8;
double numTemps = 1e3;
double tempSpacing = (tempEnd - tempStart) / numTemps;

// Create arrays for results storage
double caseAResults[(int) numTemps];
double caseBResults[(int) numTemps];

// Set value of my_chemistry.CaseBRecombination
for (int caseB = 0; caseB < 2; caseB++) {
    my_chemistry.CaseBRecombination = caseB;
    // Iterate over temperature range
    for (int i = 0; i < numTemps; i++) {
        double temp = tempStart + i*tempSpacing;
        // Store results in appropriate array
        if (caseB == 0) {
            caseAResults[i] = reHII_rate(temp, 1., &my_chemistry);
        } else {
            caseBResults[i] = reHII_rate(temp, 1., &my_chemistry);
        }
    }
}

where we have created an array of reHII coefficients for both settings of chemistry_data.CaseBRecombination over the same temperature range.

The k13dd Rate Function

Structure

The k13dd rate function, which describes the density-dependent dissociation of molecular hydrogen, is similar in form to the general rate functions, the only difference being its additional input parameter. This is a pointer to an array of length 14 * sizeof(double), which will hold the outputs of the function. The function always calculates fourteen rate parameters, the first seven of which correspond to direct collisional dissociation, whilst the latter seven correspond to dissociative tunneling – please see Martin, Schwarz & Mandy, 1996 for further details on how these are calculated. The structure of the function is then:

void k13dd_rate(double T, double units, double *results_array, chemistry_data *my_chemistry);

Calculates the density-dependent molecular hydrogen dissociation rate. The calculated rate will be divided through by the value given in the units argument. If set to 1, the units will be CGS.

Parameters:

  • T (double) – gas temperature

  • units (double) – units

  • *results_array (double) – pointer to array of length 14 * sizeof(double) in which the calculated rate coefficients will be stored.

  • *my_chemistry (chemistry_data) – pointer to chemistry_data struct containing parameters.

Return type:

void

Returns:

Nothing is returned, but values are set for the array pointed to by the results_array pointer.

Examples

Suppose we would like to print the rate coefficients for the dissociation of molecular hydrogen via the tunneling process at a temperature of 1e5 K. Given we have already configured our chemistry parameters within a chemistry_data struct named my_chemistry, we can obtain the coefficients by the following:

#include "grackle_rate_functions.h"

// Create an array of the correct size for result storage.
double results[14];

// Call the function at the desired temperature, getting results in cgs units.
k13dd_rate(1e5, 1., &results, &my_chemistry);

// Print the results corresponding to dissociative tunneling.
for (int i = 7; i < 14; i++) {
    printf("k13dd %e", results[i]);
}

The h2dust Rate Function

Structure

The h2dust rate function, which describes the formation of molecular hydrogen on dust grains, is similar in form to the general rate functions, the only difference being its additional input parameter; a double which represents the dust temperature. The function returns a double just as the general rate function, its structure is then:

double h2dust_rate(double T, double T_dust, double units, chemistry_data *my_chemistry);
Parameters:

  • T (double) – gas temperature

  • T_gas (double) – dust temperature

  • units (double) – units

  • *my_chemistry (chemistry_data) – pointer to chemistry_data struct containing parameters.

Return type:

double

Returns:

value of rate coefficient

Examples

Suppose we would like to calculate the h2dust rate coefficients for a gas temperature of 1e4 K, with a varying dust temperature. Given we have already configured our chemistry parameters within a :c:data`chemistry_data` struct named my_chemistry, we can obtain the coefficients by the following:

#include "grackle_rate_functions.h"

// Define dust temperature range to calculate coefficients over.
double tempStart_dust = 1e1;
double tempEnd_dust = 1e6;
double numTemps_dust = 1e3;
double tempSpacing_dust = (tempEnd_dust - tempStart_dust) / numTemps_dust;

// Create array for results storage.
double h2dust_results[(int) numTemps_dust];

// Loop over dust temperatures.
for (int i=0; i < numTemps_dust; i++){
    double temp_dust = tempStart_dust + i*tempSpacing_dust;
    h2dust_results[i] = h2dust_rate(1e4, temp_dust, 1., &my_chemistry);
}

The Scalar Rate Functions

Structure

The scalar rate functions (comp, gammah, gamma_isrf) are simpler than the general rate functions due to their temperature independence. They require only two inputs and return a single double, their structure is as below, where SCALAR_NAME is the name of the scalar rate coefficient you wish to calculate. These are called in the same way as the general rate functions, ignoring the temperature dependancy – please see their documentation for basic examples.

double SCALAR_NAME_rate(double units, chemistry_data *my_chemistry);

This returns the value of the rate. The calculated rate will be divided through by the value given in the units argument. If set to 1, the rate will be returned in CGS units.

Parameters:

  • units (double) – units

  • *my_chemistry (chemistry_data) – pointer to chemistry_data struct containing parameters.

Return type:

double

Returns:

value of rate coefficient