Network Comparison: small vs small_parthenope

This notebook compares the primordial abundances and time evolution predicted by the small network and the small_parthenope network.

import importlib.util, subprocess, sys

# In Colab (or any environment without a local primat checkout), pip-install
# the released package; local dev runs already have it importable (editable
# install / repo on sys.path via the cell below), so this is a no-op there.
if importlib.util.find_spec("primat") is None:
    subprocess.run([sys.executable, "-m", "pip", "install", "-q", "primat"], check=True)

import sys
import os
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

# Ensure the repo root is in path
sys.path.append(os.getcwd())

from primat.backend import run_bbn
from primat.evolution import Y_interpolator, T_gamma_interpolator
from primat.constants import CONST

# Uses primat.backend.run_bbn() only (works on either backend); the curve data
# below comes from primat.evolution's backend-agnostic interpolators built on
# each run's EvolutionResult (output_time_evolution=True).

# Configuration
params = {
    "output_time_evolution": True,
    "output_file": None,
    "verbose": False
}

# Solve for both networks
print("Solving small network...")
res_small = run_bbn({**params, "network": "small"})

print("Solving small_parthenope network...")
res_parthenope = run_bbn({**params, "network": "small_parthenope"})

print("Done.")

Solving small network...
Solving small_parthenope network...

Done.

[primat]  HT.  MT.  LT.  done.
[primat]  HT.  MT.  LT.  done.

# Compare final abundances using Y_final
nuclides = [n for n in res_small['Y_final'] if n in res_parthenope['Y_final']]
data = {
    "small": [res_small['Y_final'][n] for n in nuclides],
    "small_parthenope": [res_parthenope['Y_final'][n] for n in nuclides],
}
df = pd.DataFrame(data, index=nuclides)
df["rel_diff"] = (df["small_parthenope"] - df["small"]) / df["small"]

print("Comparison of final abundances:")
print(df)

Comparison of final abundances:
            small  small_parthenope  rel_diff
n    3.997469e-16      2.727988e-16 -0.317571
p    7.529428e-01      7.530503e-01  0.000143
H2   1.834100e-05      1.882624e-05  0.026457
H3   5.852036e-08      6.045781e-08  0.033107
He3  7.771620e-06      7.731663e-06 -0.005141
He4  6.174925e-02      6.172216e-02 -0.000439
Li7  2.181379e-11      2.075999e-11 -0.048309
Be7  3.966364e-10      3.416178e-10 -0.138713

# Plot relative difference of nuclide time evolution
plt.figure(figsize=(10, 6))

ev_small = res_small['evolution']
ev_parthenope = res_parthenope['evolution']

# Use each run's own time grid and the backend-agnostic T_gamma interpolator.
t_small = ev_small.t
t_parthenope = ev_parthenope.t
# T_gamma_interpolator returns MeV; convert to Kelvin for the x-axis label.

T_small = T_gamma_interpolator(ev_small)(t_small) * CONST.MeV_to_Kelvin
T_parthenope = T_gamma_interpolator(ev_parthenope)(t_parthenope) * CONST.MeV_to_Kelvin

for n in nuclides:
    Y_s = Y_interpolator(ev_small, n)(t_small)
    Y_p = Y_interpolator(ev_parthenope, n)(t_parthenope)

    # Interpolate parthenope evolution onto small grid (T decreases with time)
    Y_p_interp = np.interp(T_small, T_parthenope[::-1], Y_p[::-1])
    rel_diff = (Y_p_interp - Y_s) / Y_s

    mask = Y_s > 1e-25
    plt.plot(T_small[mask], rel_diff[mask], label=n)

plt.xscale("log")
plt.gca().invert_xaxis()
plt.xlabel("Photon Temperature (K)")
plt.ylabel("Relative Difference (Parthenope - Small) / Small")
plt.legend()
plt.grid(True)
plt.show()

/var/folders/x_/t8nzcr8d7_91h7fdrxdr8b180000gp/T/ipykernel_96394/1532134381.py:21: RuntimeWarning: invalid value encountered in divide
  rel_diff = (Y_p_interp - Y_s) / Y_s
/var/folders/x_/t8nzcr8d7_91h7fdrxdr8b180000gp/T/ipykernel_96394/1532134381.py:32: UserWarning: FigureCanvasAgg is non-interactive, and thus cannot be shown
  plt.show()