Merge pull request #3036 from SciML/myb/const_unit

We cannot fold unitful constants

src/structural_transformation/utils.jl

@@ -203,7 +203,9 @@ function find_eq_solvables!(state::TearingState, ieq, to_rm = Int[], coeffs = no
             all_int_vars = false
             if !allow_symbolic
                 if allow_parameter
-                    all(ModelingToolkit.isparameter, vars(a)) || continue
+                    all(
+                        x -> ModelingToolkit.isparameter(x) || ModelingToolkit.isconstant(x),
+                        vars(a)) || continue
                 else
                     continue
                 end

src/utils.jl

@@ -836,7 +836,11 @@ function fold_constants(ex)
         maketerm(typeof(ex), operation(ex), map(fold_constants, arguments(ex)),
             metadata(ex))
     elseif issym(ex) && isconstant(ex)
-        getdefault(ex)
+        if (unit = getmetadata(ex, VariableUnit, nothing); unit !== nothing)
+            ex # we cannot fold constant with units
+        else
+            getdefault(ex)
+        end
     else
         ex
     end

test/dq_units.jl

@@ -226,6 +226,16 @@ DD = Differential(tt)
 eqs = [DD(X) ~ p - d * X + d * X]
 @test ModelingToolkit.validate(eqs)
 
+@constants begin
+    to_m = 1, [unit = u"m"]
+end
+@variables begin
+    L(t), [unit = u"m"]
+    L_out(t), [unit = u"1"]
+end
+@test to_m in ModelingToolkit.vars(ModelingToolkit.fold_constants(Symbolics.unwrap(L_out *
+                                                                                   -to_m)))
+
 # test units for registered functions
 let
     mm(X, v, K) = v * X / (X + K)