How Does Dengo Generate Solvers?¶

This tutorial tells you how the solver are generated with the help of dengo.ChemicalNetwork, and Jinja2. In short, dengo.ChemicalNetwork carries the full information of the chemical reactions and cooling actions of interest. It internally generates the symbolic representation of the dynamics of each chemical species. This can be exported as C++ or python code with a pre-written templates which can be found under dengo/templates. In this example we will be demonstrating how to generate rhs and solve the initial value problem withscipy.odeint

Import Libraries and Create the Network¶

Primordial rates and cooling for the 9-species network are included in the default dengo library in dengo.primordial_rates and dengo.primordial_cooling. The reactions and cooling are added automatically to the reaction_registry, cooling_registry and species_registry with the call to dengo.primordial_rates.setup_primordial.

Here we setup the same sample network we demonstrated in the last chapter with k01, k02 and reHII.

import dengo
from dengo.chemical_network import \
 ChemicalNetwork, \
 reaction_registry, \
 cooling_registry, species_registry
import dengo.primordial_rates
import dengo.primordial_cooling

dengo.primordial_rates.setup_primordial()

simpleNetwork = ChemicalNetwork()
simpleNetwork.add_reaction("k01")
simpleNetwork.add_reaction("k02")
simpleNetwork.add_cooling("reHII")
simpleNetwork.init_temperature((1e0, 1e8))

Adding reaction: k01 : 1*H_1 + 1*de => 1*H_2 + 2*de
Adding reaction: k02 : 1*H_2 + 1*de => 1*H_1

Building the solver¶

In this example, we will walk you through how to write a template from scratch that can be fed into a scipy solver. This can be done with ChemicalNetwork.write_solver and the combination of templates available under dengo/templates.

Evaluate the reaction rates¶

Reaction rates usually have a temperature dependence. For example, for reactions following the Arrhenius equation usually have the forms of \(k(T) = A e^{-\frac{E_a}{RT}}\), where \(k\) is the reaction rate, \(E_a\) is the activation energy of the reaction, \(T\) is the temperature, \(A\), \(R\) are the pre-exponential factor, and the universal gas constant respectively. \(A\) is sometimes depends further on temperature in Modified Arrhenius equation.

Evaluating these rates on the fly would be computationally expensive. One possible way of reducing the computational time is to interpolate from a pre-calculated reaction rates table. The rates are specified when the reactions rxn are first created with the @reaction decorator. They can be evaluated handily with rxn.coeff_fn(chemicalnetwork). The range of temperature of interest for example \(T = \rm (1, 10^8) K\) can be first specified with ChemicalNetwork.init_temperature(T_bounds=(1e0, 1e8), n_bins=1024). The added reaction objects can be accessed with ChemicalNetwork.reactions. For example, the reaction rates of k01 can the accessed with the snippet below

rxn_rate = cn.reactions['k01'].coeff_fn(ChemicalNetwork)

The output rxn_rate is an numpy array with a length of [n_bins].

A reaction rate table is generated and exported to a hdf5 file below.

with h5py.File(ofn, 'r') as f:
    print(f.keys())
>>> <KeysViewHDF5 ['T', 'k01', 'k02', 'reHII_reHII']>

import os
import h5py

solver_name = 'simple'
output_dir  = "."

ofn = os.path.join(output_dir, f"{solver_name}_tables.h5")

if os.path.exists(ofn):
    os.system(f"rm {ofn}")

with h5py.File(ofn, 'w') as f:
    for rxn in sorted(simpleNetwork.reactions.values()):
        f.create_dataset(
            f"/{rxn.name}", data=rxn.coeff_fn(simpleNetwork).astype("float64")
        )

    for rxn in sorted(simpleNetwork.cooling_actions.values()):
        if hasattr(rxn, "tables"):
            for tab in rxn.tables:
                f.create_dataset(
                    f"/{rxn.name}_{tab}",
                    data=rxn.tables[tab](simpleNetwork).astype("float64"),
                )
    f.create_dataset(f"/T", data=simpleNetwork.T.astype("float64"))

Using `Jinja2`¶

Jinja2 fills out variables dynamically with the user-given inputs. It is however not limited to variables subsitution. You can also use conditions (if/else), for loops (for), filter, and tests in Jinja2 template. To construct our solver template for scipy, we will be using for-loop, if-else, and variable substitution. And we will go through these usages briefly below their respective syntax. We refer our readers for more details in Jinja2 Documentation.

For example if the template is Dengo {{token}}, Jinja2 would look for the user-fed input for the variable token, and output the file with the word in curly bracket replaced with the token = "works".

from jinja2 import Template
t = Template("Dengo {{ token }}!")
print(t.render(token="works"))

Dengo works!

Jinja2 also provided an straightforward way to iterate over iterables like lists, with the syntax {% for item in list %} {% endfor %}. In the example below, we iterate over the required species in our chemical network, and had Jinja2 print out their name.

t = Template("""
{% for s in network.required_species | sort %}
{{s.name}},
{% endfor %}
""")
print(t.render(network = simpleNetwork))

H_1,

H_2,

de,

ge,

Preparing the rates interpolation template¶

In the template below, we demonstrate a simple implementation for interpolating reaction and cooling rates with Jinja2 and our ChemicalNetwork object. In brief, the template below first reads in the rates from the hdf5 we created above, and two functions interpolate_rates and interpolate_cooling_rates are created dynamically based on the reactions avaiable in our ChemicalNetwork.

reaction_rates_template = Template("""
import numpy as np
import h5py

# read rates in as global variables
rates_table = "{{solver_name}}_tables.h5"
ratef = h5py.File(rates_table, 'r')

# Reaction Rates
{% for k in network.reactions.keys()%}
out{{k}}dev = ratef['{{k}}'][:]
{%- endfor %} 

# Cooling Rates
{%- for name, rate in network.cooling_actions | dictsort %}
{%- for name2 in rate.tables | sort %}
out_{{name}}_{{name2}} = ratef["{{name}}_{{name2}}"][:]
{%- endfor %}
{%- endfor %}
tdev = ratef['T'][:]
ratef.close()

def interpolate_rates(T):
    {% for k in network.reactions.keys()%}
    {{k}} = np.interp(T, tdev, out{{k}}dev)
    {%- endfor %} 
    return (
    {%- for k in network.reactions.keys() -%}
    {{k}}, 
    {%- endfor -%}
    )
def interpolate_cooling_rates(T):
    {%- for name, rate in network.cooling_actions | dictsort %}
    {%- for name2 in rate.tables | sort %}
    {{name}}_{{name2}} = np.interp(T, tdev, out_{{name}}_{{name2}})
    {%- endfor %}
    {%- endfor %}
    return (
    {%- for name, rate in network.cooling_actions | dictsort -%}
    {%- for name2 in rate.tables | sort -%}
    {{name}}_{{name2}}, 
    {%- endfor -%}
    {%- endfor -%}
    )
    """)

template_vars = dict(
    network=simpleNetwork, 
    solver_name=solver_name
)
output_file = f'{solver_name}_reaction_rate.py'
with open(output_file, 'w') as f:
    f.write(reaction_rates_template.render(template_vars))

Import the Rendered `simple_reaction_rate.py` and Check if it works!¶

from simple_reaction_rate import interpolate_rates, interpolate_cooling_rates
import matplotlib.pyplot as plt

temperature = simpleNetwork.T
rxnk01_rate, rxnk02_rate = interpolate_rates(temperature)
coolreHII_rate,          = interpolate_cooling_rates(temperature)

plt.loglog(temperature, rxnk01_rate, label='k01')
plt.loglog(temperature, rxnk02_rate, label='k02')
plt.xlabel(r'Temperature $(\rm K)$')
plt.ylabel(r'Reaction Rate $(\rm g^{-1} cm^{3} s^{-1} )$')
plt.legend()
ax2 = plt.twinx()
ax2.loglog(temperature, coolreHII_rate,color='C2', label='reHII')
plt.ylabel(r'Cooling Rate reHII')
plt.legend()

<matplotlib.legend.Legend at 0x7efd3c629550>

Dengo Cookbook

How Does Dengo Generate Solvers?

Contents

How Does Dengo Generate Solvers?¶

Import Libraries and Create the Network¶

Building the solver¶

Evaluate the reaction rates¶

Using `Jinja2`¶

Preparing the rates interpolation template¶

Import the Rendered `simple_reaction_rate.py` and Check if it works!¶

Implement the temperature function¶

Implement the RHS function¶

Yay! you have learnt how to write a Jinja2 template from scratch¶

Dengo Cookbook

How Does Dengo Generate Solvers?

Contents

How Does Dengo Generate Solvers?¶

Import Libraries and Create the Network¶

Building the solver¶

Evaluate the reaction rates¶

Using Jinja2¶

Preparing the rates interpolation template¶

Import the Rendered simple_reaction_rate.py and Check if it works!¶

Implement the temperature function¶

Implement the RHS function¶

Yay! you have learnt how to write a Jinja2 template from scratch¶

Using `Jinja2`¶

Import the Rendered `simple_reaction_rate.py` and Check if it works!¶