Add DiffEqBase.stripunits helper and use it for tTypeNoUnits

Adds a narrow "strip units only" primitive that is orthogonal to the existing
`value` (strips everything, including AD) and `unitfulvalue` (strips AD, keeps
units) helpers. Defaults to the identity; Unitful, FlexUnits, and
DynamicQuantities extensions override it to return the underlying numeric value.

Switches the `tTypeNoUnits` derivation in `OrdinaryDiffEqCore/src/solve.jl`
from `DiffEqBase.value(oneunit(first(tspan)))` to
`DiffEqBase.stripunits(oneunit(first(tspan)))` so that ForwardDiff Duals in t0
are preserved through `tTypeNoUnits`. The previous `value`-based form collapsed
Duals to Float64, which then made `NLSolver{true, tTypeNoUnits}` use Float64
fields for ηold etc., and SDIRK solvers (KenCarp4/47/5/58, TRBDF2) that run
with a Dual t0 errored with `MethodError: Float64(::Dual)`.

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

Author: ChrisRackauckas

lib/DiffEqBase/ext/DiffEqBaseDynamicQuantitiesExt.jl

@@ -10,6 +10,8 @@ import DiffEqBase: default_factorize
     return real(abs2(ustrip(x)))
 end
 
+DiffEqBase.stripunits(x::UnionAbstractQuantity) = ustrip(x)
+
 DiffEqBase._rate_prototype(u, t::UnionAbstractQuantity, onet) = u / oneunit(t)
 DiffEqBase.timedepentdtmin(t::UnionAbstractQuantity, dtmin) =
     abs(ustrip(dtmin / oneunit(t)) * oneunit(t))

lib/DiffEqBase/ext/DiffEqBaseFlexUnitsExt.jl

@@ -11,6 +11,8 @@ value(x::Quantity{T, U}) where {T, U} = dstrip(x)
 unitfulvalue(::Type{T}) where {T <: Quantity} = T
 unitfulvalue(x::Quantity) = x
 
+DiffEqBase.stripunits(x::Quantity) = dstrip(x)
+
 @inline function DiffEqBase.ODE_DEFAULT_NORM(
         u::AbstractArray{
             <:Quantity,

lib/DiffEqBase/ext/DiffEqBaseUnitfulExt.jl

@@ -11,6 +11,8 @@ value(x::Unitful.AbstractQuantity) = x.val
 unitfulvalue(x::Type{T}) where {T <: Unitful.AbstractQuantity} = T
 unitfulvalue(x::Unitful.AbstractQuantity) = x
 
+DiffEqBase.stripunits(x::Unitful.AbstractQuantity) = Unitful.ustrip(x)
+
 @inline function DiffEqBase.ODE_DEFAULT_NORM(
         u::AbstractArray{
             <:Unitful.AbstractQuantity,

lib/DiffEqBase/src/utils.jl

@@ -52,6 +52,12 @@ end
 # for the non-unitful case the correct type is just u
 _rate_prototype(u, t::T, onet::T) where {T} = u
 
+# Strip only the unit wrapper, leaving AD/uncertainty wrappers (Dual, Measurement,
+# Tracker, etc.) intact. Extensions for Unitful, DynamicQuantities, and FlexUnits
+# override this to return the underlying numeric value.
+# Complementary to `value` (strips everything) and `unitfulvalue` (strips AD, keeps units).
+stripunits(x) = x
+
 # Nonlinear Solve functionality
 @inline __fast_scalar_indexing(args...) = all(ArrayInterface.fast_scalar_indexing, args)
 

lib/OrdinaryDiffEqCore/src/solve.jl

@@ -317,7 +317,7 @@ Base.@constprop :aggressive function _ode_init(
     uBottomEltypeNoUnits = recursive_unitless_bottom_eltype(u)
 
     uEltypeNoUnits = recursive_unitless_eltype(u)
-    tTypeNoUnits = typeof(DiffEqBase.value(oneunit(first(tspan))))
+    tTypeNoUnits = typeof(DiffEqBase.stripunits(oneunit(first(tspan))))
 
     scalar_type_tol =
         uBottomEltypeNoUnits == uBottomEltype &&