Integration Example

We can use pynucastro to generate the righthand side function for an astrophysical reaction network.

We’ll create a H=burning network that has the rates for all p+p and the basic CNO cycle

import numpy as np
import matplotlib.pyplot as plt
import pynucastro as pyna

rate_names = ["p(p,)d",
              "d(p,g)he3",
              "he3(he3,pp)he4",
              "c12(p,g)n13",
              "c13(p,g)n14",
              "n13(,)c13",
              "n14(p,g)o15",
              "n15(p,a)c12",
              "o15(,)n15"]

rl = pyna.ReacLibLibrary()
rates = rl.get_rate_by_name(rate_names)
rc = pyna.RateCollection(rates=rates)

We can visualize the network and rates linking the nuclei

fig = rc.plot()

pynucastro can write out the python code needed to evaluate the reaction rates

pynet = pyna.PythonNetwork(rates=rates)
pynet.write_network("cno_integration_example.py")

Show code cell content

Hide code cell content

%cat cno_integration_example.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
jd = 1
jhe3 = 2
jhe4 = 3
jc12 = 4
jc13 = 5
jn13 = 6
jn14 = 7
jn15 = 8
jo15 = 9
nnuc = 10

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

A[jp] = 1
A[jd] = 2
A[jhe3] = 3
A[jhe4] = 4
A[jc12] = 12
A[jc13] = 13
A[jn13] = 13
A[jn14] = 14
A[jn15] = 15
A[jo15] = 15

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

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

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

mass[jp] = 0.0015040963047307696
mass[jd] = 0.0030058819195053215
mass[jhe3] = 0.004501176706825056
mass[jhe4] = 0.0059735574859708365
mass[jc12] = 0.017909017027273523
mass[jc13] = 0.01940644192976114
mass[jn13] = 0.01940999951603316
mass[jn14] = 0.020898440897976135
mass[jn15] = 0.022386433805845516
mass[jo15] = 0.02239084645968795

names = []
names.append("H1")
names.append("H2")
names.append("He3")
names.append("He4")
names.append("C12")
names.append("C13")
names.append("N13")
names.append("N14")
names.append("N15")
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_p__d__weak__bet_pos_", numba.float64),
    ("p_p__d__weak__electron_capture", numba.float64),
    ("p_d__He3", numba.float64),
    ("He3_He3__p_p_He4", numba.float64),
    ("p_C12__N13", numba.float64),
    ("p_C13__N14", numba.float64),
    ("N13__C13__weak__wc12", numba.float64),
    ("p_N14__O15", numba.float64),
    ("p_N15__He4_C12", numba.float64),
    ("O15__N15__weak__wc12", numba.float64),
])
class RateEval:
    def __init__(self):
        self.p_p__d__weak__bet_pos_ = np.nan
        self.p_p__d__weak__electron_capture = np.nan
        self.p_d__He3 = np.nan
        self.He3_He3__p_p_He4 = np.nan
        self.p_C12__N13 = np.nan
        self.p_C13__N14 = np.nan
        self.N13__C13__weak__wc12 = np.nan
        self.p_N14__O15 = np.nan
        self.p_N15__He4_C12 = 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_p__d__weak__bet_pos_(rate_eval, tf):
    # p + p --> d
    rate = 0.0

    # bet+w
    rate += np.exp(  -34.7863 + -3.51193*tf.T913i + 3.10086*tf.T913
                  + -0.198314*tf.T9 + 0.0126251*tf.T953 + -1.02517*tf.lnT9)

    rate_eval.p_p__d__weak__bet_pos_ = rate

@numba.njit()
def p_p__d__weak__electron_capture(rate_eval, tf):
    # p + p --> d
    rate = 0.0

    #   ecw
    rate += np.exp(  -43.6499 + -0.00246064*tf.T9i + -2.7507*tf.T913i + -0.424877*tf.T913
                  + 0.015987*tf.T9 + -0.000690875*tf.T953 + -0.207625*tf.lnT9)

    rate_eval.p_p__d__weak__electron_capture = rate

@numba.njit()
def p_d__He3(rate_eval, tf):
    # d + p --> He3
    rate = 0.0

    # de04 
    rate += np.exp(  8.93525 + -3.7208*tf.T913i + 0.198654*tf.T913
                  + 0.333333*tf.lnT9)
    # de04n
    rate += np.exp(  7.52898 + -3.7208*tf.T913i + 0.871782*tf.T913
                  + -0.666667*tf.lnT9)

    rate_eval.p_d__He3 = rate

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

    # nacrn
    rate += np.exp(  24.7788 + -12.277*tf.T913i + -0.103699*tf.T913
                  + -0.0649967*tf.T9 + 0.0168191*tf.T953 + -0.666667*tf.lnT9)

    rate_eval.He3_He3__p_p_He4 = rate

@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

    # 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)
    # nacrr
    rate += np.exp(  15.1825 + -13.5543*tf.T9i
                  + -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_N14__O15(rate_eval, tf):
    # N14 + p --> O15
    rate = 0.0

    # 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)
    # 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)

    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 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_p__d__weak__bet_pos_(rate_eval, tf)
    p_p__d__weak__electron_capture(rate_eval, tf)
    p_d__He3(rate_eval, tf)
    He3_He3__p_p_He4(rate_eval, tf)
    p_C12__N13(rate_eval, tf)
    p_C13__N14(rate_eval, tf)
    N13__C13__weak__wc12(rate_eval, tf)
    p_N14__O15(rate_eval, tf)
    p_N15__He4_C12(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, 1, 1)
        scor = screen_func(plasma_state, scn_fac)
        rate_eval.p_p__d__weak__bet_pos_ *= scor
        rate_eval.p_p__d__weak__electron_capture *= scor

        scn_fac = ScreenFactors(1, 1, 1, 2)
        scor = screen_func(plasma_state, scn_fac)
        rate_eval.p_d__He3 *= scor

        scn_fac = ScreenFactors(2, 3, 2, 3)
        scor = screen_func(plasma_state, scn_fac)
        rate_eval.He3_He3__p_p_He4 *= scor

        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, 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] = (
          + -2*5.00000000000000e-01*rho*Y[jp]**2*rate_eval.p_p__d__weak__bet_pos_  +
          + -2*5.00000000000000e-01*rho**2*ye(Y)*Y[jp]**2*rate_eval.p_p__d__weak__electron_capture  +
          -rho*Y[jp]*Y[jd]*rate_eval.p_d__He3  +
          + 2*5.00000000000000e-01*rho*Y[jhe3]**2*rate_eval.He3_He3__p_p_He4  +
          -rho*Y[jp]*Y[jc12]*rate_eval.p_C12__N13  +
          -rho*Y[jp]*Y[jc13]*rate_eval.p_C13__N14  +
          -rho*Y[jp]*Y[jn14]*rate_eval.p_N14__O15  +
          -rho*Y[jp]*Y[jn15]*rate_eval.p_N15__He4_C12
       )

    dYdt[jd] = (
          +5.00000000000000e-01*rho*Y[jp]**2*rate_eval.p_p__d__weak__bet_pos_  +
          +5.00000000000000e-01*rho**2*ye(Y)*Y[jp]**2*rate_eval.p_p__d__weak__electron_capture  +
          -rho*Y[jp]*Y[jd]*rate_eval.p_d__He3
       )

    dYdt[jhe3] = (
          +rho*Y[jp]*Y[jd]*rate_eval.p_d__He3  +
          + -2*5.00000000000000e-01*rho*Y[jhe3]**2*rate_eval.He3_He3__p_p_He4
       )

    dYdt[jhe4] = (
          +5.00000000000000e-01*rho*Y[jhe3]**2*rate_eval.He3_He3__p_p_He4  +
          +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] = (
          +rho*Y[jp]*Y[jc12]*rate_eval.p_C12__N13  +
          -Y[jn13]*rate_eval.N13__C13__weak__wc12
       )

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

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

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

    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_p__d__weak__bet_pos_(rate_eval, tf)
    p_p__d__weak__electron_capture(rate_eval, tf)
    p_d__He3(rate_eval, tf)
    He3_He3__p_p_He4(rate_eval, tf)
    p_C12__N13(rate_eval, tf)
    p_C13__N14(rate_eval, tf)
    N13__C13__weak__wc12(rate_eval, tf)
    p_N14__O15(rate_eval, tf)
    p_N15__He4_C12(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, 1, 1)
        scor = screen_func(plasma_state, scn_fac)
        rate_eval.p_p__d__weak__bet_pos_ *= scor
        rate_eval.p_p__d__weak__electron_capture *= scor

        scn_fac = ScreenFactors(1, 1, 1, 2)
        scor = screen_func(plasma_state, scn_fac)
        rate_eval.p_d__He3 *= scor

        scn_fac = ScreenFactors(2, 3, 2, 3)
        scor = screen_func(plasma_state, scn_fac)
        rate_eval.He3_He3__p_p_He4 *= scor

        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, 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] = (
       -2*5.00000000000000e-01*rho*2*Y[jp]*rate_eval.p_p__d__weak__bet_pos_
       -2*5.00000000000000e-01*rho**2*ye(Y)*2*Y[jp]*rate_eval.p_p__d__weak__electron_capture
       -rho*Y[jd]*rate_eval.p_d__He3
       -rho*Y[jc12]*rate_eval.p_C12__N13
       -rho*Y[jc13]*rate_eval.p_C13__N14
       -rho*Y[jn14]*rate_eval.p_N14__O15
       -rho*Y[jn15]*rate_eval.p_N15__He4_C12
       )

    jac[jp, jd] = (
       -rho*Y[jp]*rate_eval.p_d__He3
       )

    jac[jp, jhe3] = (
       +2*5.00000000000000e-01*rho*2*Y[jhe3]*rate_eval.He3_He3__p_p_He4
       )

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

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

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

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

    jac[jd, jp] = (
       -rho*Y[jd]*rate_eval.p_d__He3
       +5.00000000000000e-01*rho*2*Y[jp]*rate_eval.p_p__d__weak__bet_pos_
       +5.00000000000000e-01*rho**2*ye(Y)*2*Y[jp]*rate_eval.p_p__d__weak__electron_capture
       )

    jac[jd, jd] = (
       -rho*Y[jp]*rate_eval.p_d__He3
       )

    jac[jhe3, jp] = (
       +rho*Y[jd]*rate_eval.p_d__He3
       )

    jac[jhe3, jd] = (
       +rho*Y[jp]*rate_eval.p_d__He3
       )

    jac[jhe3, jhe3] = (
       -2*5.00000000000000e-01*rho*2*Y[jhe3]*rate_eval.He3_He3__p_p_He4
       )

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

    jac[jhe4, jhe3] = (
       +5.00000000000000e-01*rho*2*Y[jhe3]*rate_eval.He3_He3__p_p_He4
       )

    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[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
       )

    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[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[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

Now we can import the network that was just created

import cno_integration_example as cno

from scipy.integrate import solve_ivp

Now we’ll set the thermodynamic conditions. We’ll pick the conditions inside the core of the Sun. We initialize mass fractions and then convert to molar fractions, since that’s what the RHS uses

rho = 150
T = 1.5e7

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

Y0 = X / cno.A
Y0

array([0.7       , 0.        , 0.        , 0.07      , 0.00166667,
       0.        , 0.        , 0.        , 0.        , 0.        ])

We’ll use the BDF solver from SciPy–this works well with “stiff” reaction networks.

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-8, atol=1.e-10)

sol

  message: The solver successfully reached the end of the integration interval.
  success: True
   status: 0
        t: [ 0.000e+00  1.000e-04 ...  5.672e+19  1.000e+20]
        y: [[ 7.000e-01  7.000e-01 ... -1.270e-19 -3.698e-19]
            [ 0.000e+00  2.992e-22 ... -5.674e-12 -2.414e-11]
            ...
            [ 0.000e+00  1.888e-71 ...  3.859e-08  3.859e-08]
            [ 0.000e+00  3.937e-65 ... -3.857e-36 -1.123e-35]]
      sol: <scipy.integrate._ivp.common.OdeSolution object at 0x7f7b8c92d4d0>
 t_events: None
 y_events: None
     nfev: 1247
     njev: 30
      nlu: 164

Now we can plot the mass fractions.

fig = plt.figure()
ax = fig.add_subplot(111)

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

#ax.set_xlim(1.e8, 1.e13)
ax.set_ylim(1.e-8, 1.0)
ax.legend(fontsize="small")
ax.set_xlabel("time (yr)")
ax.set_ylabel("X")
ax.grid(ls=":")