Clean up controllers

LLNL · Aug 1, 2024 · 2caec18 · 2caec18
1 parent b75e78e
commit 2caec18
Show file tree

Hide file tree

Showing 2 changed files with 19 additions and 23 deletions.
diff --git a/src/sunadaptcontroller/imexgus/sunadaptcontroller_imexgus.c b/src/sunadaptcontroller/imexgus/sunadaptcontroller_imexgus.c
@@ -132,31 +132,32 @@ 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 = 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   = 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 */
-  if (SACIMEXGUS_FIRSTSTEP(C)) { *hnew = h * SUNRpowerR(e, k); }
+  if (SACIMEXGUS_FIRSTSTEP(C)) {
+    // TODO: could this case be handled by setting e1 = e2?
+    /* 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 (isnan(*hnew))
+  if (!isfinite(*hnew))
   {
-    /* hnew can be NAN if multiple e's are 0. In that case, make hnew inf with
-     * the same sign as h */
+    /* 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 = h / SUN_RCONST(0.0);
   }
 

diff --git a/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c b/src/sunadaptcontroller/soderlind/sunadaptcontroller_soderlind.c
@@ -322,14 +322,9 @@ SUNErrCode SUNAdaptController_EstimateStep_Soderlind(SUNAdaptController C,
   const sunrealtype e[] = {SODERLIND_BIAS(C) * dsm, SODERLIND_EP(C), SODERLIND_EPP(C)};
   const sunrealtype hrat[] = {h / SODERLIND_HP(C), SODERLIND_HP(C) / SODERLIND_HPP(C)};
 
-  // This assumes ke[0] is nonzero. Otherwise we could have pow(0, 0) which is undefined behavior in C
   *hnew = h * SUNRpowerR(e[0], ke[0]);
   for (int i = 0; i < SUNMIN(2, SODERLIND_FIRSTSTEPS(C)); i++) {
-    *hnew *= SUNRpowerR(hrat[i], kh[i]);
-    if (ke[i + 1] != 0) {
-      // This check avoids pow(0, 0)
-      *hnew *= SUNRpowerR(e[i + 1], ke[i + 1]);
-    }
+    *hnew *= SUNRpowerR(hrat[i], kh[i]) * SUNRpowerR(e[i + 1], ke[i + 1]);
   }
 
   if (!isfinite(*hnew))