Merge remote-tracking branch 'origin/feature/sunstepper' into feature…

…/sunadjoint

CHANGELOG.md

@@ -12,6 +12,14 @@ SUNStepper can be created from an ARKODE memory block with the new function
 `ARKodeCreateSUNStepper`. To enable interoperability with `MRIStepInnerStepper`,
 the function `MRIStepInnerStepper_CreateFromSUNStepper` was added.
 
+The following DIRK schemes now have coefficients accurate to quad precision:
+* `ARKODE_BILLINGTON_3_3_2`
+* `ARKODE_KVAERNO_4_2_3`
+* `ARKODE_CASH_5_2_4`
+* `ARKODE_CASH_5_3_4`
+* `ARKODE_KVAERNO_5_3_4`
+* `ARKODE_KVAERNO_7_4_5`
+
 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
@@ -54,6 +62,9 @@ using `ENABLE_KLU=ON` in combination with a static-only build of SuiteSparse.
 
 ### Deprecation Notices
 
+The ARKODE stepper specific functions to retrieve the number of right-hand side
+function evaluations have been deprecated. Use `ARKodeGetNumRhsEvals` instead.
+
 ## Changes to SUNDIALS in release 7.1.1
 
 ### Bug Fixes

benchmarks/advection_reaction_3D/kokkos/arkode_driver.cpp

@@ -189,8 +189,10 @@ int EvolveProblemDIRK(N_Vector y, UserData* udata, UserOptions* uopt)
   check_retval(&retval, "ARKodeGetNumSteps", 1, udata->myid);
   retval = ARKodeGetNumStepAttempts(arkode_mem, &nst_a);
   check_retval(&retval, "ARKodeGetNumStepAttempts", 1, udata->myid);
-  retval = ARKStepGetNumRhsEvals(arkode_mem, &nfe, &nfi);
-  check_retval(&retval, "ARKStepGetNumRhsEvals", 1, udata->myid);
+  retval = ARKodeGetNumRhsEvals(arkode_mem, 0, &nfe);
+  check_retval(&retval, "ARKodeGetNumRhsEvals", 1, udata->myid);
+  retval = ARKodeGetNumRhsEvals(arkode_mem, 1, &nfi);
+  check_retval(&retval, "ARKodeGetNumRhsEvals", 1, udata->myid);
   retval = ARKodeGetNumErrTestFails(arkode_mem, &netf);
   check_retval(&retval, "ARKodeGetNumErrTestFails", 1, udata->myid);
   retval = ARKodeGetNumNonlinSolvIters(arkode_mem, &nni);
@@ -392,8 +394,10 @@ int EvolveProblemIMEX(N_Vector y, UserData* udata, UserOptions* uopt)
   check_retval(&retval, "ARKodeGetNumSteps", 1, udata->myid);
   retval = ARKodeGetNumStepAttempts(arkode_mem, &nst_a);
   check_retval(&retval, "ARKodeGetNumStepAttempts", 1, udata->myid);
-  retval = ARKStepGetNumRhsEvals(arkode_mem, &nfe, &nfi);
-  check_retval(&retval, "ARKStepGetNumRhsEvals", 1, udata->myid);
+  retval = ARKodeGetNumRhsEvals(arkode_mem, 0, &nfe);
+  check_retval(&retval, "ARKodeGetNumRhsEvals", 1, udata->myid);
+  retval = ARKodeGetNumRhsEvals(arkode_mem, 1, &nfi);
+  check_retval(&retval, "ARKodeGetNumRhsEvals", 1, udata->myid);
   retval = ARKodeGetNumErrTestFails(arkode_mem, &netf);
   check_retval(&retval, "ARKodeGetNumErrTestFails", 1, udata->myid);
   retval = ARKodeGetNumNonlinSolvIters(arkode_mem, &nni);
@@ -516,8 +520,8 @@ int EvolveProblemExplicit(N_Vector y, UserData* udata, UserOptions* uopt)
   check_retval(&retval, "ARKodeGetNumSteps", 1, udata->myid);
   retval = ARKodeGetNumStepAttempts(arkode_mem, &nst_a);
   check_retval(&retval, "ARKodeGetNumStepAttempts", 1, udata->myid);
-  retval = ERKStepGetNumRhsEvals(arkode_mem, &nfe);
-  check_retval(&retval, "ERKStepGetNumRhsEvals", 1, udata->myid);
+  retval = ARKodeGetNumRhsEvals(arkode_mem, 0, &nfe);
+  check_retval(&retval, "ARKodeGetNumRhsEvals", 1, udata->myid);
   retval = ARKodeGetNumErrTestFails(arkode_mem, &netf);
   check_retval(&retval, "ARKodeGetNumErrTestFails", 1, udata->myid);
 

benchmarks/advection_reaction_3D/raja/arkode_driver.cpp

@@ -190,8 +190,10 @@ int EvolveProblemDIRK(N_Vector y, UserData* udata, UserOptions* uopt)
   check_retval(&retval, "ARKodeGetNumSteps", 1, udata->myid);
   retval = ARKodeGetNumStepAttempts(arkode_mem, &nst_a);
   check_retval(&retval, "ARKodeGetNumStepAttempts", 1, udata->myid);
-  retval = ARKStepGetNumRhsEvals(arkode_mem, &nfe, &nfi);
-  check_retval(&retval, "ARKStepGetNumRhsEvals", 1, udata->myid);
+  retval = ARKodeGetNumRhsEvals(arkode_mem, 0, &nfe);
+  check_retval(&retval, "ARKodeGetNumRhsEvals", 1, udata->myid);
+  retval = ARKodeGetNumRhsEvals(arkode_mem, 1, &nfi);
+  check_retval(&retval, "ARKodeGetNumRhsEvals", 1, udata->myid);
   retval = ARKodeGetNumErrTestFails(arkode_mem, &netf);
   check_retval(&retval, "ARKodeGetNumErrTestFails", 1, udata->myid);
   retval = ARKodeGetNumNonlinSolvIters(arkode_mem, &nni);
@@ -394,8 +396,10 @@ int EvolveProblemIMEX(N_Vector y, UserData* udata, UserOptions* uopt)
   check_retval(&retval, "ARKodeGetNumSteps", 1, udata->myid);
   retval = ARKodeGetNumStepAttempts(arkode_mem, &nst_a);
   check_retval(&retval, "ARKodeGetNumStepAttempts", 1, udata->myid);
-  retval = ARKStepGetNumRhsEvals(arkode_mem, &nfe, &nfi);
-  check_retval(&retval, "ARKStepGetNumRhsEvals", 1, udata->myid);
+  retval = ARKodeGetNumRhsEvals(arkode_mem, 0, &nfe);
+  check_retval(&retval, "ARKodeGetNumRhsEvals", 1, udata->myid);
+  retval = ARKodeGetNumRhsEvals(arkode_mem, 1, &nfi);
+  check_retval(&retval, "ARKodeGetNumRhsEvals", 1, udata->myid);
   retval = ARKodeGetNumErrTestFails(arkode_mem, &netf);
   check_retval(&retval, "ARKodeGetNumErrTestFails", 1, udata->myid);
   retval = ARKodeGetNumNonlinSolvIters(arkode_mem, &nni);
@@ -519,8 +523,8 @@ int EvolveProblemExplicit(N_Vector y, UserData* udata, UserOptions* uopt)
   check_retval(&retval, "ARKodeGetNumSteps", 1, udata->myid);
   retval = ARKodeGetNumStepAttempts(arkode_mem, &nst_a);
   check_retval(&retval, "ARKodeGetNumStepAttempts", 1, udata->myid);
-  retval = ERKStepGetNumRhsEvals(arkode_mem, &nfe);
-  check_retval(&retval, "ERKStepGetNumRhsEvals", 1, udata->myid);
+  retval = ARKodeGetNumRhsEvals(arkode_mem, 0, &nfe);
+  check_retval(&retval, "ARKodeGetNumRhsEvals", 1, udata->myid);
   retval = ARKodeGetNumErrTestFails(arkode_mem, &netf);
   check_retval(&retval, "ARKodeGetNumErrTestFails", 1, udata->myid);
 

doc/arkode/guide/source/Butcher.rst

@@ -1407,12 +1407,12 @@ Here, the higher-order embedding is less stable than the lower-order method
 
    \renewcommand{\arraystretch}{1.5}
    \begin{array}{r|ccc}
-     0.292893218813 & 0.292893218813 & 0 & 0 \\
-     1.091883092037 & 0.798989873223 & 0.292893218813 & 0 \\
-     1.292893218813 & 0.740789228841 & 0.259210771159 & 0.292893218813 \\
+     1 - \frac{1}{\sqrt{2}} & 1 - \frac{1}{\sqrt{2}} & 0 & 0 \\
+     \frac{27}{\sqrt{2}} - 18 & 14 \sqrt{2} - 19 & 1 - \frac{1}{\sqrt{2}} & 0 \\
+     2 - \frac{1}{\sqrt{2}} & \frac{53 - 5 \sqrt{2}}{62} & \frac{9 + 5 \sqrt{2}}{62} & 1 - \frac{1}{\sqrt{2}} \\
      \hline
-     2 & 0.740789228840 & 0.259210771159 & 0 \\
-     3 & 0.691665115992 & 0.503597029883 & -0.195262145876
+     2 & \frac{53 - 5 \sqrt{2}}{62} & \frac{9 + 5 \sqrt{2}}{62} & 0 \\
+     3 & \frac{263 - 95 \sqrt{2}}{186} & \frac{47 + 33 \sqrt{2}}{186} & \frac{\sqrt{2} - 2}{3}
    \end{array}
 
 .. figure:: /figs/arkode/stab_region_12.png

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

@@ -3160,6 +3160,9 @@ Main solver optional output functions
       The *nfi_evals* value does not account for calls made to
       :math:`f^I` by a linear solver or preconditioner module.
 
+   .. deprecated:: x.y.z
+
+      Use :c:func:`ARKodeGetNumRhsEvals` instead.
 
 
 .. c:function:: int ARKStepGetNumErrTestFails(void* arkode_mem, long int* netfails)

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

@@ -1743,6 +1743,9 @@ Main solver optional output functions
       * *ARK_SUCCESS* if successful
       * *ARK_MEM_NULL* if the ERKStep memory was ``NULL``
 
+   .. deprecated:: x.y.z
+
+      Use :c:func:`ARKodeGetNumRhsEvals` instead.
 
 
 .. c:function:: int ERKStepGetNumErrTestFails(void* arkode_mem, long int* netfails)

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

@@ -1899,6 +1899,9 @@ Main solver optional output functions
 
    * *ARK_MEM_NULL* if the MRIStep memory was ``NULL``
 
+   .. deprecated:: x.y.z
+
+      Use :c:func:`ARKodeGetNumRhsEvals` instead.
 
 
 .. c:function:: int MRIStepGetNumStepSolveFails(void* arkode_mem, long int* ncnf)

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

@@ -789,6 +789,9 @@ Main solver optional output functions
    :retval ARK_SUCCESS: if successful
    :retval ARK_MEM_NULL: if the SPRKStep memory was ``NULL``
 
+   .. deprecated:: x.y.z
+
+      Use :c:func:`ARKodeGetNumRhsEvals` instead.
 
 
 .. c:function:: int SPRKStepGetCurrentMethod(void* arkode_mem, ARKodeSPRKTable *sprk_table)

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

@@ -3194,6 +3194,7 @@ Name of constant associated with a return flag         :c:func:`ARKodeGetReturnF
 No. of explicit stability-limited steps                :c:func:`ARKodeGetNumExpSteps`
 No. of accuracy-limited steps                          :c:func:`ARKodeGetNumAccSteps`
 No. of attempted steps                                 :c:func:`ARKodeGetNumStepAttempts`
+No. of RHS evaluations                                 :c:func:`ARKodeGetNumRhsEvals`
 No. of local error test failures that have occurred    :c:func:`ARKodeGetNumErrTestFails`
 No. of failed steps due to a nonlinear solver failure  :c:func:`ARKodeGetNumStepSolveFails`
 Estimated local truncation error vector                :c:func:`ARKodeGetEstLocalErrors`
@@ -3499,6 +3500,39 @@ Retrieve a pointer for user data                       :c:func:`ARKodeGetUserDat
    .. versionadded:: 6.1.0
 
 
+.. c:function:: int ARKodeGetNumRhsEvals(void* arkode_mem, int partition_index, long int* num_rhs_evals)
+
+   Returns the number of calls to the user's right-hand side function (so far).
+   For implicit methods or methods with an implicit partition, the count does
+   not include calls made by a linear solver or preconditioner.
+
+   :param arkode_mem: pointer to the ARKODE memory block.
+   :param num_partition: the right-hand side partition index:
+
+                   * For ERKStep, ``0`` corresponds to :math:`f(t,y)`
+
+                   * For ARKStep, ``0`` corresponds to :math:`f^E(t,y)` and
+                     ``1`` to :math:`f^I(t,y)`
+
+                   * For MRIStep, ``0`` corresponds to :math:`f^E(t,y)` and
+                     ``1`` to :math:`f^I(t,y)`
+
+                   * For SPRKStep, ``0`` corresponds to :math:`f_1(t,p)` and
+                     ``1`` to :math:`f_2(t,q)`
+
+                   A negative index will return the sum of the evaluations for
+                   each partition.
+
+   :param num_rhs_evals: the number of right-hand side evaluations.
+
+   :retval ARK_SUCCESS: the function exited successfully.
+   :retval ARK_MEM_NULL: if ``arkode_mem`` was ``NULL``.
+   :retval ARK_ILL_INPUT: if ``num_partiton`` was invalid for the stepper or
+                          ``num_rhs_evals`` was ``NULL``
+
+   .. versionadded:: x.y.z
+
+
 .. c:function:: int ARKodeGetNumErrTestFails(void* arkode_mem, long int* netfails)
 
    Returns the number of local error test failures that

doc/arkode/guide/source/sunstepper/index.rst

@@ -12,31 +12,7 @@
    SUNDIALS Copyright End
    ----------------------------------------------------------------
 
-.. _SUNStepper:
-
-#####################################
-Stepper Data Structure
-#####################################
-
-This section presents the :c:type:`SUNStepper` base class which represents a
-generic solution procedure for IVPs of the form
-
-.. math::
-   \dot{v}(t) = f(t, v) + r(t), \qquad v(t_0) = v_0,
-   :label: SUNStepper_IVP
-
-on an interval :math:`t \in [t_0, t_f]`. The time dependent forcing term,
-:math:`r_i(t)`, is given by
-
-.. math::
-   r(t) = \sum_{k = 0}^{n_{\text{forcing}}}
-   \frac{t - t_{\text{shift}}}{t_{\text{scale}}} \widehat{f}_k.
-   :label: SUNStepper_forcing
-
-:c:type:`SUNStepper` provides an abstraction over SUNDIALS integrators, custom
-integrators, exact solution procedures, or other approaches for solving
-:eq:`SUNStepper_IVP`. These are used, for example, in operator splitting and
-forcing methods to solve inner IVPs in a flexible way.
+.. include:: ../../../../shared/sunstepper/SUNStepper_Structure.rst
 
 .. toctree::
    :maxdepth: 1

doc/shared/RecentChanges.rst

@@ -15,6 +15,20 @@ can be created from an ARKODE memory block with the new function
 :c:type:`MRIStepInnerStepper`, the function
 :c:func:`MRIStepInnerStepper_CreateFromSUNStepper` was added.
 
+The following DIRK schemes now have coefficients accurate to quad precision:
+
+* ``ARKODE_BILLINGTON_3_3_2``
+
+* ``ARKODE_KVAERNO_4_2_3``
+
+* ``ARKODE_CASH_5_2_4``
+
+* ``ARKODE_CASH_5_3_4``
+
+* ``ARKODE_KVAERNO_5_3_4``
+
+* ``ARKODE_KVAERNO_7_4_5``
+
 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
 only valid for Volta GPUs while the automatically selected value will vary
@@ -56,3 +70,7 @@ Fixed a CMake configuration issue related to aliasing an ``ALIAS`` target when
 using ``ENABLE_KLU=ON`` in combination with a static-only build of SuiteSparse.
 
 **Deprecation Notices**
+
+The ARKODE stepper specific functions to retrieve the number of right-hand side
+function evaluations have been deprecated. Use :c:func:`ARKodeGetNumRhsEvals`
+instead.

doc/shared/sunstepper/SUNStepper_Description.rst

@@ -22,7 +22,7 @@ The SUNStepper API
 
 As with other SUNDIALS classes, the :c:type:`SUNStepper` abstract base class is
 implemented using a C structure containing a ``content`` pointer to the derived
-class member data and a structure of function pointers the derived class
+class member data and a structure of function pointers to the derived class
 implementations of the virtual methods.
 
 .. c:type:: SUNStepper
@@ -70,7 +70,7 @@ from an ARKODE integrator.
 
       /* create an instance of the base class */
       SUNStepper stepper = NULL;
-      SUNErrCode err = SUNStepper_Create(&stepper);
+      SUNErrCode err = SUNStepper_Create(sunctx, &stepper);
 
    .. note::
 
@@ -98,29 +98,29 @@ from an ARKODE integrator.
 Stepping Functions
 ^^^^^^^^^^^^^^^^^^
 
-.. c:function:: SUNErrCode SUNStepper_Evolve(SUNStepper stepper, sunrealtype t0, sunrealtype tout, N_Vector vout, sunrealtype* tret)
+.. c:function:: SUNErrCode SUNStepper_Evolve(SUNStepper stepper, sunrealtype tout, N_Vector vret, sunrealtype* tret)
 
-   This function evolves the ODE :eq:`SUNStepper_IVP` from time *t0* to time
-   *tout* and stores the solution in *vout*.
+   This function evolves the ODE :eq:`SUNStepper_IVP` towards the time ``tout``
+   and stores the solution at time ``tret`` in ``vret``.
 
    :param stepper: the stepper object.
-   :param t0: the initial time for the integration.
-   :param tout: the final time for the integration.
-   :param vout: on output, the state at time *tout*.
-   :param tret: the final time corresponding to the output value *vout*.
+   :param tout: the time to evolve towards.
+   :param vret: on output, the state at time ``tret``.
+   :param tret: the time corresponding to the output value ``vret``.
    :return: A :c:type:`SUNErrCode` indicating success or failure.
 
 
-.. c:function:: SUNErrCode SUNStepper_OneStep(SUNStepper stepper, sunrealtype t0, sunrealtype tout, N_Vector vout, sunrealtype* tret)
+.. c:function:: SUNErrCode SUNStepper_FullRhs(SUNStepper stepper, sunrealtype t, N_Vector v, N_Vector f)
 
-   This function takes a single internal step starting at *t0* towards *tout*
-   and stores the next step solution in *vout*.
+   This function computes the full right-hand side function of the ODE,
+   :math:`f(t, v) + r(t)` in :eq:`SUNStepper_IVP` for a given value of the
+   independent variable ``t`` and state vector ``v``.
 
    :param stepper: the stepper object.
-   :param t0: the initial time for the integration.
-   :param tout: the final time for the integration.
-   :param vout: on output, the state after a single internal step.
-   :param tret: the final time corresponding to the output value *vout*.
+   :param t: the current value of the independent variable.
+   :param v: the current value of the dependent variable vector.
+   :param f: the output vector for the ODE right-hand side,
+      :math:`f(t, v) + r(t)`, in :eq:`SUNStepper_IVP`.
    :return: A :c:type:`SUNErrCode` indicating success or failure.
 
 
@@ -205,7 +205,7 @@ function to get and set a separate flag associated with a stepper.
 
 .. c:function:: SUNErrCode SUNStepper_SetLastFlag(SUNStepper stepper, int last_flag)
 
-   This function sets a flag than can be used by sunstepper implementations to
+   This function sets a flag that can be used by :c:type:`SUNStepper` implementations to
    indicate warnings or errors that occurred during an operation, e.g.,
    :c:func:`SUNStepper_Evolve`.
 
@@ -215,7 +215,7 @@ function to get and set a separate flag associated with a stepper.
 
 .. c:function:: SUNErrCode SUNStepper_GetLastFlag(SUNStepper stepper, int *last_flag)
 
-   This function provides the last value of the flag used by the sunstepper
+   This function provides the last value of the flag used by the :c:type:`SUNStepper`
    implementation to indicate warnings or errors that occurred during an
    operation, e.g., :c:func:`SUNStepper_Evolve`.
 
@@ -244,16 +244,6 @@ determined by the "consumer" of the :c:type:`SUNStepper`.
    :return: A :c:type:`SUNErrCode` indicating success or failure.
 
 
-.. c:function:: SUNErrCode SUNStepper_SetOneStepFn(SUNStepper stepper, SUNStepperOneStepFn fn)
-
-   This function attaches a :c:type:`SUNStepperOneStepFn` function to a
-   :c:type:`SUNStepper` object.
-
-   :param stepper: a stepper object.
-   :param fn: the :c:type:`SUNStepperOneStepFn` function to attach.
-   :return: A :c:type:`SUNErrCode` indicating success or failure.
-
-
 .. c:function:: SUNErrCode SUNStepper_SetFullRhsFn(SUNStepper stepper, SUNStepperFullRhsFn fn)
 
    This function attaches a :c:type:`SUNStepperFullRhsFn` function to a
@@ -303,30 +293,16 @@ This section describes the virtual methods defined by the :c:type:`SUNStepper`
 abstract base class.
 
 
-.. c:type:: SUNErrCode (*SUNStepperEvolveFn)(SUNStepper stepper, sunrealtype t0, sunrealtype tout, N_Vector v, sunrealtype* tret, int* stop_reason)
+.. c:type:: SUNErrCode (*SUNStepperEvolveFn)(SUNStepper stepper, sunrealtype tout, N_Vector v, sunrealtype* tret, int* stop_reason)
 
    This type represents a function with the signature of
    :c:func:`SUNStepper_Evolve`.
 
 
-.. c:type:: SUNErrCode (*SUNStepperOneStepFn)(SUNStepper stepper, sunrealtype t0, sunrealtype tout, N_Vector v, sunrealtype* tret, int* stop_reason)
+.. c:type:: SUNErrCode (*SUNStepperFullRhsFn)(SUNStepper stepper, sunrealtype t, N_Vector v, N_Vector f)
 
    This type represents a function with the signature of
-   :c:func:`SUNStepper_OneStep`.
-
-
-.. c:type:: SUNErrCode (*SUNStepperFullRhsFn)(SUNStepper stepper, sunrealtype t, N_Vector v, N_Vector f, int mode)
-
-   This type represents a function to compute the full right-hand side function
-   of the ODE, :math:`f(t, v) + r(t)` in :eq:`SUNStepper_IVP` for a given value
-   of the independent variable *t* and state vector *v*.
-
-   :param stepper: the stepper object.
-   :param t: the current value of the independent variable.
-   :param v: the current value of the dependent variable vector.
-   :param f: the output vector for the ODE right-hand side, :math:`f(t, v)`,
-      in :eq:`SUNStepper_IVP`.
-   :return: A :c:type:`SUNErrCode` indicating success or failure.
+   :c:func:`SUNStepper_FullRhs`.
 
 
 .. c:type:: SUNErrCode (*SUNStepperResetFn)(SUNStepper stepper, sunrealtype tR, N_Vector vR)