Merge branch 'development' into ROCK4

AMReX-Astro · Oct 3, 2023 · b3d7629 · b3d7629
2 parents 4657fe4 + e3e866c
commit b3d7629
Show file tree

Hide file tree

Showing 24 changed files with 63 additions and 61 deletions.
diff --git a/.github/workflows/good_defines.txt b/.github/workflows/good_defines.txt
@@ -18,6 +18,7 @@ RATES
 REACTIONS
 SCREENING
 SCREEN_METHOD
+SDC
 SIMPLIFIED_SDC
 STRANG
 TRUE_SDC

diff --git a/CHANGES.md b/CHANGES.md
@@ -1,3 +1,16 @@
+# 23.10
+
+  * The simplified-SDC and true-SDC code paths for integration
+    have been merged (#1338,  #1340, #1341).
+
+  * All pynucastro networks have been updated with the latest
+    version of pynucastro (2.1.0) (#1342)
+
+  * The neutrino cooling terms now use templating on derivatives
+    (#1329)
+
+  * `NUM_EXTRA_SPECIES` was removed (#1321)
+
 # 23.09
 
   * The file NETWORK_PROPERTIES has been removed from each network,
@@ -349,7 +362,7 @@
 # 20.08
 
    * Several of the unit tests had separate C++ and Fortran
-     implementions.  These have been unified (#343, #344, #345)
+     implementations.  These have been unified (#343, #344, #345)
 
    * The VBDF integrator was removed (#348)
 

diff --git a/Make.Microphysics_extern b/Make.Microphysics_extern
@@ -1,10 +1,17 @@
 # This is the main include makefile for applications that want to use Microphysics
 # You should set NETWORK_OUTPUT_PATH before including this file
 
+# for common code paths between simplified-SDC and true-SDC, we'll
+# use the SDC define
+
 ifeq ($(USE_SIMPLIFIED_SDC), TRUE)
   DEFINES += -DSIMPLIFIED_SDC
+  DEFINES += -DSDC
+  USE_ALL_SDC := TRUE
 else ifeq ($(USE_TRUE_SDC), TRUE)
   DEFINES += -DTRUE_SDC
+  DEFINES += -DSDC
+  USE_ALL_SDC := TRUE
 else
   DEFINES += -DSTRANG
 endif

diff --git a/integration/BackwardEuler/Make.package b/integration/BackwardEuler/Make.package
@@ -1,9 +1,7 @@
-ifeq ($(USE_SIMPLIFIED_SDC), TRUE)
+ifeq ($(USE_ALL_SDC), TRUE)
   CEXE_headers += actual_integrator_simplified_sdc.H
 else
-  ifneq ($(USE_TRUE_SDC), TRUE)
-    CEXE_headers += actual_integrator.H
-  endif
+  CEXE_headers += actual_integrator.H
 endif
 
 CEXE_headers += be_integrator.H

diff --git a/integration/BackwardEuler/be_integrator.H b/integration/BackwardEuler/be_integrator.H
@@ -13,7 +13,7 @@
 #ifdef STRANG
 #include <integrator_rhs_strang.H>
 #endif
-#ifdef SIMPLIFIED_SDC
+#ifdef SDC
 #include <integrator_rhs_simplified_sdc.H>
 #endif
 #include <integrator_data.H>
@@ -173,7 +173,7 @@ int be_integrator (BurnT& state, BeT& be)
 
     be.n_rhs += 1;
 
-    Real dt_sub = initial_dt(state, be, ydot);
+    Real dt_sub = initial_react_dt(state, be, ydot);
 
     // main timestepping loop
 

diff --git a/integration/BackwardEuler/be_type.H b/integration/BackwardEuler/be_type.H
@@ -10,7 +10,7 @@
 #ifdef STRANG
 #include <integrator_type_strang.H>
 #endif
-#ifdef SIMPLIFIED_SDC
+#ifdef SDC
 #include <integrator_type_simplified_sdc.H>
 #endif
 #include <network.H>

diff --git a/integration/Make.package b/integration/Make.package
@@ -16,23 +16,19 @@ CEXE_headers += integrator.H
 CEXE_headers += integrator_data.H
 CEXE_headers += integrator_type.H
 
-ifeq ($(USE_SIMPLIFIED_SDC), TRUE)
+ifeq ($(USE_ALL_SDC), TRUE)
   CEXE_headers += integrator_rhs_simplified_sdc.H
   CEXE_headers += integrator_type_simplified_sdc.H
 else
-  ifneq ($(USE_TRUE_SDC), TRUE)
-    CEXE_headers += integrator_type_strang.H
-    CEXE_headers += integrator_rhs_strang.H
-  endif
+  CEXE_headers += integrator_type_strang.H
+  CEXE_headers += integrator_rhs_strang.H
 endif
 
 ifeq ($(USE_NSE_TABLE), TRUE)
-  ifeq ($(USE_SIMPLIFIED_SDC), TRUE)
+  ifeq ($(USE_ALL_SDC), TRUE)
     CEXE_headers += nse_update_simplified_sdc.H
   else
-    ifneq ($(USE_TRUE_SDC), TRUE)
-      CEXE_headers += nse_update_strang.H
-    endif
+    CEXE_headers += nse_update_strang.H
   endif
 endif
 

diff --git a/integration/RKC/Make.package b/integration/RKC/Make.package
@@ -1,9 +1,7 @@
-ifeq ($(USE_SIMPLIFIED_SDC), TRUE)
+ifeq ($(USE_ALL_SDC), TRUE)
   CEXE_headers += actual_integrator_simplified_sdc.H
 else
-  ifneq ($(USE_TRUE_SDC), TRUE)
-    CEXE_headers += actual_integrator.H
-  endif
+  CEXE_headers += actual_integrator.H
 endif
 
 CEXE_headers += rkc_type.H

diff --git a/integration/RKC/rkc.H b/integration/RKC/rkc.H
@@ -8,7 +8,7 @@
 #include <integrator_type_strang.H>
 #include <integrator_rhs_strang.H>
 #endif
-#ifdef SIMPLIFIED_SDC
+#ifdef SDC
 #include <integrator_type_simplified_sdc.H>
 #include <integrator_rhs_simplified_sdc.H>
 #endif

diff --git a/integration/RKC/rkc_type.H b/integration/RKC/rkc_type.H
@@ -94,7 +94,7 @@ struct rkc_t {
 
 };
 
-#ifdef SIMPLIFIED_SDC
+#ifdef SDC
 #include <integrator_type_simplified_sdc.H>
 #endif
 #ifdef STRANG

diff --git a/integration/VODE/Make.package b/integration/VODE/Make.package
@@ -1,9 +1,7 @@
-ifeq ($(USE_SIMPLIFIED_SDC), TRUE)
+ifeq ($(USE_ALL_SDC), TRUE)
   CEXE_headers += actual_integrator_simplified_sdc.H
 else
-  ifneq ($(USE_TRUE_SDC), TRUE)
-    CEXE_headers += actual_integrator.H
-  endif
+  CEXE_headers += actual_integrator.H
 endif
 
 # by default we do not enable Jacobian caching on GPUs to save memory

diff --git a/integration/VODE/vode_dvhin.H b/integration/VODE/vode_dvhin.H
@@ -4,12 +4,9 @@
 #ifdef STRANG
 #include <integrator_rhs_strang.H>
 #endif
-#ifdef SIMPLIFIED_SDC
+#ifdef SDC
 #include <integrator_rhs_simplified_sdc.H>
 #endif
-#ifdef TRUE_SDC
-#include <vode_rhs_true_sdc.H>
-#endif
 
 template <typename BurnT, typename DvodeT>
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE

diff --git a/integration/VODE/vode_dvjac.H b/integration/VODE/vode_dvjac.H
@@ -8,12 +8,9 @@
 #ifdef STRANG
 #include <integrator_rhs_strang.H>
 #endif
-#ifdef SIMPLIFIED_SDC
+#ifdef SDC
 #include <integrator_rhs_simplified_sdc.H>
 #endif
-#ifdef TRUE_SDC
-#include <vode_rhs_true_sdc.H>
-#endif
 
 template <typename IArray1D, typename BurnT, typename DvodeT>
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE

diff --git a/integration/VODE/vode_dvode.H b/integration/VODE/vode_dvode.H
@@ -7,12 +7,9 @@
 #ifdef STRANG
 #include <integrator_rhs_strang.H>
 #endif
-#ifdef SIMPLIFIED_SDC
+#ifdef SDC
 #include <integrator_rhs_simplified_sdc.H>
 #endif
-#ifdef TRUE_SDC
-#include <vode_rhs_true_sdc.H>
-#endif
 #ifdef NSE_TABLE
 #include <nse_table_check.H>
 #endif
@@ -236,7 +233,7 @@ int dvode (BurnT& state, DvodeT& vstate)
 #ifdef STRANG
            update_thermodynamics(state, vstate);
 #endif
-#ifdef SIMPLIFIED_SDC
+#ifdef SDC
            int_to_burn(vstate.tn, vstate, state);
 #endif
 

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
-#if defined(SIMPLIFIED_SDC) || defined(TRUE_SDC)
+#ifdef 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;
 
-#if defined(SIMPLIFIED_SDC) || defined(TRUE_SDC)
+#ifdef SDC
 	if (vstate.y(SEINT+1) < 0.0_rt) {
 	    valid_update = false;
 	}
@@ -259,19 +259,21 @@ int dvstep (BurnT& state, DvodeT& vstate)
             // or greater than one (within a tolerance).
             if (vstate.y(i) < -species_failure_tolerance) {
                 valid_update = false;
+                break;
             }
 
             // Don't enforce the condition below if
             // vstate.y contains number densities
             if (!use_number_densities) {
                 if (vstate.y(i) > 1.0_rt + species_failure_tolerance) {
                     valid_update = false;
+                    break;
                 }
             }
 
 #endif
 
-#if defined(SIMPLIFIED_SDC) || defined(TRUE_SDC)
+#ifdef SDC
 
             // these are basically the same checks as with Strang
             // above, except now we are evolving rhoX instead of X.
@@ -295,10 +297,12 @@ int dvstep (BurnT& state, DvodeT& vstate)
 
             if (vstate.y(SFS+i) < -species_failure_tolerance * rho_current) {
                 valid_update = false;
+                break;
             }
 
             if (vstate.y(SFS+i) > (1.0_rt + species_failure_tolerance) * rho_current) {
                 valid_update = false;
+                break;
             }
 
 #endif

diff --git a/integration/VODE/vode_type.H b/integration/VODE/vode_type.H
@@ -238,7 +238,7 @@ void print_state(dvode_t<int_neqs>& dvode_state)
 }
 #endif
 
-#ifdef SIMPLIFIED_SDC
+#ifdef SDC
 #include <integrator_type_simplified_sdc.H>
 #endif
 #ifdef STRANG

diff --git a/integration/integrator.H b/integration/integrator.H
@@ -1,7 +1,7 @@
 #ifndef INTEGRATOR_H
 #define INTEGRATOR_H
 
-#ifdef SIMPLIFIED_SDC
+#ifdef SDC
 #include <actual_integrator_simplified_sdc.H>
 #else
 #include <actual_integrator.H>

diff --git a/integration/utils/circle_theorem.H b/integration/utils/circle_theorem.H
@@ -6,7 +6,7 @@
 #include <integrator_type_strang.H>
 #include <integrator_rhs_strang.H>
 #endif
-#ifdef SIMPLIFIED_SDC
+#ifdef SDC
 #include <integrator_type_simplified_sdc.H>
 #include <integrator_rhs_simplified_sdc.H>
 #endif

diff --git a/integration/utils/initial_timestep.H b/integration/utils/initial_timestep.H
@@ -12,8 +12,8 @@
 
 template <typename BurnT, typename IntT>
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-Real initial_dt (BurnT& burn_state, IntT& int_state,
-                 Array1D<Real, 1, INT_NEQS>& ydot)
+Real initial_react_dt (BurnT& burn_state, IntT& int_state,
+                     Array1D<Real, 1, INT_NEQS>& ydot)
 {
     // this is a version of the initial timestep estimator from VODE
     // we come in with int_state..y() storing the initial solution

diff --git a/integration/utils/nonaka_plot.H b/integration/utils/nonaka_plot.H
@@ -25,7 +25,7 @@ void nonaka_init() {
     std::stringstream ss;
     ss << std::setw(FIELD_WIDTH) << "time" << " ";
 
-#ifdef SIMPLIFIED_SDC
+#ifdef SDC
     ss << std::setw(FIELD_WIDTH) << "sdc-iter" << " ";
 #endif
 
@@ -80,7 +80,7 @@ void nonaka_rhs(const Real time, const BurnT& state, const YdotNetArray1D& ydot_
 
         nf << std::setw(FIELD_WIDTH) << simulation_time << " ";
 
-#ifdef SIMPLIFIED_SDC
+#ifdef SDC
         nf << std::setw(FIELD_WIDTH) << state.sdc_iter << " ";
 #endif
 

diff --git a/interfaces/burner.H b/interfaces/burner.H
@@ -16,7 +16,7 @@
 #ifdef STRANG
 #include <nse_update_strang.H>
 #endif
-#ifdef SIMPLIFIED_SDC
+#ifdef SDC
 #include <nse_update_simplified_sdc.H>
 #endif
 #endif
@@ -26,16 +26,14 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
 void burner (BurnT& state, Real dt)
 {
 
-#ifndef TRUE_SDC
-
 #ifdef NSE
 
     bool nse_check = in_nse(state);
 
     if (nse_check) {
 
         // call the table
-#ifdef SIMPLIFIED_SDC
+#ifdef SDC
         sdc_nse_burn(state, dt);
 #else
         nse_burn(state, dt);
@@ -45,7 +43,7 @@ void burner (BurnT& state, Real dt)
         // burn as usual
         integrator(state, dt);
 
-#if !defined(SIMPLIFIED_SDC) && defined(AUX_THERMO)
+#if !defined(SDC) && defined(AUX_THERMO)
         // update the aux from the new X's this is not needed with the
         // SDC method, since we took care of that internally.
         set_aux_comp_from_X(state);
@@ -74,7 +72,7 @@ void burner (BurnT& state, Real dt)
 
             // This will append to state.e the amount additional
             // energy released from adjusting to the new NSE state
-#ifdef SIMPLIFIED_SDC
+#ifdef SDC
             sdc_nse_burn(state, dt_remaining);
 #else
             nse_burn(state, dt_remaining);
@@ -94,8 +92,6 @@ void burner (BurnT& state, Real dt)
     integrator(state, dt);
 #endif
 
-#endif
-
 }
 
 #endif
diff --git a/util/gcem/README.md b/util/gcem/README.md
@@ -129,7 +129,7 @@ constexpr
 return_t<T>
 erf(const T x) noexcept;
 ```
-A set of internal templated `constexpr` functions will implement a continued fraction expansion and return a value of type `return_t<T>`. The output type ('`return_t<T>`') is generally determined by the input type, e.g., `int`, `float`, `double`, `long double`, etc.; when `T` is an intergral type, the output will be upgraded to `return_t<T> = double`, otherwise `return_t<T> = T`. For types not covered by `std::is_integral`, recasts should be used.
+A set of internal templated `constexpr` functions will implement a continued fraction expansion and return a value of type `return_t<T>`. The output type ('`return_t<T>`') is generally determined by the input type, e.g., `int`, `float`, `double`, `long double`, etc.; when `T` is an integral type, the output will be upgraded to `return_t<T> = double`, otherwise `return_t<T> = T`. For types not covered by `std::is_integral`, recasts should be used.
 
 ## Examples
 

diff --git a/util/gcem/docs/source/index.rst b/util/gcem/docs/source/index.rst
@@ -54,7 +54,7 @@ GCE-Math functions are written as C++ templates with ``constexpr`` specifiers. F
     return_t<T>
     erf(const T x) noexcept;
 
-A set of internal templated ``constexpr`` functions will implement a continued fraction expansion and return a value of type ``return_t<T>``. The output type ('``return_t<T>``') is generally determined by the input type, e.g., ``int``, ``float``, ``double``, ``long double``, etc; when ``T`` is an intergral type, the output will be upgraded to ``return_t<T> = double``, otherwise ``return_t<T> = T``. For types not covered by ``std::is_integral``, recasts should be used.
+A set of internal templated ``constexpr`` functions will implement a continued fraction expansion and return a value of type ``return_t<T>``. The output type ('``return_t<T>``') is generally determined by the input type, e.g., ``int``, ``float``, ``double``, ``long double``, etc; when ``T`` is an integral type, the output will be upgraded to ``return_t<T> = double``, otherwise ``return_t<T> = T``. For types not covered by ``std::is_integral``, recasts should be used.
 
 
 Contents