Simplify/Fix Adaptivity for Small Errors or no CFL

CHANGELOG.md

@@ -14,6 +14,12 @@ architecture for their system.
 
 ### 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 omitting unavailable terms.
+
 Fixed the loading of ARKStep's default first order explicit method.
 
 Fixed a CMake bug regarding usage of missing "print_warning" macro

doc/shared/RecentChanges.rst

@@ -10,6 +10,13 @@ override this value with the architecture for their system.
 
 **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 omitting unavailable terms.
+
 Fixed the loading of ARKStep's default first order explicit method.
 
 Fixed a CMake bug regarding usage of missing "print_warning" macro

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

@@ -158,6 +158,30 @@ extern "C" {
 #endif
 #endif
 
+/*
+ * -----------------------------------------------------------------
+ * Function : SUNRcopysign
+ * -----------------------------------------------------------------
+ * Usage : sunrealtype z;
+ *         z = SUNRcopysign(x, y);
+ * -----------------------------------------------------------------
+ * SUNRcopysign(x, y) returns x with the sign of y.
+ * -----------------------------------------------------------------
+ */
+
+#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
+
 /*
  * -----------------------------------------------------------------
  * Function : SUNRpowerI

src/arkode/arkode_adapt.c

@@ -95,7 +95,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;
 
   /* Request error-based step size from adaptivity controller */
@@ -120,9 +120,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)
@@ -131,7 +128,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",
@@ -141,34 +138,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
@@ -133,27 +132,35 @@ SUNErrCode SUNAdaptController_EstimateStep_ImExGus(SUNAdaptController C,
   /* 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;
 }

src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c

@@ -62,7 +62,6 @@
 #define DEFAULT_IMPGUS_K1 SUN_RCONST(0.98) /* Implicit Gustafsson parameters */
 #define DEFAULT_IMPGUS_K2 SUN_RCONST(0.95)
 #define DEFAULT_BIAS      SUN_RCONST(1.5)
-#define TINY              SUN_RCONST(1.0e-10)
 
 /* -----------------------------------------------------------------
  * exported functions
@@ -76,12 +75,8 @@ SUNAdaptController SUNAdaptController_Soderlind(SUNContext sunctx)
 {
   SUNFunctionBegin(sunctx);
 
-  SUNAdaptController C;
-  SUNAdaptControllerContent_Soderlind content;
-
   /* Create an empty controller object */
-  C = NULL;
-  C = SUNAdaptController_NewEmpty(sunctx);
+  SUNAdaptController C = SUNAdaptController_NewEmpty(sunctx);
   SUNCheckLastErrNull();
 
   /* Attach operations */
@@ -95,12 +90,9 @@ SUNAdaptController SUNAdaptController_Soderlind(SUNContext sunctx)
   C->ops->space        = SUNAdaptController_Space_Soderlind;
 
   /* Create content */
-  content = NULL;
-  content = (SUNAdaptControllerContent_Soderlind)malloc(sizeof *content);
-  SUNAssertNull(content, SUN_ERR_MALLOC_FAIL);
-
-  /* Attach content */
-  C->content = content;
+  C->content = (SUNAdaptControllerContent_Soderlind)malloc(
+    sizeof(struct _SUNAdaptControllerContent_Soderlind));
+  SUNAssertNull(C->content, SUN_ERR_MALLOC_FAIL);
 
   /* Fill content with default/reset values */
   SUNCheckCallNull(SUNAdaptController_SetDefaults_Soderlind(C));
@@ -169,8 +161,7 @@ SUNAdaptController SUNAdaptController_PI(SUNContext sunctx)
 {
   SUNFunctionBegin(sunctx);
 
-  SUNAdaptController C;
-  C = SUNAdaptController_Soderlind(sunctx);
+  SUNAdaptController C = SUNAdaptController_Soderlind(sunctx);
   SUNCheckLastErrNull();
 
   SUNCheckCallNull(
@@ -203,8 +194,7 @@ SUNAdaptController SUNAdaptController_I(SUNContext sunctx)
 {
   SUNFunctionBegin(sunctx);
 
-  SUNAdaptController C;
-  C = SUNAdaptController_Soderlind(sunctx);
+  SUNAdaptController C = SUNAdaptController_Soderlind(sunctx);
   SUNCheckLastErrNull();
 
   SUNCheckCallNull(SUNAdaptController_SetParams_I(C, DEFAULT_I_K1));
@@ -235,8 +225,7 @@ SUNAdaptController SUNAdaptController_ExpGus(SUNContext sunctx)
 {
   SUNFunctionBegin(sunctx);
 
-  SUNAdaptController C;
-  C = SUNAdaptController_Soderlind(sunctx);
+  SUNAdaptController C = SUNAdaptController_Soderlind(sunctx);
   SUNCheckLastErrNull();
 
   SUNCheckCallNull(SUNAdaptController_SetParams_ExpGus(C, DEFAULT_EXPGUS_K1,
@@ -269,8 +258,7 @@ SUNAdaptController SUNAdaptController_ImpGus(SUNContext sunctx)
 {
   SUNFunctionBegin(sunctx);
 
-  SUNAdaptController C;
-  C = SUNAdaptController_Soderlind(sunctx);
+  SUNAdaptController C = SUNAdaptController_Soderlind(sunctx);
   SUNCheckLastErrNull();
 
   SUNCheckCallNull(SUNAdaptController_SetParams_ImpGus(C, DEFAULT_IMPGUS_K1,
@@ -316,28 +304,36 @@ SUNErrCode SUNAdaptController_EstimateStep_Soderlind(SUNAdaptController C,
   /* order parameter to use */
   const int ord = p + 1;
 
-  /* set usable time-step adaptivity parameters */
-  const sunrealtype k1    = -SODERLIND_K1(C) / ord;
-  const sunrealtype k2    = -SODERLIND_K2(C) / ord;
-  const sunrealtype k3    = -SODERLIND_K3(C) / ord;
-  const sunrealtype k4    = SODERLIND_K4(C);
-  const sunrealtype k5    = SODERLIND_K5(C);
-  const sunrealtype e1    = SUNMAX(SODERLIND_BIAS(C) * dsm, TINY);
-  const sunrealtype e2    = SUNMAX(SODERLIND_EP(C), TINY);
-  const sunrealtype e3    = SUNMAX(SODERLIND_EPP(C), TINY);
-  const sunrealtype hrat  = h / SODERLIND_HP(C);
-  const sunrealtype hrat2 = SODERLIND_HP(C) / SODERLIND_HPP(C);
-
-  /* compute estimated optimal time step size */
-  if (SODERLIND_FIRSTSTEPS(C) < 1) { *hnew = h * SUNRpowerR(e1, k1); }
-  else if (SODERLIND_FIRSTSTEPS(C) < 2)
+  if (SODERLIND_FIRSTSTEPS(C) > 1)
   {
-    *hnew = h * SUNRpowerR(e1, k1) * SUNRpowerR(e2, k2) * SUNRpowerR(hrat, k4);
+    /* After the first 2 steps, there is sufficient history */
+    const sunrealtype k1    = -SODERLIND_K1(C) / ord;
+    const sunrealtype k2    = -SODERLIND_K2(C) / ord;
+    const sunrealtype k3    = -SODERLIND_K3(C) / ord;
+    const sunrealtype k4    = SODERLIND_K4(C);
+    const sunrealtype k5    = SODERLIND_K5(C);
+    const sunrealtype e1    = SODERLIND_BIAS(C) * dsm;
+    const sunrealtype e2    = SODERLIND_EP(C);
+    const sunrealtype e3    = SODERLIND_EPP(C);
+    const sunrealtype hrat  = h / SODERLIND_HP(C);
+    const sunrealtype hrat2 = SODERLIND_HP(C) / SODERLIND_HPP(C);
+
+    *hnew = h * SUNRpowerR(e1, k1) * SUNRpowerR(e2, k2) * SUNRpowerR(e3, k3) *
+            SUNRpowerR(hrat, k4) * SUNRpowerR(hrat2, k5);
   }
   else
   {
-    *hnew = h * SUNRpowerR(e1, k1) * SUNRpowerR(e2, k2) * SUNRpowerR(e3, k3) *
-            SUNRpowerR(hrat, k4) * SUNRpowerR(hrat2, k5);
+    /* Use an I controller on the first two steps */
+    const sunrealtype k = -SUN_RCONST(1.0) / ord;
+    const sunrealtype e = SODERLIND_BIAS(C) * dsm;
+    *hnew               = h * SUNRpowerR(e, k);
+  }
+
+  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 */

src/sundials/sundials_math.c

@@ -36,8 +36,6 @@ sunrealtype SUNRpowerI(sunrealtype base, int exponent)
 
 sunrealtype SUNRpowerR(sunrealtype base, sunrealtype exponent)
 {
-  if (base <= SUN_RCONST(0.0)) { return (SUN_RCONST(0.0)); }
-
 #if defined(SUNDIALS_DOUBLE_PRECISION)
   return (pow(base, exponent));
 #elif defined(SUNDIALS_SINGLE_PRECISION)

test/unit_tests/arkode/C_serial/CMakeLists.txt

@@ -16,6 +16,7 @@
 
 # List of test tuples of the form "name\;args"
 set(ARKODE_unit_tests
+    "ark_test_adapt\;"
     "ark_test_arkstepsetforcing\;1 0"
     "ark_test_arkstepsetforcing\;1 1"
     "ark_test_arkstepsetforcing\;1 2"