Correct OrdinaryDiffEq DAE init integration

src/common_interface/integrator_utils.jl

@@ -99,8 +99,20 @@ function handle_callback_modifiers!(integrator::ARKODEIntegrator)
                   integrator.t, integrator.u)
 end
 
-function handle_callback_modifiers!(integrator::IDAIntegrator)
+"""
+    IDAReinit!(integrator)
+
+The `integrator` object keeps a shadow copy of `u`, `du` and `t`. If these are
+modified, this function needs to be called in order to update the solver's
+internal datastructures to re-gain consistency.
+"""
+function IDAReinit!(integrator::IDAIntegrator)
     IDAReInit(integrator.mem, integrator.t, integrator.u, integrator.du)
+    integrator.u_modified = false
+end
+
+function handle_callback_modifiers!(integrator::IDAIntegrator)
+    # Implicitly does IDAReinit!
     DiffEqBase.initialize_dae!(integrator, IDADefaultInit())
 end
 
@@ -175,6 +187,9 @@ end
 
 function DiffEqBase.initialize_dae!(integrator::IDAIntegrator,
                                     initializealg::IDADefaultInit)
+    if integrator.u_modified
+        IDAReinit!(integrator)
+    end
     integrator.f(integrator.tmp, integrator.du, integrator.u, integrator.p, integrator.t)
     tstart, tend = integrator.sol.prob.tspan
     if any(abs.(integrator.tmp) .>= integrator.opts.reltol)
@@ -190,6 +205,10 @@ function DiffEqBase.initialize_dae!(integrator::IDAIntegrator,
         end
         dt = integrator.dt == tstart ? tend : integrator.dt
         integrator.flag = IDACalcIC(integrator.mem, init_type, dt)
+
+        # Reflect consistent initial conditions back into the integrator's
+        # shadow copy. N.B.: ({du, u}_nvec are aliased to {du, u}).
+        IDAGetConsistentIC(integrator.mem, integrator.u_nvec, integrator.du_nvec)
     end
     if integrator.t == tstart && integrator.flag < 0
         integrator.sol = SciMLBase.solution_new_retcode(integrator.sol,

src/common_interface/solve.jl

@@ -1223,26 +1223,9 @@ function DiffEqBase.__init(prob::DiffEqBase.AbstractDAEProblem{uType, duType, tu
         jac = nothing
     end
 
-    tout = [tspan[1]]
-    if save_start
-        if save_idxs === nothing
-            ures = Vector{uType}()
-            dures = Vector{uType}()
-            save_value!(ures, u0, uType, save_idxs)
-            if dense
-                save_value!(dures, du0, uType, save_idxs)
-            end
-        else
-            ures = [u0[save_idxs]]
-            if dense
-                dures = [du0[save_idxs]]
-            end
-        end
-    else
-        ures = Vector{uType}()
-        dures = Vector{uType}()
-    end
-
+    tout = Float64[]
+    ures = Vector{uType}()
+    dures = Vector{uType}()
     tmp = isnothing(callbacks_internal) ? u0 : similar(u0)
     uprev = isnothing(callbacks_internal) ? u0 : similar(u0)
     retcode = flag >= 0 ? ReturnCode.Default : ReturnCode.InitialFailure
@@ -1311,6 +1294,14 @@ function DiffEqBase.__init(prob::DiffEqBase.AbstractDAEProblem{uType, duType, tu
                                initializealg)
 
     DiffEqBase.initialize_dae!(integrator, initializealg)
+    integrator.u_modified && IDAReinit!(integrator)
+
+    if save_start
+        push!(tout, integrator.t)
+        save_value!(ures, integrator.u, uType, save_idxs)
+        save_value!(dures, integrator.u, uType, save_idxs)
+    end
+
     initialize_callbacks!(integrator)
     integrator
 end # function solve

test/common_interface/ida.jl

@@ -93,3 +93,18 @@ tspan = (0.0, 10.0)
 dae_prob = DAEProblem(f!, du0, u0, tspan; differential_vars = [false])
 #Currently broken because initialize_dae! in OrdiniaryDiffEq specializes on ODEIntegrator
 init(dae_prob, IDA(), initializealg = NoInit()).u == [0.0]
+
+# test that initializers which modify states actually modify the states
+struct DumbInit <: DiffEqBase.DAEInitializationAlgorithm end
+function DiffEqBase.initialize_dae!(integrator::Sundials.IDAIntegrator, initializealg::DumbInit)
+    integrator.u .= 1
+    integrator.u_modified = true
+    DiffEqBase.initialize_dae!(integrator, Sundials.IDADefaultInit())
+end
+f(du, u, p, t) = du - u # u(t) = exp(t)
+prob = DAEProblem(f, zeros(1), zeros(1), (0, 1), differential_vars = trues(1))
+sol = solve(prob, IDA(), initializealg = DumbInit())
+# test that initializealg is refected in the sol
+isapprox(only(sol.u[begin]), 1, rtol = 1e-3)
+# test that solve produced the right answer.
+isapprox(only(sol.u[end]), exp(1), rtol = 1e-3)