Perform bi-directional sync of integrator u/du values in DAE init

Sundials.jl got an implementation of the OrdinaryDiffEq DAE initialization
interface in SciML#401. Unfortunately, it turns out that this doesn't quite work,
because the u/du in Sundials' integrator struct are not the actual memory
that IDA uses (it has an internal copy). In [1], I propose re-using `u_modified`
to inform us that the shadow copy was modified during initialization. This
is the Sundials side of that interface.

[1] SciML/OrdinaryDiffEq.jl#1969

src/common_interface/integrator_utils.jl

@@ -122,8 +122,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
 
@@ -198,6 +210,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)
@@ -213,6 +228,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_nvce)
     end
     if integrator.t == tstart && integrator.flag < 0
         integrator.sol = SciMLBase.solution_new_retcode(integrator.sol,

src/common_interface/solve.jl

@@ -1311,6 +1311,7 @@ function DiffEqBase.__init(prob::DiffEqBase.AbstractDAEProblem{uType, duType, tu
                                initializealg)
 
     DiffEqBase.initialize_dae!(integrator, initializealg)
+    integrator.u_modified && IDAReinit!(integrator)
     initialize_callbacks!(integrator)
     integrator
 end # function solve