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::
+
+      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 ellipse or the exact stability region for 
+      stability checks. The stability region check for RKC and RKL methods is performed with the dominant
+      eigenvalue and the current step size to ensure stability. While this is sufficient for a stability
+      region of disk, the stability region of RKC and RKL methods is not a disk but rather a more complicated 
+      region that can be approximated by an inscribed ellipse. By default, the ellipse is used for stability 
+      checks, which is a conservative approximation of the stability region that possibly reduces the step 
+      sizes to ensure stability. Setting use_analytic_stab_region to SUNTRUE allows the use of the exact 
+      stability region, which can potentially allow for larger step sizes but possibly cause stability failures 
+      for the second dominant eigenvalue since it might be outside of the stability region even if the dominant 
+      eigenvalue is inside the stability region. Using ellipse for stability checks can be beneficial when two 
+      dominant eigenvalues are close to the stability boundary. Nevertheless, unless the full spectrum is used 
+      for stability checks, there is always a risk of stability failures one way or another. Thus, users 
+      have the option to choose their preference for a more conservative or more aggressive 
+      approach to stability. This input is only used for RKC and RKL methods.
+
+      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,19 @@ Changelog
 Changes to SUNDIALS in release 7.7.0
 ====================================
 
+.. SED_REPLACEMENT_KEY
+
+.. _Changelog.7.7.0:
+
+Changes to SUNDIALS in release 7.7.0
+====================================
+
+**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.6.0:

doc/shared/RecentChanges.rst

@@ -62,6 +62,12 @@ an eigenvalue estimate is requested on the first step in a subsequent call to
 :c:func:`ARKodeEvolve` unless the output vector passed contained the most
 recently returned solution.
 
+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.
+
+Fixed a bug in logging output from ARKODE, where for some time stepping modules,
+the current "time" output in the logger was incorrect.
 Fixed a potential bug in LSRKStep's :c:enumerator:`ARKODE_LSRK_SSP_S_3` method,
 where a real number was used instead of an integer, potentially resulting in a
 rounding error.

doc/shared/sundomeigest/SUNDomEigEst_API.rst

@@ -209,6 +209,44 @@ instead of supplying a dummy routine.
       A :c:type:`SUNErrCode`.
 
 
+.. c:function:: SUNErrCode SUNDomEigEstimator_SetRhs(SUNDomEigEstimator DEE, void* rhs_data, SUNRhsFn RHSfn)
+
+   This *optional* function provides a :c:type:`SUNRhsFn` function for performing
+   evaluations of the right-hand side function, as well as a ``void*`` pointer to a data structure
+   used by this routine, to the dominant eigenvalue estimator. This function allows the estimator to
+   perform a discrete Jacobian-vector product using quotient approximations of the Jacobian and estimate
+   dominant eigenvalues of the Jacobian without requiring the user to provide a matrix-vector product 
+   function. This can be useful when the user does not have access to the Jacobian or its matrix-vector 
+   product, but can evaluate the right-hand side function.
+
+   **Arguments:**
+
+      * *DEE* -- a SUNDomEigEstimator object.
+      * *rhs_data* -- pointer to structure for ``RHSfn``.
+      * *RHSfn* -- function pointer to perform right-hand side evaluations.
+
+   **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`. This allows the estimator to perform a discrete Jacobian-vector product using 
+   quotient approximations of the Jacobian at a specified linearization point.
+
+   **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 +282,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 +298,59 @@ instead of supplying a dummy routine.
 
    .. note::
 
+      The relative tolerance is used as a stopping criterion for the Power Iteration method. Specifically,
+      it defines the acceptable relative change between successive dominant eigenvalue estimates. It also
+      serves as a threshold for determining whether the dominant eigenvalue is real or complex.
+
+      A threshold, with :math:`\varepsilon` representing the machine precision,
+
+      .. math::
+         \mathtt{gram\_det\_tol} = 10 \cdot \max\left(\varepsilon,\; \mathtt{rel\_tol}\right)
+
+      is used to assess the numerical rank of the 2×2 Gram matrix formed by the current and previous
+      iterates in the Power Iteration method. If the determinant of this Gram matrix is less than or equal to
+      this threshold, the iterates are considered nearly linearly dependent, and the dominant eigenvalue is
+      treated as real.
+
+      To ensure that the Power Iteration method captures small imaginary parts of the dominant eigenvalue,
+      the relative tolerance 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.
+
+      When used in combination with Arnoldi Iteration, this routine sets the stopping criterion for 
+      the preprocessing Power Iteration phase. In this workflow, Power Iteration is first applied to 
+      perform an initial convergence check for the dominant eigenvalue estimate. The convergence
+      criterion is based on the change in magnitude between successive estimates and is defined as
+
+      .. math::
+
+         \left|\,|\lambda_{k}| - |\lambda_{k-1}|\,\right|
+         \le \mathtt{rel\_tol} \cdot |\lambda_{k}|.
+
+      The implementation performs this inexpensive preprocessing check using only the magnitude of 
+      the eigenvalue estimates. Arnoldi Iteration is executed only after this convergence criterion is 
+      satisfied. This approach ensures that Arnoldi is performed only once, rather than repeatedly 
+      running Arnoldi and checking convergence based on its complex-valued results. By relying on 
+      the cheaper magnitude-based preprocessing step, the routine avoids multiple Arnoldi runs that would
+      yield only marginal improvements in accuracy while incurring significantly higher computational cost.
+
+      When this routine is used in combination with Power Iteration, for ``rel_tol <= 0`` or 
+      ``rel_tol >= (1 - \varepsilon)``, a default value of ``rel_tol = 0.005`` is applied.
+
+      In the case of Arnoldi Iteration, this routine sets ``rel_tol = 0.005`` for the preprocessing 
+      Power Iteration phase for ``rel_tol = 0`` or ``rel_tol >= (1 - \varepsilon)``, and disables 
+      preprocessing when ``rel_tol < 0``.
+
       This routine will be called by :c:func:`SUNDomEigEstimator_SetOptions`
       when using the key "Did.rel_tol".
 
@@ -391,13 +483,23 @@ dominant eigenvalue estimator.  *All routines are optional.*
 Functions provided by SUNDIALS packages
 ---------------------------------------------
 
-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
+To interface with SUNDomEigEst modules, the SUNDIALS packages supply
+:c:type:`SUNATimesFn` and :c:type:`SUNRHSFn` functions 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:: int (*SUNRHSFn)(sunrealtype t, N_Vector y, N_Vector ydot, void* rhs_data)
+
+   Used to compute the right-hand side of an ODE or DAE 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:
 
 The generic SUNDomEigEstimator module

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`` - type of 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 (default is 0.0),
+
+* ``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,