Simplify/Fix Adaptivity for Small Errors or no CFL

CHANGELOG.md

@@ -25,6 +25,12 @@ specific linker flags e.g., MKL.
 
 ### Bug Fixes
 
+Removed error floors from the SUNAdaptController implementations which could
+unnecessarily limit the time size growth, particularly after the first step.
+
+On the first two time steps, the Soderlind controller uses an I controller
+instead of omitting unavailable terms.
+
 Fixed `ARKodeResize` not using the default `hscale` when an argument of `0` was
 provided.
 

doc/shared/RecentChanges.rst

@@ -21,6 +21,13 @@ specific linker flags e.g., MKL.
 
 **Bug Fixes**
 
+Removed error floors from the SUNAdaptController implementations which could
+unnecessarily limit the time size growth, particularly after the first step.
+
+On the first two time steps, the
+:ref:`Soderlind controller <SUNAdaptController.Soderlind>` uses an I controller
+instead of omitting unavailable terms.
+
 Fixed c:func:`ARKodeResize` not using the default ``hscale`` when an argument of
 ``0`` was provided.
 

doc/shared/sunadaptcontroller/SUNAdaptController_ImExGus.rst

@@ -45,9 +45,7 @@ with the form
 
 In the above formulas, the default values of :math:`k_1^E`, :math:`k_2^E`,
 :math:`k_1^I`, and :math:`k_2^I` are 0.367, 0.268, 0.98, and 0.95, respectively,
-and :math:`p` is the global order of the time integration method.  In these
-estimates, a floor of :math:`\varepsilon_* > 10^{-10}` is enforced to avoid
-division-by-zero errors.
+and :math:`p` is the global order of the time integration method.
 
 The SUNAdaptController_ImExGus controller implements both formulas
 :eq:`expGusController` and :eq:`impGusController`, and sets its recommended step

doc/shared/sunadaptcontroller/SUNAdaptController_Soderlind.rst

@@ -27,11 +27,10 @@ and :cite:p:`Sod:06`.  This controller has the form
 
 with default parameter values :math:`k_1 = 1.25`, :math:`k_2 = 0.5`,
 :math:`k_3 = -0.75`, :math:`k_4 = 0.25`, and :math:`k_5 = 0.75`, where
-:math:`p` is the global order of the time integration method.  In this estimate,
-a floor of :math:`\varepsilon_* > 10^{-10}` is enforced to avoid division-by-zero
-errors.  During the first two steps (when :math:`\varepsilon_{n-2}`,
-:math:`\varepsilon_{n-1}`, :math:`h_{n-2}`, and :math:`h_{n-2}` may be unavailable),
-the corresponding terms are merely omitted during estimation of :math:`h'`.
+:math:`p` is the global order of the time integration method.  During the first
+two steps (when :math:`\varepsilon_{n-2}`, :math:`\varepsilon_{n-1}`,
+:math:`h_{n-2}`, and :math:`h_{n-2}` may be unavailable), the I controller
+:math:`h' = h_n \varepsilon_1^{-1/(p+1)}` is used.
 
 The SUNAdaptController_Soderlind controller is implemented as a derived
 SUNAdaptController class, and defines its *content* field as:

include/sundials/sundials_math.h

@@ -160,18 +160,27 @@ extern "C" {
 
 /*
  * -----------------------------------------------------------------
- * Function : SUNRpowerI
+ * Function : SUNRcopysign
  * -----------------------------------------------------------------
- * Usage : int exponent;
- *         sunrealtype base, ans;
- *         ans = SUNRpowerI(base,exponent);
+ * Usage : sunrealtype z;
+ *         z = SUNRcopysign(x, y);
  * -----------------------------------------------------------------
- * SUNRpowerI returns the value of base^exponent, where base is of type
- * sunrealtype and exponent is of type int.
+ * SUNRcopysign(x, y) returns x with the sign of y.
  * -----------------------------------------------------------------
  */
 
-SUNDIALS_EXPORT sunrealtype SUNRpowerI(sunrealtype base, int exponent);
+#ifndef SUNRcopysign
+#if defined(SUNDIALS_DOUBLE_PRECISION)
+#define SUNRcopysign(x, y) (copysign((x), (y)))
+#elif defined(SUNDIALS_SINGLE_PRECISION)
+#define SUNRcopysign(x, y) (copysignf((x), (y)))
+#elif defined(SUNDIALS_EXTENDED_PRECISION)
+#define SUNRcopysign(x, y) (copysignl((x), (y)))
+#else
+#error \
+  "SUNDIALS precision not defined, report to github.com/LLNL/sundials/issues"
+#endif
+#endif
 
 /*
  * -----------------------------------------------------------------
@@ -181,12 +190,36 @@ SUNDIALS_EXPORT sunrealtype SUNRpowerI(sunrealtype base, int exponent);
  *         ans = SUNRpowerR(base,exponent);
  * -----------------------------------------------------------------
  * SUNRpowerR returns the value of base^exponent, where both base and
- * exponent are of type sunrealtype. If base < ZERO, then SUNRpowerR
- * returns ZERO.
+ * exponent are of type sunrealtype.
+ * -----------------------------------------------------------------
+ */
+#ifndef SUNRpowerR
+#if defined(SUNDIALS_DOUBLE_PRECISION)
+#define SUNRpowerR(base, exponent) (pow(base, exponent))
+#elif defined(SUNDIALS_SINGLE_PRECISION)
+#define SUNRpowerR(base, exponent) (powf(base, exponent))
+#elif defined(SUNDIALS_EXTENDED_PRECISION)
+#define SUNRpowerR(base, exponent) (powl(base, exponent))
+#else
+#error \
+  "SUNDIALS precision not defined, report to github.com/LLNL/sundials/issues"
+#endif
+#endif
+
+/*
+ * -----------------------------------------------------------------
+ * Function : SUNRpowerI
+ * -----------------------------------------------------------------
+ * Usage : int exponent;
+ *         sunrealtype base, ans;
+ *         ans = SUNRpowerI(base,exponent);
+ * -----------------------------------------------------------------
+ * SUNRpowerI returns the value of base^exponent, where base is of type
+ * sunrealtype and exponent is of type int.
  * -----------------------------------------------------------------
  */
 
-SUNDIALS_EXPORT sunrealtype SUNRpowerR(sunrealtype base, sunrealtype exponent);
+SUNDIALS_EXPORT sunrealtype SUNRpowerI(sunrealtype base, int exponent);
 
 /*
  * -----------------------------------------------------------------

src/arkode/arkode_adapt.c

@@ -98,7 +98,7 @@ int arkAdapt(ARKodeMem ark_mem, ARKodeHAdaptMem hadapt_mem, N_Vector ycur,
              sunrealtype tcur, sunrealtype hcur, sunrealtype dsm)
 {
   int retval;
-  sunrealtype h_acc, h_cfl, int_dir;
+  sunrealtype h_acc, h_cfl;
   int controller_order;
 
   /* Return with no stepsize adjustment if the controller is NULL */
@@ -130,9 +130,6 @@ int arkAdapt(ARKodeMem ark_mem, ARKodeHAdaptMem hadapt_mem, N_Vector ycur,
     return (ARK_CONTROLLER_ERR);
   }
 
-  /* determine direction of integration */
-  int_dir = hcur / SUNRabs(hcur);
-
   /* Call explicit stability function */
   retval = hadapt_mem->expstab(ycur, tcur, &h_cfl, hadapt_mem->estab_data);
   if (retval != ARK_SUCCESS)
@@ -141,7 +138,7 @@ int arkAdapt(ARKodeMem ark_mem, ARKodeHAdaptMem hadapt_mem, N_Vector ycur,
                     "Error in explicit stability function.");
     return (ARK_ILL_INPUT);
   }
-  if (h_cfl <= ZERO) { h_cfl = SUN_RCONST(1.0e30) * SUNRabs(hcur); }
+  if (h_cfl <= ZERO) { h_cfl = INFINITY; }
 
 #if SUNDIALS_LOGGING_LEVEL >= SUNDIALS_LOGGING_INFO
   SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_INFO, "ARKODE::arkAdapt",
@@ -151,34 +148,39 @@ int arkAdapt(ARKodeMem ark_mem, ARKodeHAdaptMem hadapt_mem, N_Vector ycur,
 
   /* enforce safety factors */
   h_acc *= hadapt_mem->safety;
-  h_cfl *= hadapt_mem->cfl * int_dir;
+  h_cfl *= hadapt_mem->cfl;
 
   /* enforce maximum bound on time step growth */
-  h_acc = int_dir * SUNMIN(SUNRabs(h_acc), SUNRabs(hadapt_mem->etamax * hcur));
+  h_acc = SUNMIN(SUNRabs(h_acc), SUNRabs(hadapt_mem->etamax * hcur));
 
   /* enforce minimum bound time step reduction */
-  h_acc = int_dir * SUNMAX(SUNRabs(h_acc), SUNRabs(hadapt_mem->etamin * hcur));
+  h_acc = SUNMAX(h_acc, SUNRabs(hadapt_mem->etamin * hcur));
 
 #if SUNDIALS_LOGGING_LEVEL >= SUNDIALS_LOGGING_INFO
   SUNLogger_QueueMsg(ARK_LOGGER, SUN_LOGLEVEL_INFO, "ARKODE::arkAdapt",
                      "new-step-after-max-min-bounds",
-                     "h_acc = %" RSYM ", h_cfl = %" RSYM, h_acc, h_cfl);
+                     "h_acc = %" RSYM ", h_cfl = %" RSYM,
+                     SUNRcopysign(h_acc, hcur), h_cfl);
 #endif
 
   /* increment the relevant step counter, set desired step */
-  if (SUNRabs(h_acc) < SUNRabs(h_cfl)) { hadapt_mem->nst_acc++; }
-  else { hadapt_mem->nst_exp++; }
-  h_acc = int_dir * SUNMIN(SUNRabs(h_acc), SUNRabs(h_cfl));
+  if (h_acc <= h_cfl) { hadapt_mem->nst_acc++; }
+  else
+  {
+    hadapt_mem->nst_exp++;
+    h_acc = h_cfl;
+  }
 
   /* enforce adaptivity bounds to retain Jacobian/preconditioner accuracy */
   if (dsm <= ONE)
   {
-    if ((SUNRabs(h_acc) > SUNRabs(hcur * hadapt_mem->lbound * ONEMSM)) &&
-        (SUNRabs(h_acc) < SUNRabs(hcur * hadapt_mem->ubound * ONEPSM)))
+    if ((h_acc > SUNRabs(hcur * hadapt_mem->lbound * ONEMSM)) &&
+        (h_acc < SUNRabs(hcur * hadapt_mem->ubound * ONEPSM)))
     {
       h_acc = hcur;
     }
   }
+  h_acc = SUNRcopysign(h_acc, hcur);
 
   /* set basic value of ark_eta */
   ark_mem->eta = h_acc / hcur;

src/sunadaptcontroller/imexgus/sunadaptcontroller_imexgus.c

@@ -48,7 +48,6 @@
 #define DEFAULT_K1I  SUN_RCONST(0.95)
 #define DEFAULT_K2I  SUN_RCONST(0.95)
 #define DEFAULT_BIAS SUN_RCONST(1.5)
-#define TINY         SUN_RCONST(1.0e-10)
 
 /* -----------------------------------------------------------------
  * exported functions
@@ -128,32 +127,38 @@ SUNErrCode SUNAdaptController_EstimateStep_ImExGus(SUNAdaptController C,
 {
   SUNFunctionBegin(C->sunctx);
 
-  SUNAssert(hnew, SUN_ERR_ARG_CORRUPT);
-
   /* order parameter to use */
   const int ord = p + 1;
 
-  /* set usable time-step adaptivity parameters -- first step */
-  const sunrealtype k = -SUN_RCONST(1.0) / ord;
-  const sunrealtype e = SUNMAX(SACIMEXGUS_BIAS(C) * dsm, TINY);
-
-  /* set usable time-step adaptivity parameters -- subsequent steps */
-  const sunrealtype k1e  = -SACIMEXGUS_K1E(C) / ord;
-  const sunrealtype k2e  = -SACIMEXGUS_K2E(C) / ord;
-  const sunrealtype k1i  = -SACIMEXGUS_K1I(C) / ord;
-  const sunrealtype k2i  = -SACIMEXGUS_K2I(C) / ord;
-  const sunrealtype e1   = SUNMAX(SACIMEXGUS_BIAS(C) * dsm, TINY);
-  const sunrealtype e2   = e1 / SUNMAX(SACIMEXGUS_EP(C), TINY);
-  const sunrealtype hrat = h / SACIMEXGUS_HP(C);
-
   /* compute estimated time step size, modifying the first step formula */
-  if (SACIMEXGUS_FIRSTSTEP(C)) { *hnew = h * SUNRpowerR(e, k); }
+  if (SACIMEXGUS_FIRSTSTEP(C))
+  {
+    /* set usable time-step adaptivity parameters -- first step */
+    const sunrealtype k = -SUN_RCONST(1.0) / ord;
+    const sunrealtype e = SACIMEXGUS_BIAS(C) * dsm;
+    *hnew               = h * SUNRpowerR(e, k);
+  }
   else
   {
+    /* set usable time-step adaptivity parameters -- subsequent steps */
+    const sunrealtype k1e  = -SACIMEXGUS_K1E(C) / ord;
+    const sunrealtype k2e  = -SACIMEXGUS_K2E(C) / ord;
+    const sunrealtype k1i  = -SACIMEXGUS_K1I(C) / ord;
+    const sunrealtype k2i  = -SACIMEXGUS_K2I(C) / ord;
+    const sunrealtype e1   = SACIMEXGUS_BIAS(C) * dsm;
+    const sunrealtype e2   = e1 / SACIMEXGUS_EP(C);
+    const sunrealtype hrat = h / SACIMEXGUS_HP(C);
     *hnew = h * SUNMIN(hrat * SUNRpowerR(e1, k1i) * SUNRpowerR(e2, k2i),
                        SUNRpowerR(e1, k1e) * SUNRpowerR(e2, k2e));
   }
 
+  if (!isfinite(*hnew))
+  {
+    /* hnew can be INFINITY or NAN if multiple e's are 0 or there are overflows.
+     * In that case, make hnew INFINITY with the same sign as h */
+    *hnew = SUNRcopysign(INFINITY, h);
+  }
+
   /* return with success */
   return SUN_SUCCESS;
 }