Fix #2959

Author: YingboMa

src/structural_transformation/StructuralTransformations.jl

@@ -11,7 +11,7 @@ using SymbolicUtils: maketerm, iscall
 
 using ModelingToolkit
 using ModelingToolkit: ODESystem, AbstractSystem, var_from_nested_derivative, Differential,
-                       unknowns, equations, vars, Symbolic, diff2term, value,
+                       unknowns, equations, vars, Symbolic, diff2term_with_unit, value,
                        operation, arguments, Sym, Term, simplify, symbolic_linear_solve,
                        isdiffeq, isdifferential, isirreducible,
                        empty_substitutions, get_substitutions,

src/structural_transformation/symbolics_tearing.jl

@@ -274,7 +274,7 @@ function tearing_reassemble(state::TearingState, var_eq_matching,
         dv === nothing && continue
         if var_eq_matching[var] !== SelectedState()
             dd = fullvars[dv]
-            v_t = setio(diff2term(unwrap(dd)), false, false)
+            v_t = setio(diff2term_with_unit(unwrap(dd), unwrap(iv)), false, false)
             for eq in 𝑑neighbors(graph, dv)
                 dummy_sub[dd] = v_t
                 neweqs[eq] = fast_substitute(neweqs[eq], dd => v_t)

src/utils.jl

@@ -859,3 +859,13 @@ function eval_or_rgf(expr::Expr; eval_expression = false, eval_module = @__MODUL
         return drop_expr(RuntimeGeneratedFunction(eval_module, eval_module, expr))
     end
 end
+
+function diff2term_with_unit(x, t)
+    x = diff2term(x)
+    if hasmetadata(x, VariableUnit) && (t isa Symbolic && hasmetadata(t, VariableUnit))
+        xu = getmetadata(x, VariableUnit)
+        tu = getmetadata(t, VariableUnit)
+        x = setmetadata(x, VariableUnit, xu / tu)
+    end
+    return x
+end