reorganize the He nets (#1687)

They are now all placed under he-burn and named according to how many nuclei they have
and what approximations they use.

This also adds a new network with 31 nuclei that has a better iron group that the old
subch_simple, and uses the (nn,gamma) approximation to reduce the nuclei, and
swaps protons out for NSE protons at A >= 48.

networks/subch_base/actual_rhs.H → networks/he-burn/he-burn-18a/actual_rhs.H

@@ -798,7 +798,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-18a/he_burn_18a.py

@@ -0,0 +1,52 @@
+# Network for He/C burning with key rates
+# to bypass the C12(a,g)O16 rate.
+
+import pynucastro as pyna
+from pynucastro.networks import AmrexAstroCxxNetwork
+
+import he_burn_core
+
+
+DO_DERIVED_RATES = True
+
+
+def doit():
+
+    subch = he_burn_core.get_core_library(include_n14_sequence=False,
+                                          include_zn=False,
+                                          do_detailed_balance=DO_DERIVED_RATES)
+
+    # these are the rates that we are going to allow to be optionally
+    # zeroed
+    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=False, 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)}")
+
+    comp = pyna.Composition(net.get_nuclei())
+    comp.set_all(0.1)
+    comp.set_nuc("he4", 0.95)
+    comp.normalize()
+
+    rho = 1.e6
+    T = 1.e9
+
+    net.plot(rho, T, comp, outfile="he-burn-18a.png",
+             rotated=True, hide_xalpha=True, curved_edges=True,
+             size=(1500, 450),
+             node_size=500, node_font_size=11, node_color="#337dff", node_shape="s",
+             Z_range=(1, 29))
+
+    net.write_network()
+
+
+if __name__ == "__main__":
+    doit()

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

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

networks/subch_base/reaclib_rates.H → networks/he-burn/he-burn-18a/reaclib_rates.H

@@ -5001,7 +5001,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/subch_simple/actual_rhs.H → networks/he-burn/he-burn-22a/actual_rhs.H

@@ -1032,7 +1032,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-22a/he_burn_22a.py

@@ -0,0 +1,52 @@
+# an approximate network for He/C burning with key rates
+# to bypass the C12(a,g)O16 rate.  This version uses some
+# (a,p)(p,g) approximations.
+
+import pynucastro as pyna
+from pynucastro.networks import AmrexAstroCxxNetwork
+
+import he_burn_core
+
+
+DO_DERIVED_RATES = False
+
+
+def doit():
+
+    subch = he_burn_core.get_core_library(include_n14_sequence=True,
+                                          include_zn=False,
+                                          do_detailed_balance=DO_DERIVED_RATES)
+
+    # these are the rates that we are going to allow to be optionally
+    # zeroed
+    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.make_ap_pg_approx(intermediate_nuclei=["cl35", "k39", "sc43", "v47", "mn51", "co55"])
+    net.remove_nuclei(["cl35", "k39", "sc43", "v47", "mn51", "co55"])
+
+    print(f"number of nuclei: {len(net.unique_nuclei)}")
+    print(f"number of rates: {len(net.rates)}")
+
+    comp = pyna.Composition(net.get_nuclei())
+    comp.set_all(0.1)
+    comp.set_nuc("he4", 0.95)
+    comp.normalize()
+
+    rho = 1.e6
+    T = 1.e9
+
+    net.plot(rho, T, comp, outfile="he-burn-22a.png",
+             rotated=True, hide_xalpha=True, curved_edges=True,
+             size=(1500, 450),
+             node_size=500, node_font_size=11, node_color="#337dff", node_shape="s",
+             Z_range=(1,29))
+
+    net.write_network()
+
+
+if __name__ == "__main__":
+    doit()

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

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

networks/subch_simple/reaclib_rates.H → networks/he-burn/he-burn-22a/reaclib_rates.H

@@ -4614,7 +4614,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{};