Simplify/Fix Adaptivity for Small Errors or no CFL

src/arkode/arkode_adapt.c

@@ -131,7 +131,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",
@@ -156,7 +156,7 @@ int arkAdapt(ARKodeMem ark_mem, ARKodeHAdaptMem hadapt_mem, N_Vector ycur,
 #endif
 
   /* increment the relevant step counter, set desired step */
-  if (SUNRabs(h_acc) < SUNRabs(h_cfl)) { hadapt_mem->nst_acc++; }
+  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));
 

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
@@ -135,15 +134,15 @@ SUNErrCode SUNAdaptController_EstimateStep_ImExGus(SUNAdaptController C,
 
   /* 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);
+  const sunrealtype e = SACIMEXGUS_BIAS(C) * dsm;
 
   /* 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 e1   = SACIMEXGUS_BIAS(C) * dsm;
+  const sunrealtype e2   = e1 / SACIMEXGUS_EP(C);
   const sunrealtype hrat = h / SACIMEXGUS_HP(C);
 
   /* compute estimated time step size, modifying the first step formula */
@@ -154,6 +153,13 @@ SUNErrCode SUNAdaptController_EstimateStep_ImExGus(SUNAdaptController C,
                        SUNRpowerR(e1, k1e) * SUNRpowerR(e2, k2e));
   }
 
+  if (isnan(*hnew))
+  {
+    /* hnew can be NAN if multiple e's are 0. In that case, make hnew inf with
+     * the same sign as h */
+    *hnew = h / SUN_RCONST(0.0);
+  }
+
   /* 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
@@ -322,9 +321,9 @@ SUNErrCode SUNAdaptController_EstimateStep_Soderlind(SUNAdaptController C,
   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 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);
 
@@ -340,6 +339,13 @@ SUNErrCode SUNAdaptController_EstimateStep_Soderlind(SUNAdaptController C,
             SUNRpowerR(hrat, k4) * SUNRpowerR(hrat2, k5);
   }
 
+  if (isnan(*hnew))
+  {
+    /* hnew can be NAN if multiple e's are 0. In that case, make hnew inf with
+     * the same sign as h */
+    *hnew = h / SUN_RCONST(0.0);
+  }
+
   /* return with success */
   return SUN_SUCCESS;
 }

src/sundials/sundials_math.c

@@ -22,6 +22,8 @@
 #include <sundials/sundials_math.h>
 #include <sundials/sundials_types.h>
 
+#include <assert.h>
+
 sunrealtype SUNRpowerI(sunrealtype base, int exponent)
 {
   int i, expt;
@@ -36,7 +38,9 @@ sunrealtype SUNRpowerI(sunrealtype base, int exponent)
 
 sunrealtype SUNRpowerR(sunrealtype base, sunrealtype exponent)
 {
-  if (base <= SUN_RCONST(0.0)) { return (SUN_RCONST(0.0)); }
+  // TODO(SBR): cleanup this and header
+  // if (base <= SUN_RCONST(0.0)) { return (SUN_RCONST(0.0)); }
+  assert(base > 0.0);
 
 #if defined(SUNDIALS_DOUBLE_PRECISION)
   return (pow(base, exponent));