Feautre/adapt opt

CHANGELOG.md

@@ -25,6 +25,8 @@ Improved the efficiency of default ARKODE methods with the following changes:
 | 4th Order ARK      | `ARKODE_ARK436L2SA_ERK_6_3_4` and `ARKODE_ARK436L2SA_DIRK_6_3_4` | `ARKODE_ARK437L2SA_ERK_7_3_4` and `ARKODE_ARK437L2SA_DIRK_7_3_4`   |
 | 5th Order ARK      | `ARKODE_ARK548L2SA_ERK_8_4_5` and `ARKODE_ARK548L2SA_DIRK_8_4_5` | `ARKODE_ARK548L2SAb_ERK_8_4_5` and `ARKODE_ARK548L2SAb_DIRK_8_4_5` |
 
+Added new temporal adaptivity controller utility constructors, `SUNAdaptController_H0211`, `SUNAdaptController_H211`, and `SUNAdaptController_H312`.
+
 The default value of `CMAKE_CUDA_ARCHITECTURES` is no longer set to `70` and is
 now determined automatically by CMake. The previous default was only valid for
 Volta GPUs while the automatically selected value will vary across compilers and
@@ -44,6 +46,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/FAQ.rst

@@ -286,6 +286,50 @@ CVODE(S) / IDA(S) / ARKODE
    integrals (very often a large number of them) need to be computed for adjoint sensitivity.
 
 
+.. collapse:: When should I select a non-default temporal adaptivity controller in ARKODE?
+
+   The default temporal adaptivity controller in ARKODE was selected due to its robust performance on test
+   problems that ranged in difficulty and stiffness, and when running with a wide range of solution tolerances
+   and method orders.  While we hope that this default runs well on most applications, it is unlikely to be
+   optimal.  A prime indicator that an alternate adaptivity controller may be useful is if the default results
+   in a large number of rejected steps.  Alternately, for higher-cost calculations where a reduction of 10%-20%
+   in the number of time steps would be important, users may want to try another controller option.
+
+   The default temporal adaptivity controller in ARKODE is the industry-standard *I* controller:
+
+   .. math::
+      h' = h_n \varepsilon_n^{-1/(p+1)}
+
+   where :math:`\varepsilon_n = \text{bias}*\text{dsm}_n`, and :math:`\text{dsm}_n` is the WRMS norm of the
+   local temporal error when using the time step :math:`h_n`, weighted by the user-requested tolerances.
+   However, the :ref:`SUNAdaptController class <SUNAdaptController>` in SUNDIALS provides a range of more
+   advanced temporal error controllers that could be applied to a given application problem, including the
+   :math:`H_{0}321`, :math:`H_{0}211`,  :math:`H211`, and :math:`H312` controllers from :cite:p:`Sod:03`,
+   as well as a variety of controllers (*PI*, *PID*, *ExpGus*, *ImpGus*, and *ImExGus*) that were included
+   in the initial ARKODE release.
+
+   The adaptive time-stepping controllers introduced by Soderlind in :cite:p:`Sod:03` can be classified into
+   two groups; *deadbeat* controllers and *non-deadbeat* controllers. A controller is known as deadbeat if
+   the roots of its characteristic equation are located at the origin. These controllers are generalized
+   forms of the *I* controller, are denoted with a zero subscript as part of the controller name (e.g.,
+   :math:`H_{0}321`), and are generally recommended for applications with smooth solutions as a function of time.
+
+   While it is impossible to exhaustively explore the question of controller optimality for all application
+   problems, we have performed tests on the range of provided controllers on a variety of stiff and non-stiff
+   problems, at varying tolerances (:math:`10^{-9} \to 10^{-3}`), and for Runge--Kutta methods at a wide range
+   of orders of accuracy (ERK orders 2-9, and DIRK orders 2-5). From our experiments, stiff problems benefitted
+   from the *I*, *PI*, :math:`H_{0}211`, :math:`H_{0}321`, and *ImpGus* controllers.  On the other hand,
+   non-stiff test problems ran most efficiently when using the *PID*, *I*, *ExpGus*, :math:`H211`, and
+   :math:`H312` controllers.
+
+   Lastly, we note that the internal controller parameters for the legacy ARKODE controllers (*PI*, *PID*,
+   *ExpGus*, *ImpGus*, and *ImExGus*) were determined via numerical optimization over a given set of test
+   problems.  As a result, although those controllers work very well for some applications, they may not
+   work well with others.  For users who are interested in exploring novel controller methods, we point out
+   that the :ref:`"Soderlind" <SUNAdaptController.Soderlind>` SUNAdaptController allows complete control
+   over all internal parameters via the :c:func:`SUNAdaptController_SetParams_Soderlind` function.
+
+
 KINSOL
 ------
 

doc/shared/RecentChanges.rst

@@ -35,6 +35,7 @@ Improved the efficiency of default ARKODE methods with the following changes:
 |                    | ``ARKODE_ARK548L2SA_DIRK_8_4_5``    | ``ARKODE_ARK548L2SAb_DIRK_8_4_5``    |
 +--------------------+-------------------------------------+--------------------------------------+
 
+Added new temporal adaptivity controller utility constructors, :c:func:`SUNAdaptController_H0211`, :c:func:`SUNAdaptController_H211`, and :c:func:`SUNAdaptController_H312`.
 
 The default value of :cmakeop:`CMAKE_CUDA_ARCHITECTURES` is no longer set to
 ``70`` and is now determined automatically by CMake. The previous default was
@@ -55,6 +56,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,12 @@ 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.  We note that
+this combination of parameters corresponds with the :math:`H_{0}321`
+controller described in :cite:p:`Sod:03`.  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_n^{-1/(p+1)}` is used.
 
 The SUNAdaptController_Soderlind controller is implemented as a derived
 SUNAdaptController class, and defines its *content* field as:
@@ -72,11 +73,11 @@ The header file to be included when using this module is
 ``sunadaptcontroller/sunadaptcontroller_soderlind.h``.
 
 We note that through appropriate selection of the parameters :math:`k_*`,
-this controller may create a wide range of proposed temporal adaptivity controllers,
-including the PID, PI, I, as well as Gustafsson's explicit and implicit controllers,
-:cite:p:`Gust:91` and :cite:p:`Gust:94`.  As a convenience, utility routines to
-create these controllers and set their parameters (as special cases of the
-SUNAdaptController_Soderlind) are provided.
+this controller may create a wide range of proposed temporal adaptivity
+controllers, including the :math:`H_{0}321`, I, PI, PID, as well as Gustafsson's
+explicit and implicit controllers, :cite:p:`Gust:91` and :cite:p:`Gust:94`,
+among others.  As a convenience, utility routines to create a variety of these
+controllers and set their parameters (as special cases of SUNAdaptController_Soderlind) are provided.
 
 The SUNAdaptController_Soderlind class provides implementations of all operations
 relevant to a ``SUN_ADAPTCONTROLLER_H`` controller listed in
@@ -243,7 +244,7 @@ also provides the following additional user-callable routines:
 .. c:function:: SUNErrCode SUNAdaptController_SetParams_ExpGus(SUNAdaptController C, sunrealtype k1_hat, sunrealtype k2_hat)
 
    This user-callable function provides control over the relevant parameters
-   above for the explicit Gustafsson controller, setting :math:`k_3 = k_4 = k_5 = 0`. 
+   above for the explicit Gustafsson controller, setting :math:`k_3 = k_4 = k_5 = 0`.
    This should be called *before* the time integrator is called to evolve the problem.
 
    .. note::
@@ -313,3 +314,54 @@ also provides the following additional user-callable routines:
    .. code-block:: c
 
       retval = SUNAdaptController_SetParams_ImpGus(C, 0.98, 0.95);
+
+.. c:function:: SUNAdaptController SUNAdaptController_H0211(SUNContext sunctx)
+
+   This constructor creates and allocates memory for a
+   SUNAdaptController_Soderlind object, set up to replicate the :math:`H_{0}211`
+   controller from :cite:p:`Sod:03`, corresponding with the parameters :math:`k_1=0.5`,
+   :math:`k_2=0.5`, :math:`k_4=-0.5`, and :math:`k_3=k_5=0`.
+
+   :param sunctx: the current :c:type:`SUNContext` object.
+   :return: if successful, a usable :c:type:`SUNAdaptController` object;
+            otherwise it will return ``NULL``.
+
+   Usage:
+
+   .. code-block:: c
+
+      SUNAdaptController C = SUNAdaptController_H0211(sunctx);
+
+.. c:function:: SUNAdaptController SUNAdaptController_H211(SUNContext sunctx)
+
+   This constructor creates and allocates memory for a
+   SUNAdaptController_Soderlind object, set up to replicate the :math:`H211`
+   controller from :cite:p:`Sod:03`, corresponding with the parameters :math:`k_1=0.25`,
+   :math:`k_2=0.25`, :math:`k_4=-0.25`, and :math:`k_3=k_5=0`.
+
+   :param sunctx: the current :c:type:`SUNContext` object.
+   :return: if successful, a usable :c:type:`SUNAdaptController` object;
+            otherwise it will return ``NULL``.
+
+   Usage:
+
+   .. code-block:: c
+
+      SUNAdaptController C = SUNAdaptController_H211(sunctx);
+
+.. c:function:: SUNAdaptController SUNAdaptController_H312(SUNContext sunctx)
+
+   This constructor creates and allocates memory for a
+   SUNAdaptController_Soderlind object, set up to replicate the :math:`H312`
+   controller from :cite:p:`Sod:03`, corresponding with the parameters :math:`k_1=0.125`,
+   :math:`k_2=0.25`, :math:`k_3=0.125`, :math:`k_4=-0.375`, and :math:`k_5=-0.125`.
+
+   :param sunctx: the current :c:type:`SUNContext` object.
+   :return: if successful, a usable :c:type:`SUNAdaptController` object;
+            otherwise it will return ``NULL``.
+
+   Usage:
+
+   .. code-block:: c
+
+      SUNAdaptController C = SUNAdaptController_H312(sunctx);

examples/arkode/CXX_serial/ark_heat2D.cpp

@@ -106,7 +106,9 @@ struct UserData
   sunrealtype atol;   // absolute tolerance
   sunrealtype hfixed; // fixed step size
   int order;          // ARKode method order
-  int controller;     // step size adaptivity method
+  int controller;     // step size adaptivity method: 0=PID, 1=PI,
+                      //    2=I, 3=ExpGus, 4=ImpGus, 5=ImExGus,
+                      //    6=H0321, 7=H0211, 8=H211, 9=H312
   int maxsteps;       // max number of steps between outputs
   bool linear;        // enable/disable linearly implicit option
   bool diagnostics;   // output diagnostics
@@ -360,12 +362,16 @@ int main(int argc, char* argv[])
   {
     switch (udata->controller)
     {
-    case (ARK_ADAPT_PID): C = SUNAdaptController_PID(ctx); break;
-    case (ARK_ADAPT_PI): C = SUNAdaptController_PI(ctx); break;
-    case (ARK_ADAPT_I): C = SUNAdaptController_I(ctx); break;
-    case (ARK_ADAPT_EXP_GUS): C = SUNAdaptController_ExpGus(ctx); break;
-    case (ARK_ADAPT_IMP_GUS): C = SUNAdaptController_ImpGus(ctx); break;
-    case (ARK_ADAPT_IMEX_GUS): C = SUNAdaptController_ImExGus(ctx); break;
+    case (0): C = SUNAdaptController_PID(ctx); break;
+    case (1): C = SUNAdaptController_PI(ctx); break;
+    case (2): C = SUNAdaptController_I(ctx); break;
+    case (3): C = SUNAdaptController_ExpGus(ctx); break;
+    case (4): C = SUNAdaptController_ImpGus(ctx); break;
+    case (5): C = SUNAdaptController_ImExGus(ctx); break;
+    case (6): C = SUNAdaptController_Soderlind(ctx); break;
+    case (7): C = SUNAdaptController_H0211(ctx); break;
+    case (8): C = SUNAdaptController_H211(ctx); break;
+    case (9): C = SUNAdaptController_H312(ctx); break;
     }
     flag = ARKodeSetAdaptController(arkode_mem, C);
     if (check_flag(&flag, "ARKodeSetAdaptController", 1)) { return 1; }

include/sunadaptcontroller/sunadaptcontroller_soderlind.h

@@ -125,6 +125,15 @@ SUNDIALS_EXPORT
 SUNErrCode SUNAdaptController_SetParams_ImpGus(SUNAdaptController C,
                                                sunrealtype k1, sunrealtype k2);
 
+SUNDIALS_EXPORT
+SUNAdaptController SUNAdaptController_H0211(SUNContext sunctx);
+
+SUNDIALS_EXPORT
+SUNAdaptController SUNAdaptController_H211(SUNContext sunctx);
+
+SUNDIALS_EXPORT
+SUNAdaptController SUNAdaptController_H312(SUNContext sunctx);
+
 #ifdef __cplusplus
 }
 #endif

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);
 
 /*
  * -----------------------------------------------------------------