Check the validity of DEs more carefully

Author: YingboMa

src/structural_transformation/codegen.jl

@@ -524,8 +524,11 @@ function ODAEProblem{iip}(
                           parammap = DiffEqBase.NullParameters();
                           callback = nothing,
                           use_union = false,
+                          check = true,
                           kwargs...
                          ) where {iip}
+    eqs = equations(sys)
+    check && ModelingToolkit.check_operator_variables(eqs, Differential)
     fun, dvs = build_torn_function(sys; kwargs...)
     ps = parameters(sys)
     defs = defaults(sys)
@@ -535,7 +538,7 @@ function ODAEProblem{iip}(
     u0 = ModelingToolkit.varmap_to_vars(u0map, dvs; defaults=defs, tofloat=true)
     p = ModelingToolkit.varmap_to_vars(parammap, ps; defaults=defs, tofloat=!use_union, use_union)
 
-    has_difference = any(isdifferenceeq, equations(sys))
+    has_difference = any(isdifferenceeq, eqs)
     if has_continuous_events(sys)
         event_cb = generate_rootfinding_callback(sys; kwargs...)
     else

src/systems/diffeqs/abstractodesystem.jl

@@ -95,8 +95,6 @@ function generate_dae_jacobian(sys::AbstractODESystem, dvs = states(sys), ps = p
 
 end
 
-check_derivative_variables(eq) = check_operator_variables(eq, Differential)
-
 function generate_function(
         sys::AbstractODESystem, dvs = states(sys), ps = parameters(sys);
         implicit_dae=false,
@@ -107,7 +105,7 @@ function generate_function(
 
     eqs = [eq for eq in equations(sys) if !isdifferenceeq(eq)]
     if !implicit_dae
-        foreach(check_derivative_variables, eqs)
+        check_operator_variables(eqs, Differential)
     end
     check_lhs(eqs, Differential, Set(dvs))
     # substitute x(t) by just x
@@ -130,7 +128,7 @@ end
 
 function generate_difference_cb(sys::ODESystem, dvs = states(sys), ps = parameters(sys); kwargs...)
     eqs = equations(sys)
-    foreach(check_difference_variables, eqs)
+    check_operator_variables(eqs, Difference)
 
     var2eq = Dict(arguments(eq.lhs)[1] => eq for eq in eqs if isdifference(eq.lhs))
 

src/systems/discrete_system/discrete_system.jl

@@ -243,14 +243,12 @@ function get_delay_val(iv, x)
     return -delay
 end
 
-check_difference_variables(eq) = check_operator_variables(eq, Difference)
-
 function generate_function(
         sys::DiscreteSystem, dvs = states(sys), ps = parameters(sys);
         kwargs...
     )
     eqs = equations(sys)
-    foreach(check_difference_variables, eqs)
+    check_operator_variables(eqs, Difference)
     rhss = [eq.rhs for eq in eqs]
 
     u = map(x->time_varying_as_func(value(x), sys), dvs)

src/utils.jl

@@ -233,12 +233,39 @@ end
 end
 
 "Check if difference/derivative operation occurs in the R.H.S. of an equation"
-function check_operator_variables(eq, op::Type, expr=eq.rhs)
+function _check_operator_variables(eq, op::T, expr=eq.rhs) where T
     istree(expr) || return nothing
     if operation(expr) isa op
         throw_invalid_operator(expr, eq, op)
     end
-    foreach(expr -> check_operator_variables(eq, op, expr), SymbolicUtils.unsorted_arguments(expr))
+    foreach(expr -> _check_operator_variables(eq, op, expr), SymbolicUtils.unsorted_arguments(expr))
+end
+"Check if all the LHS are unique"
+function check_operator_variables(eqs, op::T) where T
+    ops = Set()
+    tmp = Set()
+    for eq in eqs
+        _check_operator_variables(eq, op)
+        vars!(tmp, eq.lhs)
+        if length(tmp) == 1
+            x = only(tmp)
+            if op === Differential
+                # Having a differece is fine for ODEs
+                is_tmp_fine = isdifferential(x) || isdifference(x)
+            else
+                is_tmp_fine = istree(x) && !(operation(x) isa op)
+            end
+        else
+            nd = count(x->istree(x) && !(operation(x) isa op), tmp)
+            is_tmp_fine = iszero(nd)
+        end
+        empty!(tmp)
+        is_tmp_fine || error("The LHS cannot contain nondifferentiated variables. Please run `structural_simplify` or use the DAE form.\nGot $eq")
+        for v in tmp
+            v in ops && error("The LHS operator must be unique. Please run `structural_simplify` or use the DAE form. $v appears in LHS more than once.")
+            push!(ops, v)
+        end
+    end
 end
 
 isoperator(expr, op) = istree(expr) && operation(expr) isa op