[WIP] implement damping in the SDC Newton solver

Source/driver/_cpp_parameters

@@ -272,6 +272,15 @@ sdc_extra                    int           0
 # which SDC nonlinear solver to use?  1 = Newton, 2 = VODE, 3 = VODE for first iter
 sdc_solver                   int           1
 
+# do we apply damping to the corrections in the SDC Newton solver?
+sdc_newton_use_damping       int           0
+
+# for the SDC Newton solver with damping, what is the mass fraction
+# threshold below which we ignore species in the damping factor
+# calculation?
+sdc_newton_damping_xcutoff   Real          1.e-4
+
+
 # Do we include geometry source terms due to local unit vectors in non-Cartesian Coord?
 # We currently support R-Z cylinderical 2D (Bernand-Champmartin) and R-THETA spherical 2D
 use_geom_source              bool          1

Source/sdc/Castro_sdc_util.H

@@ -1,6 +1,8 @@
 #ifndef CASTRO_SDC_UTIL_H
 #define CASTRO_SDC_UTIL_H
 
+#include <format>
+
 #include <Castro.H>
 #ifdef REACTIONS
 #ifdef NEW_NETWORK_IMPLEMENTATION
@@ -56,11 +58,17 @@ sdc_solve(const Real dt_m,
         // We are going to assume we already have a good guess
         // for the solve in U_new and just pass the solve onto
         // the main Newton solve
-        sdc_newton_subdivide(dt_m, U_old, U_new, C, sdc_iteration, err_out, ierr);
+        std::vector<amrex::Real> err_history{};
+        sdc_newton_subdivide(dt_m, U_old, U_new, C, sdc_iteration, err_history, ierr);
 
         // failing?
         if (ierr != newton::NEWTON_SUCCESS) {
-            amrex::Abort("Newton subcycling failed in sdc_solve");
+            std::cout << std::format("Newton subcycling failed in sdc_solve with ierr = {}\n",
+                                     ierr);
+            for (auto e : err_history) {
+                std::cout << e << std::endl;
+            }
+            amrex::Abort("Newton failure");
         }
     } else if (sdc_solver == VODE_SOLVE) {
         // Use VODE to do the solution
@@ -75,7 +83,8 @@ sdc_solve(const Real dt_m,
 
         // Now U_new is the update that VODE predicts, so we
         // will use that as the initial guess to the Newton solve
-        sdc_newton_subdivide(dt_m, U_old, U_new, C, sdc_iteration, err_out, ierr);
+        std::vector<amrex::Real> err_history{};
+        sdc_newton_subdivide(dt_m, U_old, U_new, C, sdc_iteration, err_history, ierr);
 
         // Failing?
         if (ierr != newton::NEWTON_SUCCESS) {

Source/sdc/sdc_newton_solve.H

@@ -11,13 +11,15 @@ namespace newton {
     constexpr int CONVERGENCE_FAILURE = -2;
     constexpr int BAD_MASS_FRACTIONS = -3;
 
-    constexpr Real species_failure_tolerance = 1.e-2_rt;
+    constexpr amrex::Real species_failure_tolerance = 1.e-2_rt;
 
-    constexpr int MAX_ITER = 100;
-    constexpr int MAX_NSUB = 64;
+    constexpr int MAX_ITER = 50;
+    constexpr int MAX_NSUB = 2048;
 
 };
 
+
+
 #ifdef REACTIONS
 
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
@@ -170,10 +172,53 @@ sdc_newton_solve(const Real dt_m,
 #endif
 
         // on output, f is the correction to the solution (dU)
+
+        // we want to limit the change by some amount k, such that
+        //
+        // 1/k < (U + dU) / U < k
+        //
+        // we'll write this as
+        //
+        // 1/k < 1 + φ < k
+        //
+        // where φ = dU/U
+        //
+        // we want to limit this, so we will introduce a term α, such that
+        //
+        // 1/k < 1 + α φ < k
+        //
+        // if φ > 0, then we need α < (k - 1) / φ
+        //
+        // if φ < 0, then we need α < |1/k - 1| / |φ|
+
+        amrex::Real alpha{1.0_rt};
+        constexpr amrex::Real k = 10.0;
+
+        if (castro::sdc_newton_use_damping) {
+            for (int n = 0; n < NumSpec; ++n) {
+                const amrex::Real rhoX = burn_state.y[SFS+n];
+                if (rhoX < burn_state.rho * castro::sdc_newton_damping_xcutoff) {
+                    continue;
+                }
+                const amrex::Real drhoX = f(n + 1);
+                const amrex::Real phi = drhoX / rhoX;
+
+                if (phi > 0.0_rt) {
+                    alpha = std::min(alpha, (k - 1.0_rt) / phi);
+                } else if (phi < 0.0_rt) {
+                    alpha = std::min(alpha, (1.0_rt / k - 1.0_rt) / phi);
+                }
+            }
+
+            // make sure we apply at least some correction
+            alpha = std::max(alpha, 5.e-2_rt);
+        }
+
         for (int n = 0; n < NumSpec; ++n) {
-            burn_state.y[SFS+n] += f(n + 1);
+            burn_state.y[SFS+n] += alpha * f(n + 1);
+
         }
-        burn_state.y[SEINT] += f(NumSpec+1);
+        burn_state.y[SEINT] += alpha * f(NumSpec+1);
 
         // compute the norm of the weighted error, where the
         // weights are 1/eps_tot
@@ -190,7 +235,7 @@ sdc_newton_solve(const Real dt_m,
                 eps = integrator_rp::rtol_enuc * std::abs(burn_state.y[SEINT]) +
                       integrator_rp::atol_enuc;
             }
-            err_sum += f(n) * f(n) / (eps * eps);
+            err_sum += (alpha * f(n)) * (alpha * f(n)) / (eps * eps);
         }
         err = std::sqrt(err_sum / static_cast<Real>(NumSpec+1));
 
@@ -231,7 +276,7 @@ sdc_newton_subdivide(const Real dt_m,
                      GpuArray<Real, NUM_STATE>& U_new,
                      GpuArray<Real, NUM_STATE> const& C,
                      const int sdc_iteration,
-                     Real& err_out,
+                     std::vector<amrex::Real>& err_out,
                      int& ierr) {
 
     // This is the driver for solving the nonlinear update for the
@@ -282,7 +327,9 @@ sdc_newton_subdivide(const Real dt_m,
                 U_begin[UFS + n] *= U_begin[URHO] / sum_rhoX;
             }
 
-            sdc_newton_solve(dt_sub, U_begin, U_new, C, sdc_iteration, err_out, ierr);
+            amrex::Real err;
+            sdc_newton_solve(dt_sub, U_begin, U_new, C, sdc_iteration, err, ierr);
+            err_out.push_back(err);
 
             // our solve may have resulted in mass fractions outside
             // of [0, 1] -- reject if this is the case