Merge branch 'development' into use_cxx-20

zingale · zingale · commit e42b7fa40d83 · 2026-01-10T09:18:52.000-05:00
diff --git a/EOS/helmholtz/actual_eos.H b/EOS/helmholtz/actual_eos.H
@@ -348,23 +348,26 @@ void apply_electrons (T& state)
 
 
     // electron chemical potential etaele
-    amrex::Real etaele = 0.0e0_rt;
-    amrex::constexpr_for<0, 16>([&] (auto K) {
-        constexpr int k = K;
-        constexpr auto m = electron_table_indexing::map[k];
-
-        etaele += ef[jat+m.dj][iat+m.di][m.n] * wdt[k];
-    });
+    [[maybe_unused]] amrex::Real etaele{};
+    if constexpr (has_eta<T>) {
+        amrex::constexpr_for<0, 16>([&] (auto K) {
+            constexpr int k = K;
+            constexpr auto m = electron_table_indexing::map[k];
 
+            etaele += ef[jat+m.dj][iat+m.di][m.n] * wdt[k];
+        });
+    }
 
     // electron + positron number densities
-    amrex::Real xnefer = 0.0e0_rt;
-    amrex::constexpr_for<0, 16>([&] (auto K) {
-        constexpr int k = K;
-        constexpr auto m = electron_table_indexing::map[k];
+    [[maybe_unused]] amrex::Real xnefer{};
+    if constexpr (has_xne_xnp<T>) {
+        amrex::constexpr_for<0, 16>([&] (auto K) {
+            constexpr int k = K;
+            constexpr auto m = electron_table_indexing::map[k];
 
-        xnefer += xf[jat+m.dj][iat+m.di][m.n] * wdt[k];
-    });
+            xnefer += xf[jat+m.dj][iat+m.di][m.n] * wdt[k];
+        });
+    }
 
     // the desired electron-positron thermodynamic quantities
 
@@ -375,22 +378,38 @@ void apply_electrons (T& state)
     amrex::Real x       = din * din;
     amrex::Real pele    = x * df_d;
     amrex::Real dpepdt  = x * df_dt;
-    [[maybe_unused]] amrex::Real s       = dpepdd/state.y_e - 2.0e0_rt * din * df_d;
-    [[maybe_unused]] amrex::Real dpepda  = -ytot1 * (2.0e0_rt * pele + s * din);
-    [[maybe_unused]] amrex::Real dpepdz  = state.rho*ytot1*(2.0e0_rt * din * df_d  +  s);
+    [[maybe_unused]] amrex::Real dpepda{};
+    [[maybe_unused]] amrex::Real dpepdz{};
+    if constexpr (has_dpdA<T> || has_dpdZ<T>) {
+        amrex::Real s       = dpepdd/state.y_e - 2.0e0_rt * din * df_d;
+        if constexpr (has_dpdA<T>) {
+            dpepda  = -ytot1 * (2.0e0_rt * pele + s * din);
+        }
+        if constexpr (has_dpdZ<T>) {
+            dpepdz  = state.rho*ytot1*(2.0e0_rt * din * df_d  +  s);
+        }
+    }
 
     x            = state.y_e * state.y_e;
     amrex::Real sele    = -df_t * state.y_e;
     amrex::Real dsepdt  = -df_tt * state.y_e;
     amrex::Real dsepdd  = -df_dt * x;
-    [[maybe_unused]] amrex::Real dsepda  = ytot1 * (state.y_e * df_dt * din - sele);
-    [[maybe_unused]] amrex::Real dsepdz  = -ytot1 * (state.y_e * df_dt * state.rho  + df_t);
 
     amrex::Real eele    = state.y_e*free + state.T * sele;
     amrex::Real deepdt  = state.T * dsepdt;
     amrex::Real deepdd  = x * df_d + state.T * dsepdd;
-    [[maybe_unused]] amrex::Real deepda  = -state.y_e * ytot1 * (free +  df_d * din) + state.T * dsepda;
-    [[maybe_unused]] amrex::Real deepdz  = ytot1* (free + state.y_e * df_d * state.rho) + state.T * dsepdz;
+
+    [[maybe_unused]] amrex::Real deepda{};
+    if constexpr (has_dedA<T>) {
+        amrex::Real dsepda = ytot1 * (state.y_e * df_dt * din - sele);
+        deepda = -state.y_e * ytot1 * (free +  df_d * din) + state.T * dsepda;
+    }
+
+    [[maybe_unused]] amrex::Real deepdz{};
+    if constexpr (has_dedZ<T>) {
+        amrex::Real dsepdz = -ytot1 * (state.y_e * df_dt * state.rho  + df_t);
+        deepdz = ytot1* (free + state.y_e * df_d * state.rho) + state.T * dsepdz;
+    }
 
     if constexpr (has_pressure<T>) {
         state.p    = state.p + pele;
@@ -549,14 +568,10 @@ void apply_radiation (T& state)
 
     amrex::Real dpraddd = 0.0e0_rt;
     amrex::Real dpraddt = 4.0e0_rt * prad * tempi;
-    [[maybe_unused]] amrex::Real dpradda = 0.0e0_rt;
-    [[maybe_unused]] amrex::Real dpraddz = 0.0e0_rt;
 
     amrex::Real erad    = 3.0e0_rt * prad*deni;
     amrex::Real deraddd = -erad * deni;
     amrex::Real deraddt = 3.0e0_rt * dpraddt * deni;
-    [[maybe_unused]] amrex::Real deradda = 0.0e0_rt;
-    [[maybe_unused]] amrex::Real deraddz = 0.0e0_rt;
 
     amrex::Real srad    = (prad * deni + erad) * tempi;
     amrex::Real dsraddd = (dpraddd * deni - prad * deni * deni + deraddd) * tempi;
@@ -571,10 +586,10 @@ void apply_radiation (T& state)
         state.dpdr = dpraddd;
         state.dpdT = dpraddt;
         if constexpr (has_dpdA<T>) {
-            state.dpdA = dpradda;
+            state.dpdA = 0.0_rt;
         }
         if constexpr (has_dpdZ<T>) {
-            state.dpdZ = dpraddz;
+            state.dpdZ = 0.0_rt;
         }
     }
 
@@ -583,10 +598,10 @@ void apply_radiation (T& state)
         state.dedr = deraddd;
         state.dedT = deraddt;
         if constexpr (has_dedA<T>) {
-            state.dedA = deradda;
+            state.dedA = 0.0_rt;
         }
         if constexpr (has_dedZ<T>) {
-            state.dedZ = deraddz;
+            state.dedZ = 0.0_rt;
         }
     }
 
diff --git a/integration/integrator_data.H b/integration/integrator_data.H
@@ -1,6 +1,8 @@
 #ifndef INTEGRATOR_DATA_H
 #define INTEGRATOR_DATA_H
 
+#include <limits>
+
 #include <burn_type.H>
 
 // Define the size of the ODE system that VODE will integrate
@@ -12,6 +14,10 @@ constexpr int INT_NEQS = NumSpec + 1;
 // -failure_tolerance <= X <= 1.0 + failure_tolerance).
 constexpr amrex::Real species_failure_tolerance = 1.e-2_rt;
 
+// for outputting the burn state failure, use enough precision
+// to be reproducible (not guaranteed for std::fixed)
+constexpr int OUTDIGITS = std::numeric_limits<amrex::Real>::max_digits10;
+
 #ifdef NSE
 constexpr int MIN_NSE_BAILOUT_STEPS = 10;
 #endif
diff --git a/integration/integrator_setup_sdc.H b/integration/integrator_setup_sdc.H
@@ -219,47 +219,47 @@ void integrator_cleanup (IntegratorT& int_state, BurnT& state,
             }
             std::cout << "zone = (" << state.i << ", " << state.j << ", " << state.k << ")" << std::endl;
             std::cout << "time = " << state.time << std::endl;
-            std::cout << "dt = " << std::setprecision(16) << dt << std::endl;
-            std::cout << "dens start = " << std::setprecision(16) << state.rho_orig << std::endl;
-            std::cout << "temp start = " << std::setprecision(16) << state_save.T_in << std::endl;
-            std::cout << "rhoe start = " << std::setprecision(16) << state_save.rhoe_in << std::endl;
+            std::cout << "dt = " << std::setprecision(OUTDIGITS) << dt << std::endl;
+            std::cout << "dens start = " << std::setprecision(OUTDIGITS) << state.rho_orig << std::endl;
+            std::cout << "temp start = " << std::setprecision(OUTDIGITS) << state_save.T_in << std::endl;
+            std::cout << "rhoe start = " << std::setprecision(OUTDIGITS) << state_save.rhoe_in << std::endl;
             std::cout << "xn start = ";
             for (const auto X : state_save.xn_in) {
-                std::cout << std::setprecision(16) << X << " ";
+                std::cout << std::setprecision(OUTDIGITS) << X << " ";
             }
             std::cout << std::endl;
 #ifdef AUX_THERMO
             std::cout << "aux start = ";
             for (const auto aux : state_save.aux_in) {
-                std::cout << std::setprecision(16) << aux << " ";
+                std::cout << std::setprecision(OUTDIGITS) << aux << " ";
             }
             std::cout << std::endl;
 #endif
-            std::cout << "dens current = " << std::setprecision(16) << state.rho << std::endl;
-            std::cout << "temp current = " << std::setprecision(16) << state.T << std::endl;
+            std::cout << "dens current = " << std::setprecision(OUTDIGITS) << state.rho << std::endl;
+            std::cout << "temp current = " << std::setprecision(OUTDIGITS) << state.T << std::endl;
             std::cout << "xn current = ";
             for (int n = 0; n < NumSpec; ++n) {
-                std::cout << std::setprecision(16) << state.xn[n] << " ";
+                std::cout << std::setprecision(OUTDIGITS) << state.xn[n] << " ";
             }
             std::cout << std::endl;
 #ifdef AUX_THERMO
             std::cout << "aux current = ";
             for (int n = 0; n < NumAux; ++n) {
-                std::cout << std::setprecision(16) << state.aux[n] << " ";
+                std::cout << std::setprecision(OUTDIGITS) << state.aux[n] << " ";
             }
             std::cout << std::endl;
 #endif
-            std::cout << "A(rho) = " << std::setprecision(16) << state.ydot_a[SRHO] << std::endl;
-            std::cout << "A(rho e) = " << std::setprecision(16) << state.ydot_a[SEINT] << std::endl;
+            std::cout << "A(rho) = " << std::setprecision(OUTDIGITS) << state.ydot_a[SRHO] << std::endl;
+            std::cout << "A(rho e) = " << std::setprecision(OUTDIGITS) << state.ydot_a[SEINT] << std::endl;
             std::cout << "A(rho X_k) = ";
             for (int n = 0; n < NumSpec; n++) {
-                std::cout << std::setprecision(16) << state.ydot_a[SFS+n] << " ";
+                std::cout << std::setprecision(OUTDIGITS) << state.ydot_a[SFS+n] << " ";
             }
             std::cout << std::endl;
 #ifdef AUX_THERMO
             std::cout << "A(rho aux_k) = ";
             for (int n = 0; n < NumAux; n++) {
-                std::cout << std::setprecision(16) << state.ydot_a[SFX+n] << " ";
+                std::cout << std::setprecision(OUTDIGITS) << state.ydot_a[SFX+n] << " ";
             }
             std::cout << std::endl;
 #endif
diff --git a/integration/integrator_setup_strang.H b/integration/integrator_setup_strang.H
@@ -190,18 +190,18 @@ void integrator_cleanup (IntegratorT& int_state, BurnT& state,
         }
         std::cout << "zone = (" << state.i << ", " << state.j << ", " << state.k << ")" << std::endl;
         std::cout << "time = " << int_state.t << std::endl;
-        std::cout << "dt = " << std::setprecision(16) << dt << std::endl;
-        std::cout << "temp start = " << std::setprecision(16) << state_save.T_in << std::endl;
+        std::cout << "dt = " << std::setprecision(OUTDIGITS) << dt << std::endl;
+        std::cout << "temp start = " << std::setprecision(OUTDIGITS) << state_save.T_in << std::endl;
         std::cout << "xn start = ";
         for (const double X : state_save.xn_in) {
-            std::cout << std::scientific << std::setprecision(16) << X << " ";
+            std::cout << std::scientific << std::setprecision(OUTDIGITS) << X << " ";
         }
         std::cout << std::endl;
-        std::cout << "dens current = " << std::setprecision(16) << state.rho << std::endl;
-        std::cout << "temp current = " << std::setprecision(16) << state.T << std::endl;
+        std::cout << "dens current = " << std::setprecision(OUTDIGITS) << state.rho << std::endl;
+        std::cout << "temp current = " << std::setprecision(OUTDIGITS) << state.T << std::endl;
         std::cout << "xn current = ";
         for (const double X : state.xn) {
-            std::cout << std::scientific << std::setprecision(16) << X << " ";
+            std::cout << std::scientific << std::setprecision(OUTDIGITS) << X << " ";
         }
         std::cout << std::endl;
         std::cout << "energy generated = " << state.e << std::endl;
