From a193152ee75957b365d3fb0eabde22615709c1b6 Mon Sep 17 00:00:00 2001
From: Keno Fischer <keno@juliacomputing.com>
Date: Wed, 25 Oct 2023 06:33:03 +0000
Subject: [PATCH] Give a warning for underconstrained variables forced to 0
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

@oxinabox noticed that the following MSL test
https://github.com/SciML/ModelingToolkitStandardLibrary.jl/blob/0f0fae138fe202140a4862eccaaac734051dc44f/test/Mechanical/translational.jl

has a test that depends on an underconstrained variable (either
`free.f` or `acc.f` may be set to an arbitrary value without
affecting the dynamics of the system). Our structural singularity
removal logic detects these situations and chooses one variable
arbitrarily to set to 0. However, it does so silently, so users
are unaware that they likely have a modeling bug. This PR adds
a warning when this happens. For example, the above mentioned
test now gives:

```
┌ Warning: The model is under-constrained. Variable acc₊flange₊f(t) was arbitrarily chosen to be set to 0. This may indicate a model bug!
└ @ ModelingToolkit ~/.julia/dev/ModelingToolkit/src/systems/alias_elimination.jl:54
```

This is styled as an overridable callback, so higher level modeling
frameworks that use MTK as a library can hook into this to give
more domain-specific errors if desired.
---
 src/systems/alias_elimination.jl | 43 +++++++++++++++++++++-----------
 1 file changed, 29 insertions(+), 14 deletions(-)

diff --git a/src/systems/alias_elimination.jl b/src/systems/alias_elimination.jl
index 0437cb72e6..067f0937ad 100644
--- a/src/systems/alias_elimination.jl
+++ b/src/systems/alias_elimination.jl
@@ -1,14 +1,15 @@
 using SymbolicUtils: Rewriters
 using Graphs.Experimental.Traversals
 
-function alias_eliminate_graph!(state::TransformationState; kwargs...)
+function alias_eliminate_graph!(state::TransformationState;
+    variable_underconstrained! = force_var_to_zero!, kwargs...)
     mm = linear_subsys_adjmat!(state; kwargs...)
     if size(mm, 1) == 0
         return mm # No linear subsystems
     end
 
     @unpack graph, var_to_diff, solvable_graph = state.structure
-    mm = alias_eliminate_graph!(state, mm)
+    mm = alias_eliminate_graph!(state, mm; variable_underconstrained!)
     s = state.structure
     for g in (s.graph, s.solvable_graph)
         g === nothing && continue
@@ -47,9 +48,17 @@ function alias_elimination!(state::TearingState; kwargs...)
     sys = state.sys
     complete!(state.structure)
     graph_orig = copy(state.structure.graph)
-    mm = alias_eliminate_graph!(state; kwargs...)
-
     fullvars = state.fullvars
+
+    function variable_underconstrained_mtk!(structure::SystemStructure,
+        ils::SparseMatrixCLIL,
+        v::Int)
+        @warn "The model is under-constrained. Variable $(fullvars[v]) was arbitrarily chosen to be set to 0. This may indicate a model bug!"
+        return force_var_to_zero!(structure, ils, v)
+    end
+    mm = alias_eliminate_graph!(state;
+        variable_underconstrained! = variable_underconstrained_mtk!, kwargs...)
+
     @unpack var_to_diff, graph, solvable_graph = state.structure
 
     subs = Dict()
@@ -348,7 +357,21 @@ function do_bareiss!(M, Mold, is_linear_variables, is_highest_diff)
     (rank1, rank2, rank3, pivots)
 end
 
-function alias_eliminate_graph!(state::TransformationState, ils::SparseMatrixCLIL)
+function force_var_to_zero!(structure::SystemStructure, ils::SparseMatrixCLIL, v::Int)
+    @unpack graph, solvable_graph, eq_to_diff = structure
+    @set! ils.nparentrows += 1
+    push!(ils.nzrows, ils.nparentrows)
+    push!(ils.row_cols, [v])
+    push!(ils.row_vals, [convert(eltype(ils), 1)])
+    add_vertex!(graph, SRC)
+    add_vertex!(solvable_graph, SRC)
+    add_edge!(graph, ils.nparentrows, v)
+    add_edge!(solvable_graph, ils.nparentrows, v)
+    add_vertex!(eq_to_diff)
+end
+
+function alias_eliminate_graph!(state::TransformationState, ils::SparseMatrixCLIL;
+    variable_underconstrained! = force_var_to_zero!)
     @unpack structure = state
     @unpack graph, solvable_graph, var_to_diff, eq_to_diff = state.structure
     # Step 1: Perform Bareiss factorization on the adjacency matrix of the linear
@@ -360,15 +383,7 @@ function alias_eliminate_graph!(state::TransformationState, ils::SparseMatrixCLI
     rk1vars = BitSet(@view pivots[1:rank1])
     for v in solvable_variables
         v in rk1vars && continue
-        @set! ils.nparentrows += 1
-        push!(ils.nzrows, ils.nparentrows)
-        push!(ils.row_cols, [v])
-        push!(ils.row_vals, [convert(eltype(ils), 1)])
-        add_vertex!(graph, SRC)
-        add_vertex!(solvable_graph, SRC)
-        add_edge!(graph, ils.nparentrows, v)
-        add_edge!(solvable_graph, ils.nparentrows, v)
-        add_vertex!(eq_to_diff)
+        variable_underconstrained!(structure, ils, v)
     end
 
     return ils