Add SUNNonlinearSolver_Auto module that automatically switches between newton and fixed-point based on stiffness

CHANGELOG.md

@@ -6,6 +6,13 @@
 
 ### New Features and Enhancements
 
+We added a new SUNNonlinearSolver implementation, `SUNNonlinearSolver_Auto`,
+which uses an algorithm described in https://doi.org/10.1007/BF01933714 to
+switch between modified Newton iteration and fixed-point iteration based on
+an estimate of stiffness. This solver may be useful to pair with the BDF method
+in CVODE/CVODES for users who are unsure about the stiffness of their problem.
+See the module documentation for more information.
+
 Updated the Kokkos N_Vector to support Kokkos 5.x versions.
 
 ### Bug Fixes

bindings/sundials4py/CMakeLists.txt

@@ -111,6 +111,7 @@ set(sundials_SOURCES
     sunmatrix/sunmatrix_dense.cpp
     sunmatrix/sunmatrix_sparse.cpp
     sunmemory/sunmemory_system.cpp
+    sunnonlinsol/sunnonlinsol_auto.cpp
     sunnonlinsol/sunnonlinsol_fixedpoint.cpp
     sunnonlinsol/sunnonlinsol_newton.cpp
     test/sundials4py_test.cpp)
@@ -140,6 +141,9 @@ target_link_libraries(
           sundials_sunlinsolspfgmr
           sundials_sunlinsolsptfqmr
           sundials_sunlinsolpcg
+          sundials_sunnonlinsolauto
+          sundials_sunnonlinsolfixedpoint
+          sundials_sunnonlinsolnewton
           sundials_sunmatrixband
           sundials_sunmatrixdense
           sundials_sunmatrixsparse

bindings/sundials4py/sundials/sundials_nonlinearsolver_generated.hpp

@@ -9,6 +9,7 @@ auto pyEnumSUNNonlinearSolver_Type =
                                      nb::is_arithmetic(), "")
     .value("SUNNONLINEARSOLVER_ROOTFIND", SUNNONLINEARSOLVER_ROOTFIND, "")
     .value("SUNNONLINEARSOLVER_FIXEDPOINT", SUNNONLINEARSOLVER_FIXEDPOINT, "")
+    .value("SUNNONLINEARSOLVER_HYBRID", SUNNONLINEARSOLVER_HYBRID, "")
     .export_values();
 // #ifndef SWIG
 //
@@ -87,8 +88,26 @@ m.def(
     return SUNNonlinSolGetNumConvFails_adapt_modifiable_immutable_to_return(NLS);
   },
   nb::arg("NLS"));
+
+m.def(
+  "SUNNonlinSolGetDeltaNorm",
+  [](SUNNonlinearSolver NLS) -> std::tuple<SUNErrCode, sunrealtype>
+  {
+    auto SUNNonlinSolGetDeltaNorm_adapt_modifiable_immutable_to_return =
+      [](SUNNonlinearSolver NLS) -> std::tuple<SUNErrCode, sunrealtype>
+    {
+      sunrealtype delnrm_adapt_modifiable;
+
+      SUNErrCode r = SUNNonlinSolGetDeltaNorm(NLS, &delnrm_adapt_modifiable);
+      return std::make_tuple(r, delnrm_adapt_modifiable);
+    };
+
+    return SUNNonlinSolGetDeltaNorm_adapt_modifiable_immutable_to_return(NLS);
+  },
+  nb::arg("NLS"));
 m.attr("SUN_NLS_CONTINUE")   = +901;
 m.attr("SUN_NLS_CONV_RECVR") = +902;
+m.attr("SUN_NLS_SWITCH")     = +903;
 // #ifdef __cplusplus
 //
 // #endif

bindings/sundials4py/sundials4py.cpp

@@ -65,6 +65,7 @@ void bind_sunmatrix_sparse(nb::module_& m);
 
 void bind_sunnonlinsol_fixedpoint(nb::module_& m);
 void bind_sunnonlinsol_newton(nb::module_& m);
+void bind_sunnonlinsol_auto(nb::module_& m);
 
 } // namespace sundials4py
 
@@ -143,4 +144,5 @@ NB_MODULE(sundials4py, m)
 
   sundials4py::bind_sunnonlinsol_fixedpoint(core_m);
   sundials4py::bind_sunnonlinsol_newton(core_m);
+  sundials4py::bind_sunnonlinsol_auto(core_m);
 }

bindings/sundials4py/sunnonlinsol/generate.yaml

@@ -53,3 +53,10 @@ modules:
     fn_exclude_by_name__regex:
     # We do not support getting a function pointer back to Python
     - "^SUNNonlinSolGetSysFn_Newton$"
+  sunnonlinsol_auto:
+      path: sunnonlinsol/sunnonlinsol_auto_generated.hpp
+      headers:
+      - ../../include/sunnonlinsol/sunnonlinsol_auto.h
+      fn_exclude_by_name__regex:
+      # We do not support getting a function pointer back to Python
+      - "^SUNNonlinSolGetSysFn_Auto$"

bindings/sundials4py/sunnonlinsol/sunnonlinsol_auto.cpp

@@ -0,0 +1,32 @@
+/* -----------------------------------------------------------------
+ * Programmer(s): Cody J. Balos @ LLNL
+ * -----------------------------------------------------------------
+ * SUNDIALS Copyright Start
+ * Copyright (c) 2025-2025, Lawrence Livermore National Security,
+ * University of Maryland Baltimore County, and the SUNDIALS contributors.
+ * Copyright (c) 2013-2025, Lawrence Livermore National Security
+ * and Southern Methodist University.
+ * Copyright (c) 2002-2013, Lawrence Livermore National Security.
+ * All rights reserved.
+ *
+ * See the top-level LICENSE and NOTICE files for details.
+ *
+ * SPDX-License-Identifier: BSD-3-Clause
+ * SUNDIALS Copyright End
+ * -----------------------------------------------------------------*/
+
+#include <sundials/sundials_nonlinearsolver.hpp>
+#include <sunnonlinsol/sunnonlinsol_auto.h>
+#include "sundials4py.hpp"
+
+namespace nb = nanobind;
+using namespace sundials::experimental;
+
+namespace sundials4py {
+
+void bind_sunnonlinsol_auto(nb::module_& m)
+{
+#include "sunnonlinsol_auto_generated.hpp"
+}
+
+} // namespace sundials4py

bindings/sundials4py/sunnonlinsol/sunnonlinsol_auto_generated.hpp

@@ -0,0 +1,99 @@
+// #ifndef SUNDIALS_NONLINSOLAUTO_H_
+//
+// #ifdef __cplusplus
+//
+// #endif
+//
+
+auto pyEnumSUNNonlinSolAutoType =
+  nb::enum_<SUNNonlinSolAutoType>(m, "SUNNonlinSolAutoType",
+                                  nb::is_arithmetic(), "")
+    .value("SUNNONLINSOL_AUTO_FIXEDPOINT", SUNNONLINSOL_AUTO_FIXEDPOINT, "")
+    .value("SUNNONLINSOL_AUTO_NEWTON", SUNNONLINSOL_AUTO_NEWTON, "")
+    .export_values();
+// #ifndef SWIG
+//
+// #endif
+//
+
+auto pyClass_SUNNonlinearSolverContent_Auto =
+  nb::class_<_SUNNonlinearSolverContent_Auto>(m,
+                                              "_SUNNonlinearSolverContent_Auto",
+                                              "")
+    .def(nb::init<>()) // implicit default constructor
+  ;
+
+m.def(
+  "SUNNonlinSol_Auto",
+  [](N_Vector y, int m, SUNNonlinSolAutoType active_solver_type,
+     SUNContext sunctx) -> std::shared_ptr<std::remove_pointer_t<SUNNonlinearSolver>>
+  {
+    auto SUNNonlinSol_Auto_adapt_return_type_to_shared_ptr =
+      [](N_Vector y, int m, SUNNonlinSolAutoType active_solver_type,
+         SUNContext sunctx)
+      -> std::shared_ptr<std::remove_pointer_t<SUNNonlinearSolver>>
+    {
+      auto lambda_result = SUNNonlinSol_Auto(y, m, active_solver_type, sunctx);
+
+      return our_make_shared<std::remove_pointer_t<SUNNonlinearSolver>,
+                             SUNNonlinearSolverDeleter>(lambda_result);
+    };
+
+    return SUNNonlinSol_Auto_adapt_return_type_to_shared_ptr(y, m,
+                                                             active_solver_type,
+                                                             sunctx);
+  },
+  nb::arg("y"), nb::arg("m"), nb::arg("active_solver_type"), nb::arg("sunctx"),
+  "nb::keep_alive<0, 4>()", nb::keep_alive<0, 4>());
+
+m.def("SUNNonlinSolSetSysFns_Auto", SUNNonlinSolSetSysFns_Auto, nb::arg("NLS"),
+      nb::arg("root_fn"), nb::arg("fixed_point_fn"));
+
+m.def("SUNNonlinSolSetSwitchingParameters_Auto",
+      SUNNonlinSolSetSwitchingParameters_Auto, nb::arg("NLS"),
+      nb::arg("newt_to_fp_threshold"), nb::arg("newt_to_fp_delay"),
+      nb::arg("fp_to_newt_threshold"), nb::arg("fp_to_newt_delay"));
+
+m.def(
+  "SUNNonlinSolGetNumItersByType_Auto",
+  [](SUNNonlinearSolver NLS) -> std::tuple<SUNErrCode, long, long>
+  {
+    auto SUNNonlinSolGetNumItersByType_Auto_adapt_modifiable_immutable_to_return =
+      [](SUNNonlinearSolver NLS) -> std::tuple<SUNErrCode, long, long>
+    {
+      long fp_iters_adapt_modifiable;
+      long newt_iters_adapt_modifiable;
+
+      SUNErrCode r =
+        SUNNonlinSolGetNumItersByType_Auto(NLS, &fp_iters_adapt_modifiable,
+                                           &newt_iters_adapt_modifiable);
+      return std::make_tuple(r, fp_iters_adapt_modifiable,
+                             newt_iters_adapt_modifiable);
+    };
+
+    return SUNNonlinSolGetNumItersByType_Auto_adapt_modifiable_immutable_to_return(
+      NLS);
+  },
+  nb::arg("NLS"));
+
+m.def(
+  "SUNNonlinSolGetDeltaNorm_Auto",
+  [](SUNNonlinearSolver NLS) -> std::tuple<SUNErrCode, sunrealtype>
+  {
+    auto SUNNonlinSolGetDeltaNorm_Auto_adapt_modifiable_immutable_to_return =
+      [](SUNNonlinearSolver NLS) -> std::tuple<SUNErrCode, sunrealtype>
+    {
+      sunrealtype delnrm_adapt_modifiable;
+
+      SUNErrCode r = SUNNonlinSolGetDeltaNorm_Auto(NLS, &delnrm_adapt_modifiable);
+      return std::make_tuple(r, delnrm_adapt_modifiable);
+    };
+
+    return SUNNonlinSolGetDeltaNorm_Auto_adapt_modifiable_immutable_to_return(NLS);
+  },
+  nb::arg("NLS"));
+// #ifdef __cplusplus
+//
+// #endif
+//
+// #endif

bindings/sundials4py/sunnonlinsol/sunnonlinsol_fixedpoint_generated.hpp

@@ -53,6 +53,48 @@ m.def(
 
 m.def("SUNNonlinSolSetDamping_FixedPoint", SUNNonlinSolSetDamping_FixedPoint,
       nb::arg("NLS"), nb::arg("beta"));
+
+m.def("SUNNonlinSolSetCrateConstant_FixedPoint",
+      SUNNonlinSolSetCrateConstant_FixedPoint, nb::arg("NLS"),
+      nb::arg("crate_const"));
+
+m.def(
+  "SUNNonlinSolGetConvRate_FixedPoint",
+  [](SUNNonlinearSolver NLS) -> std::tuple<SUNErrCode, sunrealtype>
+  {
+    auto SUNNonlinSolGetConvRate_FixedPoint_adapt_modifiable_immutable_to_return =
+      [](SUNNonlinearSolver NLS) -> std::tuple<SUNErrCode, sunrealtype>
+    {
+      sunrealtype crate_adapt_modifiable;
+
+      SUNErrCode r =
+        SUNNonlinSolGetConvRate_FixedPoint(NLS, &crate_adapt_modifiable);
+      return std::make_tuple(r, crate_adapt_modifiable);
+    };
+
+    return SUNNonlinSolGetConvRate_FixedPoint_adapt_modifiable_immutable_to_return(
+      NLS);
+  },
+  nb::arg("NLS"));
+
+m.def(
+  "SUNNonlinSolGetDeltaNorm_FixedPoint",
+  [](SUNNonlinearSolver NLS) -> std::tuple<SUNErrCode, sunrealtype>
+  {
+    auto SUNNonlinSolGetDeltaNorm_FixedPoint_adapt_modifiable_immutable_to_return =
+      [](SUNNonlinearSolver NLS) -> std::tuple<SUNErrCode, sunrealtype>
+    {
+      sunrealtype delnrm_adapt_modifiable;
+
+      SUNErrCode r =
+        SUNNonlinSolGetDeltaNorm_FixedPoint(NLS, &delnrm_adapt_modifiable);
+      return std::make_tuple(r, delnrm_adapt_modifiable);
+    };
+
+    return SUNNonlinSolGetDeltaNorm_FixedPoint_adapt_modifiable_immutable_to_return(
+      NLS);
+  },
+  nb::arg("NLS"));
 // #ifdef __cplusplus
 //
 // #endif

bindings/sundials4py/sunnonlinsol/sunnonlinsol_newton_generated.hpp

@@ -50,6 +50,25 @@ m.def(
   },
   nb::arg("count"), nb::arg("y"), nb::arg("sunctx"), "nb::keep_alive<0, 3>()",
   nb::keep_alive<0, 3>());
+
+m.def(
+  "SUNNonlinSolGetDeltaNorm_Newton",
+  [](SUNNonlinearSolver NLS) -> std::tuple<SUNErrCode, sunrealtype>
+  {
+    auto SUNNonlinSolGetDeltaNorm_Newton_adapt_modifiable_immutable_to_return =
+      [](SUNNonlinearSolver NLS) -> std::tuple<SUNErrCode, sunrealtype>
+    {
+      sunrealtype delnrm_adapt_modifiable;
+
+      SUNErrCode r = SUNNonlinSolGetDeltaNorm_Newton(NLS,
+                                                     &delnrm_adapt_modifiable);
+      return std::make_tuple(r, delnrm_adapt_modifiable);
+    };
+
+    return SUNNonlinSolGetDeltaNorm_Newton_adapt_modifiable_immutable_to_return(
+      NLS);
+  },
+  nb::arg("NLS"));
 // #ifdef __cplusplus
 //
 // #endif

doc/arkode/guide/source/sunnonlinsol/SUNNonlinSol_links.rst

@@ -15,4 +15,5 @@
 
 .. include:: ../../../../shared/sunnonlinsol/SUNNonlinSol_Newton.rst
 .. include:: ../../../../shared/sunnonlinsol/SUNNonlinSol_FixedPoint.rst
+.. include:: ../../../../shared/sunnonlinsol/SUNNonlinSol_Auto.rst
 .. include:: ../../../../shared/sunnonlinsol/SUNNonlinSol_PetscSNES.rst

doc/cvode/guide/source/sunnonlinsol/SUNNonlinSol_links.rst

@@ -15,4 +15,5 @@
 
 .. include:: ../../../../shared/sunnonlinsol/SUNNonlinSol_Newton.rst
 .. include:: ../../../../shared/sunnonlinsol/SUNNonlinSol_FixedPoint.rst
+.. include:: ../../../../shared/sunnonlinsol/SUNNonlinSol_Auto.rst
 .. include:: ../../../../shared/sunnonlinsol/SUNNonlinSol_PetscSNES.rst

doc/cvodes/guide/source/sunnonlinsol/SUNNonlinSol_links.rst

@@ -15,4 +15,5 @@
 
 .. include:: ../../../../shared/sunnonlinsol/SUNNonlinSol_Newton.rst
 .. include:: ../../../../shared/sunnonlinsol/SUNNonlinSol_FixedPoint.rst
+.. include:: ../../../../shared/sunnonlinsol/SUNNonlinSol_Auto.rst
 .. include:: ../../../../shared/sunnonlinsol/SUNNonlinSol_PetscSNES.rst

doc/ida/guide/source/sunnonlinsol/SUNNonlinSol_links.rst

@@ -15,4 +15,5 @@
 
 .. include:: ../../../../shared/sunnonlinsol/SUNNonlinSol_Newton.rst
 .. include:: ../../../../shared/sunnonlinsol/SUNNonlinSol_FixedPoint.rst
+.. include:: ../../../../shared/sunnonlinsol/SUNNonlinSol_Auto.rst
 .. include:: ../../../../shared/sunnonlinsol/SUNNonlinSol_PetscSNES.rst

doc/idas/guide/source/sunnonlinsol/SUNNonlinSol_links.rst

@@ -15,4 +15,5 @@
 
 .. include:: ../../../../shared/sunnonlinsol/SUNNonlinSol_Newton.rst
 .. include:: ../../../../shared/sunnonlinsol/SUNNonlinSol_FixedPoint.rst
+.. include:: ../../../../shared/sunnonlinsol/SUNNonlinSol_Auto.rst
 .. include:: ../../../../shared/sunnonlinsol/SUNNonlinSol_PetscSNES.rst

doc/shared/RecentChanges.rst

@@ -5,6 +5,13 @@
 
 **New Features and Enhancements**
 
+We added a new SUNNonlinearSolver implementation,
+:ref:`SUNNonlinearSolver_Auto <SUNNonlinSol.Auto>`, which uses an algorithm described in
+:cite:p:`norsett1986switching` to switch between modified Newton
+iteration and fixed-point iteration based on an estimate of stiffness. This
+solver may be useful to pair with the BDF method in CVODE/CVODES for users who
+are unsure about the stiffness of their problem.
+
 Updated the Kokkos N_Vector to support Kokkos 5.x versions.
 
 **Bug Fixes**

doc/shared/sundials.bib

@@ -2677,3 +2677,19 @@ @article{vonmises29
   doi = {10.1002/zamm.19290090206},
   year={1929}
 }
+
+%
+% Nonlinear solver switching
+%
+
+@article{norsett1986switching,
+  title={Switching between modified Newton and fix-point iteration for implicit {ODE}-solvers},
+  author={N{\o}rsett, Syvert P and Thomsen, Per G},
+  journal={BIT Numerical Mathematics},
+  volume={26},
+  number={3},
+  pages={339--348},
+  year={1986},
+  publisher={Springer},
+  doi={10.1007/BF01933714}
+}

doc/shared/sunnonlinsol/SUNNonlinSol_API.rst

@@ -381,6 +381,25 @@ linear solver module; otherwise :c:func:`SUNNonlinSolGetCurIter` is optional.
       * A :c:type:`SUNErrCode`
 
 
+.. c:function:: SUNErrCode SUNNonlinSolGetDeltaNorm(SUNNonlinearSolver NLS, sunrealtype *delnrm)
+
+   This *optional* function returns the norm of the most recent nonlinear
+   solver update (often denoted :math:`\|\delta\|}`) computed by the
+   nonlinear solver. The norm used is left up to the nonlinear solver implementation,
+   but typically it is the WRMS norm.
+
+   **Arguments:**
+      * *NLS* -- a SUNNonlinSol object.
+      * *delnrm* -- the update norm.
+
+   **Return value:**
+      * A :c:type:`SUNErrCode`
+
+   **Notes:**
+      If the nonlinear solver does not implement this operation, the return
+      value will be :c:macro:`SUN_ERR_NOT_IMPLEMENTED`.
+
+
 .. _SUNNonlinSol.API.SUNSuppliedFn:
 
 Functions provided by SUNDIALS integrators