NSE vs. Direct Integration

import pynucastro as pyna

To further demonstrate the validity of the NSE equation, we will compare the equilibrium mass fractions computed via direct integration and the mass fractions obtained from NSE equations.

We will create the PythonNetwork needed for direct integration using the same library as our previous example.

First we recreate our library, but we make sure that we use reverse rates calculated via detailed balance. In pynucastro, this is the DerivedRate.

rl = pyna.ReacLibLibrary()

all_nuclei = ["p", "he4",
              "c12", "n13",
              "o16", "ne20", "na23",
              "mg24", "al27", "si28",
              "p31", "s32", "cl35",
              "ar36", "k39", "ca40",
              "sc43", "ti44", "v47",
              "cr48", "mn51",
              "fe52","co55","ni56"]

lib = rl.linking_nuclei(all_nuclei)

Here we search through the existing ReacLib reverse rates and remove them. So that Library only contains forward rates. We use the backward method for this.

reverse_rates = lib.backward().get_rates()
for r in reverse_rates:
    lib.remove_rate(r)

Now create the corresponding DerivedRate using the forward rate and add them to Library.

for r in lib.get_rates():
    d = pyna.DerivedRate(rate=r, compute_Q=True, use_pf=True)
    lib.add_rate(d)

We can also use ModifiedRate where we modify the reactants and/or the products of a reaction rate but still use the original rate. This is compatible with NSE as long as a corresponding DerivedRate is added.

Warning

Currently stoichiometry does NOT work with NSE.

other_rates = [("c12(c12,n)mg23", "mg24"),
               ("o16(o16,n)s31", "s32"),
               ("o16(c12,n)si27", "si28")]
for r, mp in other_rates:
    _r = rl.get_rate_by_name(r)
    forward_rate = pyna.ModifiedRate(_r, new_products=[mp])
    derived_rate = pyna.DerivedRate(rate=forward_rate, compute_Q=True, use_pf=True)
    lib.add_rates([forward_rate, derived_rate])

Now create the PythonNetwork using the modified Library object.

net = pyna.PythonNetwork(libraries=lib)

We can also do some approximations, such as the \((\alpha,p)(p,\gamma)\) approximation. This is fully compatible with NSE.

net.make_ap_pg_approx(intermediate_nuclei=["cl35", "k39", "sc43", "v47", "mn51", "co55"])
net.remove_nuclei(["cl35", "k39", "sc43", "v47", "mn51", "co55"])

Finally, we write out python module that contains essential functions needed for integration, such as the rhs and jacobian.

net.write_network(outfile="network.py")

/opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages/pynucastro/rates/derived_rate.py:126: UserWarning: C12 partition function is not supported by tables: set pf = 1.0 by default
  warnings.warn(UserWarning(f'{nuc} partition function is not supported by tables: set pf = 1.0 by default'))
/opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages/pynucastro/rates/derived_rate.py:126: UserWarning: N13 partition function is not supported by tables: set pf = 1.0 by default
  warnings.warn(UserWarning(f'{nuc} partition function is not supported by tables: set pf = 1.0 by default'))

Important

In order for the integrated equilibrium result to be fully compatible with NSE equations, watch out for a number of things:

1. Make sure that all forward rates in the network has a corresponding DerivedRate. Otherwise there is no detailed balance.

2. All species are connected.

3. NSE only guarantees equilibrium of strong reactions. So if there are weak rates in the network, \(Y_e\) will evolve in time.

Direct integration

Now we want to integrate. Since we will compare to NSE, we want to use the same screening that NSE uses–potekhin_1998

from pynucastro.screening import potekhin_1998
import network

and we’ll integrate using solve_ivp from scipy.integrate

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

Set up an initial condition and let it evolve until we reach equilibrium.

rho = 1e7
T = 6e9

X0 = np.ones(network.nnuc)
X0 /= np.sum(X0)
Y0 = X0/network.A

ye = network.ye(Y0)

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

Let’s plot the evolution to see how mass fraction evolve with time. Equilibrium is indicated by flat lines.

nucs = list(map(pyna.Nucleus, network.names))
nuc_labels = [f"${n.pretty}$" for n in nucs]

fig, ax = plt.subplots()
for i in range(network.nnuc):
    ax.loglog(sol.t, sol.y[i,:]*network.A[i], label=rf"X(${nucs[i].pretty}$)")

ax.set_xlim(1e-16, tmax)
ax.set_xlabel(f"time [s]", fontsize=12)
ax.set_ylabel(f"mass fractions", fontsize=12)
ax.legend(loc="best", fontsize="small", ncol=3)
ax.grid(ls=":")
fig.tight_layout()

NSE calculation

We can create an NSENetwork by simply grabbing the rates from our current network and using them to construct the NSE network

nse = pyna.NSENetwork(rates=net.get_rates())
nse.summary()

Network summary
---------------
  explicitly carried nuclei: 18
  approximated-out nuclei: 6
  inert nuclei (included in carried): 0

  total number of rates: 105

  rates explicitly connecting nuclei: 66
  hidden rates: 39

  reaclib rates: 45
  tabular rates: 0
  approximate rates: 12
  derived rates: 45
  modified rates: 3
  custom rates: 0

Now let’s get the directly-computed NSE state

nse_comp = nse.get_comp_nse(rho, T, ye, init_guess=(-3.5, -14.0), use_coulomb_corr=True)

Finally, we can make some plots to see the direct comparison between the two results.

X_net = sol.y[:,-1]*network.A
X_nse = np.array(list(nse_comp.X.values()))

x = np.arange(len(network.names))
width = 0.35

fig, ax = plt.subplots()
ax.bar(x - width/2.0, X_nse, width, label = 'NSE')
ax.bar(x + width/2.0, X_net, width, label = 'Integration')
ax.set_xlabel("nucleus", fontsize=14)
ax.set_xticks(x, labels=nuc_labels, rotation=90, fontsize=14)
ax.set_ylabel("mass fractions", fontsize=14)
ax.legend(fontsize=14)
ax.set_yscale("log")
ax.set_ylim(ymax = 1)
ax.legend()
fig.tight_layout()

Now a plot of absolute and relative difference shows that the mass fraction obtained from direct integration matches with NSE equations to machine precision. The differences between the two are largely set by the tolerances of both the NSE solver and the ODE integrator.

diff = np.abs(X_net - X_nse) 
diff_rel = diff / X_net

fig, ax = plt.subplots()
ax.scatter(x, diff, label='Absolute Error', marker='o', color='r')
ax.scatter(x, diff_rel, label='Relative Error', marker ='x', color='b')
ax.set_xticks(x, labels=nucs, rotation=90, fontsize=14)
ax.set_xlabel("nucleus", fontsize=14)
ax.set_yscale("log")
ax.legend(fontsize=14)
ax.set_ylim(ymax = 1)
ax.legend()
fig.tight_layout()

Visualizing equilibrium

At equilibrium, the forward and reverse rates linking nuclei should come into balance. We can visualize this graphically using the network plot function, and plotting the difference between the forward and reverse rates (use_net_rate=True). We will also normalize by the total activity of the rate (forward + reverse) by doing normalize_net_rate=True.

First we’ll convert the final composition from the integration into a pynucastro Composition, keeping in mind that the network evolved molar fractions (\(Y = X / A\)).

comp = pyna.Composition(net.unique_nuclei)
comp.set_array(sol.y[:, -1] * network.A)

fig = net.plot(rho, T, comp,
               screen_func=potekhin_1998,
               rotated=True, size=(1200, 500),
               node_size=650, node_font_size="9",
               use_net_rate=True, normalize_net_rate=True,
               color_nodes_by_abundance=True,
               hide_xalpha=True)

/opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages/pynucastro/rates/derived_rate.py:126: UserWarning: C12 partition function is not supported by tables: set pf = 1.0 by default
  warnings.warn(UserWarning(f'{nuc} partition function is not supported by tables: set pf = 1.0 by default'))
/opt/hostedtoolcache/Python/3.11.14/x64/lib/python3.11/site-packages/pynucastro/rates/derived_rate.py:126: UserWarning: N13 partition function is not supported by tables: set pf = 1.0 by default
  warnings.warn(UserWarning(f'{nuc} partition function is not supported by tables: set pf = 1.0 by default'))

Here we see that the normalized net rates are all < \(10^{-12}\), so we are effectively seeing the integration tolerance, and are very close to machine precision.