Handle dummy derivatives and more in initialize_equations

Also throws a better error for the case where a differential is not caught in the initialization system handling. Fixes #3029

Author: ChrisRackauckas

src/structural_transformation/symbolics_tearing.jl

@@ -453,6 +453,7 @@ function tearing_reassemble(state::TearingState, var_eq_matching,
             # convert it into the mass matrix form.
             # We cannot solve the differential variable like D(x)
             if isdervar(iv)
+                isnothing(D) && error("Differential found in a non-differential system. Likely this is a bug in the construction of an initialization system. Please report this issue with a reproducible example. Offending equation: $(equations(sys)[ieq])")
                 order, lv = var_order(iv)
                 dx = D(simplify_shifts(lower_varname_withshift(
                     fullvars[lv], idep, order - 1)))

src/systems/nonlinear/initializesystem.jl

@@ -91,13 +91,16 @@ function generate_initializesystem(sys::ODESystem;
         end
     end
 
-    pars = [parameters(sys); get_iv(sys)]
-    nleqs = if algebraic_only
-        [eqs_ics; observed(sys)]
-    else
-        [eqs_ics; get_initialization_eqs(sys); initialization_eqs; observed(sys)]
+    if !algebraic_only
+        for eq in [get_initialization_eqs(sys); initialization_eqs]
+            _eq = ModelingToolkit.fixpoint_sub(eq, full_diffmap)
+            push!(eqs_ics,_eq)
+        end
     end
 
+    pars = [parameters(sys); get_iv(sys)]
+    nleqs = [eqs_ics; observed(sys)]
+
     sys_nl = NonlinearSystem(nleqs,
         full_states,
         pars;

test/initializationsystem.jl

@@ -508,3 +508,13 @@ end
     sol2 = solve(prob2, Tsit5())
     @test all(sol2[x] .== 2) && all(sol2[y] .== 2)
 end
+
+# https://github.com/SciML/ModelingToolkit.jl/issues/3029
+@testset "Derivatives in Initialization Equations" begin
+    @variables x(t)
+    sys = ODESystem([D(D(x)) ~ 0], t; initialization_eqs = [x ~ 0, D(x) ~ 1], name = :sys) |> structural_simplify
+    @test_nowarn ODEProblem(sys, [], (0.0, 1.0), []);
+
+    sys = ODESystem([D(D(x)) ~ 0], t; initialization_eqs = [x ~ 0, D(D(x)) ~ 0], name = :sys) |> structural_simplify
+    @test_nowarn ODEProblem(sys, [D(x) => 1.0], (0.0, 1.0), []);
+end