Feature/oo adaptivity2

CHANGELOG.md

@@ -3,6 +3,21 @@
 
 ## Changes to SUNDIALS in release X.X.X
 
+Added the `SUNAdaptController` base class, ported ARKODE's internal
+implementations of time step controllers into implementations of this class,
+and updated ARKODE to use these objects instead of its own implementations.
+Added `ARKStepSetAdaptController` and `ERKStepSetAdaptController` routines
+so that users can modify controller parameters, or even provide custom
+implementations.
+
+Added the routines `ARKStepSetAdaptivityAdjustment` and
+`ERKStepSetAdaptivityAdjustment`, that allow users to adjust the
+value for the method order supplied to the temporal adaptivity controllers.
+The ARKODE default for this adjustment has been -1 since its initial
+release, but for some applications a value of 0 is more appropriate.
+Users who notice that their simulations encounter a large number of
+temporal error test failures may want to experiment with adjusting this value.
+
 Fixed a regression introduced by the stop time bug fix in v6.6.1 where ARKODE,
 CVODE, CVODES, IDA, and IDAS would return at the stop time rather than the
 requested output time if the stop time was reached in the same step in which the

cmake/macros/SundialsAddLibrary.cmake

@@ -176,6 +176,15 @@ macro(sundials_add_library target)
       target_link_libraries(${obj_target} ${_all_libs})
     endif()
 
+    if(SUNDIALS_BUILD_WITH_PROFILING)
+      if(ENABLE_CALIPER)
+        target_link_libraries(${obj_target} PUBLIC caliper)
+      endif()
+      if(ENABLE_ADIAK)
+        target_link_libraries(${obj_target} PUBLIC adiak::adiak ${CMAKE_DL_LIBS})
+      endif()
+    endif()
+
     # add includes to object library
     target_include_directories(${obj_target}
       PUBLIC

doc/arkode/guide/source/Constants.rst

@@ -468,6 +468,8 @@ contains the ARKODE output constants.
    | :index:`ARK_RELAX_JAC_FAIL`         | -46  | The relaxation Jacobian function returned an unrecoverable |
    |                                     |      | error                                                      |
    +-------------------------------------+------+------------------------------------------------------------+
+   | :index:`ARK_CONTROLLER_ERR`         | -47  | An error with a SUNAdaptController object was encountered. |
+   +-------------------------------------+------+------------------------------------------------------------+
    | :index:`ARK_UNRECOGNIZED_ERROR`     | -99  | An unknown error was encountered.                          |
    +-------------------------------------+------+------------------------------------------------------------+
    |                                                                                                         |

doc/arkode/guide/source/Introduction.rst

@@ -133,6 +133,21 @@ Changes from previous versions
 Changes in vX.X.X
 -----------------
 
+Added the :c:type:`SUNAdaptController` base class, ported ARKODE's internal
+implementations of time step controllers into implementations of this class,
+and updated ARKODE to use these objects instead of its own implementations.  Added
+:c:func:`ARKStepSetAdaptController` and :c:func:`ERKStepSetAdaptController`
+routines so that users can modify controller parameters, or even provide custom
+implementations.
+
+Added the routines :c:func:`ARKStepSetAdaptivityAdjustment` and
+:c:func:`ERKStepSetAdaptivityAdjustment`, that allow users to adjust the
+value for the method order supplied to the temporal adaptivity controllers.
+The ARKODE default for this adjustment has been :math:`-1` since its initial
+release, but for some applications a value of :math:`0` is more appropriate.
+Users who notice that their simulations encounter a large number of
+temporal error test failures may want to experiment with adjusting this value.
+
 Fixed the build system support for MAGMA when using a NVIDIA HPC SDK installation of CUDA
 and fixed the targets used for rocBLAS and rocSPARSE.
 

doc/arkode/guide/source/Mathematics.rst

@@ -906,147 +906,24 @@ steps, :math:`t_{n-3} \to t_{n-2} \to t_{n-1} \to t_n`.  These local
 error history values are all initialized to 1 upon program
 initialization, to accommodate the few initial time steps of a
 calculation where some of these error estimates have not yet been
-computed.  With these estimates, ARKODE supports a variety of error
-control algorithms, as specified in the subsections below.
+computed.  With these estimates, ARKODE supports one of two approaches
+for temporal error control.
 
+First, any valid implementation of the SUNAdaptController class
+:numref:`SUNAdaptController.Description` may be used by ARKODE's adaptive
+time-stepping modules to provide a candidate error-based prospective step
+size :math:`h'`.
 
-.. _ARKODE.Mathematics.Adaptivity.ErrorControl.PID:
-
-PID controller
------------------
-
-This is the default time adaptivity controller used by the ARKStep and
-ERKStep modules.  It derives from those found in :cite:p:`KenCarp:03`, :cite:p:`Sod:98`, :cite:p:`Sod:03` and
-:cite:p:`Sod:06`, and uses all three of the local error estimates
-:math:`\varepsilon_n`, :math:`\varepsilon_{n-1}` and
-:math:`\varepsilon_{n-2}` in determination of a prospective step size,
-
-.. math::
-   h' \;=\; h_n\; \varepsilon_n^{-k_1/p}\; \varepsilon_{n-1}^{k_2/p}\;
-        \varepsilon_{n-2}^{-k_3/p},
-
-where the constants :math:`k_1`, :math:`k_2` and :math:`k_3` default
-to 0.58, 0.21 and 0.1, respectively, and may be modified by the user.
-In this estimate, a floor of :math:`\varepsilon > 10^{-10}` is
-enforced to avoid division-by-zero errors.
-
-
-
-.. _ARKODE.Mathematics.Adaptivity.ErrorControl.PI:
-
-PI controller
-----------------------
-
-Like with the previous method, the PI controller derives from those
-found in :cite:p:`KenCarp:03`, :cite:p:`Sod:98`, :cite:p:`Sod:03` and :cite:p:`Sod:06`, but it differs in
-that it only uses the two most recent step sizes in its adaptivity
-algorithm,
-
-.. math::
-   h' \;=\; h_n\; \varepsilon_n^{-k_1/p}\; \varepsilon_{n-1}^{k_2/p}.
-
-Here, the default values of :math:`k_1` and :math:`k_2` default
-to 0.8 and 0.31, respectively, though they may be changed by the user.
-
-
-
-.. _ARKODE.Mathematics.Adaptivity.ErrorControl.I:
-
-I controller
-----------------------
-
-This is the standard time adaptivity control algorithm in use by most
-publicly-available ODE solver codes.  It bases the prospective time step
-estimate entirely off of the current local error estimate,
-
-.. math::
-   h' \;=\; h_n\; \varepsilon_n^{-k_1/p}.
-
-By default, :math:`k_1=1`, but that may be modified by the user.
-
-
-
-
-.. _ARKODE.Mathematics.Adaptivity.ErrorControl.eGus:
-
-Explicit Gustafsson controller
----------------------------------
-
-This step adaptivity algorithm was proposed in :cite:p:`Gust:91`, and
-is primarily useful with explicit Runge--Kutta methods.
-In the notation of our earlier controllers, it has the form
-
-.. math::
-   h' \;=\; \begin{cases}
-      h_1\; \varepsilon_1^{-1/p}, &\quad\text{on the first step}, \\
-      h_n\; \varepsilon_n^{-k_1/p}\;
-        \left(\dfrac{\varepsilon_n}{\varepsilon_{n-1}}\right)^{k_2/p}, &
-      \quad\text{on subsequent steps}.
-   \end{cases}
-   :label: ARKODE_expGus
-
-The default values of :math:`k_1` and :math:`k_2` are 0.367 and 0.268,
-respectively, and may be modified by the user.
-
-
-
-
-.. _ARKODE.Mathematics.Adaptivity.ErrorControl.iGus:
-
-Implicit Gustafsson controller
----------------------------------
-
-A version of the above controller suitable for implicit Runge--Kutta
-methods was introduced in :cite:p:`Gust:94`, and has the form
-
-.. math::
-   h' = \begin{cases}
-      h_1 \varepsilon_1^{-1/p}, &\quad\text{on the first step}, \\
-      h_n \left(\dfrac{h_n}{h_{n-1}}\right) \varepsilon_n^{-k_1/p}
-        \left(\dfrac{\varepsilon_n}{\varepsilon_{n-1}}\right)^{-k_2/p}, &
-      \quad\text{on subsequent steps}.
-   \end{cases}
-   :label: ARKODE_impGus
-
-The algorithm parameters default to :math:`k_1 = 0.98` and
-:math:`k_2 = 0.95`, but may be modified by the user.
-
-
-
-
-.. _ARKODE.Mathematics.Adaptivity.ErrorControl.ieGus:
-
-ImEx Gustafsson controller
----------------------------------
-
-An ImEx version of these two preceding controllers is also available.
-This approach computes the estimates :math:`h'_1` arising from
-equation :eq:`ARKODE_expGus` and the estimate :math:`h'_2` arising from
-equation :eq:`ARKODE_impGus`, and selects
-
-.. math::
-   h' = \frac{h}{|h|}\min\left\{|h'_1|, |h'_2|\right\}.
-
-Here, equation :eq:`ARKODE_expGus` uses :math:`k_1` and
-:math:`k_2` with default values of 0.367 and 0.268, while equation
-:eq:`ARKODE_impGus` sets both parameters to the input :math:`k_3` that
-defaults to 0.95.  All of these values may be modified by the user.
-
-
-
-.. _ARKODE.Mathematics.Adaptivity.ErrorControl.User:
-
-User-supplied controller
----------------------------------
-
-Finally, ARKODE's time-stepping modules allow the user to define their
-own time step adaptivity function,
+Second, ARKODE's adaptive time-stepping modules currently still allow the
+user to define their own time step adaptivity function,
 
 .. math::
    h' = H(y, t, h_n, h_{n-1}, h_{n-2}, \varepsilon_n, \varepsilon_{n-1}, \varepsilon_{n-2}, q, p),
 
-to allow for problem-specific choices, or for continued
-experimentation with temporal error controllers.
+allowing for problem-specific choices, or for continued
+experimentation with temporal error controllers.  We note that this
+support has been deprecated in favor of the SUNAdaptController class,
+and will be removed in a future release.
 
 
 

doc/arkode/guide/source/Usage/ARKStep_c_interface/User_callable.rst

@@ -1107,8 +1107,9 @@ Set max number of constraint failures             :c:func:`ARKStepSetMaxNumConst
       :c:func:`ARKStepSetMaxGrowth()`,
       :c:func:`ARKStepSetMinReduction()`,
       :c:func:`ARKStepSetSafetyFactor()`,
-      :c:func:`ARKStepSetSmallNumEFails()` and
-      :c:func:`ARKStepSetStabilityFn()`
+      :c:func:`ARKStepSetSmallNumEFails()`,
+      :c:func:`ARKStepSetStabilityFn()`, and
+      :c:func:`ARKStepSetAdaptController()`
       will be ignored, since temporal adaptivity is disabled.
 
       If both :c:func:`ARKStepSetFixedStep()` and
@@ -1688,26 +1689,44 @@ the code, is provided in :numref:`ARKODE.Mathematics.Adaptivity`.
 
 .. cssclass:: table-bordered
 
-========================================================   ======================================  ========
-Optional input                                             Function name                           Default
-========================================================   ======================================  ========
-Set a custom time step adaptivity function                 :c:func:`ARKStepSetAdaptivityFn()`      internal
-Choose an existing time step adaptivity method             :c:func:`ARKStepSetAdaptivityMethod()`  0
-Explicit stability safety factor                           :c:func:`ARKStepSetCFLFraction()`       0.5
-Time step error bias factor                                :c:func:`ARKStepSetErrorBias()`         1.5
-Bounds determining no change in step size                  :c:func:`ARKStepSetFixedStepBounds()`   1.0  1.5
-Maximum step growth factor on convergence fail             :c:func:`ARKStepSetMaxCFailGrowth()`    0.25
-Maximum step growth factor on error test fail              :c:func:`ARKStepSetMaxEFailGrowth()`    0.3
-Maximum first step growth factor                           :c:func:`ARKStepSetMaxFirstGrowth()`    10000.0
-Maximum allowed general step growth factor                 :c:func:`ARKStepSetMaxGrowth()`         20.0
-Minimum allowed step reduction factor on error test fail   :c:func:`ARKStepSetMinReduction()`      0.1
-Time step safety factor                                    :c:func:`ARKStepSetSafetyFactor()`      0.96
-Error fails before MaxEFailGrowth takes effect             :c:func:`ARKStepSetSmallNumEFails()`    2
-Explicit stability function                                :c:func:`ARKStepSetStabilityFn()`       none
-========================================================   ======================================  ========
+=========================================================   ==========================================  ========
+Optional input                                              Function name                               Default
+=========================================================   ==========================================  ========
+Provide a :c:type:`SUNAdaptController` for ARKStep to use   :c:func:`ARKStepSetAdaptController()`       PID
+Set a custom time step adaptivity function                  :c:func:`ARKStepSetAdaptivityFn()`          internal
+Choose an existing time step adaptivity method              :c:func:`ARKStepSetAdaptivityMethod()`      0
+Adjust the method order used in the controller              :c:func:`ERKStepSetAdaptivityAdjustment()`  -1
+Explicit stability safety factor                            :c:func:`ARKStepSetCFLFraction()`           0.5
+Time step error bias factor                                 :c:func:`ARKStepSetErrorBias()`             1.5
+Bounds determining no change in step size                   :c:func:`ARKStepSetFixedStepBounds()`       1.0  1.5
+Maximum step growth factor on convergence fail              :c:func:`ARKStepSetMaxCFailGrowth()`        0.25
+Maximum step growth factor on error test fail               :c:func:`ARKStepSetMaxEFailGrowth()`        0.3
+Maximum first step growth factor                            :c:func:`ARKStepSetMaxFirstGrowth()`        10000.0
+Maximum allowed general step growth factor                  :c:func:`ARKStepSetMaxGrowth()`             20.0
+Minimum allowed step reduction factor on error test fail    :c:func:`ARKStepSetMinReduction()`          0.1
+Time step safety factor                                     :c:func:`ARKStepSetSafetyFactor()`          0.96
+Error fails before MaxEFailGrowth takes effect              :c:func:`ARKStepSetSmallNumEFails()`        2
+Explicit stability function                                 :c:func:`ARKStepSetStabilityFn()`           none
+=========================================================   ==========================================  ========
 
 
 
+.. c:function:: int ARKStepSetAdaptController(void* arkode_mem, SUNAdaptController C)
+
+   Sets a user-supplied time-step controller object.
+
+   **Arguments:**
+      * *arkode_mem* -- pointer to the ARKStep memory block.
+      * *C* -- user-supplied time adaptivity controller.  If ``NULL`` then the PID controller will be created (see :numref:`SUNAdaptController.Soderlind`).
+
+   **Return value:**
+      * *ARK_SUCCESS* if successful
+      * *ARK_MEM_NULL* if the ARKStep memory is ``NULL``
+      * *ARK_MEM_FAIL* if *C* was ``NULL`` and the PID controller could not be allocated.
+
+   .. versionadded:: x.x.x
+
+
 .. c:function:: int ARKStepSetAdaptivityFn(void* arkode_mem, ARKAdaptFn hfun, void* h_data)
 
    Sets a user-supplied time-step adaptivity function.
@@ -1729,6 +1748,11 @@ Explicit stability function                                :c:func:`ARKStepSetSt
       :c:func:`ARKStepSetStabilityFn()` should be used instead.
 
 
+   .. deprecated:: x.x.x
+
+      Use the SUNAdaptController infrastructure instead (see :numref:`SUNAdaptController.Description`).
+
+
 
 .. c:function:: int ARKStepSetAdaptivityMethod(void* arkode_mem, int imethod, int idefault, int pq, realtype* adapt_params)
 
@@ -1744,7 +1768,8 @@ Explicit stability function                                :c:func:`ARKStepSetSt
         parameters (1), or that they will be supplied in the
         *adapt_params* argument (0).
       * *pq* -- flag denoting whether to use the embedding order of
-        accuracy *p* (0) or the method order of accuracy *q* (1)
+        accuracy *p* (0), the method order of accuracy *q* (1), or the
+        minimum of the two (any input not equal to 0 or 1)
         within the adaptivity algorithm.  *p* is the default.
       * *adapt_params[0]* -- :math:`k_1` parameter within accuracy-based adaptivity algorithms.
       * *adapt_params[1]* -- :math:`k_2` parameter within accuracy-based adaptivity algorithms.
@@ -1761,8 +1786,40 @@ Explicit stability function                                :c:func:`ARKStepSetSt
       parameter values are desired, it is recommended to instead provide
       a custom function through a call to :c:func:`ARKStepSetAdaptivityFn()`.
 
+      .. versionchanged:: x.x.x
+
+         Prior to version x.x.x, any nonzero value for *pq* would result in use of the
+         embedding order of accuracy.
+
+
+   .. deprecated:: x.x.x
+
+      Use the SUNAdaptController infrastructure instead (see :numref:`SUNAdaptController.Description`).
 
 
+
+.. c:function:: int ARKStepSetAdaptivityAdjustment(void* arkode_mem, int adjust)
+
+   Called by a user to adjust the method order supplied to the temporal adaptivity
+   controller.  For example, if the user expects order reduction due to problem stiffness,
+   they may request that the controller assume a reduced order of accuracy for the method
+   by specifying a value :math:`adjust < 0`.
+
+   **Arguments:**
+      * *arkode_mem* -- pointer to the ARKStep memory block.
+      * *adjust* -- adjustment factor (default is -1).
+
+   **Return value:**
+      * *ARK_SUCCESS* if successful
+      * *ARK_MEM_NULL* if the ARKStep memory is ``NULL``
+      * *ARK_ILL_INPUT* if an argument has an illegal value
+
+   **Notes:**
+      This should be called prior to calling :c:func:`ARKStepEvolve()`, and can only be
+      reset following a call to :c:func:`ARKStepReInit()`.
+
+   .. versionadded:: x.x.x
+
 .. c:function:: int ARKStepSetCFLFraction(void* arkode_mem, realtype cfl_frac)
 
    Specifies the fraction of the estimated explicitly stable step to use.
@@ -1800,6 +1857,13 @@ Explicit stability function                                :c:func:`ARKStepSetSt
    **Notes:**
       Any value below 1.0 will imply a reset to the default value.
 
+      If both this and one of :c:func:`ARKStepSetAdaptivityMethod` or
+      :c:func:`ARKStepSetAdaptController` will be called, then this routine must be called
+      *second*.
+
+   .. deprecated:: x.x.x
+
+      Use the SUNAdaptController infrastructure instead (see :numref:`SUNAdaptController.Description`).
 
 
 .. c:function:: int ARKStepSetFixedStepBounds(void* arkode_mem, realtype lb, realtype ub)