get the BackwardEuler solver working with simplified-SDC (#1042)

integration/BackwardEuler/Make.package

@@ -1,4 +1,10 @@
-CEXE_headers += actual_integrator.H
+ifeq ($(USE_SIMPLIFIED_SDC), TRUE)
+  CEXE_headers += actual_integrator_simplified_sdc.H
+else
+  ifneq ($(USE_TRUE_SDC), TRUE)
+    CEXE_headers += actual_integrator.H
+  endif
+endif
+
 CEXE_headers += be_integrator.H
 CEXE_headers += be_type.H
-CEXE_headers += be_type_strang.H

integration/BackwardEuler/actual_integrator_simplified_sdc.H

@@ -0,0 +1,181 @@
+#ifndef actual_integrator_H
+#define actual_integrator_H
+
+#include <iomanip>
+
+#include <be_type.H>
+#include <be_integrator.H>
+#include <extern_parameters.H>
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+void actual_integrator (burn_t& state, const Real dt)
+{
+
+    be_t be;
+
+    // set the Jacobian type
+    be.jacobian_type = jacobian;
+
+    // Start off by assuming a successful burn.
+
+    state.success = true;
+
+    // Initialize the integration time.
+
+    be.t = 0.0;
+    be.tout = dt;
+
+    // Fill in the initial integration state.
+
+    burn_to_int(state, be);
+
+    // Save the initial composition and temperature for our later diagnostics.
+
+#ifndef AMREX_USE_GPU
+    Real xn_in[NumSpec];
+    for (int n = 0; n < NumSpec; ++n) {
+        xn_in[n] = state.y[SFS+n] / state.y[SRHO];
+    }
+    // we are assuming that the temperature was valid on input
+    Real T_in = state.T;
+#ifdef AUX_THERMO
+    Real aux_in[NumAux];
+    for (int n = 0; n < NumAux; ++n) {
+        aux_in[n] = state.y[SFX+n] / state.y[SRHO];
+    }
+#endif
+    Real rhoe_in = state.y[SEINT];
+#endif
+
+
+    // Set the tolerances.
+
+    Real sdc_tol_fac = std::pow(sdc_burn_tol_factor, state.num_sdc_iters - state.sdc_iter - 1);
+
+    // we use 1-based indexing inside of BackwardEuler, so we need to shift the
+    // indices SRHO, SFS, etc by 1
+
+    Real sdc_min_density = amrex::min(state.rho, state.rho_orig + state.ydot_a[SRHO] * dt);
+
+    be.atol_enuc = sdc_min_density * atol_enuc * sdc_tol_fac;
+    be.rtol_enuc = rtol_enuc * sdc_tol_fac;
+
+    // Note: we define the input atol for species to refer only to the
+    // mass fraction part, and we multiply by a representative density
+    // so that atol becomes an absolutely tolerance on (rho X)
+
+    be.atol_spec = sdc_min_density * atol_spec * sdc_tol_fac;
+    be.rtol_spec = rtol_spec * sdc_tol_fac;
+
+    // Call the integration routine.
+
+    int istate = be_integrator(state, be);
+
+    // Get the number of RHS and Jacobian evaluations.
+
+    state.n_rhs = be.n_rhs;
+    state.n_jac = be.n_jac;
+    state.n_step = be.n_step;
+
+    // Copy the integration data back to the burn state.
+    // This will also update the aux state from X if we are using NSE
+
+    int_to_burn(be.t, be, state);
+
+    // we only evolved (rho e), not (rho E), so we need to update the
+    // total energy now to ensure we are conservative
+
+    Real rho_Sdot = 0.0_rt;
+    if (state.time > 0) {
+        rho_Sdot = (be.y(SEINT+1) - state.rhoe_orig) / state.time - state.ydot_a[SEINT];
+    }
+
+    state.y[SEDEN] += state.time * (state.ydot_a[SEDEN] + rho_Sdot);
+
+    // also momentum
+
+    state.y[SMX] += state.time * state.ydot_a[SMX];
+    state.y[SMY] += state.time * state.ydot_a[SMY];
+    state.y[SMZ] += state.time * state.ydot_a[SMZ];
+
+    // normalize the abundances on exit.  We'll assume that the driver
+    // calling this is making use of the conserved state (state.y[]),
+    // so that is what will be normalized.
+
+    normalize_abundances_sdc_burn(state);
+
+    if (istate < 0) {
+        state.success = false;
+    }
+
+
+#ifndef AMREX_USE_GPU
+    if (burner_verbose) {
+        // Print out some integration statistics, if desired.
+        std::cout <<  "integration summary: " << std::endl;
+        std::cout <<  "dens: " << state.rho << " temp: " << state.T << std::endl;
+        std::cout <<  "energy released: " << state.e << std::endl;
+        std::cout <<  "number of steps taken: " << state.n_step << std::endl;
+        std::cout <<  "number of f evaluations: " << state.n_rhs << std::endl;
+    }
+#endif
+
+    // If we failed, print out the current state of the integration.
+
+    if (!state.success) {
+#ifndef AMREX_USE_GPU
+        std::cout << Font::Bold << FGColor::Red << "[ERROR] integration failed in net" << ResetDisplay << std::endl;
+        std::cout << "istate = " << istate << std::endl;
+        std::cout << "zone = (" << state.i << ", " << state.j << ", " << state.k << ")" << std::endl;
+        std::cout << "time = " << be.t << 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) << T_in << std::endl;
+        std::cout << "rhoe start = " << std::setprecision(16) << rhoe_in << std::endl;
+        std::cout << "xn start = ";
+        for (int n = 0; n < NumSpec; ++n) {
+            std::cout << std::setprecision(16) << xn_in[n] << " ";
+        }
+        std::cout << std::endl;
+#ifdef AUX_THERMO
+        std::cout << "aux start = ";
+        for (int n = 0; n < NumAux; ++n) {
+            std::cout << std::setprecision(16) << aux_in[n] << " ";
+        }
+        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 << "xn current = ";
+        for (int n = 0; n < NumSpec; ++n) {
+            std::cout << std::setprecision(16) << 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::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 X_k) = ";
+        for (int n = 0; n < NumSpec; n++) {
+            std::cout << std::setprecision(16) << 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::endl;
+#endif
+#endif
+    }
+
+
+}
+
+#endif

integration/BackwardEuler/be_integrator.H

@@ -2,14 +2,18 @@
 #define BE_INTEGRATOR_H
 
 #include <be_type.H>
-#include <integrator_rhs_strang.H>
 #include <network.H>
 #include <actual_network.H>
 #include <actual_rhs.H>
 #include <burn_type.H>
 #include <linpack.H>
 #include <numerical_jacobian.H>
-
+#ifdef STRANG
+#include <integrator_rhs_strang.H>
+#endif
+#ifdef SIMPLIFIED_SDC
+#include <integrator_rhs_simplified_sdc.H>
+#endif
 
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
 Real initial_dt (burn_t& state, be_t& be,
@@ -48,9 +52,9 @@ Real initial_dt (burn_t& state, be_t& be,
        // routine
 
        for (int ii = 1; ii <= NumSpec; ii++) {
-           ewt(ii) = rtol_spec * std::abs(y_old(ii)) + atol_spec;
+           ewt(ii) = be.rtol_spec * std::abs(y_old(ii)) + be.atol_spec;
        }
-       ewt(net_ienuc) = rtol_enuc * std::abs(y_old(net_ienuc)) + atol_enuc;
+       ewt(net_ienuc) = be.rtol_enuc * std::abs(y_old(net_ienuc)) + be.atol_enuc;
 
        // Construct the trial point.
 
@@ -143,13 +147,15 @@ int single_step (burn_t& state, be_t& be, const Real dt)
 
         if (be.jacobian_type == 1) {
             jac(be.t, state, be, be.jac);
-
         } else {
             jac_info_t jac_info;
             jac_info.h = dt;
             numerical_jac(state, jac_info, be.jac);
+            be.n_rhs += (NumSpec+1);
         }
 
+        be.n_jac++;
+
         if (!integrate_energy) {
             for (int j = 1; j <= BE_NEQS; j++) {
                 be.jac(net_ienuc, j) = 0.0_rt;
@@ -242,6 +248,7 @@ int be_integrator (burn_t& state, be_t& be)
 {
 
     be.n_rhs = 0;
+    be.n_jac = 0;
     be.n_step = 0;
 
     int ierr;
@@ -304,9 +311,9 @@ int be_integrator (burn_t& state, be_t& be)
 
         Array1D<Real, 1, BE_NEQS> w;
         for (int n = 1; n <= NumSpec; n++) {
-            w(n) = 1.0_rt / (rtol_spec * std::abs(y_fine(n)) + atol_spec);
+            w(n) = 1.0_rt / (be.rtol_spec * std::abs(y_fine(n)) + be.atol_spec);
         }
-        w(net_ienuc) = 1.0_rt / (rtol_enuc * std::abs(y_fine(net_ienuc)) + atol_enuc);
+        w(net_ienuc) = 1.0_rt / (be.rtol_enuc * std::abs(y_fine(net_ienuc)) + be.atol_enuc);
 
         // now look for w |y_fine - y_coarse| < 1
 

integration/BackwardEuler/be_type.H

@@ -9,7 +9,9 @@
 #include <integrator_data.H>
 #ifdef STRANG
 #include <integrator_type_strang.H>
-#include <integrator_rhs_strang.H>
+#endif
+#ifdef SIMPLIFIED_SDC
+#include <integrator_type_simplified_sdc.H>
 #endif
 #include <network.H>
 
@@ -50,7 +52,6 @@ struct be_t {
     JacNetArray2D jac;
 
     short jacobian_type;
-
 };
 
 #endif

integration/Make.package

@@ -14,13 +14,15 @@ endif
 
 CEXE_headers += integrator.H
 CEXE_headers += integrator_data.H
-ifneq ($(USE_TRUE_SDC), TRUE)
-  CEXE_headers += integrator_type_strang.H
-  CEXE_headers += integrator_rhs_strang.H
-endif
+
 ifeq ($(USE_SIMPLIFIED_SDC), TRUE)
-  CEXE_headers += integrator_type_simplified_sdc.H
   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
 endif
 
 INCLUDE_LOCATIONS += $(MICROPHYSICS_HOME)/integration/utils

integration/VODE/_parameters

@@ -1,8 +1,5 @@
 @namespace: integrator
 
-# SDC iteration tolerance adjustment factor
-sdc_burn_tol_factor          real         1.d0
-
 # for the step rejection logic on mass fractions, we only consider
 # species that are > X_reject_buffer * atol_spec
 X_reject_buffer              real         1.0

integration/_parameters

@@ -72,3 +72,7 @@ use_number_densities     logical          .false.
 
 # flag for tuning on the subtraction of internal energy
 subtract_internal_energy   logical       .true.
+
+# SDC iteration tolerance adjustment factor
+sdc_burn_tol_factor          real         1.d0
+

integration/integrator_rhs_simplified_sdc.H

@@ -16,7 +16,7 @@
 #include <nonaka_plot.H>
 #endif
 
-// The f_rhs routine provides the right-hand-side for the integrator.
+// The f_rhs routine provides the right-hand-side for the integration solver.
 // This is a generic interface that calls the specific RHS routine in the
 // network you're actually using.
 

integration/integrator_type_simplified_sdc.H

@@ -95,7 +95,7 @@ void clean_state(const Real time, burn_t& state, T& int_state)
         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);
+                                            state.rho * SMALL_X_SAFE);
         }
     }