Give a warning for underconstrained variables forced to 0

@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

@@ -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