Finding other important rates#
import pynucastro as pyna
import mesa_reader as mr
from pynucastro import mesa_utils
from pynucastro.screening import chugunov_2007
Read in the MESA model#
model = mr.MesaData("profile2.data")
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
Cell In[3], line 1
----> 1 model = mr.MesaData("profile2.data")
File /opt/hostedtoolcache/Python/3.14.5/x64/lib/python3.14/site-packages/mesa_reader/__init__.py:134, in MesaData.__init__(self, file_name, file_type)
132 self.header_data = None
133 self.header_names = None
--> 134 self.read_data()
File /opt/hostedtoolcache/Python/3.14.5/x64/lib/python3.14/site-packages/mesa_reader/__init__.py:183, in MesaData.read_data(self)
181 self.read_model_data()
182 elif self.file_type == "log":
--> 183 self.read_log_data()
184 else:
185 raise UnknownFileTypeError("Unknown file type {}".format(self.file_type))
File /opt/hostedtoolcache/Python/3.14.5/x64/lib/python3.14/site-packages/mesa_reader/__init__.py:199, in MesaData.read_log_data(self)
187 def read_log_data(self):
188 """Reads in or update data from the original log (.data or .log) file.
189
190 This re-reads the data from the originally-provided file name. Mostly
(...) 197 None
198 """
--> 199 self.bulk_data = np.genfromtxt(
200 self.file_name,
201 skip_header=MesaData.bulk_names_line - 1,
202 names=True,
203 ndmin=1, # Make sure a single entry is still a 1D array
204 dtype=None,
205 )
206 self.bulk_names = self.bulk_data.dtype.names
207 header_data = []
File /opt/hostedtoolcache/Python/3.14.5/x64/lib/python3.14/site-packages/numpy/lib/_npyio_impl.py:1978, in genfromtxt(fname, dtype, comments, delimiter, skip_header, skip_footer, converters, missing_values, filling_values, usecols, names, excludelist, deletechars, replace_space, autostrip, case_sensitive, defaultfmt, unpack, usemask, loose, invalid_raise, max_rows, encoding, ndmin, like)
1976 fname = os.fspath(fname)
1977 if isinstance(fname, str):
-> 1978 fid = np.lib._datasource.open(fname, 'rt', encoding=encoding)
1979 fid_ctx = contextlib.closing(fid)
1980 else:
File /opt/hostedtoolcache/Python/3.14.5/x64/lib/python3.14/site-packages/numpy/lib/_datasource.py:192, in open(path, mode, destpath, encoding, newline)
155 """
156 Open `path` with `mode` and return the file object.
157
(...) 188
189 """
191 ds = DataSource(destpath)
--> 192 return ds.open(path, mode, encoding=encoding, newline=newline)
File /opt/hostedtoolcache/Python/3.14.5/x64/lib/python3.14/site-packages/numpy/lib/_datasource.py:529, in DataSource.open(self, path, mode, encoding, newline)
526 return _file_openers[ext](found, mode=mode,
527 encoding=encoding, newline=newline)
528 else:
--> 529 raise FileNotFoundError(f"{path} not found.")
FileNotFoundError: profile2.data not found.
model.bulk_names
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[4], line 1
----> 1 model.bulk_names
NameError: name 'model' is not defined
nuclei = mesa_utils.get_nuclei(model)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[5], line 1
----> 1 nuclei = mesa_utils.get_nuclei(model)
NameError: name 'model' is not defined
nuclei
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[6], line 1
----> 1 nuclei
NameError: name 'nuclei' is not defined
mesa_zones = mesa_utils.get_all_data(model)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[7], line 1
----> 1 mesa_zones = mesa_utils.get_all_data(model)
NameError: name 'model' is not defined
Create the network using the same species as the MESA model#
Full network#
By default, network_helper will find all the rates connecting the nuclei, some of which were likely not used in the MESA simulation.
net_full = pyna.network_helper(nuclei)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[8], line 1
----> 1 net_full = pyna.network_helper(nuclei)
NameError: name 'nuclei' is not defined
net_full.summary()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[9], line 1
----> 1 net_full.summary()
NameError: name 'net_full' is not defined
fig = net_full.plot()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[10], line 1
----> 1 fig = net_full.plot()
NameError: name 'net_full' is not defined
reduced network#
separately create a network that has just the rates that we used in MESA.
rl = pyna.ReacLibLibrary()
tl = pyna.TabularLibrary()
all_lib = rl + tl
The MESA run uses a rate that approximates the \({}^{16}\mathrm{O} + {}^{16}\mathrm{O}\) burning through \({}^{31}\mathrm{P}\). We’ll create that rate first.
from pynucastro.rates.aprox_family_rates import make_CO_approx_rates
o_rates = make_CO_approx_rates(rl.get_rates(), "O")
o_lib = pyna.Library(rates=o_rates)
r1616 = o_lib.get_rate_by_name("o16(o16,a)si28")
type(r1616)
pynucastro.rates.approximate_rates.ApproximateRate
We see that r1616 is an ApproximateRate that will carry multiple rates to do the approximation.
Now get all the weak rates we used in the MESA simulation
rate_names = ["ne20(,)f20", "f20(,)ne20",
"f20(,)o20", "o20(,)f20",
"mg24(,)na24", "na24(,)mg24",
"na24(,)ne24", "ne24(,)na24",
"na23(,)ne23", "ne23(,)na23",
"mg25(,)na25", "na25(,)mg25",
"na25(,)ne25", "ne25(,)na25"]
rates = all_lib.get_rate_by_name(rate_names)
and put it all together, eliminating any duplicates (e.g. ReacLib and Suzuki both provide the link—we prefer the Suzuki)
lib = pyna.Library(rates=[r1616] + rates)
lib.eliminate_duplicates()
net = pyna.PythonNetwork(libraries=lib)
net.summary()
Network summary
---------------
explicitly carried nuclei: 14
approximated-out nuclei: 1
inert nuclei (included in carried): 0
NSE compatible? False
total number of rates: 22
rates explicitly connecting nuclei: 15
hidden rates: 7
reaclib rates: 7
starlib rates: 0
temperature tabular rates: 0
weak tabular rates: 14
approximate rates: 1
derived rates: 0
modified rates: 0
custom rates: 0
fig = net.plot()
Check nuclei for other branches#
full_library = pyna.full_library()
_, missing_rates = net.validate(full_library, return_dict=True)
validation: Si28 produced in O16 + O16 ⟶ Si28 + He4 never consumed.
validation: He4 produced in O16 + O16 ⟶ Si28 + He4 never consumed.
len(missing_rates)
82
This finds 82 other potential rates
lib = pyna.Library(rates=list(missing_rates.keys()))
lib
O16 + He4 ⟶ Ne20 + 𝛾 [Q = 4.73 MeV] (O16 + He4 --> Ne20 <reaclib_co10>)
O16 + O16 ⟶ p + P31 [Q = 7.68 MeV] (O16 + O16 --> p + P31 <reaclib_cf88>)
O16 + O16 ⟶ n + S31 [Q = 1.50 MeV] (O16 + O16 --> n + S31 <reaclib_cf88>)
O20 + He4 ⟶ n + Ne23 [Q = 3.30 MeV] (O20 + He4 --> n + Ne23 <reaclib_rpsm>)
O20 + He4 ⟶ Ne24 + 𝛾 [Q = 12.17 MeV] (O20 + He4 --> Ne24 <reaclib_rpsm>)
F20 + He4 ⟶ p + Ne23 [Q = 0.27 MeV] (F20 + He4 --> p + Ne23 <reaclib_rpsm>)
F20 + He4 ⟶ n + Na23 [Q = 3.87 MeV] (F20 + He4 --> n + Na23 <reaclib_rpsm>)
F20 + He4 ⟶ Na24 + 𝛾 [Q = 10.82 MeV] (F20 + He4 --> Na24 <reaclib_rpsm>)
Ne20 + He4 ⟶ Mg24 + 𝛾 [Q = 9.32 MeV] (Ne20 + He4 --> Mg24 <reaclib_il10>)
Ne23 + He4 ⟶ n + Mg26 [Q = 5.42 MeV] (Ne23 + He4 --> n + Mg26 <reaclib_rath>)
Ne23 + He4 ⟶ Mg27 + 𝛾 [Q = 11.86 MeV] (Ne23 + He4 --> Mg27 <reaclib_rath>)
Ne24 + He4 ⟶ n + Mg27 [Q = 2.99 MeV] (Ne24 + He4 --> n + Mg27 <reaclib_rath>)
Ne24 + He4 ⟶ Mg28 + 𝛾 [Q = 11.50 MeV] (Ne24 + He4 --> Mg28 <reaclib_rath>)
Ne25 + He4 ⟶ n + Mg28 [Q = 7.31 MeV] (Ne25 + He4 --> n + Mg28 <reaclib_rath>)
Ne25 + He4 ⟶ Mg29 + 𝛾 [Q = 11.03 MeV] (Ne25 + He4 --> Mg29 <reaclib_rath>)
Na23 + He4 ⟶ p + Mg26 [Q = 1.82 MeV] (Na23 + He4 --> p + Mg26 <reaclib_ths8>)
Na23 + He4 ⟶ Al27 + 𝛾 [Q = 10.09 MeV] (Na23 + He4 --> Al27 <reaclib_ths8>)
Na24 + He4 ⟶ p + Mg27 [Q = 1.30 MeV] (Na24 + He4 --> p + Mg27 <reaclib_rath>)
Na24 + He4 ⟶ n + Al27 [Q = 3.13 MeV] (Na24 + He4 --> n + Al27 <reaclib_rath>)
Na24 + He4 ⟶ Al28 + 𝛾 [Q = 10.86 MeV] (Na24 + He4 --> Al28 <reaclib_rath>)
Na25 + He4 ⟶ p + Mg28 [Q = 0.80 MeV] (Na25 + He4 --> p + Mg28 <reaclib_rath>)
Na25 + He4 ⟶ n + Al28 [Q = 1.85 MeV] (Na25 + He4 --> n + Al28 <reaclib_rath>)
Na25 + He4 ⟶ Al29 + 𝛾 [Q = 11.28 MeV] (Na25 + He4 --> Al29 <reaclib_rath>)
Mg24 + He4 ⟶ Si28 + 𝛾 [Q = 9.98 MeV] (Mg24 + He4 --> Si28 <reaclib_st08>)
Mg25 + He4 ⟶ n + Si28 [Q = 2.65 MeV] (Mg25 + He4 --> n + Si28 <reaclib_nacr>)
Mg25 + He4 ⟶ Si29 + 𝛾 [Q = 11.13 MeV] (Mg25 + He4 --> Si29 <reaclib_cf88>)
Si28 + He4 ⟶ S32 + 𝛾 [Q = 6.95 MeV] (Si28 + He4 --> S32 <reaclib_ths8>)
He4 + He4 ⟶ p + Li7 [Q = -17.35 MeV] (He4 + He4 --> p + Li7 <reaclib_de04>)
He4 + He4 ⟶ n + Be7 [Q = -18.99 MeV] (He4 + He4 --> n + Be7 <reaclib_wag>)
O16 + He4 ⟶ p + F19 [Q = -8.11 MeV] (O16 + He4 --> p + F19 <reaclib_nacr>)
O16 + He4 ⟶ n + Ne19 [Q = -12.13 MeV] (O16 + He4 --> n + Ne19 <reaclib_ths8>)
O20 + He4 ⟶ p + F23 [Q = -4.40 MeV] (O20 + He4 --> p + F23 <reaclib_rpsm>)
O20 ⟶ n + O19 [Q = -7.61 MeV] (O20 --> n + O19 <reaclib_wies>)
F20 ⟶ p + O19 [Q = -10.64 MeV] (F20 --> p + O19 <reaclib_rpsm>)
F20 ⟶ n + F19 [Q = -6.60 MeV] (F20 --> n + F19 <reaclib_ks03>)
Ne20 + He4 ⟶ C12 + C12 [Q = -4.62 MeV] (Ne20 + He4 --> C12 + C12 <reaclib_cf88>)
Ne20 + He4 ⟶ p + Na23 [Q = -2.38 MeV] (Ne20 + He4 --> p + Na23 <reaclib_il10>)
Ne20 + He4 ⟶ n + Mg23 [Q = -7.21 MeV] (Ne20 + He4 --> n + Mg23 <reaclib_ths8>)
Ne20 ⟶ He4 + O16 [Q = -4.73 MeV] (Ne20 --> He4 + O16 <reaclib_co10>)
Ne20 ⟶ p + F19 [Q = -12.84 MeV] (Ne20 --> p + F19 <reaclib_nacr>)
Ne20 ⟶ n + Ne19 [Q = -16.86 MeV] (Ne20 --> n + Ne19 <reaclib_ths8>)
Ne20 ⟶ Na20 + e⁻ + 𝜈 [Q = -13.89 MeV] (Ne20 --> Na20 <weaktab_oda>)
Ne23 + He4 ⟶ p + Na26 [Q = -3.12 MeV] (Ne23 + He4 --> p + Na26 <reaclib_rath>)
Ne23 ⟶ He4 + O19 [Q = -10.91 MeV] (Ne23 --> He4 + O19 <reaclib_rpsm>)
Ne23 ⟶ p + F22 [Q = -15.24 MeV] (Ne23 --> p + F22 <reaclib_rpsm>)
Ne23 + e⁻ ⟶ F23 + 𝜈 [Q = -8.44 MeV] (Ne23 --> F23 <weaktab_oda>)
Ne23 ⟶ n + Ne22 [Q = -5.20 MeV] (Ne23 --> n + Ne22 <reaclib_ks03>)
Ne24 + He4 ⟶ p + Na27 [Q = -5.23 MeV] (Ne24 + He4 --> p + Na27 <reaclib_rath>)
Ne24 ⟶ He4 + O20 [Q = -12.17 MeV] (Ne24 --> He4 + O20 <reaclib_rpsm>)
Ne24 ⟶ p + F23 [Q = -16.57 MeV] (Ne24 --> p + F23 <reaclib_rpsm>)
Ne24 ⟶ n + Ne23 [Q = -8.87 MeV] (Ne24 --> n + Ne23 <reaclib_wies>)
Ne25 + He4 ⟶ p + Na28 [Q = -5.89 MeV] (Ne25 + He4 --> p + Na28 <reaclib_rath>)
Ne25 ⟶ He4 + O21 [Q = -12.54 MeV] (Ne25 --> He4 + O21 <reaclib_rpsm>)
Ne25 ⟶ p + F24 [Q = -16.89 MeV] (Ne25 --> p + F24 <reaclib_rpsm>)
Ne25 ⟶ n + Ne24 [Q = -4.28 MeV] (Ne25 --> n + Ne24 <reaclib_wies>)
Na23 + He4 ⟶ n + Al26 [Q = -2.97 MeV] (Na23 + He4 --> n + Al26 <reaclib_ol11>)
Na23 ⟶ He4 + F19 [Q = -10.47 MeV] (Na23 --> He4 + F19 <reaclib_rpsm>)
Na23 ⟶ p + Ne22 [Q = -8.79 MeV] (Na23 --> p + Ne22 <reaclib_ke17>)
Na23 ⟶ n + Na22 [Q = -12.42 MeV] (Na23 --> n + Na22 <reaclib_ths8>)
Na23 ⟶ Mg23 + e⁻ + 𝜈 [Q = -4.06 MeV] (Na23 --> Mg23 <weaktab_oda>)
Na24 ⟶ He4 + F20 [Q = -10.82 MeV] (Na24 --> He4 + F20 <reaclib_rpsm>)
Na24 ⟶ p + Ne23 [Q = -10.55 MeV] (Na24 --> p + Ne23 <reaclib_rath>)
Na24 ⟶ n + Na23 [Q = -6.96 MeV] (Na24 --> n + Na23 <reaclib_ks03>)
Na25 ⟶ He4 + F21 [Q = -11.73 MeV] (Na25 --> He4 + F21 <reaclib_rpsm>)
Na25 ⟶ p + Ne24 [Q = -10.70 MeV] (Na25 --> p + Ne24 <reaclib_rath>)
Na25 ⟶ n + Na24 [Q = -9.01 MeV] (Na25 --> n + Na24 <reaclib_rath>)
Mg24 + He4 ⟶ C12 + O16 [Q = -6.77 MeV] (Mg24 + He4 --> C12 + O16 <reaclib_cf88>)
Mg24 + He4 ⟶ p + Al27 [Q = -1.60 MeV] (Mg24 + He4 --> p + Al27 <reaclib_il10>)
Mg24 + He4 ⟶ n + Si27 [Q = -7.20 MeV] (Mg24 + He4 --> n + Si27 <reaclib_ths8>)
Mg24 ⟶ He4 + Ne20 [Q = -9.32 MeV] (Mg24 --> He4 + Ne20 <reaclib_il10>)
Mg24 ⟶ p + Na23 [Q = -11.69 MeV] (Mg24 --> p + Na23 <reaclib_il10>)
Mg24 ⟶ n + Mg23 [Q = -16.53 MeV] (Mg24 --> n + Mg23 <reaclib_ths8>)
Mg24 ⟶ Al24 + e⁻ + 𝜈 [Q = -13.88 MeV] (Mg24 --> Al24 <weaktab_oda>)
Mg25 + He4 ⟶ p + Al28 [Q = -1.21 MeV] (Mg25 + He4 --> p + Al28 <reaclib_cf88>)
Mg25 ⟶ He4 + Ne21 [Q = -9.88 MeV] (Mg25 --> He4 + Ne21 <reaclib_cf88>)
Mg25 ⟶ p + Na24 [Q = -12.06 MeV] (Mg25 --> p + Na24 <reaclib_rath>)
Mg25 ⟶ n + Mg24 [Q = -7.33 MeV] (Mg25 --> n + Mg24 <reaclib_ks03>)
Mg25 ⟶ Al25 + e⁻ + 𝜈 [Q = -4.28 MeV] (Mg25 --> Al25 <weaktab_oda>)
Si28 + He4 ⟶ O16 + O16 [Q = -9.59 MeV] (Si28 + He4 --> O16 + O16 <reaclib_cf88>)
Si28 + He4 ⟶ C12 + Ne20 [Q = -12.02 MeV] (Si28 + He4 --> C12 + Ne20 <reaclib_rolf>)
Si28 + He4 ⟶ p + P31 [Q = -1.92 MeV] (Si28 + He4 --> p + P31 <reaclib_il10>)
Si28 + He4 ⟶ n + S31 [Q = -8.09 MeV] (Si28 + He4 --> n + S31 <reaclib_ths8>)
Determine if any missing rates are important#
We’ll pick some zones from the MESA model to consider
len(mesa_zones)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[24], line 1
----> 1 len(mesa_zones)
NameError: name 'mesa_zones' is not defined
import numpy as np
zones = np.linspace(385, 785, 11, dtype=int)
zones
array([385, 425, 465, 505, 545, 585, 625, 665, 705, 745, 785])
We only consider a rate important if it is within this threshold for the fastest capture on the same nucleus. Since some existing rates may be very small, we also put a floor, ydot_cutoff on the # of reactions / sec we care about
thresh = 0.2
ydot_cutoff = 1.e-15
possibly_important = []
for izone in zones:
rho, T, comp = mesa_zones[izone]
for nuc in net.unique_nuclei:
# evaluate all existing rates working with this nucleus
r_current = {r : r.eval_full_rate(rho, T, comp, screen_func=chugunov_2007)
for r in net.get_rates() if nuc == max(r.reactants)}
# evaluate all missing rates working with this nucleus
r_missing = {r : r.eval_full_rate(rho, T, comp, screen_func=chugunov_2007)
for r in missing_rates if nuc == max(r.reactants)}
# compute the fastest (abs value) current rate consuming this nucleus
if not (r_current and r_missing):
continue
r_fastest = max(r_current, key=r_current.get)
lambda_fastest = r_current[r_fastest]
# now check if any missing rate is within the threshold
possibly_important += [r for r, v in r_missing.items()
if v > max(ydot_cutoff, thresh * lambda_fastest)]
possibly_important = set(possibly_important)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[28], line 4
1 possibly_important = []
2
3 for izone in zones:
----> 4 rho, T, comp = mesa_zones[izone]
5 for nuc in net.unique_nuclei:
6 # evaluate all existing rates working with this nucleus
7 r_current = {r : r.eval_full_rate(rho, T, comp, screen_func=chugunov_2007)
NameError: name 'mesa_zones' is not defined
comp
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[29], line 1
----> 1 comp
NameError: name 'comp' is not defined
possibly_important
[]
which of these involve only the nuclei we are already carrying? (also allow protons)
okay_nuc = net.unique_nuclei + [pyna.Nucleus("p")]
filtered_rates = []
for r in possibly_important:
okay = True
for nuc in r.reactants + r.products:
if nuc not in okay_nuc:
okay = False
break
if okay:
filtered_rates.append(r)
filtered_rates
[]
Network with all missing rates#
rates = net.get_rates() + list(possibly_important)
big_net = pyna.PythonNetwork(rates=rates)
big_net.summary()
Network summary
---------------
explicitly carried nuclei: 14
approximated-out nuclei: 1
inert nuclei (included in carried): 0
NSE compatible? False
total number of rates: 22
rates explicitly connecting nuclei: 15
hidden rates: 7
reaclib rates: 7
starlib rates: 0
temperature tabular rates: 0
weak tabular rates: 14
approximate rates: 1
derived rates: 0
modified rates: 0
custom rates: 0
fig = big_net.plot()
what if we restrict it to just the rates that use the same nuclei?
rates = net.get_rates() + list(filtered_rates)
new_net = pyna.PythonNetwork(rates=rates)
new_net.summary()
Network summary
---------------
explicitly carried nuclei: 14
approximated-out nuclei: 1
inert nuclei (included in carried): 0
NSE compatible? False
total number of rates: 22
rates explicitly connecting nuclei: 15
hidden rates: 7
reaclib rates: 7
starlib rates: 0
temperature tabular rates: 0
weak tabular rates: 14
approximate rates: 1
derived rates: 0
modified rates: 0
custom rates: 0
fig = new_net.plot()