more updates

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)}")
 

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<do_T_derivatives, T>(tfactors, rate_eval);
+    fill_approx_rates<do_T_derivatives, T>(tfactors, state.rho, Y, rate_eval);
 
     // Calculate tabular rates
 

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()

networks/he-burn/he-burn-36a/he_burn_core.py

@@ -0,0 +1 @@
+../he_burn_core.py

networks/he-burn/he-burn-36a/reaclib_rates.H

@@ -9208,7 +9208,10 @@ fill_reaclib_rates(const tf_t& tfactors, T& rate_eval)
 template <int do_T_derivatives, typename T>
 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<amrex::Real, 1, NumSpec>& Y,
+                  [[maybe_unused]] T& rate_eval)
 {
 
     [[maybe_unused]] amrex::Real rate{};

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