Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

add a method to generate an NSE table #676

Merged
merged 27 commits into from
Nov 6, 2023
Merged

add a method to generate an NSE table #676

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
27 commits
Select commit Hold shift + click to select a range
ace086b
separate NSE stuff into its own NSENetwork
zingale Oct 29, 2023
329a7d7
add missing file
zingale Oct 29, 2023
b517733
some missing
zingale Oct 29, 2023
9ad0199
update docs
zingale Oct 29, 2023
bc62055
fix tests
zingale Oct 29, 2023
e898c36
fix docs
zingale Oct 29, 2023
4bf28cc
update source API
zingale Oct 29, 2023
bb362b4
update notebook
zingale Oct 29, 2023
8961bb8
try again
zingale Oct 29, 2023
f219271
add a method to generate an NSE table
zingale Oct 29, 2023
74a5da9
this seems to work now
zingale Oct 29, 2023
7cfde65
flake8
zingale Oct 29, 2023
25c8972
fix verbose
zingale Oct 29, 2023
0131217
Merge branch 'main' into nse_table
zingale Oct 31, 2023
da6d6d4
fix missing newline
zingale Oct 31, 2023
d308dc3
update Microphysics benchmarks due to sneut5 changes
zingale Oct 31, 2023
568dc80
a little doc reorganization
zingale Oct 31, 2023
9c178ab
Merge branch 'some_docs_reorg' into nse_table
zingale Oct 31, 2023
a175cf5
Merge branch 'update_microphysics_bench' into nse_table
zingale Oct 31, 2023
2ab9cca
add a unit test
zingale Oct 31, 2023
4b6fc19
fix flake8
zingale Oct 31, 2023
b4975df
Merge branch 'main' into nse_table
zingale Oct 31, 2023
541a735
Merge branch 'main' into nse_table
zingale Nov 2, 2023
37e0497
add d<B/A>/dt
zingale Nov 2, 2023
777aa8b
fix
zingale Nov 2, 2023
04fb95c
Merge branch 'main' into nse_table
zingale Nov 4, 2023
bd60131
update
zingale Nov 4, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
137 changes: 134 additions & 3 deletions pynucastro/networks/nse_network.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,9 +4,69 @@
from scipy import constants
from scipy.optimize import fsolve

from pynucastro._version import version
from pynucastro.networks.rate_collection import Composition, RateCollection
from pynucastro.nucdata import Nucleus
from pynucastro.screening import NseState
from pynucastro.rates import TabularRate
from pynucastro.screening import NseState, potekhin_1998


class NSETableEntry:
def __init__(self, rho, T, Ye, *,
comp=None, ydots=None, enu=None,
comp_reduction_func=None):
"""a simple container to hold a single entry in the NSE table.

Here, comp_reduction_func(comp) is a function that converts
the NSE composition into a smaller set of nuclei. It takes a
Composition object and returns a dictionary with the nucleus
name (like "Ni56") as the key and the corresponding mass fraction
as the value. It should be ordered in the way you want the nuclei
output into the NSE table file.

"""

self.rho = rho
self.T = T
self.Ye = Ye

self.comp = comp
self.ydots = ydots
self.enu = enu

# compute the bits we need for the table

if comp:
# mean molecular weight of the full NSE state
self.abar = comp.eval_abar()

# average binding energy / nucleon for the full NSE state
self.bea = sum(q.nucbind * self.comp.X[q] for q in self.comp.X)

# evolution of the electron fraction from weak rates alone
self.dYedt = sum(q.Z * self.ydots[q] for q in self.comp.X)

# evolution of Abar from weak rates alone
self.dabardt = -self.abar**2 * sum(self.ydots[q] for q in self.comp.X)

# evolution of B/A from weak rates alone
self.dbeadt = sum(self.ydots[q] * q.nucbind for q in self.comp.X)

self.X = None
if comp_reduction_func:
self.X = comp_reduction_func(self.comp)

def __str__(self):
return f"({self.rho:12.6g}, {self.T:12.6g}, {self.Ye:6.4f}): {self.abar:6.4f} {self.bea:6.4f} {self.dYedt:12.6g} {self.enu:12.6g}"

def value(self):
"""a simple integer used for sorting. This has the format
(logrho)(logT)(1-Ye)"""

return int(f"{np.log10(self.rho):08.5f}{np.log10(self.T):08.5f}{1-self.Ye:08.5f}".replace(".", ""))

def __lt__(self, other):
return self.value() < other.value()


class NSENetwork(RateCollection):
Expand Down Expand Up @@ -101,7 +161,8 @@ def _constraint_eq(self, u, u_c, state):

return [eq1, eq2]

def get_comp_nse(self, rho, T, ye, init_guess=(-3.5, -15), tol=1.0e-11, use_coulomb_corr=False, return_sol=False):
def get_comp_nse(self, rho, T, ye, init_guess=(-3.5, -15),
tol=1.0e-11, use_coulomb_corr=False, return_sol=False):
"""
Returns the NSE composition given density, temperature and prescribed electron fraction
using scipy.fsolve.
Expand All @@ -123,7 +184,7 @@ def get_comp_nse(self, rho, T, ye, init_guess=(-3.5, -15), tol=1.0e-11, use_coul
return_sol: Whether to return the solution of the proton and neutron chemical potential.
"""

#From here we convert the init_guess list into a np.array object:
# From here we convert the init_guess list into a np.array object:

init_guess = np.array(init_guess)
state = NseState(T, rho, ye)
Expand Down Expand Up @@ -182,3 +243,73 @@ def get_comp_nse(self, rho, T, ye, init_guess=(-3.5, -15), tol=1.0e-11, use_coul
init_guess[0] -= 0.5

raise ValueError("Unable to find a solution, try to adjust initial guess manually")

def generate_table(self, rho_values=None, T_values=None, Ye_values=None,
comp_reduction_func=None,
verbose=False, outfile="nse.tbl"):

# initial guess
mu_p0 = -3.5
mu_n0 = -15.0

# arrays to cache the chemical potentials as mu_p(rho, Ye)
mu_p = np.ones((len(rho_values), len(Ye_values)), dtype=np.float64) * mu_p0
mu_n = np.ones((len(rho_values), len(Ye_values)), dtype=np.float64) * mu_n0

nse_states = []
for T in reversed(T_values):
for irho, rho in enumerate(reversed(rho_values)):
for iye, ye in enumerate(reversed(Ye_values)):
initial_guess = (mu_p[irho, iye], mu_n[irho, iye])
try:
comp, sol = self.get_comp_nse(rho, T, ye, use_coulomb_corr=True,
init_guess=initial_guess,
return_sol=True)
except ValueError:
initial_guess = (-3.5, -15)
comp, sol = self.get_comp_nse(rho, T, ye, use_coulomb_corr=True,
init_guess=initial_guess,
return_sol=True)

mu_p[irho, iye] = sol[0]
mu_n[irho, iye] = sol[1]

# get the dY/dt for just the weak rates
ydots = self.evaluate_ydots(rho, T, comp,
screen_func=potekhin_1998,
rate_filter=lambda r: isinstance(r, TabularRate))

_, enu = self.evaluate_energy_generation(rho, T, comp,
screen_func=potekhin_1998,
return_enu=True)

nse_states.append(NSETableEntry(rho, T, ye,
comp=comp, ydots=ydots, enu=enu,
comp_reduction_func=comp_reduction_func))
if verbose:
print(nse_states[-1])

with open(outfile, "w") as of:

# write the header
of.write(f"# NSE table generated by pynucastro {version}\n")
of.write(f"# original NSENetwork had {len(self.unique_nuclei)} nuclei\n")
of.write("#\n")
of.write(f"# {'log10(rho)':^15} {'log10(T)':^15} {'Ye':^15} ")
of.write(f"{'Abar':^15} {'<B/A>':^15} {'dYe/dt':^15} {'dAbar/dt':^15} {'d<B/A>/dt':^15} {'e_nu':^15} ")

if nse_states[0].X:
for nuc, _ in nse_states[0].X:
_tmp = f"X({nuc})"
of.write(f"{_tmp:^15} ")

of.write("\n")

for entry in sorted(nse_states):
of.write(f"{np.log10(entry.rho):15.10f} {np.log10(entry.T):15.10f} {entry.Ye:15.10f} ")
of.write(f"{entry.abar:15.10f} {entry.bea:15.10f} {entry.dYedt:15.8g} {entry.dabardt:15.8g} {entry.dbeadt:15.8g} {entry.enu:15.8g} ")

if entry.X:
for _, val in entry.X:
of.write(f"{val:15.10g} ")
of.write("\n")
40 changes: 40 additions & 0 deletions pynucastro/networks/tests/_nse_table/nse.tbl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,40 @@
# NSE table generated by pynucastro 2.0.1.post2+g184debe8
# original NSENetwork had 13 nuclei
#
# log10(rho) log10(T) Ye Abar <B/A> dYe/dt dAbar/dt d<B/A>/dt e_nu X(n) X(p) X(He4) X(Fe52) X(Fe54) X(Ni56)
7.0000000000 9.6000000000 0.5000000000 50.2057480859 8.6280577146 -4.3803693e-05 1.7080366e-17 2.3287883e-06 7.4773253e+13 8.545469735e-12 0.001886132128 0.0006595854585 0.02734776519 0.08349445037 0.8866120669
7.0000000000 9.6000000000 0.4650000000 55.7461136880 8.7870609269 -5.5255921e-09 -1.8576389e-21 3.9319171e-10 1.39808e+10 9.276985777e-09 1.02578239e-06 0.0002299104701 3.490397709e-06 0.999765231 3.330322786e-07
7.0000000000 9.6000000000 0.4300000000 11.0637644849 8.1412160175 0.00017699678 -2.9553491e-26 -1.8453762e-12 2.7707299e+14 0.07384628728 1.129016861e-14 1.733537583e-06 1e-16 0.9261519792 6.132411849e-37
7.0000000000 10.0000000000 0.5000000000 1.2148718034 1.6681978229 0.084056626 -4.1382279e-56 5.7341747e-41 7.8635859e+17 0.3820880783 0.3820880783 0.2358238434 1e-16 1e-16 5.808870622e-40
7.0000000000 10.0000000000 0.4650000000 1.2129151675 1.6556736890 0.10066497 -2.9530597e-56 3.521478e-41 8.3544486e+17 0.4179733119 0.3479733119 0.2340533762 1e-16 1e-16 5.158080357e-40
7.0000000000 10.0000000000 0.4300000000 1.2071212811 1.6183497719 0.11786424 -1.1606804e-56 5.9185559e-42 8.8944373e+17 0.4556114491 0.3156114491 0.2287771018 1e-16 1e-16 3.696650891e-40
7.0000000000 10.4000000000 0.5000000000 1.0000000000 0.0000000002 5.7463273 2.7835262e-218 6.7437141e-203 2.3050257e+20 0.5 0.5 2.860399672e-11 1e-16 1e-16 5.077180341e-210
7.0000000000 10.4000000000 0.4650000000 1.0000000000 0.0000000002 7.2542055 -3.8299905e-219 5.8186031e-203 2.360043e+20 0.535 0.465 2.832423608e-11 1e-16 1e-16 4.40642192e-210
7.0000000000 10.4000000000 0.4300000000 1.0000000000 0.0000000002 8.7650901 2.0886837e-218 3.8274857e-203 2.4155562e+20 0.57 0.43 2.749345673e-11 1e-16 1e-16 2.892199116e-210
8.0000000000 9.6000000000 0.5000000000 54.3195129673 8.6383992184 -0.00084875851 -0 4.3700394e-05 1.4743336e+15 3.164505989e-13 0.0005147240754 6.760472336e-05 0.02041760492 0.02585276493 0.9531473014
8.0000000000 9.6000000000 0.4650000000 55.8970081563 8.7873915908 -4.9040515e-08 -2.5845056e-21 4.1959348e-09 1.2982252e+11 9.216342753e-10 1.055115869e-07 2.410540533e-05 3.55855299e-06 0.9999717635 4.660594303e-07
8.0000000000 9.6000000000 0.4300000000 11.0638264093 8.1412199515 9.2047376e-06 -3.3885784e-29 -5.6069156e-13 1.1254068e+13 0.07384616075 1.803180722e-16 8.965651972e-08 1e-16 0.9261537496 4.095837587e-41
8.0000000000 10.0000000000 0.5000000000 2.5712542408 5.7636820538 -0.019589298 -1.8701585e-33 1.7211947e-19 1.6311e+17 0.09261017033 0.09261017033 0.8147796593 3.169318178e-16 4.009614583e-14 4.062724573e-19
8.0000000000 10.0000000000 0.4650000000 2.5124362034 5.6778066014 -0.0041256975 5.0812947e-34 2.1065809e-19 1.5065038e+17 0.1336800378 0.06368003776 0.8026399245 3.471532475e-16 1.203035445e-13 3.184715791e-19
8.0000000000 10.0000000000 0.4300000000 2.3664067080 5.4461446407 0.0093041487 -6.5089027e-34 1.1442051e-19 1.6169818e+17 0.1850544189 0.0450544189 0.7698911622 2.916205557e-16 2.279128e-13 1.716449901e-19
8.0000000000 10.4000000000 0.5000000000 1.0000000214 0.0000002024 1.213841 7.2035365e-165 9.8460077e-150 2.1894106e+20 0.4999999857 0.4999999857 2.861178959e-08 1e-16 1e-16 6.391509301e-155
8.0000000000 10.4000000000 0.4650000000 1.0000000212 0.0000002004 2.8346072 1.8145634e-165 8.4158915e-150 2.2127881e+20 0.5349999858 0.4649999858 2.833167452e-08 1e-16 1e-16 5.508122145e-155
8.0000000000 10.4000000000 0.4300000000 1.0000000206 0.0000001945 4.4726034 2.4538451e-165 5.481869e-150 2.2397331e+20 0.5699999863 0.4299999863 2.750039456e-08 1e-16 1e-16 3.588734491e-155
9.0000000000 9.6000000000 0.5000000000 55.5136288840 8.6413421175 -0.060154436 -0 0.0031211006 1.5552556e+17 1.081519091e-14 0.000138934617 5.722712952e-06 0.01322003813 0.007498457684 0.9791368469
9.0000000000 9.6000000000 0.4650000000 55.9133553598 8.7874247178 -3.9839589e-06 7.6546118e-19 3.6046977e-07 1.0369536e+13 9.223625557e-11 9.830431177e-09 2.121537165e-06 3.550122966e-06 0.9999935846 7.338313681e-07
9.0000000000 9.6000000000 0.4300000000 11.0638296395 8.1412201568 1.0225013e-08 -9.8142004e-34 -5.2463937e-15 1.1664591e+10 0.07384615415 1.005278937e-16 3.918883194e-09 1e-16 0.9261538419 3.163731285e-45
9.0000000000 10.0000000000 0.5000000000 4.3121615772 7.1515898299 -0.10742116 -7.8751776e-19 2.9757777e-05 4.7264232e+17 0.01126674095 0.02343959399 0.7751073898 0.0006022183527 0.189572885 1.117195931e-05
9.0000000000 10.0000000000 0.4650000000 5.1247240232 7.3995844381 -0.019430675 -5.218123e-20 3.3567416e-06 8.7381876e+16 0.04663429738 0.004754114617 0.5461854789 2.394198789e-05 0.4024020903 7.681641422e-08
9.0000000000 10.0000000000 0.4300000000 4.4162729378 7.0777579412 -0.0059568824 1.4777727e-20 5.1854108e-07 3.2607638e+16 0.1092190107 0.001798710239 0.4287763004 2.233180113e-06 0.4602037431 2.387093977e-09
9.0000000000 10.4000000000 0.5000000000 1.0000214391 0.0002025494 -17.715252 1.5106478e-112 3.7266649e-96 2.7735801e+20 0.4999857015 0.4999857015 2.863327916e-05 1e-16 1e-16 1.16193348e-99
9.0000000000 10.4000000000 0.4650000000 1.0000212298 0.0002005604 -15.005795 -5.5593438e-112 3.0826278e-96 2.5962747e+20 0.5349858413 0.4649858413 2.83521049e-05 1e-16 1e-16 9.841859964e-100
9.0000000000 10.4000000000 0.4300000000 1.0000206076 0.0001946699 -12.327315 -6.3674509e-112 1.9400005e-96 2.4324142e+20 0.5699862565 0.4299862565 2.751939962e-05 1e-16 1e-16 6.297337464e-100
10.0000000000 9.6000000000 0.5000000000 55.8692197296 8.6423237057 -5.2276094 -0 0.2720124 2.9376254e+19 4.270621897e-16 3.352682994e-05 3.219778383e-07 0.006375782936 0.001914399324 0.9916759689
10.0000000000 9.6000000000 0.4650000000 55.9152992729 8.7874226105 -0.00092940692 -2.648266e-16 8.5387004e-05 4.3026908e+15 9.257520104e-12 7.48047395e-10 1.2818476e-07 3.511642285e-06 0.9999947952 1.564231185e-06
10.0000000000 9.6000000000 0.4300000000 11.0638297827 8.1412201659 5.0971987e-15 -9.3107493e-42 -6.8646484e-21 5801.4456 0.07384615386 1.00002858e-16 1.18839261e-10 1e-16 0.926153846 3.322511754e-49
10.0000000000 10.0000000000 0.5000000000 13.9370868136 8.2743459363 -5.5785841 1.0782619e-14 0.10961074 3.9381535e+19 0.0001306586039 0.02534088625 0.1230842434 0.1129118461 0.7005437567 0.03798860903
10.0000000000 10.0000000000 0.4650000000 28.2438967467 8.6492551831 -0.075078082 2.0321472e-16 0.00050736054 5.2201774e+17 0.004701104652 0.0004649937736 0.05362316347 5.467021651e-05 0.9411557864 2.814928867e-07
10.0000000000 10.0000000000 0.4300000000 10.2419420841 8.0844779277 -0.0025107223 -1.1694167e-19 1.1073417e-06 1.7236419e+16 0.07573604143 1.89744648e-05 0.02316161685 1.349398617e-08 0.9010833538 1.832667784e-12
10.0000000000 10.4000000000 0.5000000000 1.0197762456 0.1829100336 -196.6445 -4.638832e-57 3.90578e-42 2.3555844e+21 0.4870715132 0.4870715132 0.02585697363 1e-16 1e-16 1.174425453e-44
10.0000000000 10.4000000000 0.4650000000 1.0195858994 0.1811833443 -167.54221 3.2294093e-57 3.0442073e-42 1.9968628e+21 0.5221935594 0.4521935594 0.02561288117 1e-16 1e-16 9.619150094e-45
10.0000000000 10.4000000000 0.4300000000 1.0190251146 0.1760925439 -140.99607 -1.1373217e-57 1.8183284e-42 1.6720803e+21 0.5575533885 0.4175533885 0.02489322305 1e-16 1e-16 6.01760526e-45
73 changes: 73 additions & 0 deletions pynucastro/networks/tests/test_nse_table.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,73 @@
import os

import numpy as np
import pytest

import pynucastro as pyna
from pynucastro import Nucleus


class TestNSETable:
@pytest.fixture(scope="class")
def nse_net(self, reaclib_library, tabular_library):

nucs = [Nucleus("p"), Nucleus("n"), Nucleus("he4"),
Nucleus("fe52"), Nucleus("fe53"), Nucleus("fe54"),
Nucleus("fe55"), Nucleus("fe56"),
Nucleus("co54"), Nucleus("co55"), Nucleus("co56"),
Nucleus("ni56"), Nucleus("ni57")]

tlib = tabular_library.linking_nuclei(nucs)
rlib = reaclib_library.linking_nuclei(nucs)

all_lib = rlib + tlib

dupes = all_lib.find_duplicate_links()

rates_to_remove = []
for d in dupes:
rates_to_remove += [r for r in d if isinstance(r, pyna.rates.ReacLibRate)]

for r in rates_to_remove:
all_lib.remove_rate(r)

nse = pyna.NSENetwork(libraries=[all_lib])

return nse

@staticmethod
def get_reduced_comp(comp):
nuc_list = [Nucleus("p"), Nucleus("n"), Nucleus("he4"),
Nucleus("fe52"), Nucleus("fe54"), Nucleus("ni56")]
reduced_comp = comp.bin_as(nuc_list, exclude=[Nucleus("ni56")])
X = []
for n in reduced_comp.X:
X.append((f"{n}", reduced_comp.X[n]))
return X

def test_generate_table(self, nse_net):

Ts = np.logspace(9.6, 10.4, 3)
rhos = np.logspace(7, 10, 4)
yes = np.linspace(0.43, 0.5, 3)

nse_net.generate_table(rho_values=rhos,
T_values=Ts,
Ye_values=yes,
comp_reduction_func=self.get_reduced_comp)

# this creates a file called `nse.tbl`, which we want to compare
# to the stored benchmark

base_path = os.path.relpath(os.path.dirname(__file__))
ref_path = os.path.join(base_path, "_nse_table")

with open("nse.tbl") as new_table, open(f"{ref_path}/nse.tbl") as ref_table:

new_lines = new_table.readlines()
ref_lines = ref_table.readlines()

for new, ref in zip(new_lines, ref_lines):
if new.startswith("#"):
continue
assert new == ref