Merge branch 'development' into ROCK4

AMReX-Astro · Sep 23, 2023 · f9203a5 · f9203a5
2 parents 9a717a9 + d55561c
commit f9203a5
Show file tree

Hide file tree

Showing 40 changed files with 1,075 additions and 466 deletions.
diff --git a/.github/workflows/test_neutrinos.yml b/.github/workflows/test_neutrinos.yml
@@ -0,0 +1,49 @@
+name: test_neutrinos
+
+on: [pull_request]
+jobs:
+  test_neutrinos:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v3
+        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: Build the fextrema tool
+        run: |
+          cd external/amrex/Tools/Plotfile
+          make programs=fextrema -j 2
+
+      - name: Compile
+        run: |
+          cd unit_test/test_neutrino_cooling
+          make realclean
+          make -j 4
+
+      - name: Run test_neutrino_cooling
+        run: |
+          cd unit_test/test_neutrino_cooling
+          ./main3d.gnu.ex inputs
+          ../../external/amrex/Tools/Plotfile/fextrema.gnu.ex test_sneut5 > test.out
+
+      - name: Compare to stored output
+        run: |
+          cd unit_test/test_neutrino_cooling
+          diff test.out ci-benchmarks/test_sneut5.out
diff --git a/CHANGES.md b/CHANGES.md
@@ -3,6 +3,11 @@
   * The file NETWORK_PROPERTIES has been removed from each network,
     as the legacy screening method is no longer used. (#1310)
 
+  * The rprox network was updated and the Jacobian was fixed (#1300)
+
+  * The primordial_chem EOS now can take density and pressure as
+    inputs (#1302)
+
 # 23.07
 
   * The preprocessor variable EXTRA_THERMO has been removed.
@@ -17,7 +22,7 @@
     for estimating the spectral radius of our ODE system (#1222)
 
   * removed SDC_EVOLVE_ENTHALPY -- this was not being used (#1204)
-  
+
 # 23.06
 
   * Added a new Runge-Kutta-Chebyshev integrator (#1191)

diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -63,7 +63,7 @@ function(setup_target_for_microphysics_compilation network_name output_dir)
                              ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/EOS ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/EOS/primordial_chem
                              ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks/primordial_chem ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks 
                              ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/constants 
-                             ${CMAKE_CURRENT_FUNCTION_LIST_DIR}/networks/general_null PARENT_SCOPE)
+                             PARENT_SCOPE)
 
     #we need to have extern_parameters.cpp be available at configure time
     #the script write_probin.py writes this .cpp file so we call it here

diff --git a/EOS/breakout/actual_eos.H b/EOS/breakout/actual_eos.H
@@ -123,7 +123,7 @@ void actual_eos (I input, T& state)
                 state.gam1 = gamma_const;
 
                 // sound speed
-                state.cs = sqrt(gamma_const * poverrho);
+                state.cs = std::sqrt(gamma_const * poverrho);
 
                 state.dpdr_e = poverrho;
                 state.dpde = (gamma_const - 1.0_rt) * state.rho;

diff --git a/integration/VODE/vode_dvhin.H b/integration/VODE/vode_dvhin.H
@@ -44,24 +44,6 @@ void dvhin (BurnT& state, DvodeT& vstate, Real& H0, int& NITER, int& IER)
     // Set an upper bound on h based on vstate.tout-vstate.t and the initial Y and YDOT.
     Real HUB = PT1 * TDIST;
 
-#ifdef TRUE_SDC
-    for (int i = 1; i <= int_neqs; ++i) {
-        Real atol;
-        if (i == 1) {
-            atol = vstate.atol_dens;
-        } else if (i == NumSpec+2) {
-            atol = vstate.atol_enuc;
-        } else {
-            atol = vstate.atol_spec;
-        }
-        Real DELYI = PT1 * std::abs(vstate.yh(i,1)) + atol;
-        Real AFI = std::abs(vstate.yh(i,2));
-        if (AFI * HUB > DELYI) {
-            HUB = DELYI / AFI;
-        }
-    }
-
-#else
     for (int i = 1; i <= int_neqs; ++i) {
         Real atol = i <= NumSpec ? vstate.atol_spec : vstate.atol_enuc;
         Real DELYI = PT1 * std::abs(vstate.yh(i,1)) + atol;
@@ -70,7 +52,6 @@ void dvhin (BurnT& state, DvodeT& vstate, Real& H0, int& NITER, int& IER)
             HUB = DELYI / AFI;
         }
     }
-#endif
 
     // Set initial guess for h as geometric mean of upper and lower bounds.
     int iter = 0;

diff --git a/integration/VODE/vode_dvode.H b/integration/VODE/vode_dvode.H
@@ -72,24 +72,12 @@ int dvode (BurnT& state, DvodeT& vstate)
     vstate.NQ = 1;
     vstate.H = 1.0_rt;
 
-#ifdef TRUE_SDC
-    vstate.ewt(1) = vstate.rtol_dens * std::abs(vstate.yh(1,1)) + vstate.atol_dens;
-    vstate.ewt(1) = 1.0_rt / vstate.ewt(1);
-
-    for (int i = 2; i <= NumSpec+1; ++i) {
-        vstate.ewt(i) = vstate.rtol_spec * std::abs(vstate.yh(i,1)) + vstate.atol_spec;
-        vstate.ewt(i) = 1.0_rt / vstate.ewt(i);
-    }
-    vstate.ewt(NumSpec+2) = vstate.rtol_enuc * std::abs(vstate.yh(NumSpec+2,1)) + vstate.atol_enuc;
-    vstate.ewt(NumSpec+2) = 1.0_rt / vstate.ewt(NumSpec+2);
-#else
     for (int i = 1; i <= NumSpec; ++i) {
         vstate.ewt(i) = vstate.rtol_spec * std::abs(vstate.yh(i,1)) + vstate.atol_spec;
         vstate.ewt(i) = 1.0_rt / vstate.ewt(i);
     }
     vstate.ewt(NumSpec+1) = vstate.rtol_enuc * std::abs(vstate.yh(NumSpec+1,1)) + vstate.atol_enuc;
     vstate.ewt(NumSpec+1) = 1.0_rt / vstate.ewt(NumSpec+1);
-#endif
 
     // Call DVHIN to set initial step size H0 to be attempted.
     H0 = 0.0_rt;
@@ -157,24 +145,12 @@ int dvode (BurnT& state, DvodeT& vstate)
 
            }
 
-#ifdef TRUE_SDC
-           vstate.ewt(1) = vstate.rtol_dens * std::abs(vstate.yh(1,1)) + vstate.atol_dens;
-           vstate.ewt(1) = 1.0_rt / vstate.ewt(1);
-
-           for (int i = 2; i <= NumSpec+1; ++i) {
-               vstate.ewt(i) = vstate.rtol_spec * std::abs(vstate.yh(i,1)) + vstate.atol_spec;
-               vstate.ewt(i) = 1.0_rt / vstate.ewt(i);
-           }
-           vstate.ewt(NumSpec+2) = vstate.rtol_enuc * std::abs(vstate.yh(NumSpec+2,1)) + vstate.atol_enuc;
-           vstate.ewt(NumSpec+2) = 1.0_rt / vstate.ewt(NumSpec+2);
-#else
            for (int i = 1; i <= NumSpec; ++i) {
                vstate.ewt(i) = vstate.rtol_spec * std::abs(vstate.yh(i,1)) + vstate.atol_spec;
                vstate.ewt(i) = 1.0_rt / vstate.ewt(i);
            }
            vstate.ewt(NumSpec+1) = vstate.rtol_enuc * std::abs(vstate.yh(NumSpec+1,1)) + vstate.atol_enuc;
            vstate.ewt(NumSpec+1) = 1.0_rt / vstate.ewt(NumSpec+1);
-#endif
 
        }
        else {

diff --git a/integration/VODE/vode_dvstep.H b/integration/VODE/vode_dvstep.H
@@ -149,7 +149,7 @@ int dvstep (BurnT& state, DvodeT& vstate)
         y_save(i) = vstate.y(i);
     }
 #endif
-#ifdef SIMPLIFIED_SDC
+#if defined(SIMPLIFIED_SDC) || defined(TRUE_SDC)
     Real rho_old = state.rho_orig + TOLD * state.ydot_a[SRHO];
     for (int i = 1; i <= int_neqs; ++i) {
         y_save(i) = vstate.y(i);
@@ -224,7 +224,7 @@ int dvstep (BurnT& state, DvodeT& vstate)
 
         bool valid_update = true;
 
-#ifdef SIMPLIFIED_SDC
+#if defined(SIMPLIFIED_SDC) || defined(TRUE_SDC)
 	if (vstate.y(SEINT+1) < 0.0_rt) {
 	    valid_update = false;
 	}
@@ -271,7 +271,7 @@ int dvstep (BurnT& state, DvodeT& vstate)
 
 #endif
 
-#ifdef SIMPLIFIED_SDC
+#if defined(SIMPLIFIED_SDC) || defined(TRUE_SDC)
 
             // these are basically the same checks as with Strang
             // above, except now we are evolving rhoX instead of X.
@@ -303,24 +303,6 @@ int dvstep (BurnT& state, DvodeT& vstate)
 
 #endif
 
-#ifdef TRUE_SDC
-            // as above, we want to make sure that the mass fractions stay bounded, but we are
-            // evolving:
-            //
-            //  y(1) = density
-            //  y(2:2+NumSpec-1) = rho X
-            //  y(NumSpec+2) = rho e
-
-            if (vstate.y(1+i) < -species_failure_tolerance * vstate.y(1)) {
-                valid_update = false;
-            }
-
-            if (vstate.y(1+i) > (1.0_rt + species_failure_tolerance) * vstate.y(1)) {
-                valid_update = false;
-            }
-#endif
-
-
         }
 
         // The corrector has converged (NFLAG = 0). The local error test is

diff --git a/integration/integrator_data.H b/integration/integrator_data.H
@@ -14,7 +14,7 @@ const int INT_NEQS = SVAR_EVOLVE;
 #endif
 
 #ifdef TRUE_SDC
-const int INT_NEQS = NumSpec + 2;
+const int INT_NEQS = NumSpec + 1;
 #endif
 
 
@@ -56,7 +56,7 @@ constexpr int integrator_neqs ()
 #endif
 
 #ifdef TRUE_SDC
-    int_neqs = NumSpec + 2;
+    int_neqs = NumSpec + 1;
 #endif
 
     return int_neqs;

diff --git a/interfaces/burn_type.H b/interfaces/burn_type.H
@@ -131,12 +131,6 @@ struct burn_t
   // dx is useful for estimating timescales for equilibriation
   Real dx;
 
-#ifdef TRUE_SDC
-  Real E_var;
-  Real mom[3];
-  Real f_source[NumSpec+2];
-#endif
-
   Real mu;
   Real mu_e;
   Real y_e;

diff --git a/networks/CNO_extras/actual_rhs.H b/networks/CNO_extras/actual_rhs.H
@@ -794,8 +794,8 @@ void actual_rhs (burn_t& state, Array1D<Real, 1, neqs>& ydot)
     // Get the thermal neutrino losses
 
     Real sneut, dsneutdt, dsneutdd, snuda, snudz;
-
-    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
+    constexpr int do_derivatives{0};
+    sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
 
     // Append the energy equation (this is erg/g/s)
 
@@ -1197,7 +1197,8 @@ void actual_jac(const burn_t& state, MatrixType& jac)
     // Account for the thermal neutrino losses
 
     Real sneut, dsneutdt, dsneutdd, snuda, snudz;
-    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
+    constexpr int do_derivatives{1};
+    sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
 
     for (int j = 1; j <= NumSpec; ++j) {
        Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz);

diff --git a/networks/ECSN/actual_rhs.H b/networks/ECSN/actual_rhs.H
@@ -412,8 +412,8 @@ void actual_rhs (burn_t& state, Array1D<Real, 1, neqs>& ydot)
     // Get the thermal neutrino losses
 
     Real sneut, dsneutdt, dsneutdd, dsnuda, dsnudz;
-
-    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, dsnuda, dsnudz);
+    constexpr int do_derivatives{0};
+    sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, dsnuda, dsnudz);
 
     // Append the energy equation (this is erg/g/s)
 
@@ -629,7 +629,8 @@ void actual_jac(const burn_t& state, MatrixType& jac)
     // Account for the thermal neutrino losses
 
     Real sneut, dsneutdt, dsneutdd, dsnuda, dsnudz;
-    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, dsnuda, dsnudz);
+    constexpr int do_derivatives{1};
+    sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, dsnuda, dsnudz);
 
     for (int j = 1; j <= NumSpec; ++j) {
        Real b1 = (-state.abar * state.abar * dsnuda + (zion[j-1] - state.zbar) * state.abar * dsnudz);

diff --git a/networks/ase/actual_rhs.H b/networks/ase/actual_rhs.H
@@ -1093,8 +1093,8 @@ void actual_rhs (burn_t& state, Array1D<Real, 1, neqs>& ydot)
     // Get the thermal neutrino losses
 
     Real sneut, dsneutdt, dsneutdd, snuda, snudz;
-
-    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
+    constexpr int do_derivatives{0};
+    sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
 
     // Append the energy equation (this is erg/g/s)
 
@@ -1646,7 +1646,8 @@ void actual_jac(const burn_t& state, MatrixType& jac)
     // Account for the thermal neutrino losses
 
     Real sneut, dsneutdt, dsneutdd, snuda, snudz;
-    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
+    constexpr int do_derivatives{1};
+    sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
 
     for (int j = 1; j <= NumSpec; ++j) {
        Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz);

diff --git a/networks/general_null/network_header.template b/networks/general_null/network_header.template
@@ -23,9 +23,7 @@
 
 constexpr int NumSpec = @@NSPEC@@;
 
-#define NUM_EXTRA_SPECIES @@NEXTRASPEC@@
-
-constexpr int NumSpecExtra = NUM_EXTRA_SPECIES;
+constexpr int NumSpecExtra = @@NEXTRASPEC@@;
 
 constexpr int NumSpecTotal = NumSpec + NumSpecExtra;
 

diff --git a/networks/ignition_reaclib/C-burn-simple/actual_rhs.H b/networks/ignition_reaclib/C-burn-simple/actual_rhs.H
@@ -212,8 +212,8 @@ void actual_rhs (burn_t& state, Array1D<Real, 1, neqs>& ydot)
     // Get the thermal neutrino losses
 
     Real sneut, dsneutdt, dsneutdd, snuda, snudz;
-
-    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
+    constexpr int do_derivatives{0};
+    sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
 
     // Append the energy equation (this is erg/g/s)
 
@@ -309,7 +309,8 @@ void actual_jac(const burn_t& state, MatrixType& jac)
     // Account for the thermal neutrino losses
 
     Real sneut, dsneutdt, dsneutdd, snuda, snudz;
-    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
+    constexpr int do_derivatives{1};
+    sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
 
     for (int j = 1; j <= NumSpec; ++j) {
        Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz);

diff --git a/networks/ignition_reaclib/URCA-simple/actual_rhs.H b/networks/ignition_reaclib/URCA-simple/actual_rhs.H
@@ -234,8 +234,8 @@ void actual_rhs (burn_t& state, Array1D<Real, 1, neqs>& ydot)
     // Get the thermal neutrino losses
 
     Real sneut, dsneutdt, dsneutdd, snuda, snudz;
-
-    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
+    constexpr int do_derivatives{0};
+    sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
 
     // Append the energy equation (this is erg/g/s)
 
@@ -343,7 +343,8 @@ void actual_jac(const burn_t& state, MatrixType& jac)
     // Account for the thermal neutrino losses
 
     Real sneut, dsneutdt, dsneutdd, snuda, snudz;
-    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
+    constexpr int do_derivatives{1};
+    sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
 
     for (int j = 1; j <= NumSpec; ++j) {
        Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz);

diff --git a/networks/nova/actual_rhs.H b/networks/nova/actual_rhs.H
@@ -487,8 +487,8 @@ void actual_rhs (burn_t& state, Array1D<Real, 1, neqs>& ydot)
     // Get the thermal neutrino losses
 
     Real sneut, dsneutdt, dsneutdd, snuda, snudz;
-
-    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
+    constexpr int do_derivatives{0};
+    sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
 
     // Append the energy equation (this is erg/g/s)
 
@@ -746,7 +746,8 @@ void actual_jac(const burn_t& state, MatrixType& jac)
     // Account for the thermal neutrino losses
 
     Real sneut, dsneutdt, dsneutdd, snuda, snudz;
-    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
+    constexpr int do_derivatives{1};
+    sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
 
     for (int j = 1; j <= NumSpec; ++j) {
        Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz);