Non-Equilibrium Chemistry

This example shows some of the inner-workings of the chemical network and solver using a simple hydrogen-only test problem in an isothermal medium.

To begin, first import a few things:

import ares
import numpy as np
import matplotlib.pyplot as pl

Let’s initialize a grid of 64 cells:

# Initialize grid object
grid = ares.static.Grid(grid_cells=64)

So far, this grid object only knows how many cells it has. To give it some physical properties, we’ll call several setter routines:

# Set initial conditions
grid.set_physics(isothermal=True)
grid.set_chemistry(include_He=False)
grid.set_density(nH=1.)
grid.set_ionization()
grid.set_temperature(np.logspace(3, 5, 64))

The above commands initialize the grid to be isothermal, composed of hydrogen only, with a density of 1 hydrogen atom per cubic centimeter, initialized to be neutral and with temperatures between \(10^3 \leq T /\ \mathrm{K} \leq 10^5\).

To see how the ion fractions evolve with time, we can pass this grid off to the chemistry solver:

# Initialize chemistry network / solver
chem = ares.solvers.Chemistry(grid, rt=False)

# Compute rate coefficients (only need to do this once; isothermal)
chem.chemnet.SourceIndependentCoefficients(grid.data['Tk'])

Now, to actually run the thing:

# To compute timestep
timestep = ares.util.RestrictTimestep(grid)

# Set initial time-step and maximum allowed change
data = grid.data
dt = ares.physics.Constants.s_per_myr / 1e3
dt_max = 1e2 * ares.physics.Constants.s_per_myr
t = 0.0
tf = ares.physics.Constants.s_per_gyr

# Initialize progress bar [optional]
pb = ares.util.ProgressBar(tf)
pb.start()

# Start calculation
while t < tf:
    pb.update(t)

    # Evolve system for time dt
    data = chem.Evolve(data, t=t, dt=dt)
    t += dt

    # Limit time-step based on maximum rate of change in grid quantities
    new_dt = timestep.Limit(chem.chemnet.q, chem.chemnet.dqdt)

    # Limit to factor of 2x increase in timestep
    dt = min(new_dt, 2 * dt)

    # Impose maximum timestep
    dt = min(dt, dt_max)

    # Make sure we end at t == tf
    dt = min(dt, tf - t)

pb.finish()

All of this work is done for you each time you call ares.simulations.Global21cm and ares.simulations.RaySegment.

To visualize the results:

ax = pl.subplot(111)
ax.loglog(T, data['h_1'], color='k', ls='-')
ax.loglog(T, data['h_2'], color='k', ls='--')
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel(r'$T \ (\mathrm{K})$')
ax.set_ylabel('Species Fraction')
ax.set_ylim(1e-4, 2)
pl.draw()