save du and k for DAEs when dense==true

src/integrators/integrator_utils.jl

@@ -108,6 +108,14 @@ function _savevalues!(integrator, force_save, reduce_size)::Tuple{Bool, Bool}
                 copyat_or_push!(integrator.sol.alg_choice, integrator.saveiter,
                     integrator.cache.current)
             end
+            if integrator.opts.save_du
+                if integrator.opts.save_idxs === nothing
+                    copyat_or_push!(integrator.sol.du, integrator.saveiter, integrator(integrator.t, Val{1}))
+                else
+                    copyat_or_push!(integrator.sol.du, integrator.saveiter,
+                        integrator(integrator.t, Val{1}, idxs=integrator.opts.save_idxs), false)
+                end
+            end
         end
     end
     if force_save || (integrator.opts.save_everystep &&
@@ -140,6 +148,14 @@ function _savevalues!(integrator, force_save, reduce_size)::Tuple{Bool, Bool}
             copyat_or_push!(integrator.sol.alg_choice, integrator.saveiter,
                 integrator.cache.current)
         end
+        if integrator.opts.save_du
+            if integrator.opts.save_idxs === nothing
+                copyat_or_push!(integrator.sol.du, integrator.saveiter, integrator(integrator.t, Val{1}))
+            else
+                copyat_or_push!(integrator.sol.du, integrator.saveiter,
+                    integrator(integrator.t, Val{1}, idxs=integrator.opts.save_idxs), false)
+            end
+        end
     end
     reduce_size && resize!(integrator.k, integrator.kshortsize)
     return saved, savedexactly
@@ -493,7 +509,7 @@ function reset_fsal!(integrator)
             integrator.fsalfirst = integrator.f(integrator.u, integrator.p, integrator.t)
         end
     end
-    
+
     # Do not set false here so it can be checked in the algorithm
     # integrator.reeval_fsal = false
 end

src/integrators/type.jl

@@ -46,6 +46,7 @@ mutable struct DEOptions{absType, relType, QT, tType, Controller, F1, F2, F3, F4
     force_dtmin::Bool
     advance_to_tstop::Bool
     stop_at_next_tstop::Bool
+    save_du::Bool
 end
 
 TruncatedStacktraces.@truncate_stacktrace DEOptions

src/solve.jl

@@ -25,7 +25,7 @@ function DiffEqBase.__init(prob::Union{DiffEqBase.AbstractODEProblem,
     save_end = nothing,
     callback = nothing,
     dense = save_everystep &&
-                !(alg isa Union{DAEAlgorithm, FunctionMap}) &&
+                !(alg isa FunctionMap) &&
                 isempty(saveat),
     calck = (callback !== nothing && callback !== CallbackSet()) ||
                 (dense) || !isempty(saveat), # and no dense output
@@ -70,6 +70,7 @@ function DiffEqBase.__init(prob::Union{DiffEqBase.AbstractODEProblem,
     alias_u0 = false,
     alias_du0 = false,
     initializealg = DefaultInit(),
+    save_du = false,
     kwargs...) where {recompile_flag}
     if prob isa DiffEqBase.AbstractDAEProblem && alg isa OrdinaryDiffEqAlgorithm
         error("You cannot use an ODE Algorithm with a DAEProblem")
@@ -407,16 +408,24 @@ function DiffEqBase.__init(prob::Union{DiffEqBase.AbstractODEProblem,
         unstable_check,
         verbose, calck, force_dtmin,
         advance_to_tstop,
-        stop_at_next_tstop)
+        stop_at_next_tstop,
+        save_du)
 
     stats = SciMLBase.DEStats(0)
     differential_vars = prob isa DAEProblem ? prob.differential_vars : get_differential_vars(f, u)
 
     id = InterpolationData(f, timeseries, ts, ks, alg_choice, dense, cache, differential_vars, false)
-    sol = DiffEqBase.build_solution(prob, _alg, ts, timeseries,
-        dense = dense, k = ks, interp = id,
-        alg_choice = alg_choice,
-        calculate_error = false, stats = stats)
+    if _alg isa DAEAlgorithm
+        sol = DiffEqBase.build_solution(prob, _alg, ts, timeseries, Vector{typeof(du)}(undef,0),
+            dense = dense, k = ks, interp = id,
+            alg_choice = alg_choice,
+            calculate_error = false, stats = stats)
+    else
+        sol = DiffEqBase.build_solution(prob, _alg, ts, timeseries,
+            dense = dense, k = ks, interp = id,
+            alg_choice = alg_choice,
+            calculate_error = false, stats = stats)
+    end
 
     if recompile_flag == true
         FType = typeof(f)