Feature: SUNDomEigEstimator upgrades

bindings/sundials4py/arkode/arkode_lsrkstep_generated.hpp

@@ -42,6 +42,9 @@ m.def("LSRKStepSetMaxNumStages", LSRKStepSetMaxNumStages, nb::arg("arkode_mem"),
 m.def("LSRKStepSetDomEigSafetyFactor", LSRKStepSetDomEigSafetyFactor,
       nb::arg("arkode_mem"), nb::arg("dom_eig_safety"));
 
+m.def("LSRKStepSetUseAnalyticStabRegion", LSRKStepSetUseAnalyticStabRegion,
+      nb::arg("arkode_mem"), nb::arg("analytic_stab_region"));
+
 m.def("LSRKStepSetNumDomEigEstInitPreprocessIters",
       LSRKStepSetNumDomEigEstInitPreprocessIters, nb::arg("arkode_mem"),
       nb::arg("num_iters"));

bindings/sundials4py/sundials/generate.yaml

@@ -63,6 +63,7 @@ modules:
     - ../../include/sundials/sundials_domeigestimator.h
     fn_exclude_by_name__regex:
     - "SUNDomEigEstimator_SetATimes" # nanobind cannot bind to functions which take a function pointer, so we do something custom
+    - "SUNDomEigEstimator_SetRhs" # nanobind cannot bind to functions which take a function pointer, so we do something custom
   sundials_errors:
     path: sundials/sundials_errors_generated.hpp
     headers:

bindings/sundials4py/sundials/sundials_domeigestimator.cpp

@@ -83,6 +83,26 @@ void bind_sundomeigestimator(nb::module_& m)
       else { return SUNDomEigEstimator_SetATimes(dee, fntable, nullptr); }
     },
     nb::arg("DEE"), nb::arg("ATimes").none());
+
+  m.def(
+    "SUNDomEigEstimator_SetRhs",
+    [](SUNDomEigEstimator DEE,
+       std::function<std::remove_pointer_t<SUNRhsFn>> RHSfn) -> SUNErrCode
+    {
+      if (!DEE->python) { DEE->python = new SUNDomEigEstimatorFunctionTable; }
+
+      auto fntable = static_cast<SUNDomEigEstimatorFunctionTable*>(DEE->python);
+
+      fntable->deerhs = nb::cast(RHSfn);
+
+      if (RHSfn)
+      {
+        return SUNDomEigEstimator_SetRhs(DEE, fntable,
+                                         sundomeigestimator_setrhs_wrapper);
+      }
+      else { return SUNDomEigEstimator_SetRhs(DEE, fntable, nullptr); }
+    },
+    nb::arg("DEE"), nb::arg("RHSfn").none());
 }
 
 } // namespace sundials4py

bindings/sundials4py/sundials/sundials_domeigestimator_usersupplied.hpp

@@ -33,6 +33,7 @@ using namespace sundials::experimental;
 struct SUNDomEigEstimatorFunctionTable
 {
   nb::object atimes;
+  nb::object deerhs;
 };
 
 template<typename... Args>
@@ -43,4 +44,12 @@ SUNErrCode sundomeigestimator_atimes_wrapper(Args... args)
     3>(&SUNDomEigEstimatorFunctionTable::atimes, std::forward<Args>(args)...);
 }
 
+template<typename... Args>
+SUNErrCode sundomeigestimator_setrhs_wrapper(Args... args)
+{
+  return sundials4py::user_supplied_fn_caller<
+    std::remove_pointer_t<SUNRhsFn>, SUNDomEigEstimatorFunctionTable,
+    1>(&SUNDomEigEstimatorFunctionTable::deerhs, std::forward<Args>(args)...);
+}
+
 #endif

bindings/sundials4py/sundomeigest/generate.yaml

@@ -26,11 +26,13 @@ modules:
     fn_exclude_by_name__regex:
       # Don't interface the implementation specific overrides of the generic routines
       - "^SUNDomEigEstimator_SetATimes_.*"
+      - "^SUNDomEigEstimator_SetRhs_.*"
       - "^SUNDomEigEstimator_SetOptions_.*"
       - "^SUNDomEigEstimator_SetMaxIters_.*"
       - "^SUNDomEigEstimator_SetNumPreprocessIters_.*"
       - "^SUNDomEigEstimator_SetRelTol_.*"
       - "^SUNDomEigEstimator_SetInitialGuess_.*"
+      - "^SUNDomEigEstimator_SetRhsLinearizationPoint_.*"
       - "^SUNDomEigEstimator_Initialize_.*"
       - "^SUNDomEigEstimator_Estimate_.*"
       - "^SUNDomEigEstimator_GetRes_.*"

bindings/sundials4py/sundomeigest/sundomeigest_power_generated.hpp

@@ -32,6 +32,9 @@ m.def(
   },
   nb::arg("q"), nb::arg("max_iters"), nb::arg("rel_tol"), nb::arg("sunctx"),
   "nb::keep_alive<0, 4>()", nb::keep_alive<0, 4>());
+
+m.def("SUNDomEigEstimator_SetIsReal_Power", SUNDomEigEstimator_SetIsReal_Power,
+      nb::arg("DEE"), nb::arg("real"));
 // #ifdef __cplusplus
 //
 // #endif

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

@@ -315,6 +315,47 @@ Allowable Method Families
       This routine will be called by :c:func:`ARKodeSetOptions`
       when using the key "arkid.dom_eig_safety_factor".
 
+.. c:function:: int LSRKStepSetUseAnalyticStabRegion(void* arkode_mem, sunbooleantype analytic_stab_region);
+
+   Specifies whether to use the analytic stability region for determining the number of stages in STS methods.
+
+   **Arguments:**
+      * *arkode_mem* -- pointer to the LSRKStep memory block.
+      * *analytic_stab_region* -- Use the analytic stability region if ``SUNTRUE``; use the inscribed ellipse stability region if ``SUNFALSE``.
+
+   **Return value:**
+      * *ARK_SUCCESS* if successful
+      * *ARK_MEM_NULL* if ``arkode_mem`` was ``NULL``.
+
+   .. note::
+
+      This input is only used for RKC and RKL methods.
+
+      If :c:func:`LSRKStepSetUseAnalyticStabRegion` is not called, then the default
+      ``analytic_stab_region`` is set to ``SUNFALSE``.  This routine will be called by
+      :c:func:`ARKodeSetOptions` when using the key "arkid.use_analytic_stab_region".
+
+      :c:func:`LSRKStepSetUseAnalyticStabRegion` sets whether to use the exact stability region or an 
+      ellipse that is fully inscribed in the stability region for determining stability in RKC and RKL 
+      methods.  Whichever region is selected, LSRKStep will ensure that the complex number 
+      :math:`z=h\lambda`, where :math:`h` is the current time step size and :math:`\lambda` is the 
+      estimated dominant eigenvalue, is in this region.  Since the ellipse is smaller than the analytical
+      stability region it provides a more conservative estimate, which may be appropriate for problems 
+      wherein sub-dominant eigenvalues may also limit stability of the method.  Thus, the ellipse may 
+      result in a smaller time step size, more internal stages, or both, when compared to using the 
+      analytical stability region. We note that in either case, since LSRKStep does not examine the 
+      full eigenvalue spectrum, there is always some risk of instability.  For applications where 
+      this is a concern, we recommend calling :c:func:`LSRKStepSetDomEigSafetyFactor` with a `dom_eig_safety` 
+      value significantly larger than 1.
+
+      If :c:func:`LSRKStepSetUseAnalyticStabRegion` is called during integration, the change will take effect 
+      at the next step attempt. Both analytic and ellipse stability regions of RKC and RKL methods with 10 stages
+      are shown in the figure below. 
+
+      .. figure:: ../../../../../shared/figs/arkode/STS2_region_s10.png
+         :alt: Stability region of RKL method with 10 stages
+         :align: center
+         :width: 50%
 
 .. c:function:: int LSRKStepSetNumDomEigEstInitPreprocessIters(void* arkode_mem, int num_iters);
 

doc/shared/Changelog.rst

@@ -29,6 +29,12 @@ Changelog
 Changes to SUNDIALS in release X.Y.Z
 ====================================
 
+**Bug Fixes**
+
+Fixed a minor bug where the number of required stages for STS methods 
+in the LSRKStep module was incorrectly computed using the spectral 
+radius instead of the real part of the Jacobian eigenvalues.
+
 .. include:: RecentChanges_link.rst
 
 .. _Changelog.7.7.0:

doc/shared/RecentChanges.rst

@@ -7,4 +7,8 @@
 
 **Bug Fixes**
 
+Fixed a minor bug where the number of required stages for STS methods 
+in the LSRKStep module was incorrectly computed using the spectral 
+radius instead of the real part of the Jacobian eigenvalues.
+
 **Deprecation Notices**

doc/shared/sundomeigest/SUNDomEigEst_API.rst

@@ -209,6 +209,46 @@ instead of supplying a dummy routine.
       A :c:type:`SUNErrCode`.
 
 
+.. c:function:: SUNErrCode SUNDomEigEstimator_SetRhs(SUNDomEigEstimator DEE, void* rhs_data, SUNRhsFn RHSfn)
+
+   For applications that do not provide a :c:type:`SUNATimesFn` function to :c:func:`SUNDomEigEstimator_SetATimes`, 
+   the action of this matrix-vector product may be approximated internally.  If the matrix corresponds with 
+   the Jacobian of a vector-valued function, :math:`A = \frac{\partial f_{RHS}}{\partial y}(t,y)`, then the 
+   function :math:`f_{RHS}` may optionally be input via this routine, and the Jacobian-vector products will 
+   be approximated as 
+
+   .. math::
+      Av \approx \frac{1}{\sigma}\left[ f_{RHS}(t,y+\sigma v) - f_{RHS}(t,y)\right]
+
+   The linearization point :math:`(t,y)` should be separately supplied by calling :c:func:`SUNDomEigEstimator_SetRhsLinearizationPoint`.
+
+   **Arguments:**
+
+      * *DEE* -- a SUNDomEigEstimator object.
+      * *rhs_data* -- pointer to structure for ``RHSfn``.
+      * *RHSfn* -- function pointer to perform right-hand side evaluations.  This is typically the same as the problem-defining function supplied to CVODE or ARKODE.
+
+   **Return value:**
+
+      A :c:type:`SUNErrCode`.
+
+
+.. c:function:: SUNErrCode SUNDomEigEstimator_SetRhsLinearizationPoint(SUNDomEigEstimator DEE, sunrealtype t, N_Vector y)
+
+   This *optional* function sets the linearization point for the right-hand side function when using
+   :c:func:`SUNDomEigEstimator_SetRhs`.
+
+   **Arguments:**
+
+      * *DEE* -- a SUNDomEigEstimator object.
+      * *t* -- the time at which the linearization point is specified.
+      * *y* -- the linearization point for the right-hand side function.
+
+   **Return value:**
+
+      A :c:type:`SUNErrCode`.
+
+
 .. c:function:: SUNErrCode SUNDomEigEstimator_SetNumPreprocessIters(SUNDomEigEstimator DEE, int num_iters)
 
    This *optional* routine sets the number of preprocessing matrix-vector
@@ -244,6 +284,7 @@ instead of supplying a dummy routine.
       This routine will be called by :c:func:`SUNDomEigEstimator_SetOptions`
       when using the key "Did.num_preprocess_iters".
 
+
 .. c:function:: SUNErrCode SUNDomEigEstimator_SetRelTol(SUNDomEigEstimator DEE, sunrealtype rel_tol)
 
    This *optional* routine sets the estimator's :ref:`relative tolerance <pi_rel_tol>`.
@@ -259,6 +300,52 @@ instead of supplying a dummy routine.
 
    .. note::
 
+      The relative tolerance is used as a stopping criterion for either the Power Iteration method, or for the 
+      preprocessing phase of the Arnoldi iteration. Specifically, it defines the acceptable relative change 
+      between successive dominant eigenvalue estimates. It also serves as a threshold for determining 
+      whether the dominant eigenvalue estimated by the Power Iteration is real- or complex-valued.
+
+      When used to check convergence, we declare the iteration converged when the change in magnitude between successive estimates satisfies
+
+      .. math::
+
+         \left|\,|\lambda_{k}| - |\lambda_{k-1}|\,\right|
+         \le \mathtt{rel\_tol} \cdot |\lambda_{k}|.
+
+      When used to assess preprocessing iterations for the Arnoldi method, once the above test is satisfied, 
+      the Arnoldi iteration commences, ensuring that Arnoldi is performed only once, as opposed to restarting 
+      Arnoldi repeatedly.
+
+      When ``rel_tol`` is used as a threshold to assess whether the Power Iteration dominant eigenvalue is 
+      real versus complex, we define 
+
+      .. math::
+         \mathtt{gram\_det\_tol} = 10 \cdot \max\left(\varepsilon,\; \mathtt{rel\_tol}\right)
+
+      where :math:`\varepsilon` represents the machine precision.  The ``gram_det_tol`` value is used to 
+      assess the numerical rank of the 2×2 Gram matrix formed by the two most recent iterates
+      in the Power Iteration method -- if the determinant of this Gram matrix is less than or equal to 
+      ``gram_det_tol``, the iterates are considered linearly dependent, and the dominant eigenvalue is 
+      treated as real.  For this use case, to ensure that small imaginary parts of the dominant eigenvalue are 
+      not ignored, ``rel_tol`` should not be chosen too large. In practice, the smallest reliably detectable
+      imaginary part is proportional to the chosen relative tolerance, i.e.,
+
+      .. math::
+         |\beta| \gtrsim \mathcal{O}(\mathtt{rel\_tol}).
+
+      Therefore, to resolve an expected imaginary part of magnitude :math:`|\beta|`, it is recommended to choose
+
+      .. math::
+         \mathtt{rel\_tol} \ll |\beta|.
+
+      Choosing a smaller relative tolerance improves the ability to detect weakly complex eigenvalues,
+      but may increase computational cost.
+
+      Acceptable inputs to this routine include :math:`0 < \mathtt{rel\_tol} < 1 - \varepsilon`.  For Power 
+      Iteration, values outside this range result in the default value of ``rel_tol = 0.005``.  For 
+      Arnoldi, ``rel_tol < 0`` disables preprocessing, while 
+      :math:`\mathtt{rel\_tol} \ge 1-\varepsilon` result in the default value of ``rel_tol = 0.005``.
+
       This routine will be called by :c:func:`SUNDomEigEstimator_SetOptions`
       when using the key "Did.rel_tol".
 
@@ -388,14 +475,25 @@ dominant eigenvalue estimator.  *All routines are optional.*
 
 .. _SUNDomEigEst.SUNSuppliedFn:
 
-Functions provided by SUNDIALS packages
+User or SUNDIALS package provided functions
 ---------------------------------------------
 
 To interface with SUNDomEigEst modules, the SUNDIALS packages supply a
-:c:type:`SUNATimesFn` function for evaluating the matrix-vector product. This
-package-provided routine translates between the user-supplied ODE or DAE systems
-and the generic dominant eigenvalue estimator API. The function types for these
-routines are defined in the header file ``sundials/sundials_iterative.h``.
+:c:type:`SUNATimesFn` function for evaluating the matrix-vector 
+product. This package-provided routine translates between the user-supplied ODE or DAE 
+systems and the generic dominant eigenvalue estimator API. The function type :c:type:`SUNATimesFn` is defined in the header file ``sundials/sundials_iterative.h``.
+
+Users who wish to use a SUNDomEigEst module in "standalone" mode, however, must  provide either a :c:type:`SUNATimesFn` or a :c:type:`SUNRhsFn`, as described below.
+
+
+.. c:type:: int (*SUNRhsFn)(sunrealtype t, N_Vector y, N_Vector ydot, void* rhs_data)
+
+   Used to compute the right-hand side of an ODE system. This function is used 
+   when the dominant eigenvalue estimator is configured to perform a discrete 
+   Jacobian-vector product using quotient approximations of the Jacobian. The parameter
+   *rhs_data* is a pointer to any information about RHS which the function needs in order 
+   to do its job. The time parameter :math:`t` and the vector :math:`y` should be left 
+   unchanged.
 
 
 .. _SUNDomEigEst.Generic:

doc/shared/sundomeigest/SUNDomEigEst_Arnoldi.rst

@@ -116,15 +116,24 @@ The SUNDomEigEstimator_Arnoldi module defines the *content* field of a
      void* ATdata;
      N_Vector* V;
      N_Vector q;
+     N_Vector rhs_linY;
+     N_Vector Fy;
+     N_Vector work;
      int kry_dim;
      int num_warmups;
      long int num_iters;
+     sunbooleantype warmup_to_tol;
+     sunrealtype tol_preprocess;
+     sunrealtype rhs_linT;
      long int num_ATimes;
+     SUNRhsFn rhsfn;
+     void* rhs_data;
+     long int nfevals;
      sunrealtype* LAPACK_A;
      sunrealtype* LAPACK_wr;
      sunrealtype* LAPACK_wi;
      sunrealtype* LAPACK_work;
-     snuindextype LAPACK_lwork;
+     sunindextype LAPACK_lwork;
      sunrealtype** LAPACK_arr;
      sunrealtype** Hes;
    };
@@ -137,7 +146,7 @@ information:
 
 * ``ATData`` - pointer to structure for ``ATimes``,
 
-* ``V, q``   - vectors used for workspace by the Arnoldi algorithm.
+* ``V, q, Fy, work``   - vectors used for workspace.
 
 * ``kry_dim`` - dimension of Krylov subspaces (default is 3),
 
@@ -146,6 +155,21 @@ information:
 * ``num_iters`` - number of iterations (preprocessing and estimation) in the
   last :c:func:`SUNDomEigEstimator_Estimate` call,
 
+* ``warmup_to_tol`` - enable warmup iterations (default is ``SUNFALSE``)
+
+* ``tol_preprocess`` - tolerance for preprocessing iterations (default is 0.005; 
+  only used if ``warmup_to_tol`` is ``SUNTRUE``),
+
+* ``rhs_linY`` - state vector for linearization point,
+
+* ``rhs_linT`` - time value for linearization point,
+
+* ``rhsfn`` - user provided RHS function,
+
+* ``rhs_data`` - pointer to the data structure for ``rhsfn``,
+
+* ``nfevals`` - number of RHS evaluations,
+
 * ``num_ATimes`` - number of calls to the ``ATimes`` function,
 
 * ``LAPACK_A, LAPACK_wr, LAPACK_wi, LAPACK_work`` - ``sunrealtype`` used for workspace by LAPACK,