Merge branch 'development' into jacobian_correction

EOS/helmholtz/actual_eos.H

@@ -156,9 +156,9 @@ void apply_electrons (T& state)
 
     // hash locate this temperature and density
     int jat = int((std::log10(state.T) - tlo) * tstpi) + 1;
-    jat = amrex::max(1, amrex::min(jat, jmax-1)) - 1;
+    jat = std::clamp(jat, 1, jmax-1) - 1;
     int iat = int((std::log10(din) - dlo) * dstpi) + 1;
-    iat = amrex::max(1, amrex::min(iat, imax-1)) - 1;
+    iat = std::clamp(iat, 1, imax-1) - 1;
 
     amrex::Real fi[36];
 
@@ -1005,7 +1005,7 @@ void single_iter_update (T& state, int var, int dvar,
     amrex::Real xnew = x - (v - v_want) / dvdx;
 
     // Don't let the temperature/density change by more than a factor of two
-    xnew = amrex::max(0.5_rt * x, amrex::min(xnew, 2.0_rt * x));
+    xnew = std::clamp(xnew, 0.5_rt * x, 2.0_rt * x);
 
     // Don't let us freeze/evacuate
     xnew = amrex::max(smallx, xnew);
@@ -1120,8 +1120,8 @@ void double_iter_update (T& state, int var1, int var2,
 
     // Don't let the temperature or density change by more
     // than a factor of two
-    tnew = amrex::max(0.5e0_rt * told, amrex::min(tnew, 2.0e0_rt * told));
-    rnew = amrex::max(0.5e0_rt * rold, amrex::min(rnew, 2.0e0_rt * rold));
+    tnew = std::clamp(tnew, 0.5e0_rt * told, 2.0e0_rt * told);
+    rnew = std::clamp(rnew, 0.5e0_rt * rold, 2.0e0_rt * rold);
 
     // Don't let us freeze or evacuate
     tnew = amrex::max(EOSData::mintemp, tnew);

integration/BackwardEuler/be_integrator.H

@@ -272,7 +272,7 @@ int be_integrator (BurnT& state, BeT& be)
             // backward-Euler has a local truncation error of dt**2
 
             amrex::Real dt_new = dt_sub * std::pow(1.0_rt / rel_error, 0.5_rt);
-            dt_sub = amrex::min(amrex::max(dt_new, dt_sub / 2.0), 2.0 * dt_sub);
+            dt_sub = std::clamp(dt_new, dt_sub / 2.0, 2.0 * dt_sub);
 
         } else {
 

integration/ForwardEuler/actual_integrator.H

@@ -86,7 +86,7 @@ void clean_state (const amrex::Real time, IntT& int_state, BurnT& state)
 
     // Ensure that the temperature always stays within reasonable limits.
 
-    state.T = amrex::min(MAX_TEMP, amrex::max(state.T, EOSData::mintemp));
+    state.T = std::clamp(state.T, EOSData::mintemp, MAX_TEMP);
 
 }
 

integration/RKC/rkc.H

@@ -308,7 +308,7 @@ int rkclow (BurnT& state, RkcT& rstate)
             }
         }
         absh = std::max(0.1_rt, fac) * absh;
-        absh = std::max(hmin, std::min(rstate.hmax, absh));
+        absh = std::clamp(absh, hmin, rstate.hmax);
         errold = err;
         hold = h;
         h = tdir * absh;

integration/VODE/vode_dvhin.H

@@ -1,6 +1,8 @@
 #ifndef VODE_DVHIN_H
 #define VODE_DVHIN_H
 
+#include <AMReX_REAL.H>
+
 #ifdef STRANG
 #include <integrator_rhs_strang.H>
 #endif
@@ -23,7 +25,7 @@ void dvhin (BurnT& state, DvodeT& vstate, amrex::Real& H0, int& NITER, int& IER)
 
     constexpr int int_neqs = integrator_neqs<BurnT>();
 
-    const amrex::Real PT1 = 0.1e0_rt;
+    constexpr amrex::Real PT1 = 0.1e0_rt;
 
     NITER = 0;
     const amrex::Real TDIST = std::abs(vstate.tout - vstate.t);

integration/VODE/vode_dvjac.H

@@ -116,7 +116,7 @@ void dvjac (int& IERPJ, BurnT& state, DvodeT& vstate)
                 R0 = 1.0_rt;
             }
 
-            const bool in_jacobian = true;
+            constexpr bool in_jacobian = true;
             for (int j = 1; j <= int_neqs; ++j) {
                 const amrex::Real yj = vstate.y(j);
 

integration/VODE/vode_dvnlsd.H

@@ -20,11 +20,11 @@ amrex::Real dvnlsd (int& NFLAG, BurnT& state, DvodeT& vstate)
     // It then handles the corrector phase of this integration package.
 
     // Parameter declarations
-    const amrex::Real CCMAX = 0.3e0_rt;
-    const amrex::Real CRDOWN = 0.3e0_rt;
-    const amrex::Real RDIV  = 2.0e0_rt;
-    const int MAXCOR = 3;
-    const int MSBP = 20;
+    constexpr amrex::Real CCMAX = 0.3e0_rt;
+    constexpr amrex::Real CRDOWN = 0.3e0_rt;
+    constexpr amrex::Real RDIV  = 2.0e0_rt;
+    constexpr int MAXCOR = 3;
+    constexpr int MSBP = 20;
 
     amrex::Real DEL{};
     int M{};

integration/VODE/vode_dvset.H

@@ -20,7 +20,7 @@ void dvset (DvodeT& vstate)
     // H*xi(i) = t sub n  -  t sub (n-i)
     //         = H + TAU(1) + TAU(2) + ... TAU(i-1).
 
-    const amrex::Real CORTES = 0.1e0_rt;
+    constexpr amrex::Real CORTES = 0.1e0_rt;
 
     const amrex::Real FLOTL = vstate.L;
     const int NQM1 = vstate.NQ - 1;

integration/VODE/vode_dvstep.H

@@ -67,20 +67,20 @@ int dvstep (BurnT& state, DvodeT& vstate)
     // On a successful return, ETAMAX is reset and ACOR is scaled.
 
     // Parameter declarations
-    const int KFC = -3;
-    const int KFH = -7;
-    const int MXNCF = 10;
-    const amrex::Real ADDON = 1.0e-6_rt;
-    const amrex::Real BIAS1 = 6.0e0_rt;
-    const amrex::Real BIAS2 = 6.0e0_rt;
-    const amrex::Real BIAS3 = 10.0e0_rt;
-    const amrex::Real ETACF = 0.25e0_rt;
-    const amrex::Real ETAMIN = 0.1e0_rt;
-    const amrex::Real ETAMXF = 0.2e0_rt;
-    const amrex::Real ETAMX2 = 10.0e0_rt;
-    const amrex::Real ETAMX3 = 10.0e0_rt;
-    const amrex::Real ONEPSM = 1.00001e0_rt;
-    const amrex::Real THRESH = 1.5e0_rt;
+    constexpr int KFC = -3;
+    constexpr int KFH = -7;
+    constexpr int MXNCF = 10;
+    constexpr amrex::Real ADDON = 1.0e-6_rt;
+    constexpr amrex::Real BIAS1 = 6.0e0_rt;
+    constexpr amrex::Real BIAS2 = 6.0e0_rt;
+    constexpr amrex::Real BIAS3 = 10.0e0_rt;
+    constexpr amrex::Real ETACF = 0.25e0_rt;
+    constexpr amrex::Real ETAMIN = 0.1e0_rt;
+    constexpr amrex::Real ETAMXF = 0.2e0_rt;
+    constexpr amrex::Real ETAMX2 = 10.0e0_rt;
+    constexpr amrex::Real ETAMX3 = 10.0e0_rt;
+    constexpr amrex::Real ONEPSM = 1.00001e0_rt;
+    constexpr amrex::Real THRESH = 1.5e0_rt;
 
     amrex::Real CNQUOT{}, DDN{}, DSM{}, DUP{}, TOLD{};
     amrex::Real FLOTL{}, R{};

integration/integrator_type_sdc.H

@@ -21,8 +21,7 @@ void clean_state(const amrex::Real time, BurnT& state, T& int_state)
     if (do_species_clip) {
         for (int n = 1; n <= NumSpec; ++n) {
             // we use 1-based indexing, so we need to offset SFS
-            int_state.y(SFS+n) = amrex::max(amrex::min(int_state.y(SFS+n), state.rho),
-                                            state.rho * SMALL_X_SAFE);
+            int_state.y(SFS+n) = std::clamp(int_state.y(SFS+n), state.rho * SMALL_X_SAFE, state.rho);
         }
     }
 

integration/integrator_type_strang.H

@@ -55,7 +55,7 @@ void clean_state (const amrex::Real time, BurnT& state, I& int_state)
 
     if (do_species_clip) {
         for (int n = 1; n <= NumSpec; ++n) {
-            int_state.y(n) = amrex::max(amrex::min(int_state.y(n), 1.0_rt), SMALL_X_SAFE);
+            int_state.y(n) = std::clamp(int_state.y(n), SMALL_X_SAFE, 1.0_rt);
         }
     }
 

integration/utils/initial_timestep.H

@@ -87,7 +87,7 @@ amrex::Real initial_react_dt (BurnT& burn_state, IntT& int_state,
     // Save the final timestep, with a bias factor.
 
     amrex::Real dt = h / 2.0_rt;
-    dt = amrex::min(amrex::max(h, hL), hU);
+    dt = std::clamp(h, hL, hU);
 
     dt = amrex::min(dt, ode_max_dt);
 

networks/CNO_He_burn/CNO_He_burn.py

@@ -0,0 +1,123 @@
+# a blend of CNO_extras and subch_simple
+
+import pynucastro as pyna
+from pynucastro.networks import AmrexAstroCxxNetwork
+
+DO_DERIVED_RATES = False
+
+def get_library():
+
+    reaclib_lib = pyna.ReacLibLibrary()
+
+    all_reactants = ["h1", "he4",
+                     "c12", "c13",
+                     "n13", "n14", "n15",
+                     "o14", "o15", "o16", "o17", "o18",
+                     "f17", "f18", "f19",
+                     "ne18", "ne19", "ne20", "ne21",
+                     "na22", "na23",
+                     "mg22", "mg24",
+                     "al27", "si28", "p31", "s32",
+                     "cl35", "ar36", "k39", "ca40",
+                     "sc43", "ti44", "v47", "cr48",
+                     "mn51", "fe52", "co55", "ni56"]
+
+    subch = 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)
+        subch += pyna.Library(rates=[_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 subch.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")]]:
+            subch.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 = subch.get_rate_by_name(r)
+        subch.remove_rate(_r)
+
+    if DO_DERIVED_RATES:
+        rates_to_derive = []
+        for r in subch.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 = subch.get_rate_by_nuclei(r.products, r.reactants)
+            if fr:
+                print(f"modifying {r} from {fr}")
+                subch.remove_rate(r)
+                d = pyna.DerivedRate(rate=fr, compute_Q=False, use_pf=True)
+                subch.add_rate(d)
+
+    return subch
+
+def doit():
+
+    subch = get_library()
+
+    # 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"])
+
+    # 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="CNO_He_burn.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/CNO_He_burn/Make.package

@@ -0,0 +1,14 @@
+CEXE_headers += network_properties.H
+
+ifeq ($(USE_REACT),TRUE)
+  CEXE_sources += actual_network_data.cpp
+  CEXE_headers += actual_network.H
+  CEXE_headers += tfactors.H
+  CEXE_headers += partition_functions.H
+  CEXE_headers += actual_rhs.H
+  CEXE_headers += reaclib_rates.H
+  CEXE_headers += table_rates.H
+  CEXE_sources += table_rates_data.cpp
+  USE_SCREENING = TRUE
+  USE_NEUTRINOS = TRUE
+endif

networks/CNO_He_burn/_parameters

@@ -0,0 +1,4 @@
+@namespace: network
+
+disable_p_C12_to_N13    int     0
+disable_He4_N13_to_p_O16    int     0