diff --git a/networks/he-burn/he-burn-22a/he_burn_22a.py b/networks/he-burn/he-burn-22a/he_burn_22a.py index 0ef6826383..07958a9b59 100644 --- a/networks/he-burn/he-burn-22a/he_burn_22a.py +++ b/networks/he-burn/he-burn-22a/he_burn_22a.py @@ -22,13 +22,12 @@ def doit(): r1 = subch.get_rate_by_name("c12(p,g)n13") r2 = subch.get_rate_by_name("n13(he4,p)o16") - net = AmrexAstroCxxNetwork(libraries=[subch], symmetric_screening=True, disable_rate_params=[r1, r2]) + net = AmrexAstroCxxNetwork(libraries=[subch], symmetric_screening=True, + disable_rate_params=[r1, r2]) + net.make_ap_pg_approx(intermediate_nuclei=["cl35", "k39", "sc43", "v47", "mn51", "co55"]) net.remove_nuclei(["cl35", "k39", "sc43", "v47", "mn51", "co55"]) - # finally, the aprox nets don't include the reverse rates for - # C12+C12, C12+O16, and O16+O16, so remove those - print(f"number of nuclei: {len(net.unique_nuclei)}") print(f"number of rates: {len(net.rates)}") diff --git a/networks/he-burn/he-burn-36a/He-C-Fe-group.png b/networks/he-burn/he-burn-36a/He-C-Fe-group.png deleted file mode 100644 index 0970da4ec4..0000000000 Binary files a/networks/he-burn/he-burn-36a/He-C-Fe-group.png and /dev/null differ diff --git a/networks/he-burn/he-burn-36a/He-C-Fe-group.py b/networks/he-burn/he-burn-36a/He-C-Fe-group.py deleted file mode 100644 index 2a9a7940dd..0000000000 --- a/networks/he-burn/he-burn-36a/He-C-Fe-group.py +++ /dev/null @@ -1,142 +0,0 @@ -import pynucastro as pyna -from pynucastro.rates import ReacLibRate, TabularRate - -DO_DERIVED_RATES = True - -reaclib_lib = pyna.ReacLibLibrary() -weak_lib = pyna.TabularLibrary() - -# these are the nuclei we have in subch_simple - -all_reactants = ["p", - "he4", "c12", "o16", "ne20", "mg24", "si28", "s32", - "ar36", "ca40", "ti44", "cr48", "fe52", "ni56", - "cu59", "zn60", - "al27", "p31", "cl35", "k39", "sc43", "v47", "mn51", "co55", - "n13", "n14", "f18", "ne21", "na22", "na23"] - -# create a library of ReacLib rates - -core_lib = reaclib_lib.linking_nuclei(all_reactants) - -# in this list, we have the reactants, the actual reactants, -# and modified products that we will use instead - -other_rates = [("c12(c12,n)mg23", "mg24"), - ("o16(o16,n)s31", "s32"), - ("o16(c12,n)si27", "si28")] - -for r, mp in other_rates: - _r = reaclib_lib.get_rate_by_name(r) - _r.modify_products(mp) - core_lib.add_rate(_r) - -# finally, the aprox nets don't include the reverse rates for -# C12+C12, C12+O16, and O16+O16, so remove those - -for r in core_lib.get_rates(): - if sorted(r.products) in [[pyna.Nucleus("c12"), pyna.Nucleus("c12")], - [pyna.Nucleus("c12"), pyna.Nucleus("o16")], - [pyna.Nucleus("o16"), pyna.Nucleus("o16")]]: - core_lib.remove_rate(r) - -# C12+Ne20 and reverse -# (a,g) links between Na23 and Al27 -# (a,g) links between Al27 and P31 - -rates_to_remove = ["p31(p,c12)ne20", - "si28(a,c12)ne20", - "ne20(c12,p)p31", - "ne20(c12,a)si28", - "na23(a,g)al27", - "al27(g,a)na23", - "al27(a,g)p31", - "p31(g,a)al27"] - -for r in rates_to_remove: - print("removing: ", r) - _r = core_lib.get_rate_by_name(r) - core_lib.remove_rate(_r) - -# now create a list of iron group nuclei and find both the -# ReacLib and weak / tabular rates linking these. - -iron_peak = ["n", "p", "he4", - "mn51", "mn55", - "fe52", "fe53", "fe54", "fe55", "fe56", - "co55", "co56", "co57", - "ni56", "ni57", "ni58", - "cu59", "zn60"] - -iron_reaclib = reaclib_lib.linking_nuclei(iron_peak) - -iron_weak_lib = weak_lib.linking_nuclei(iron_peak) - -# add the libraries - -all_lib = core_lib + iron_reaclib + iron_weak_lib - -if DO_DERIVED_RATES: - rates_to_derive = [] - for r in all_lib.get_rates(): - if r.reverse: - # this rate was computed using detailed balance, regardless - # of whether Q < 0 or not. We want to remove it and then - # recompute it - rates_to_derive.append(r) - - # now for each of those derived rates, look to see if the pair exists - - for r in rates_to_derive: - fr = all_lib.get_rate_by_nuclei(r.products, r.reactants) - if fr: - print(f"modifying {r} from {fr}") - all_lib.remove_rate(r) - d = pyna.DerivedRate(rate=fr, compute_Q=False, use_pf=True) - all_lib.add_rate(d) - -# we will have duplicate rates -- we want to remove any ReacLib rates -# that we have tabular rates for - -dupes = all_lib.find_duplicate_links() - -rates_to_remove = [] -for d in dupes: - for r in d: - if isinstance(r, ReacLibRate): - rates_to_remove.append(r) - -for r in rates_to_remove: - all_lib.remove_rate(r) - -# combine all three libraries into a single network - -net = pyna.AmrexAstroCxxNetwork(libraries=[all_lib], - symmetric_screening=False) - - -# now we approximate some (alpha, p)(p, gamma) links - -net.make_ap_pg_approx(intermediate_nuclei=["cl35", "k39", "sc43", "v47"]) -net.remove_nuclei(["cl35", "k39", "sc43", "v47"]) - -# let's make a figure - -T = 6.e9 -rho = 9.e7 -comp = pyna.Composition(net.unique_nuclei) -comp.set_all(1.0) -comp.normalize() - -fig = net.plot(rho, T, comp, - rotated=True, curved_edges=True, hide_xalpha=True, - size=(1800, 900), - node_size=500, node_shape="s", node_font_size=10) - -fig.savefig("He-C-Fe-group.png") - -print(f"number of nuclei = {len(net.unique_nuclei)}") -print(f"number of ReacLib rates = {len(net.reaclib_rates)}") -print(f"number of tabular rates = {len(net.tabular_rates)}") - -net.write_network() diff --git a/networks/he-burn/he-burn-36a/actual_rhs.H b/networks/he-burn/he-burn-36a/actual_rhs.H index 21aff67909..9b2aaf46f1 100644 --- a/networks/he-burn/he-burn-36a/actual_rhs.H +++ b/networks/he-burn/he-burn-36a/actual_rhs.H @@ -1270,7 +1270,7 @@ void evaluate_rates(const burn_t& state, T& rate_eval) { // Fill approximate rates - fill_approx_rates(tfactors, rate_eval); + fill_approx_rates(tfactors, state.rho, Y, rate_eval); // Calculate tabular rates diff --git a/networks/he-burn/he-burn-36a/he-burn-36a.png b/networks/he-burn/he-burn-36a/he-burn-36a.png new file mode 100644 index 0000000000..47578149b8 Binary files /dev/null and b/networks/he-burn/he-burn-36a/he-burn-36a.png differ diff --git a/networks/he-burn/he-burn-36a/he_burn_36a.py b/networks/he-burn/he-burn-36a/he_burn_36a.py new file mode 100644 index 0000000000..4c04aede76 --- /dev/null +++ b/networks/he-burn/he-burn-36a/he_burn_36a.py @@ -0,0 +1,48 @@ +import pynucastro as pyna +from pynucastro.networks import AmrexAstroCxxNetwork + +import he_burn_core + + +DO_DERIVED_RATES = True + + +def doit(): + + lib = he_burn_core.get_core_library(include_n14_sequence=True, + include_zn=True, + include_iron_peak=True, + do_detailed_balance=DO_DERIVED_RATES) + + net = pyna.AmrexAstroCxxNetwork(libraries=[lib], + symmetric_screening=False) + + # now we approximate some (alpha, p)(p, gamma) links + + net.make_ap_pg_approx(intermediate_nuclei=["cl35", "k39", "sc43", "v47"]) + net.remove_nuclei(["cl35", "k39", "sc43", "v47"]) + + print(f"number of nuclei = {len(net.unique_nuclei)}") + print(f"number of ReacLib rates = {len(net.reaclib_rates)}") + print(f"number of tabular rates = {len(net.tabular_rates)}") + + # let's make a figure + + comp = pyna.Composition(net.unique_nuclei) + comp.set_equal() + + rho = 9.e7 + T = 6.e9 + + fig = net.plot(rho, T, comp, + rotated=True, curved_edges=True, hide_xalpha=True, + size=(1800, 900), + node_size=500, node_shape="s", node_color="#337dff", node_font_size=10) + + fig.savefig("he-burn-36a.png") + + net.write_network() + + +if __name__ == "__main__": + doit() diff --git a/networks/he-burn/he-burn-36a/he_burn_core.py b/networks/he-burn/he-burn-36a/he_burn_core.py new file mode 120000 index 0000000000..bc9cd39283 --- /dev/null +++ b/networks/he-burn/he-burn-36a/he_burn_core.py @@ -0,0 +1 @@ +../he_burn_core.py \ No newline at end of file diff --git a/networks/he-burn/he-burn-36a/reaclib_rates.H b/networks/he-burn/he-burn-36a/reaclib_rates.H index d423ce14dc..58ef65c6a4 100644 --- a/networks/he-burn/he-burn-36a/reaclib_rates.H +++ b/networks/he-burn/he-burn-36a/reaclib_rates.H @@ -9208,7 +9208,10 @@ fill_reaclib_rates(const tf_t& tfactors, T& rate_eval) template AMREX_GPU_HOST_DEVICE AMREX_INLINE void -fill_approx_rates([[maybe_unused]] const tf_t& tfactors, [[maybe_unused]] T& rate_eval) +fill_approx_rates([[maybe_unused]] const tf_t& tfactors, + [[maybe_unused]] const amrex::Real rho, + [[maybe_unused]] const amrex::Array1D& Y, + [[maybe_unused]] T& rate_eval) { [[maybe_unused]] amrex::Real rate{}; diff --git a/networks/he-burn/he_burn_core.py b/networks/he-burn/he_burn_core.py index d685ca4210..7a47e24cac 100644 --- a/networks/he-burn/he_burn_core.py +++ b/networks/he-burn/he_burn_core.py @@ -2,8 +2,12 @@ # they can then adjust these via various approximations import pynucastro as pyna +from pynucastro.rates import ReacLibRate, TabularRate -def get_core_library(*, include_n14_sequence=False, include_zn=False, +def get_core_library(*, + include_n14_sequence=False, + include_zn=False, + include_iron_peak=False, do_detailed_balance=False): reaclib_lib = pyna.ReacLibLibrary() @@ -62,6 +66,27 @@ def get_core_library(*, include_n14_sequence=False, include_zn=False, _r = core_lib.get_rate_by_name(r) core_lib.remove_rate(_r) + if include_iron_peak: + # now create a list of iron group nuclei and find both the + # ReacLib and weak / tabular rates linking these. + + iron_peak = ["n", "p", "he4", + "mn51", "mn55", + "fe52", "fe53", "fe54", "fe55", "fe56", + "co55", "co56", "co57", + "ni56", "ni57", "ni58", + "cu59", "zn60"] + + iron_reaclib = reaclib_lib.linking_nuclei(iron_peak) + + weak_lib = pyna.TabularLibrary() + iron_weak_lib = weak_lib.linking_nuclei(iron_peak) + + all_lib = core_lib + iron_reaclib + iron_weak_lib + + else: + all_lib = core_lib + if do_detailed_balance: rates_to_derive = [] for r in core_lib.get_rates(): @@ -81,4 +106,18 @@ def get_core_library(*, include_n14_sequence=False, include_zn=False, d = pyna.DerivedRate(rate=fr, compute_Q=False, use_pf=True) core_lib.add_rate(d) - return core_lib + # we may have duplicate rates -- we want to remove any ReacLib rates + # that we have tabular rates for + + dupes = all_lib.find_duplicate_links() + + rates_to_remove = [] + for d in dupes: + for r in d: + if isinstance(r, ReacLibRate): + rates_to_remove.append(r) + + for r in rates_to_remove: + all_lib.remove_rate(r) + + return all_lib