Merge pull request #1855 from SciML/myb/lowering_fix

Handle differentiation chains that are only partly dummy derivatives

Author: YingboMa

src/structural_transformation/pantelides.jl

@@ -75,7 +75,7 @@ end
 Perform Pantelides algorithm.
 """
 function pantelides!(state::TransformationState, ag::Union{AliasGraph, Nothing} = nothing;
-                     maxiters = 8000)
+                     finalize = true, maxiters = 8000)
     @unpack graph, solvable_graph, var_to_diff, eq_to_diff = state.structure
     neqs = nsrcs(graph)
     nvars = nv(var_to_diff)
@@ -163,6 +163,11 @@ function pantelides!(state::TransformationState, ag::Union{AliasGraph, Nothing}
         pathfound ||
             error("maxiters=$maxiters reached! File a bug report if your system has a reasonable index (<100), and you are using the default `maxiters`. Try to increase the maxiters by `pantelides(sys::ODESystem; maxiters=1_000_000)` if your system has an incredibly high index and it is truly extremely large.")
     end # for k in 1:neqs′
+
+    finalize && for var in 1:ndsts(graph)
+        varwhitelist[var] && continue
+        var_eq_matching[var] = unassigned
+    end
     return var_eq_matching
 end
 
@@ -175,6 +180,6 @@ instead, which calls this function internally.
 """
 function dae_index_lowering(sys::ODESystem; kwargs...)
     state = TearingState(sys)
-    var_eq_matching = pantelides!(state; kwargs...)
+    var_eq_matching = pantelides!(state; finalize = false, kwargs...)
     return invalidate_cache!(pantelides_reassemble(state, var_eq_matching))
 end

src/structural_transformation/partial_state_selection.jl

@@ -143,13 +143,6 @@ function partial_state_selection_graph!(structure::SystemStructure, var_eq_match
         level
     end
 
-    # TODO: Should pantelides just return this?
-    for var in 1:ndsts(graph)
-        if varlevel[var] !== 0
-            var_eq_matching[var] = unassigned
-        end
-    end
-
     var_eq_matching = pss_graph_modia!(structure,
                                        complete(var_eq_matching), varlevel, inv_varlevel,
                                        inv_eqlevel)
@@ -267,6 +260,10 @@ function dummy_derivative_graph!(structure::SystemStructure, var_eq_matching, ja
         end
     end
 
+    if (n_diff_eqs = count(!isnothing, diff_to_eq)) !=
+       (n_dummys = length(dummy_derivatives))
+        @warn "The number of dummy derivatives ($n_dummys) does not match the number of differentiated equations ($n_diff_eqs)."
+    end
     dummy_derivatives_set = BitSet(dummy_derivatives)
     # We can eliminate variables that are not a selected state (differential
     # variables). Selected states are differentiated variables that are not

src/structural_transformation/symbolics_tearing.jl

@@ -195,6 +195,23 @@ function to_mass_matrix_form(neweqs, ieq, graph, fullvars, isdervar::F,
     end
 end
 
+#=
+function check_diff_graph(var_to_diff, fullvars)
+    diff_to_var = invview(var_to_diff)
+    for (iv, v) in enumerate(fullvars)
+        ov, order = var_from_nested_derivative(v)
+        graph_order = 0
+        vv = iv
+        while true
+            vv = diff_to_var[vv]
+            vv === nothing && break
+            graph_order += 1
+        end
+        @assert graph_order==order "graph_order: $graph_order, order: $order for variable $v"
+    end
+end
+=#
+
 function tearing_reassemble(state::TearingState, var_eq_matching; simplify = false)
     @unpack fullvars, sys = state
     @unpack solvable_graph, var_to_diff, eq_to_diff, graph = state.structure
@@ -237,8 +254,12 @@ function tearing_reassemble(state::TearingState, var_eq_matching; simplify = fal
         possible_x_t[rhs] = i, lhs
     end
 
-    #removed_eqs = Int[]
-    #removed_vars = Int[]
+    if ModelingToolkit.has_iv(state.sys)
+        iv = get_iv(state.sys)
+        D = Differential(iv)
+    else
+        iv = D = nothing
+    end
     removed_obs = Int[]
     diff_to_var = invview(var_to_diff)
     dummy_sub = Dict()
@@ -258,7 +279,22 @@ function tearing_reassemble(state::TearingState, var_eq_matching; simplify = fal
                 neweqs[eq] = substitute(neweqs[eq], dd => v_t)
             end
             fullvars[dv] = v_t
+            # If we have:
+            # x -> D(x) -> D(D(x))
+            # We need to to transform it to:
+            # x   x_t -> D(x_t)
             # update the structural information
+            dx = dv
+            x_t = v_t
+            while (ddx = var_to_diff[dx]) !== nothing
+                dx_t = D(x_t)
+                for eq in 𝑑neighbors(graph, ddx)
+                    neweqs[eq] = substitute(neweqs[eq], fullvars[ddx] => dx_t)
+                end
+                fullvars[ddx] = dx_t
+                dx = ddx
+                x_t = dx_t
+            end
             diff_to_var[dv] = nothing
         end
     end
@@ -323,12 +359,6 @@ function tearing_reassemble(state::TearingState, var_eq_matching; simplify = fal
     # As a final note, in all the above cases where we need to introduce new
     # variables and equations, don't add them when they already exist.
 
-    if ModelingToolkit.has_iv(state.sys)
-        iv = get_iv(state.sys)
-        D = Differential(iv)
-    else
-        iv = D = nothing
-    end
     nvars = ndsts(graph)
     processed = falses(nvars)
     subinfo = NTuple{3, Int}[]
@@ -488,7 +518,7 @@ function tearing_reassemble(state::TearingState, var_eq_matching; simplify = fal
             # 0 ~ a * var + b
             # var ~ -b/a
             if ModelingToolkit._iszero(a)
-                @warn "Tearing: $eq is a singular equation!"
+                @warn "Tearing: solving $eq for $var is singular!"
                 #push!(removed_eqs, ieq)
                 #push!(removed_vars, iv)
             else

src/systems/diffeqs/abstractodesystem.jl

@@ -189,7 +189,7 @@ function generate_difference_cb(sys::ODESystem, dvs = states(sys), ps = paramete
 end
 
 function calculate_massmatrix(sys::AbstractODESystem; simplify = false)
-    eqs = [eq for eq in full_equations(sys) if !isdifferenceeq(eq)]
+    eqs = [eq for eq in equations(sys) if !isdifferenceeq(eq)]
     dvs = states(sys)
     M = zeros(length(eqs), length(eqs))
     state2idx = Dict(s => i for (i, s) in enumerate(dvs))

src/systems/systemstructure.jl

@@ -315,17 +315,16 @@ function TearingState(sys; quick_cancel = false, check = true)
 
     nvars = length(fullvars)
     diffvars = []
-    vartype = fill(DIFFERENTIAL_VARIABLE, nvars)
     var_to_diff = DiffGraph(nvars, true)
     for dervaridx in dervaridxs
-        vartype[dervaridx] = DERIVATIVE_VARIABLE
         dervar = fullvars[dervaridx]
         diffvar = arguments(dervar)[1]
         diffvaridx = var2idx[diffvar]
         push!(diffvars, diffvar)
         var_to_diff[diffvaridx] = dervaridx
     end
 
+    #=
     algvars = setdiff(states(sys), diffvars)
     for algvar in algvars
         # it could be that a variable appeared in the states, but never appeared
@@ -334,10 +333,8 @@ function TearingState(sys; quick_cancel = false, check = true)
         #if algvaridx == 0
         #    check ? throw(InvalidSystemException("The system is missing an equation for $algvar.")) : return nothing
         #end
-        if algvaridx != 0
-            vartype[algvaridx] = ALGEBRAIC_VARIABLE
-        end
     end
+    =#
 
     graph = BipartiteGraph(neqs, nvars, Val(false))
     for (ie, vars) in enumerate(symbolic_incidence), v in vars