Refactor & Generalize: System for opting alg into linear interpolation

We haven't really needed this much before since most methods had a special interpolation or default Hermite, so we had a bespoke way of opting just a few algorithms out (FunctionMap and the DAE solvers for now). But there are cases which need to be able to opt into saying they only want a linear interpolation as the Hermite would take extra function evaluations or would not satisfy the right properties, and thus for these equations this gives them an overload they can use to opt out on a per-algorithm basis via a trait function.

lib/OrdinaryDiffEqCore/src/alg_utils.jl

@@ -444,3 +444,6 @@ function Base.show(io::IO, ::MIME"text/plain", alg::OrdinaryDiffEqAlgorithm)
     end
     print(io, ")")
 end
+
+# Defaults in the current system: currently opt out DAEAlgorithms until complete
+default_linear_interpolation(alg, prob) = alg isa DAEAlgorithm || prob isa DiscreteProblem

lib/OrdinaryDiffEqCore/src/solve.jl

@@ -26,9 +26,7 @@ function DiffEqBase.__init(
                          saveat isa Number || prob.tspan[1] in saveat,
         save_end = nothing,
         callback = nothing,
-        dense = save_everystep &&
-                    !(alg isa DAEAlgorithm) && !(prob isa DiscreteProblem) &&
-                    isempty(saveat),
+        dense = save_everystep && isempty(saveat) && !default_linear_interpolation(prob, alg),
         calck = (callback !== nothing && callback !== CallbackSet()) ||
                     (dense) || !isempty(saveat), # and no dense output
         dt = isdiscretealg(alg) && isempty(tstops) ?