<jupyter_start><jupyter_text>Example of integrating a networkThis notebook illustrates how to create a python network and integrateit with the scipy library.<jupyter_code>import pynucastro as pyna<jupyter_output><empty_output><jupyter_text>We'll start again with the basic CNO network explored earlier.<jupyter_code>files = ["c12-pg-n13-ls09", 
         "c13-pg-n14-nacr",
         "n13--c13-wc12",
         "n13-pg-o14-lg06",
         "n14-pg-o15-im05",
         "n15-pa-c12-nacr",
         "o14--n14-wc12",
         "o15--n15-wc12"]<jupyter_output><empty_output><jupyter_text>A `PythonNetwork` is based on a `RateCollection` but has methods to write the RHS of the system of ODEs.<jupyter_code>pynet = pyna.PythonNetwork(files)<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>r = pynet.rates[1]
print(r)
print(pynet.ydot_string(r))<jupyter_output>rho*Y[jp]*Y[jc13]*lambda_p_c13__n14<jupyter_text>and the code needed to evaluate that rate (the T-dependent part) as:<jupyter_code>print(pynet.function_string(r))<jupyter_output>@numba.njit()
def p_c13__n14(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)
    
    return rate<jupyter_text>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.<jupyter_code>pynet.write_network("cno_test_integrate.py")
%cat cno_test_integrate.py<jupyter_output>import numpy as np
from pynucastro.rates import Tfactors
import numba

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

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

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

@numba.njit()
def p_c12__n13(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.ln[...]<jupyter_text>We can now import the network that was just created and integrate it using the SciPy ODE solvers<jupyter_code>import cno_test_integrate as cno<jupyter_output><empty_output><jupyter_text>Integrating the networkWe can use the stiff ODE integration solvers that are part of Scipy to integrate this system now<jupyter_code>from scipy.integrate import solve_ivp
import numpy as np<jupyter_output><empty_output><jupyter_text>Initialize the thermodynamic conditions and initial composition.  We express the composition as molar fractions, `Y0`.<jupyter_code>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<jupyter_output><empty_output><jupyter_text>Now we integrate.  We use the `BDF` method, since reaction networks are in general stiff<jupyter_code>tmax = 1.e20

sol = solve_ivp(cno.rhs, [0, tmax], Y0, method="BDF",
                dense_output=True, args=(rho, T), rtol=1.e-6, atol=1.e-6)<jupyter_output><empty_output><jupyter_text>A network plot<jupyter_code>import matplotlib.pyplot as plt
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")

fig.set_size_inches((10, 8))<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup