Simplify/Fix Adaptivity for Small Errors or no CFL

CHANGELOG.md

@@ -27,6 +27,12 @@ Improved the efficiency of default ARKODE methods with the following changes:
 
 ### 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 bug in `ARKodeSetFixedStep` where it could return `ARK_SUCCESS` despite
 an error occurring.
 

doc/shared/RecentChanges.rst

@@ -37,6 +37,13 @@ Improved the efficiency of default ARKODE methods with the following changes:
 
 **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 bug in :c:func:`ARKodeSetFixedStep` where it could return ``ARK_SUCCESS``
 despite an error occurring.
 

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:

examples/arkode/F2003_serial/ark_bruss_f2003.f90

@@ -305,6 +305,13 @@ program main
     stop 1
   end if
 
+  ! TODO(SBR): temporary fix until https://github.com/LLNL/sundials/pull/565 merged
+  ierr = FARKodeSetAdaptivityAdjustment(arkode_mem, 0)
+  if (ierr /= 0) then
+    print *, 'Error in FARKodeSetJacFn'
+    stop 1
+  end if
+
   sunls => FSUNLinSol_Dense(sunvec_y, sunmat_A, sunctx)
   if (.not. associated(sunls)) then
     print *, 'ERROR: sunls = NULL'

include/sundials/sundials_math.h

@@ -160,6 +160,54 @@ 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 : SUNRpowerR
+ * -----------------------------------------------------------------
+ * Usage : sunrealtype base, exponent, ans;
+ *         ans = SUNRpowerR(base,exponent);
+ * -----------------------------------------------------------------
+ * SUNRpowerR returns the value of base^exponent, where both base and
+ * 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 : SUNRround
@@ -213,21 +261,6 @@ SUNDIALS_EXPORT int SUNIpowerI(int base, int exponent);
 
 SUNDIALS_EXPORT sunrealtype SUNRpowerI(sunrealtype base, int exponent);
 
-/*
- * -----------------------------------------------------------------
- * Function : SUNRpowerR
- * -----------------------------------------------------------------
- * Usage : sunrealtype base, exponent, ans;
- *         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.
- * -----------------------------------------------------------------
- */
-
-SUNDIALS_EXPORT sunrealtype SUNRpowerR(sunrealtype base, sunrealtype exponent);
-
 /*
  * -----------------------------------------------------------------
  * Function : SUNRCompare

src/arkode/arkode.c

@@ -3010,20 +3010,6 @@ int arkRwtSetSV(ARKodeMem ark_mem, N_Vector My, N_Vector weight)
   return (0);
 }
 
-/*---------------------------------------------------------------
-  arkExpStab is the default explicit stability estimation function
-  ---------------------------------------------------------------*/
-int arkExpStab(SUNDIALS_MAYBE_UNUSED N_Vector y,
-               SUNDIALS_MAYBE_UNUSED sunrealtype t, sunrealtype* hstab,
-               SUNDIALS_MAYBE_UNUSED void* data)
-{
-  /* explicit stability not used by default,
-     set to zero to disable */
-  *hstab = SUN_RCONST(0.0);
-
-  return (ARK_SUCCESS);
-}
-
 /*---------------------------------------------------------------
   arkPredict_MaximumOrder
 

src/arkode/arkode_adapt.c

@@ -79,7 +79,7 @@ void arkPrintAdaptMem(ARKodeHAdaptMem hadapt_mem, FILE* outfile)
     fprintf(outfile, "ark_hadapt: p = %i\n", hadapt_mem->p);
     fprintf(outfile, "ark_hadapt: q = %i\n", hadapt_mem->q);
     fprintf(outfile, "ark_hadapt: adjust = %i\n", hadapt_mem->adjust);
-    if (hadapt_mem->expstab == arkExpStab)
+    if (hadapt_mem->expstab == NULL)
     {
       fprintf(outfile, "  ark_hadapt: Default explicit stability function\n");
     }
@@ -106,7 +106,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;
   int controller_order;
 
   /* Return with no stepsize adjustment if the controller is NULL */
@@ -138,49 +138,54 @@ 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)
-  {
-    arkProcessError(ark_mem, ARK_ILL_INPUT, __LINE__, __func__, __FILE__,
-                    "Error in explicit stability function.");
-    return (ARK_ILL_INPUT);
-  }
-  if (h_cfl <= ZERO) { h_cfl = SUN_RCONST(1.0e30) * SUNRabs(hcur); }
-
-  SUNLogDebug(ARK_LOGGER, "new-step-before-bounds",
-              "h_acc = " SUN_FORMAT_G ", h_cfl = " SUN_FORMAT_G, h_acc, h_cfl);
+  SUNLogDebug(ARK_LOGGER, "new-step-before-bounds", "h_acc = " SUN_FORMAT_G,
+              h_acc);
 
   /* enforce safety factors */
   h_acc *= hadapt_mem->safety;
-  h_cfl *= hadapt_mem->cfl * int_dir;
 
   /* 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));
 
-  SUNLogDebug(ARK_LOGGER, "new-step-after-max-min-bounds",
-              "h_acc = " SUN_FORMAT_G ", h_cfl = " SUN_FORMAT_G, h_acc, h_cfl);
+  if (hadapt_mem->expstab != NULL)
+  {
+    sunrealtype h_cfl = ZERO;
+    retval = hadapt_mem->expstab(ycur, tcur, &h_cfl, hadapt_mem->estab_data);
+    if (retval != ARK_SUCCESS)
+    {
+      arkProcessError(ark_mem, ARK_ILL_INPUT, __LINE__, __func__, __FILE__,
+                      "Error in explicit stability function.");
+      return (ARK_ILL_INPUT);
+    }
 
-  /* 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));
+    h_cfl *= hadapt_mem->cfl;
+    SUNLogDebug(ARK_LOGGER, "new-step-cfl", "h_cfl = " SUN_FORMAT_G, h_cfl);
+
+    if (h_cfl > ZERO && h_cfl < h_acc)
+    {
+      hadapt_mem->nst_exp++;
+      h_acc = h_cfl;
+    }
+    else { hadapt_mem->nst_acc++; }
+  }
+  else { hadapt_mem->nst_acc++; }
 
   /* 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);
+
+  SUNLogDebug(ARK_LOGGER, "new-step-after-max-min-bounds",
+              "h_acc = " SUN_FORMAT_G, h_acc);
 
   /* set basic value of ark_eta */
   ark_mem->eta = h_acc / hcur;

src/arkode/arkode_io.c

@@ -93,12 +93,12 @@ int ARKodeSetDefaults(void* arkode_mem)
   ark_mem->hadapt_mem->growth = GROWTH; /* step adaptivity growth factor */
   ark_mem->hadapt_mem->lbound = HFIXED_LB; /* step adaptivity no-change lower bound */
   ark_mem->hadapt_mem->ubound = HFIXED_UB; /* step adaptivity no-change upper bound */
-  ark_mem->hadapt_mem->expstab = arkExpStab; /* internal explicit stability fn */
-  ark_mem->hadapt_mem->estab_data = NULL;    /* no explicit stability fn data */
-  ark_mem->hadapt_mem->pq         = PQ;      /* embedding order */
-  ark_mem->hadapt_mem->p          = 0;       /* no default embedding order */
-  ark_mem->hadapt_mem->q          = 0;       /* no default method order */
-  ark_mem->hadapt_mem->adjust     = ADJUST;  /* controller order adjustment */
+  ark_mem->hadapt_mem->expstab    = NULL;   /* internal explicit stability fn */
+  ark_mem->hadapt_mem->estab_data = NULL;   /* no explicit stability fn data */
+  ark_mem->hadapt_mem->pq         = PQ;     /* embedding order */
+  ark_mem->hadapt_mem->p          = 0;      /* no default embedding order */
+  ark_mem->hadapt_mem->q          = 0;      /* no default method order */
+  ark_mem->hadapt_mem->adjust     = ADJUST; /* controller order adjustment */
 
   /* Set stepper defaults (if provided) */
   if (ark_mem->step_setdefaults)
@@ -1290,7 +1290,7 @@ int ARKodeSetStepDirection(void* arkode_mem, sunrealtype stepdir)
       return retval;
     }
 
-    if (h != ZERO && ((h > ZERO) != (stepdir > ZERO)))
+    if (h != SUNRcopysign(h, stepdir))
     {
       /* Reverse the sign of h. If adaptive, h will be overwritten anyway by the
        * initial step estimation since ARKodeReset must be called before this.
@@ -1961,16 +1961,8 @@ int ARKodeSetStabilityFn(void* arkode_mem, ARKExpStabFn EStab, void* estab_data)
   }
 
   /* NULL argument sets default, otherwise set inputs */
-  if (EStab == NULL)
-  {
-    hadapt_mem->expstab    = arkExpStab;
-    hadapt_mem->estab_data = ark_mem;
-  }
-  else
-  {
-    hadapt_mem->expstab    = EStab;
-    hadapt_mem->estab_data = estab_data;
-  }
+  hadapt_mem->expstab    = EStab;
+  hadapt_mem->estab_data = estab_data;
 
   return (ARK_SUCCESS);
 }
@@ -3143,7 +3135,7 @@ int ARKodeWriteParameters(void* arkode_mem, FILE* fp)
           ark_mem->hadapt_mem->lbound);
   fprintf(fp, "  Step growth upper bound = " SUN_FORMAT_G "\n",
           ark_mem->hadapt_mem->ubound);
-  if (ark_mem->hadapt_mem->expstab == arkExpStab)
+  if (ark_mem->hadapt_mem->expstab == NULL)
   {
     fprintf(fp, "  Default explicit stability function\n");
   }

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,28 +127,27 @@ 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));
   }