feat: allow simple newton steps

Author: avik-pal

Project.toml

@@ -1,7 +1,7 @@
 name = "Sundials"
 uuid = "c3572dad-4567-51f8-b174-8c6c989267f4"
 authors = ["Chris Rackauckas <[email protected]>"]
-version = "4.27.0"
+version = "4.28.0"
 
 [deps]
 CEnum = "fa961155-64e5-5f13-b03f-caf6b980ea82"

src/common_interface/algorithms.jl

@@ -726,6 +726,7 @@ KINSOL(;
     prec_side = 0,
     krylov_dim = 0,
     globalization_strategy = :None
+    maxsetupcalls=0
 )
 ```
 
@@ -754,6 +755,11 @@ The choices for globalization strategy are:
 
 - `:None`: No globalization strategy
 - `:LineSearch`: A line search globalization strategy
+
+Other options:
+
+- `maxsetupcalls`: Maximum number of nonlinear iterations that can be performed between
+  calls to the preconditioner or Jacobian setup function.
 """
 struct KINSOL{LinearSolver} <: SundialsNonlinearSolveAlgorithm{LinearSolver}
     jac_upper::Int
@@ -762,6 +768,7 @@ struct KINSOL{LinearSolver} <: SundialsNonlinearSolveAlgorithm{LinearSolver}
     prec_side::Int
     krylov_dim::Int
     globalization_strategy::Symbol
+    maxsetupcalls::Int
 end
 
 Base.@pure function KINSOL(;
@@ -771,7 +778,8 @@ Base.@pure function KINSOL(;
     userdata = nothing,
     prec_side = 0,
     krylov_dim = 0,
-    globalization_strategy = :None)
+    globalization_strategy = :None,
+    maxsetupcalls = 0)
     if !(linear_solver in (:None,
         :Dense,
         :LapackDense,
@@ -789,7 +797,7 @@ Base.@pure function KINSOL(;
         error("Globalization strategy not accepted.")
     end
     KINSOL{linear_solver}(jac_upper, jac_lower, userdata, prec_side, krylov_dim,
-    globalization_strategy)
+    globalization_strategy, maxsetupcalls)
 end
 
 method_choice(alg::SundialsODEAlgorithm{Method}) where {Method} = Method

src/common_interface/solve.jl

@@ -76,7 +76,8 @@ function DiffEqBase.__solve(prob::Union{
         alg.prec_side,
         alg.krylov_dim,
         maxiters,
-        strategy = alg.globalization_strategy)
+        strategy = alg.globalization_strategy,
+        alg.maxsetupcalls)
 
     f!(resid, u)
     retcode = interpret_sundials_retcode(flag)

src/simple.jl

@@ -55,7 +55,8 @@ function ___kinsol(f,
     krylov_dim::Int = 0,
     jac_prototype = nothing,
     maxiters = 1000,
-    strategy = :None)
+    strategy = :None,
+    maxsetupcalls=0)
     # f, Function to be optimized of the form f(y::Vector{Float64}, fy::Vector{Float64})
     #    where `y` is the input vector, and `fy` is the result of the function
     # y0, Vector of initial values
@@ -111,6 +112,7 @@ function ___kinsol(f,
     flag = @checkflag KINSetLinearSolver(kmem, LS, A) true
     flag = @checkflag KINSetUserData(kmem, userfun) true
     flag = @checkflag KINSetNumMaxIters(kmem, maxiters) true
+    flag = @checkflag KINSetMaxSetupCalls(kmem, maxsetupcalls) true
     ## Solve problem
     scale = ones(length(y0))
     if strategy == :None

test/kinsol_nonlinear_solve.jl

@@ -36,6 +36,15 @@ prob_oop = NonlinearProblem{false}(f_oop, u0)
     du = zeros(2)
     f_oop(sol.u, nothing)
     @test maximum(abs, du) < 1e-6
+
+    # Pure Newton Steps
+    alg = KINSOL(; linear_solver, globalization_strategy, maxsetupcalls = 1)
+    sol = solve(prob_oop, alg; abstol)
+    @test SciMLBase.successful_retcode(sol.retcode)
+
+    du = zeros(2)
+    f_oop(sol.u, nothing)
+    @test maximum(abs, du) < 1e-6
 end
 
 # Scalar