Example of integrating a network

This notebook illustrates how to create a python network and integrate it with the scipy library.

[1]:
import pynucastro as pyna

We’ll start again with the basic CNO network explored earlier. Again, we’ll read in the entire ReacLib library and pass in the names of the rates in the form of a string A(x,y)B to filter out just the rates we are interested in.

[2]:
rl = pyna.ReacLibLibrary()
[3]:
rate_names = ["c12(p,g)n13",
              "c13(p,g)n14",
              "n13(,)c13",
              "n13(p,g)o14",
              "n14(p,g)o15",
              "n15(p,a)c12",
              "o14(,)n14",
              "o15(,)n15"]

rates = rl.get_rate_by_name(rate_names)
rates
[3]:
[C12 + p ⟶ N13 + 𝛾,
 C13 + p ⟶ N14 + 𝛾,
 N13 ⟶ C13 + e⁺ + 𝜈,
 N13 + p ⟶ O14 + 𝛾,
 N14 + p ⟶ O15 + 𝛾,
 N15 + p ⟶ He4 + C12,
 O14 ⟶ N14 + e⁺ + 𝜈,
 O15 ⟶ N15 + e⁺ + 𝜈]

A PythonNetwork is based on a RateCollection but has methods to write the RHS of the system of ODEs.

[4]:
pynet = pyna.PythonNetwork(rates=rates)
[5]:
fig = pynet.plot()
_images/pynucastro-integration_7_0.png

For example, this network knows how to write the full term for a reaction that goes into the \(dY/dt\) equation of the ODE system.

Here we pick one of the rates that is part of the network an explore it.

[6]:
r = pynet.rates[1]
print(r)
C13 + p ⟶ N14 + 𝛾

a rate also knows what its contribution is to the \(dY/dt\) equation is:

[7]:
print(r.ydot_string_py())
rho*Y[jp]*Y[jc13]*rate_eval.p_C13__N14

and the code needed to evaluate that rate (the T-dependent part) is output by the Rate class:

..
[8]:
print(r.function_string_py())
@numba.njit()
def p_C13__N14(rate_eval, tf):
    # C13 + p --> N14
    rate = 0.0

    # nacrr
    rate += np.exp(  15.1825 + -13.5543*tf.T9i
                  + -1.5*tf.lnT9)
    # nacrn
    rate += np.exp(  18.5155 + -13.72*tf.T913i + -0.450018*tf.T913
                  + 3.70823*tf.T9 + -1.70545*tf.T953 + -0.666667*tf.lnT9)
    # nacrr
    rate += np.exp(  13.9637 + -5.78147*tf.T9i + -0.196703*tf.T913
                  + 0.142126*tf.T9 + -0.0238912*tf.T953 + -1.5*tf.lnT9)

    rate_eval.p_C13__N14 = rate


The temperature-dependent rate evaluation functions take a Tfactor object, which precomputes most of the commonly-used temperature factors in the rates.

The write_network() method will output the python code needed to define the RHS of a network for integration with the SciPy integrators.

Since python code can be slow, we use Numba to do just-in-time compilation of the functions to speed things up.

[9]:
pynet.write_network("cno_test_integrate.py")
[10]:
%cat cno_test_integrate.py
import numba
import numpy as np
from scipy import constants
from numba.experimental import jitclass

from pynucastro.rates import TableIndex, TableInterpolator, TabularRate, Tfactors
from pynucastro.screening import PlasmaState, ScreenFactors

jp = 0
jhe4 = 1
jc12 = 2
jc13 = 3
jn13 = 4
jn14 = 5
jn15 = 6
jo14 = 7
jo15 = 8
nnuc = 9

A = np.zeros((nnuc), dtype=np.int32)

A[jp] = 1
A[jhe4] = 4
A[jc12] = 12
A[jc13] = 13
A[jn13] = 13
A[jn14] = 14
A[jn15] = 15
A[jo14] = 14
A[jo15] = 15

Z = np.zeros((nnuc), dtype=np.int32)

Z[jp] = 1
Z[jhe4] = 2
Z[jc12] = 6
Z[jc13] = 6
Z[jn13] = 7
Z[jn14] = 7
Z[jn15] = 7
Z[jo14] = 8
Z[jo15] = 8

# masses in ergs
mass = np.zeros((nnuc), dtype=np.float64)

mass[jp] = 0.0015040963030260536
mass[jhe4] = 0.0059735574925878256
mass[jc12] = 0.017909017027273523
mass[jc13] = 0.019406441930882663
mass[jn13] = 0.01940999951603316
mass[jn14] = 0.020898440903103103
mass[jn15] = 0.0223864338056853
mass[jo14] = 0.020906683076491985
mass[jo15] = 0.02239084645968795

names = []
names.append("H1")
names.append("He4")
names.append("C12")
names.append("C13")
names.append("N13")
names.append("N14")
names.append("N15")
names.append("O14")
names.append("O15")

def to_composition(Y):
    """Convert an array of molar fractions to a Composition object."""
    from pynucastro import Composition, Nucleus
    nuclei = [Nucleus.from_cache(name) for name in names]
    comp = Composition(nuclei)
    for i, nuc in enumerate(nuclei):
        comp.X[nuc] = Y[i] * A[i]
    return comp


def energy_release(dY):
    """return the energy release in erg/g (/s if dY is actually dY/dt)"""
    enuc = 0.0
    for i, y in enumerate(dY):
        enuc += y * mass[i]
    enuc *= -1*constants.Avogadro
    return enuc

@jitclass([
    ("p_C12__N13", numba.float64),
    ("p_C13__N14", numba.float64),
    ("N13__C13__weak__wc12", numba.float64),
    ("p_N13__O14", numba.float64),
    ("p_N14__O15", numba.float64),
    ("p_N15__He4_C12", numba.float64),
    ("O14__N14__weak__wc12", numba.float64),
    ("O15__N15__weak__wc12", numba.float64),
])
class RateEval:
    def __init__(self):
        self.p_C12__N13 = np.nan
        self.p_C13__N14 = np.nan
        self.N13__C13__weak__wc12 = np.nan
        self.p_N13__O14 = np.nan
        self.p_N14__O15 = np.nan
        self.p_N15__He4_C12 = np.nan
        self.O14__N14__weak__wc12 = np.nan
        self.O15__N15__weak__wc12 = np.nan

@numba.njit()
def ye(Y):
    return np.sum(Z * Y)/np.sum(A * Y)

@numba.njit()
def p_C12__N13(rate_eval, tf):
    # C12 + p --> N13
    rate = 0.0

    # ls09n
    rate += np.exp(  17.1482 + -13.692*tf.T913i + -0.230881*tf.T913
                  + 4.44362*tf.T9 + -3.15898*tf.T953 + -0.666667*tf.lnT9)
    # ls09r
    rate += np.exp(  17.5428 + -3.77849*tf.T9i + -5.10735*tf.T913i + -2.24111*tf.T913
                  + 0.148883*tf.T9 + -1.5*tf.lnT9)

    rate_eval.p_C12__N13 = rate

@numba.njit()
def p_C13__N14(rate_eval, tf):
    # C13 + p --> N14
    rate = 0.0

    # nacrr
    rate += np.exp(  15.1825 + -13.5543*tf.T9i
                  + -1.5*tf.lnT9)
    # nacrn
    rate += np.exp(  18.5155 + -13.72*tf.T913i + -0.450018*tf.T913
                  + 3.70823*tf.T9 + -1.70545*tf.T953 + -0.666667*tf.lnT9)
    # nacrr
    rate += np.exp(  13.9637 + -5.78147*tf.T9i + -0.196703*tf.T913
                  + 0.142126*tf.T9 + -0.0238912*tf.T953 + -1.5*tf.lnT9)

    rate_eval.p_C13__N14 = rate

@numba.njit()
def N13__C13__weak__wc12(rate_eval, tf):
    # N13 --> C13
    rate = 0.0

    # wc12w
    rate += np.exp(  -6.7601)

    rate_eval.N13__C13__weak__wc12 = rate

@numba.njit()
def p_N13__O14(rate_eval, tf):
    # N13 + p --> O14
    rate = 0.0

    # lg06r
    rate += np.exp(  10.9971 + -6.12602*tf.T9i + 1.57122*tf.T913i
                  + -1.5*tf.lnT9)
    # lg06n
    rate += np.exp(  18.1356 + -15.1676*tf.T913i + 0.0955166*tf.T913
                  + 3.0659*tf.T9 + -0.507339*tf.T953 + -0.666667*tf.lnT9)

    rate_eval.p_N13__O14 = rate

@numba.njit()
def p_N14__O15(rate_eval, tf):
    # N14 + p --> O15
    rate = 0.0

    # im05n
    rate += np.exp(  17.01 + -15.193*tf.T913i + -0.161954*tf.T913
                  + -7.52123*tf.T9 + -0.987565*tf.T953 + -0.666667*tf.lnT9)
    # im05r
    rate += np.exp(  6.73578 + -4.891*tf.T9i
                  + 0.0682*tf.lnT9)
    # im05r
    rate += np.exp(  7.65444 + -2.998*tf.T9i
                  + -1.5*tf.lnT9)
    # im05n
    rate += np.exp(  20.1169 + -15.193*tf.T913i + -4.63975*tf.T913
                  + 9.73458*tf.T9 + -9.55051*tf.T953 + 0.333333*tf.lnT9)

    rate_eval.p_N14__O15 = rate

@numba.njit()
def p_N15__He4_C12(rate_eval, tf):
    # N15 + p --> He4 + C12
    rate = 0.0

    # nacrn
    rate += np.exp(  27.4764 + -15.253*tf.T913i + 1.59318*tf.T913
                  + 2.4479*tf.T9 + -2.19708*tf.T953 + -0.666667*tf.lnT9)
    # nacrr
    rate += np.exp(  -6.57522 + -1.1638*tf.T9i + 22.7105*tf.T913
                  + -2.90707*tf.T9 + 0.205754*tf.T953 + -1.5*tf.lnT9)
    # nacrr
    rate += np.exp(  20.8972 + -7.406*tf.T9i
                  + -1.5*tf.lnT9)
    # nacrr
    rate += np.exp(  -4.87347 + -2.02117*tf.T9i + 30.8497*tf.T913
                  + -8.50433*tf.T9 + -1.54426*tf.T953 + -1.5*tf.lnT9)

    rate_eval.p_N15__He4_C12 = rate

@numba.njit()
def O14__N14__weak__wc12(rate_eval, tf):
    # O14 --> N14
    rate = 0.0

    # wc12w
    rate += np.exp(  -4.62354)

    rate_eval.O14__N14__weak__wc12 = rate

@numba.njit()
def O15__N15__weak__wc12(rate_eval, tf):
    # O15 --> N15
    rate = 0.0

    # wc12w
    rate += np.exp(  -5.17053)

    rate_eval.O15__N15__weak__wc12 = rate

def rhs(t, Y, rho, T, screen_func=None):
    return rhs_eq(t, Y, rho, T, screen_func)

@numba.njit()
def rhs_eq(t, Y, rho, T, screen_func):

    tf = Tfactors(T)
    rate_eval = RateEval()

    # reaclib rates
    p_C12__N13(rate_eval, tf)
    p_C13__N14(rate_eval, tf)
    N13__C13__weak__wc12(rate_eval, tf)
    p_N13__O14(rate_eval, tf)
    p_N14__O15(rate_eval, tf)
    p_N15__He4_C12(rate_eval, tf)
    O14__N14__weak__wc12(rate_eval, tf)
    O15__N15__weak__wc12(rate_eval, tf)

    if screen_func is not None:
        plasma_state = PlasmaState(T, rho, Y, Z)

        scn_fac = ScreenFactors(1, 1, 6, 12)
        scor = screen_func(plasma_state, scn_fac)
        rate_eval.p_C12__N13 *= scor

        scn_fac = ScreenFactors(1, 1, 6, 13)
        scor = screen_func(plasma_state, scn_fac)
        rate_eval.p_C13__N14 *= scor

        scn_fac = ScreenFactors(1, 1, 7, 13)
        scor = screen_func(plasma_state, scn_fac)
        rate_eval.p_N13__O14 *= scor

        scn_fac = ScreenFactors(1, 1, 7, 14)
        scor = screen_func(plasma_state, scn_fac)
        rate_eval.p_N14__O15 *= scor

        scn_fac = ScreenFactors(1, 1, 7, 15)
        scor = screen_func(plasma_state, scn_fac)
        rate_eval.p_N15__He4_C12 *= scor

    dYdt = np.zeros((nnuc), dtype=np.float64)

    dYdt[jp] = (
       -rho*Y[jp]*Y[jc12]*rate_eval.p_C12__N13
       -rho*Y[jp]*Y[jc13]*rate_eval.p_C13__N14
       -rho*Y[jp]*Y[jn13]*rate_eval.p_N13__O14
       -rho*Y[jp]*Y[jn14]*rate_eval.p_N14__O15
       -rho*Y[jp]*Y[jn15]*rate_eval.p_N15__He4_C12
       )

    dYdt[jhe4] = (
       +rho*Y[jp]*Y[jn15]*rate_eval.p_N15__He4_C12
       )

    dYdt[jc12] = (
       -rho*Y[jp]*Y[jc12]*rate_eval.p_C12__N13
       +rho*Y[jp]*Y[jn15]*rate_eval.p_N15__He4_C12
       )

    dYdt[jc13] = (
       -rho*Y[jp]*Y[jc13]*rate_eval.p_C13__N14
       +Y[jn13]*rate_eval.N13__C13__weak__wc12
       )

    dYdt[jn13] = (
       -Y[jn13]*rate_eval.N13__C13__weak__wc12
       -rho*Y[jp]*Y[jn13]*rate_eval.p_N13__O14
       +rho*Y[jp]*Y[jc12]*rate_eval.p_C12__N13
       )

    dYdt[jn14] = (
       -rho*Y[jp]*Y[jn14]*rate_eval.p_N14__O15
       +rho*Y[jp]*Y[jc13]*rate_eval.p_C13__N14
       +Y[jo14]*rate_eval.O14__N14__weak__wc12
       )

    dYdt[jn15] = (
       -rho*Y[jp]*Y[jn15]*rate_eval.p_N15__He4_C12
       +Y[jo15]*rate_eval.O15__N15__weak__wc12
       )

    dYdt[jo14] = (
       -Y[jo14]*rate_eval.O14__N14__weak__wc12
       +rho*Y[jp]*Y[jn13]*rate_eval.p_N13__O14
       )

    dYdt[jo15] = (
       -Y[jo15]*rate_eval.O15__N15__weak__wc12
       +rho*Y[jp]*Y[jn14]*rate_eval.p_N14__O15
       )

    return dYdt

def jacobian(t, Y, rho, T, screen_func=None):
    return jacobian_eq(t, Y, rho, T, screen_func)

@numba.njit()
def jacobian_eq(t, Y, rho, T, screen_func):

    tf = Tfactors(T)
    rate_eval = RateEval()

    # reaclib rates
    p_C12__N13(rate_eval, tf)
    p_C13__N14(rate_eval, tf)
    N13__C13__weak__wc12(rate_eval, tf)
    p_N13__O14(rate_eval, tf)
    p_N14__O15(rate_eval, tf)
    p_N15__He4_C12(rate_eval, tf)
    O14__N14__weak__wc12(rate_eval, tf)
    O15__N15__weak__wc12(rate_eval, tf)

    if screen_func is not None:
        plasma_state = PlasmaState(T, rho, Y, Z)

        scn_fac = ScreenFactors(1, 1, 6, 12)
        scor = screen_func(plasma_state, scn_fac)
        rate_eval.p_C12__N13 *= scor

        scn_fac = ScreenFactors(1, 1, 6, 13)
        scor = screen_func(plasma_state, scn_fac)
        rate_eval.p_C13__N14 *= scor

        scn_fac = ScreenFactors(1, 1, 7, 13)
        scor = screen_func(plasma_state, scn_fac)
        rate_eval.p_N13__O14 *= scor

        scn_fac = ScreenFactors(1, 1, 7, 14)
        scor = screen_func(plasma_state, scn_fac)
        rate_eval.p_N14__O15 *= scor

        scn_fac = ScreenFactors(1, 1, 7, 15)
        scor = screen_func(plasma_state, scn_fac)
        rate_eval.p_N15__He4_C12 *= scor

    jac = np.zeros((nnuc, nnuc), dtype=np.float64)

    jac[jp, jp] = (
       -rho*Y[jc12]*rate_eval.p_C12__N13
       -rho*Y[jc13]*rate_eval.p_C13__N14
       -rho*Y[jn13]*rate_eval.p_N13__O14
       -rho*Y[jn14]*rate_eval.p_N14__O15
       -rho*Y[jn15]*rate_eval.p_N15__He4_C12
       )

    jac[jp, jc12] = (
       -rho*Y[jp]*rate_eval.p_C12__N13
       )

    jac[jp, jc13] = (
       -rho*Y[jp]*rate_eval.p_C13__N14
       )

    jac[jp, jn13] = (
       -rho*Y[jp]*rate_eval.p_N13__O14
       )

    jac[jp, jn14] = (
       -rho*Y[jp]*rate_eval.p_N14__O15
       )

    jac[jp, jn15] = (
       -rho*Y[jp]*rate_eval.p_N15__He4_C12
       )

    jac[jhe4, jp] = (
       +rho*Y[jn15]*rate_eval.p_N15__He4_C12
       )

    jac[jhe4, jn15] = (
       +rho*Y[jp]*rate_eval.p_N15__He4_C12
       )

    jac[jc12, jp] = (
       -rho*Y[jc12]*rate_eval.p_C12__N13
       +rho*Y[jn15]*rate_eval.p_N15__He4_C12
       )

    jac[jc12, jc12] = (
       -rho*Y[jp]*rate_eval.p_C12__N13
       )

    jac[jc12, jn15] = (
       +rho*Y[jp]*rate_eval.p_N15__He4_C12
       )

    jac[jc13, jp] = (
       -rho*Y[jc13]*rate_eval.p_C13__N14
       )

    jac[jc13, jc13] = (
       -rho*Y[jp]*rate_eval.p_C13__N14
       )

    jac[jc13, jn13] = (
       +rate_eval.N13__C13__weak__wc12
       )

    jac[jn13, jp] = (
       -rho*Y[jn13]*rate_eval.p_N13__O14
       +rho*Y[jc12]*rate_eval.p_C12__N13
       )

    jac[jn13, jc12] = (
       +rho*Y[jp]*rate_eval.p_C12__N13
       )

    jac[jn13, jn13] = (
       -rate_eval.N13__C13__weak__wc12
       -rho*Y[jp]*rate_eval.p_N13__O14
       )

    jac[jn14, jp] = (
       -rho*Y[jn14]*rate_eval.p_N14__O15
       +rho*Y[jc13]*rate_eval.p_C13__N14
       )

    jac[jn14, jc13] = (
       +rho*Y[jp]*rate_eval.p_C13__N14
       )

    jac[jn14, jn14] = (
       -rho*Y[jp]*rate_eval.p_N14__O15
       )

    jac[jn14, jo14] = (
       +rate_eval.O14__N14__weak__wc12
       )

    jac[jn15, jp] = (
       -rho*Y[jn15]*rate_eval.p_N15__He4_C12
       )

    jac[jn15, jn15] = (
       -rho*Y[jp]*rate_eval.p_N15__He4_C12
       )

    jac[jn15, jo15] = (
       +rate_eval.O15__N15__weak__wc12
       )

    jac[jo14, jp] = (
       +rho*Y[jn13]*rate_eval.p_N13__O14
       )

    jac[jo14, jn13] = (
       +rho*Y[jp]*rate_eval.p_N13__O14
       )

    jac[jo14, jo14] = (
       -rate_eval.O14__N14__weak__wc12
       )

    jac[jo15, jp] = (
       +rho*Y[jn14]*rate_eval.p_N14__O15
       )

    jac[jo15, jn14] = (
       +rho*Y[jp]*rate_eval.p_N14__O15
       )

    jac[jo15, jo15] = (
       -rate_eval.O15__N15__weak__wc12
       )

    return jac

We can now import the network that was just created and integrate it using the SciPy ODE solvers

[11]:
import cno_test_integrate as cno

Integrating the network

We can use the stiff ODE integration solvers that are part of SciPy to integrate this system now

[12]:
from scipy.integrate import solve_ivp
import numpy as np

Initialize the thermodynamic conditions and initial composition. We express the composition as molar fractions, Y0.

[13]:
rho = 150
T = 1.5e7

X0 = np.zeros(cno.nnuc)
X0[cno.jp] = 0.7
X0[cno.jhe4] = 0.28
X0[cno.jc12] = 0.02

Y0 = X0/cno.A

Now we integrate. We use the BDF method, since reaction networks are in general stiff

[14]:
tmax = 1.e20

sol = solve_ivp(cno.rhs, [0, tmax], Y0, method="BDF", jac=cno.jacobian,
                dense_output=True, args=(rho, T), rtol=1.e-6, atol=1.e-6)

Plotting the results

[15]:
import matplotlib.pyplot as plt
[16]:
fig = plt.figure()
ax = fig.add_subplot(111)

for i in range(cno.nnuc):
    ax.loglog(sol.t, sol.y[i,:] * cno.A[i], label=f"X({cno.names[i].capitalize()})")

ax.set_xlim(1.e10, 1.e20)
ax.set_ylim(1.e-8, 1.0)
ax.legend(fontsize="small")
ax.set_xlabel("t (s)")
ax.set_ylabel("X")

fig.set_size_inches((8, 6))
_images/pynucastro-integration_28_0.png

We can get the energy release from the change in molar abundances, \(\Delta Y\). This will be computed as:

\[\Delta E = -N_A \sum_i \Delta Y_i M_i\]

Here \(M_i\) is the mass of the nucleus \(i\) (this is tabulated in the network as mass[]).

[17]:
E = cno.energy_release(sol.y[:,-1] - Y0)
print(f"E = {E:20.10g}")
E =      4.510854175e+18

This result is in erg/g

We can compute the instantaneous energy generation rate as well,

\[\epsilon = -N_A \sum_i \frac{dY_i}{dt} M_i\]
[18]:
epsilon = cno.energy_release(cno.rhs(0.0, Y0, rho, T))
print(f"epsilon = {epsilon:20.10g}")
epsilon =          117.7737713