Merge branch 'master' into uprev

Author: ChrisRackauckas

src/alg_utils.jl

@@ -285,8 +285,11 @@ function DiffEqBase.prepare_alg(alg::CompositeAlgorithm, u0, p, prob)
 end
 
 has_autodiff(alg::OrdinaryDiffEqAlgorithm) = false
-has_autodiff(alg::Union{OrdinaryDiffEqAdaptiveImplicitAlgorithm, OrdinaryDiffEqImplicitAlgorithm, 
-                        CompositeAlgorithm, OrdinaryDiffEqExponentialAlgorithm, DAEAlgorithm}) = true
+function has_autodiff(alg::Union{
+        OrdinaryDiffEqAdaptiveImplicitAlgorithm, OrdinaryDiffEqImplicitAlgorithm,
+        CompositeAlgorithm, OrdinaryDiffEqExponentialAlgorithm, DAEAlgorithm})
+    true
+end
 
 # Extract AD type parameter from algorithm, returning as Val to ensure type stability for boolean options.
 function _alg_autodiff(alg::OrdinaryDiffEqAlgorithm)

src/caches/bdf_caches.jl

@@ -47,8 +47,11 @@ function alg_cache(alg::ABDF2, u, rate_prototype, ::Type{uEltypeNoUnits},
     fsalfirstprev = zero(rate_prototype)
     atmp = similar(u, uEltypeNoUnits)
     recursivefill!(atmp, false)
+    algebraic_vars = f.mass_matrix === I ? nothing :
+                     [all(iszero, x) for x in eachcol(f.mass_matrix)]
 
-    eulercache = ImplicitEulerCache(u, uprev, uprev2, fsalfirst, atmp, nlsolver)
+    eulercache = ImplicitEulerCache(
+        u, uprev, uprev2, fsalfirst, atmp, nlsolver, algebraic_vars)
 
     dtₙ₋₁ = one(dt)
     zₙ₋₁ = zero(u)

src/caches/sdirk_caches.jl

@@ -1,13 +1,14 @@
 abstract type SDIRKMutableCache <: OrdinaryDiffEqMutableCache end
 
-@cache mutable struct ImplicitEulerCache{uType, rateType, uNoUnitsType, N} <:
+@cache mutable struct ImplicitEulerCache{uType, rateType, uNoUnitsType, N, AV} <:
                       SDIRKMutableCache
     u::uType
     uprev::uType
     uprev2::uType
     fsalfirst::rateType
     atmp::uNoUnitsType
     nlsolver::N
+    algebraic_vars::AV
 end
 
 function alg_cache(alg::ImplicitEuler, u, rate_prototype, ::Type{uEltypeNoUnits},
@@ -22,7 +23,10 @@ function alg_cache(alg::ImplicitEuler, u, rate_prototype, ::Type{uEltypeNoUnits}
     atmp = similar(u, uEltypeNoUnits)
     recursivefill!(atmp, false)
 
-    ImplicitEulerCache(u, uprev, uprev2, fsalfirst, atmp, nlsolver)
+    algebraic_vars = f.mass_matrix === I ? nothing :
+                     [all(iszero, x) for x in eachcol(f.mass_matrix)]
+
+    ImplicitEulerCache(u, uprev, uprev2, fsalfirst, atmp, nlsolver, algebraic_vars)
 end
 
 mutable struct ImplicitEulerConstantCache{N} <: OrdinaryDiffEqConstantCache

src/initialize_dae.jl

@@ -134,11 +134,12 @@ end
 function _initialize_dae!(integrator, prob::Union{ODEProblem, DAEProblem},
         alg::OverrideInit, isinplace::Union{Val{true}, Val{false}})
     initializeprob = prob.f.initializeprob
-    
+
     # If it doesn't have autodiff, assume it comes from symbolic system like ModelingToolkit
     # Since then it's the case of not a DAE but has initializeprob
     # In which case, it should be differentiable
-    isAD = has_autodiff(integrator.alg) ? alg_autodiff(integrator.alg) isa AutoForwardDiff : true
+    isAD = has_autodiff(integrator.alg) ? alg_autodiff(integrator.alg) isa AutoForwardDiff :
+           true
 
     alg = default_nlsolve(alg.nlsolve, isinplace, initializeprob.u0, initializeprob, isAD)
     nlsol = solve(initializeprob, alg)

src/perform_step/sdirk_perform_step.jl

@@ -102,6 +102,13 @@ end
     end
 
     integrator.fsallast = f(u, p, t + dt)
+
+    if integrator.opts.adaptive && integrator.f.mass_matrix !== I
+        atmp = @. ifelse(!integrator.differential_vars, integrator.fsallast, false) ./
+                  integrator.opts.abstol
+        integrator.EEst += integrator.opts.internalnorm(atmp, t)
+    end
+
     integrator.stats.nf += 1
     integrator.k[1] = integrator.fsalfirst
     integrator.k[2] = integrator.fsallast
@@ -152,6 +159,12 @@ end
     end
     integrator.stats.nf += 1
     f(integrator.fsallast, u, p, t + dt)
+
+    if integrator.opts.adaptive && integrator.f.mass_matrix !== I
+        @.. broadcast=false atmp=ifelse(cache.algebraic_vars, integrator.fsallast, false) /
+                                 integrator.opts.abstol
+        integrator.EEst += integrator.opts.internalnorm(atmp, t)
+    end
 end
 
 @muladd function perform_step!(integrator, cache::ImplicitMidpointConstantCache,

test/interface/dae_initialize_integration.jl

@@ -44,17 +44,22 @@ sol = solve(prob, Rodas5(), initializealg = ShampineCollocationInit())
 # Initialize on ODEs
 # https://github.com/SciML/ModelingToolkit.jl/issues/2508
 
-function testsys(du,u,p,t)
+function testsys(du, u, p, t)
     du[1] = -2
 end
-function initsys(du,u,p)
+function initsys(du, u, p)
     du[1] = -1 + u[1]
 end
 nlprob = NonlinearProblem(initsys, [0.0])
 initprobmap(nlprob) = nlprob.u
 sol = solve(nlprob)
 
 _f = ODEFunction(testsys; initializeprob = nlprob, initializeprobmap = initprobmap)
+prob = ODEProblem(_f, [0.0], (0.0, 1.0))
+sol = solve(prob, Tsit5())
+@test SciMLBase.successful_retcode(sol)
+@test sol[1] == [1.0]
+
 prob = ODEProblem(_f, [0.0], (0.0,1.0))
 sol = solve(prob, Tsit5(), dt = 1e-10)
 @test SciMLBase.successful_retcode(sol)
@@ -66,4 +71,4 @@ sol = solve(prob, Rodas5P(), dt = 1e-10)
 @test SciMLBase.successful_retcode(sol)
 @test sol[1] == [1.0]
 @test sol[2] ≈ [0.9999999998]
-@test sol[end] ≈ [-1.0]
+@test sol[end] ≈ [-1.0]