Merge branch 'master' into min_timestep_fix

Author: ChrisRackauckas

.github/workflows/CI.yml

@@ -15,7 +15,6 @@ jobs:
           - Core
         version:
           - '1'
-          - '1.6'
     steps:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v1

Project.toml

@@ -1,7 +1,7 @@
 name = "Sundials"
 uuid = "c3572dad-4567-51f8-b174-8c6c989267f4"
 authors = ["Chris Rackauckas <accounts@chrisrackauckas.com>"]
-version = "4.20.0"
+version = "4.23.1"
 
 [deps]
 CEnum = "fa961155-64e5-5f13-b03f-caf6b980ea82"
@@ -17,14 +17,14 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Sundials_jll = "fb77eaff-e24c-56d4-86b1-d163f2edb164"
 
 [compat]
-CEnum = "0.2, 0.3, 0.4"
+CEnum = "0.5"
 DataStructures = "0.18"
 DiffEqBase = "6.122"
 PrecompileTools = "1"
-Reexport = "0.2, 1.0"
-SciMLBase = "1.92, 2"
+Reexport = "1.0"
+SciMLBase = "2.9"
 Sundials_jll = "5.2"
-julia = "1.6"
+julia = "1.9"
 
 [extras]
 AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c"

gen/generate.jl

@@ -86,7 +86,7 @@ function wrap_sundials_api(expr::Expr)
             if occursin(r"UserDataB?$", func_name)
                 # replace Ptr{Void} with Any to allow passing Julia objects through user data
                 for (i, arg_expr) in enumerate(expr.args[2].args[1].args)
-                    if !(typeof(arg_expr) <: Symbol) &&
+                    if !(arg_expr isa Symbol) &&
                        arg_expr.args[1] in values(ctor_return_type)
                         if arg_expr.args[2] == :(Ptr{Cvoid})
                             arg_expr.args[2] = Any
@@ -96,7 +96,7 @@ function wrap_sundials_api(expr::Expr)
                     end
                 end
             end
-            if !(typeof(expr) <: Symbol) && length(expr.args) > 1 &&
+            if !(expr isa Symbol) && length(expr.args) > 1 &&
                (expr.args[2].args[1].args[2].args[2] == :libsundials_sunlinsol ||
                 expr.args[2].args[1].args[2].args[2] == :libsundials_sunmatrix ||
                 expr.args[2].args[1].args[2].args[2] == :libsundials_sunnonlinsol)
@@ -124,7 +124,7 @@ function wrap_sundials_api(expr::Expr)
                 #   2) expr for local var definition, nothing if not required
                 #   3) expr for low-level wrapper call
                 # if 1)==3), then no wrapping is required
-                if typeof(expr) <: Symbol
+                if expr isa Symbol
                     arg_name_expr = expr
                     arg_type_expr = Any
                 else

src/Sundials.jl

@@ -14,8 +14,7 @@ import LinearAlgebra
 import Libdl
 using CEnum
 
-const warnkeywords = (:d_discontinuities,
-    :isoutofdomain,
+const warnkeywords = (:isoutofdomain,
     :unstable_check,
     :calck,
     :internalnorm,

src/common_interface/algorithms.jl

@@ -3,7 +3,7 @@
 # Abstract Types
 abstract type SundialsODEAlgorithm{Method, LinearSolver} <: DiffEqBase.AbstractODEAlgorithm end
 abstract type SundialsDAEAlgorithm{LinearSolver} <: DiffEqBase.AbstractDAEAlgorithm end
-abstract type SundialsNonlinearSolveAlgorithm{LinearSolver} end
+abstract type SundialsNonlinearSolveAlgorithm{LinearSolver} <: SciMLBase.AbstractNonlinearAlgorithm end
 
 SciMLBase.alg_order(alg::Union{SundialsODEAlgorithm, SundialsDAEAlgorithm}) = alg.max_order
 
@@ -723,6 +723,9 @@ KINSOL(;
     jac_upper = 0,
     jac_lower = 0,
     userdata = nothing,
+    prec_side = 0,
+    krylov_dim = 0,
+    globalization_strategy = :None
 )
 ```
 
@@ -738,7 +741,6 @@ The choices for the linear solver are:
 - `:LapackBand`: A version of the banded linear solver that uses the Julia-provided
   OpenBLAS-linked LAPACK for multithreaded operations. This will be faster than
   `:Band` on larger systems but has noticeable overhead on smaller (<100 ODE) systems.
-- `:Diagonal`: This method is specialized for diagonal Jacobians.
 - `:GMRES`: A GMRES method. Recommended first choice Krylov method.
 - `:BCG`: A biconjugate gradient method
 - `:PCG`: A preconditioned conjugate gradient method. Only for symmetric
@@ -747,19 +749,30 @@ The choices for the linear solver are:
 - `:KLU`: A sparse factorization method. Requires that the user specify a
   Jacobian. The Jacobian must be set as a sparse matrix in the `ODEProblem`
   type.
+
+The choices for globalization strategy are:
+
+- `:None`: No globalization strategy
+- `:LineSearch`: A line search globalization strategy
 """
 struct KINSOL{LinearSolver} <: SundialsNonlinearSolveAlgorithm{LinearSolver}
     jac_upper::Int
     jac_lower::Int
     userdata::Any
+    prec_side::Int
+    krylov_dim::Int
+    globalization_strategy::Symbol
 end
+
 Base.@pure function KINSOL(;
     linear_solver = :Dense,
     jac_upper = 0,
     jac_lower = 0,
-    userdata = nothing)
+    userdata = nothing,
+    prec_side = 0,
+    krylov_dim = 0,
+    globalization_strategy = :None)
     if !(linear_solver in (:None,
-        :Diagonal,
         :Dense,
         :LapackDense,
         :Band,
@@ -772,9 +785,11 @@ Base.@pure function KINSOL(;
         :KLU))
         error("Linear solver choice not accepted.")
     end
-    KINSOL{linear_solver}(jac_upper,
-        jac_lower,
-        userdata)
+    if !(globalization_strategy in (:LineSearch, :None))
+        error("Globalization strategy not accepted.")
+    end
+    KINSOL{linear_solver}(jac_upper, jac_lower, userdata, prec_side, krylov_dim,
+    globalization_strategy)
 end
 
 method_choice(alg::SundialsODEAlgorithm{Method}) where {Method} = Method

src/common_interface/function_types.jl

@@ -146,7 +146,7 @@ function massmat(t::Float64,
     tmp1::N_Vector,
     tmp2::N_Vector,
     tmp3::N_Vector)
-    if typeof(mmf.mass_matrix) <: Array
+    if mmf.mass_matrix isa Array
         M = convert(Matrix, _M)
         M .= mmf.mass_matrix
     else

src/common_interface/integrator_types.jl

@@ -20,8 +20,10 @@ mutable struct DEOptions{SType, TstopType, SType2, TstopType2, SIX, CType, relto
     advance_to_tstop::Bool
     stop_at_next_tstop::Bool
     progress::Bool
+    progress_steps::Int
     progress_name::String
     progress_message::F5
+    progress_id::Symbol
     maxiters::Int
 end
 
@@ -172,6 +174,7 @@ mutable struct IDAIntegrator{N,
     tmp::Array{Float64, N}
     uprev::Array{Float64, N}
     flag::Cint
+    iter::Int
     just_hit_tstop::Bool
     event_last_time::Int
     vector_event_last_time::Int
@@ -197,6 +200,13 @@ function (integrator::IDAIntegrator)(out,
     integrator.flag = @checkflag IDAGetDky(integrator.mem, t, Cint(T), vec(out))
     return idxs === nothing ? out : @view out[idxs]
 end
+function (integrator::IDAIntegrator)(out::SubArray,
+    t::Number,
+    deriv::Type{Val{T}} = Val{0};
+    idxs = nothing) where {T}
+    throw(ArgumentError("Views are not supported with IDA!"))
+end
+
 
 ###  Error check (retcode)
 
@@ -217,7 +227,7 @@ DiffEqBase.postamble!(integrator::AbstractSundialsIntegrator) = nothing
             tstop = first(integrator.opts.tstops)
             set_stop_time(integrator, tstop)
             integrator.tprev = integrator.t
-            if !(typeof(integrator.opts.callback.continuous_callbacks) <: Tuple{})
+            if !(integrator.opts.callback.continuous_callbacks isa Tuple{})
                 integrator.uprev .= integrator.u
             end
             solver_step(integrator, tstop)
@@ -228,7 +238,7 @@ DiffEqBase.postamble!(integrator::AbstractSundialsIntegrator) = nothing
         end
     else
         integrator.tprev = integrator.t
-        if !(typeof(integrator.opts.callback.continuous_callbacks) <: Tuple{})
+        if !(integrator.opts.callback.continuous_callbacks isa Tuple{})
             integrator.uprev .= integrator.u
         end
         if !isempty(integrator.opts.tstops)

src/common_interface/integrator_utils.jl

@@ -6,7 +6,7 @@ function handle_callbacks!(integrator)
     continuous_modified = false
     discrete_modified = false
     saved_in_cb = false
-    if !(typeof(continuous_callbacks) <: Tuple{})
+    if !(continuous_callbacks isa Tuple{})
         time, upcrossing, event_occured, event_idx, idx, counter = DiffEqBase.find_first_continuous_callback(integrator,
             continuous_callbacks...)
         if event_occured
@@ -22,7 +22,7 @@ function handle_callbacks!(integrator)
             integrator.vector_event_last_time = 1
         end
     end
-    if !(typeof(discrete_callbacks) <: Tuple{})
+    if !(discrete_callbacks isa Tuple{})
         discrete_modified, saved_in_cb = DiffEqBase.apply_discrete_callback!(integrator,
             discrete_callbacks...)
     end
@@ -46,10 +46,9 @@ function DiffEqBase.savevalues!(integrator::AbstractSundialsIntegrator,
     uType = typeof(integrator.sol.prob.u0)
     # The call to first is an overload of Base.first implemented in DataStructures
     while !isempty(integrator.opts.saveat) &&
-        integrator.tdir * first(integrator.opts.saveat) < integrator.tdir * integrator.t
+        first(integrator.opts.saveat) <= integrator.tdir * integrator.t
         saved = true
-        curt = pop!(integrator.opts.saveat)
-
+        curt = integrator.tdir * pop!(integrator.opts.saveat)
         tmp = integrator(curt)
         save_value!(integrator.sol.u, tmp, uType,
             integrator.opts.save_idxs, false)
@@ -123,13 +122,13 @@ end
 function DiffEqBase.add_tstop!(integrator::AbstractSundialsIntegrator, t)
     integrator.tdir * (t - integrator.t) < 0 &&
         error("Tried to add a tstop that is behind the current time. This is strictly forbidden")
-    push!(integrator.opts.tstops, t)
+    push!(integrator.opts.tstops, integrator.tdir * t)
 end
 
 function DiffEqBase.add_saveat!(integrator::AbstractSundialsIntegrator, t)
     integrator.tdir * (t - integrator.t) < 0 &&
         error("Tried to add a saveat that is behind the current time. This is strictly forbidden")
-    push!(integrator.opts.saveat, t)
+    push!(integrator.opts.saveat, integrator.tdir * t)
 end
 
 DiffEqBase.get_tmp_cache(integrator::AbstractSundialsIntegrator) = (integrator.tmp,)