Merge branch 'development' into fix_integrator_rp

.github/workflows/jac_cell.yml

@@ -0,0 +1,45 @@
+name: jac_cell
+
+on: [pull_request]
+jobs:
+ nse_table:
+ runs-on: ubuntu-latest
+ steps:
+ - uses: actions/checkout@v4
+ with:
+ fetch-depth: 0
+
+ - name: Get AMReX
+ run: |
+ mkdir external
+ cd external
+ git clone https://github.com/AMReX-Codes/amrex.git
+ cd amrex
+ git checkout development
+ echo 'AMREX_HOME=$(GITHUB_WORKSPACE)/external/amrex' >> $GITHUB_ENV
+ echo $AMREX_HOME
+ if [[ -n "${AMREX_HOME}" ]]; then exit 1; fi
+ cd ../..
+
+ - name: Install dependencies
+ run: |
+ sudo apt-get update -y -qq
+ sudo apt-get -qq -y install curl cmake jq clang g++>=9.3.0
+
+ - name: Compile jac_cell
+ run: |
+ cd unit_test/jac_cell
+ make -j 4
+
+ - name: Run jac_cell
+ run: |
+ cd unit_test/jac_cell
+ ./main3d.gnu.ex inputs_aprox13 > test.out
+
+ - name: Compare to stored output
+ run: |
+ cd unit_test/jac_cell
+ diff -I "^Initializing AMReX" -I "^AMReX" test.out ci-benchmarks/jac_cell_aprox13.out
+
+
+

EOS/helmholtz/actual_eos.H

@@ -154,9 +154,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];
 
@@ -1003,7 +1003,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);
@@ -1118,8 +1118,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

@@ -84,7 +84,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(integrator_rp::MAX_TEMP, amrex::max(state.T, EOSData::mintemp));
+ state.T = std::clamp(state.T, EOSData::mintemp, integrator_rp::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/integrator_type_sdc.H

@@ -19,8 +19,8 @@ void clean_state(const amrex::Real time, BurnT& state, T& int_state)
  if (integrator_rp::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 * integrator_rp::SMALL_X_SAFE);
+ int_state.y(SFS+n) = std::clamp(int_state.y(SFS+n),
+ state.rho * integrator_rp::SMALL_X_SAFE, state.rho);
  }
  }
 

integration/integrator_type_strang.H

@@ -53,7 +53,7 @@ void clean_state (const amrex::Real time, BurnT& state, I& int_state)
 
  if (integrator_rp::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), integrator_rp::SMALL_X_SAFE);
+ int_state.y(n) = std::clamp(int_state.y(n), integrator_rp::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, integrator_rp::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