Support SciMLBase v3 and DiffEqBase v7

- Widen compat for SciMLBase to "2, 3" and DiffEqBase to "6.62, 7".
- Bump version to 1.2.0.
- Accept `DEVerbosity` objects for the `verbose` kwarg by treating
  any non-`false` value as "show warnings". DiffEqBase v7 now passes
  a structured verbosity object instead of a plain `Bool`.
- Unwrap `SciMLBase.unwrapped_f(prob.f.f)` (and `prob.f.f1.f` for
  second-order problems) before invoking the user function.
  SciMLBase v3 changed the default specialization of `ODEFunction{iip}`
  to `AutoSpecialize`, which wraps the function in a
  FunctionWrappersWrapper that only accepts the exact argument types
  captured at problem construction (e.g. `Matrix{Float64}`). Passing
  a `Base.ReshapedArray` through the wrapper fails; using the unwrapped
  function avoids the type-specialization trap.

Supersedes #60 and #61.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>

Author: ChrisRackauckas

Project.toml

@@ -1,6 +1,6 @@
 name = "GeometricIntegratorsDiffEq"
 uuid = "5a33fad7-5ce4-5983-9f5d-5f26ceab5c96"
-version = "1.1.0"
+version = "1.2.0"
 authors = ["Chris Rackauckas <accounts@chrisrackauckas.com>"]
 
 [deps]
@@ -10,11 +10,11 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 
 [compat]
-DiffEqBase = "6.62"
+DiffEqBase = "6.62, 7"
 ExplicitImports = "1.14.0"
 GeometricIntegrators = "0.15, 0.16"
 Reexport = "0.2, 1"
-SciMLBase = "2"
+SciMLBase = "2, 3"
 julia = "1.10"
 
 [extras]

src/solve.jl

@@ -15,14 +15,18 @@ function DiffEqBase.__solve(
         error("dt required for fixed timestep methods.")
     end
 
+    # DiffEqBase v7 passes `verbose` as a `DEVerbosity` object instead of a Bool.
+    # Treat anything that is not literally `false` as "show warnings".
+    verbose_bool = verbose !== false
+
     isstiff = !(
         alg isa Union{
             GIImplicitEuler, GIImplicitMidpoint,
             GISRK3, GIGLRK, GIRadauIA, GIRadauIIA,
         }
     )
 
-    if verbose
+    if verbose_bool
         warned = !isempty(kwargs) && check_keywords(alg, kwargs, warnlist)
         if !(prob.f isa DiffEqBase.AbstractParameterizedFunction) && isstiff
             if DiffEqBase.has_tgrad(prob.f)
@@ -68,15 +72,20 @@ function DiffEqBase.__solve(
         # Create function wrapper for GeometricIntegrators API
         # GeometricIntegrators expects: v(v, t, q, params)
         # DiffEqBase provides: f(du, u, p, t) for inplace or f(u, p, t) for out-of-place
+        # SciMLBase v3's default AutoSpecialize wraps f in a FunctionWrappersWrapper that
+        # only accepts the exact argument types captured at problem construction (e.g.
+        # `Matrix{Float64}`), so passing a `Base.ReshapedArray` to it errors. Unwrap the
+        # underlying user-defined function before invoking it.
+        raw_f = SciMLBase.unwrapped_f(prob.f.f)
         if !isinplace && u isa AbstractArray
-            v! = (v, t, q, params) -> (v .= vec(prob.f(reshape(q, sizeu), p, t)); nothing)
+            v! = (v, t, q, params) -> (v .= vec(raw_f(reshape(q, sizeu), p, t)); nothing)
         elseif !(u isa Vector{Float64})
             v! = (
                 v, t, q,
                 params,
-            ) -> (prob.f(reshape(v, sizeu), reshape(q, sizeu), p, t); nothing)
+            ) -> (raw_f(reshape(v, sizeu), reshape(q, sizeu), p, t); nothing)
         else
-            v! = (v, t, q, params) -> prob.f(v, q, p, t)
+            v! = (v, t, q, params) -> raw_f(v, q, p, t)
         end
 
         ode = GeometricIntegrators.ODEProblem(v!, prob.tspan, dt, vec(prob.u0))
@@ -111,17 +120,21 @@ function DiffEqBase.__solve(
 
         v! = (v, t, q, p_state, params) -> (v .= p_state)  # dq/dt = p
 
+        # Unwrap past SciMLBase v3 AutoSpecialize FunctionWrappers for the acceleration
+        # function so it accepts argument types other than those captured at construction.
+        raw_f1 = SciMLBase.unwrapped_f(prob.f.f1.f)
+
         # Handle both inplace and out-of-place problems
         if isinplace
             f! = (
                 f_out, t, q, p_state,
                 params,
-            ) -> (prob.f.f1.f(f_out, p_state, q, p, t); nothing)  # dp/dt = f1(p, q)
+            ) -> (raw_f1(f_out, p_state, q, p, t); nothing)  # dp/dt = f1(p, q)
         else
             f! = (
                 f_out, t, q, p_state,
                 params,
-            ) -> (f_out .= prob.f.f1.f(p_state, q, p, t); nothing)
+            ) -> (f_out .= raw_f1(p_state, q, p, t); nothing)
         end
 
         pode = GeometricIntegrators.PODEProblem(