WIP: Initialization on non-DAE models

#2508 pointed out that the heuristic of `implicit_dae || calculate_massmatrix(sys) !== I`, i.e. "only do initialization on DAEs", while sensible for the DAE model, doesn't quite fit in the MTK context. Instead, what we need to do is check whether the varmap is correct and consistent for the chosen ODE. If it's not consistent for the selected ODE, then we still need to run an initialization step for the resulting ODE solve.

This needs a few changes in the ODE solver as though, since currently only DAE solvers will run the initialization step. Thus the test case:

```julia
using ModelingToolkit, OrdinaryDiffEq, Test
using ModelingToolkit: t_nounits as t, D_nounits as D

function System(;name)
    vars = @variables begin
        dx(t), [guess=0]
        ddx(t), [guess=0]
    end
    eqs = [
        D(dx) ~ ddx
        0 ~ ddx + dx + 1
    ]
    return ODESystem(eqs, t, vars, []; name)
end

@mtkbuild sys = System()
prob = ODEProblem(sys, [sys.dx => 1], (0,1)) # OK
prob = ODEProblem(sys, [sys.ddx => -2], (0,1), guesses = [sys.dx => 1])
sol = solve(prob, Rodas5P())
sol = solve(prob, Tsit5())
```

gives an erroneous success as it skips initialization, using `prob.u0 == [0.0]`, the sentinel value since initialization is only run on DAEs. We will need to override that so that initialization is always run on any system that provides an initializeprob.

Author: ChrisRackauckas

src/systems/diffeqs/abstractodesystem.jl

@@ -862,14 +862,18 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap;
     ps = full_parameters(sys)
     iv = get_iv(sys)
 
+    varmap = merge(defaults, todict(u0map))
+    varlist = collect(map(unwrap, dvs))
+    missingvars = setdiff(varlist, collect(keys(varmap)))
+
     # Append zeros to the variables which are determined by the initialization system
     # This essentially bypasses the check for if initial conditions are defined for DAEs
     # since they will be checked in the initialization problem's construction
     # TODO: make check for if a DAE cheaper than calculating the mass matrix a second time!
     ci = infer_clocks!(ClockInference(TearingState(sys)))
     # TODO: make it work with clocks
     # ModelingToolkit.get_tearing_state(sys) !== nothing => Requires structural_simplify first
-    if (implicit_dae || calculate_massmatrix(sys) !== I) &&
+    if (implicit_dae || !isempty(missingvars)) &&
        all(isequal(Continuous()), ci.var_domain) &&
        ModelingToolkit.get_tearing_state(sys) !== nothing
         if eltype(u0map) <: Number