Merge pull request SciML#2621 from SciML/myb/nonSI

Add support for symbolic units by normalizing them to SI units

Author: ChrisRackauckas

src/systems/unit_check.jl

@@ -48,14 +48,11 @@ function screen_unit(result)
     if result isa DQ.AbstractQuantity
         d = DQ.dimension(result)
         if d isa DQ.Dimensions
-            if result != oneunit(result)
-                throw(ValidationError("$result uses non SI unit. Please use SI unit only."))
-            end
             return result
         elseif d isa DQ.SymbolicDimensions
-            throw(ValidationError("$result uses SymbolicDimensions, please use `u\"m\"` to instantiate SI unit only."))
+            return DQ.uexpand(oneunit(result))
         else
-            throw(ValidationError("$result doesn't use SI unit, please use `u\"m\"` to instantiate SI unit only."))
+            throw(ValidationError("$result doesn't have a recognized unit"))
         end
     else
         throw(ValidationError("$result doesn't have any unit."))
@@ -69,7 +66,7 @@ get_literal_unit(x) = screen_unit(something(__get_literal_unit(x), unitless))
 Find the unit of a symbolic item.
 """
 get_unit(x::Real) = unitless
-get_unit(x::DQ.AbstractQuantity) = screen_unit(oneunit(x))
+get_unit(x::DQ.AbstractQuantity) = screen_unit(x)
 get_unit(x::AbstractArray) = map(get_unit, x)
 get_unit(x::Num) = get_unit(unwrap(x))
 get_unit(op::Differential, args) = get_unit(args[1]) / get_unit(op.x)
@@ -81,12 +78,19 @@ get_unit(op::typeof(instream), args) = get_unit(args[1])
 function get_unit(op, args) # Fallback
     result = op(get_unit.(args)...)
     try
-        oneunit(result)
+        result
     catch
         throw(ValidationError("Unable to get unit for operation $op with arguments $args."))
     end
 end
 
+function get_unit(::Union{typeof(+), typeof(-)}, args)
+    u = get_unit(args[1])
+    if all(i -> get_unit(args[i]) == u, 2:length(args))
+        return u
+    end
+end
+
 function get_unit(op::Integral, args)
     unit = 1
     if op.domain.variables isa Vector
@@ -96,7 +100,7 @@ function get_unit(op::Integral, args)
     else
         unit *= get_unit(op.domain.variables)
     end
-    return oneunit(get_unit(args[1]) * unit)
+    return get_unit(args[1]) * unit
 end
 
 equivalent(x, y) = isequal(x, y)
@@ -197,7 +201,11 @@ function _validate(terms::Vector, labels::Vector{String}; info::String = "")
                 first_label = label
             elseif !equivalent(first_unit, equnit)
                 valid = false
-                @warn("$info: units [$(first_unit)] for $(first_label) and [$(equnit)] for $(label) do not match.")
+                str = "$info: units [$(first_unit)] for $(first_label) and [$(equnit)] for $(label) do not match."
+                if oneunit(first_unit) == oneunit(equnit)
+                    str *= " If there are non-SI units in the system, please use symbolic units like `us\"ms\"`"
+                end
+                @warn(str)
             end
         end
     end
@@ -227,7 +235,11 @@ function _validate(conn::Connection; info::String = "")
                 bunit = safe_get_unit(sst[j], info * string(nameof(s)) * "#$j")
                 if !equivalent(aunit, bunit)
                     valid = false
-                    @warn("$info: connected system unknowns $x and $(sst[j]) have mismatched units.")
+                    str = "$info: connected system unknowns $x ($aunit) and $(sst[j]) ($bunit) have mismatched units."
+                    if oneunit(aunit) == oneunit(bunit)
+                        str *= " If there are non-SI units in the system, please use symbolic units like `us\"ms\"`"
+                    end
+                    @warn(str)
                 end
             end
         end

test/dq_units.jl

@@ -13,10 +13,6 @@ using ModelingToolkit: t, D
 @test MT.get_unit(0.5) == MT.unitless
 @test MT.get_unit(MT.SciMLBase.NullParameters()) == MT.unitless
 
-# Prohibited unit types
-@parameters γ [unit = 1u"ms"]
-@test_throws MT.ValidationError MT.get_unit(γ)
-
 eqs = [D(E) ~ P - E / τ
        0 ~ P]
 @test MT.validate(eqs)
@@ -59,11 +55,18 @@ end
 end
 @named p1 = Pin()
 @named p2 = Pin()
-@test_throws MT.ValidationError @named op = OtherPin()
 @named lp = LongPin()
 good_eqs = [connect(p1, p2)]
 @test MT.validate(good_eqs)
 @named sys = ODESystem(good_eqs, t, [], [])
+@named op = OtherPin()
+bad_eqs = [connect(p1, op)]
+@test !MT.validate(bad_eqs)
+@test_throws MT.ValidationError @named sys = ODESystem(bad_eqs, t, [], [])
+@named op2 = OtherPin()
+good_eqs = [connect(op, op2)]
+@test MT.validate(good_eqs)
+@named sys = ODESystem(good_eqs, t, [], [])
 
 # Array variables
 @variables x(t)[1:3] [unit = u"m"]
@@ -85,18 +88,22 @@ eqs = [
 eqs = [D(E) ~ P - E / τ
        P ~ Q]
 
+noiseeqs = [0.1us"W",
+    0.1us"W"]
+@named sys = SDESystem(eqs, noiseeqs, t, [P, E], [τ, Q])
+
 noiseeqs = [0.1u"W",
     0.1u"W"]
-@named sys = SDESystem(eqs, noiseeqs, t, [P, E], [τ, Q])
+@test_throws MT.ValidationError @named sys = SDESystem(eqs, noiseeqs, t, [P, E], [τ, Q])
 
 # With noise matrix
-noiseeqs = [0.1u"W" 0.1u"W"
-            0.1u"W" 0.1u"W"]
+noiseeqs = [0.1us"W" 0.1us"W"
+            0.1us"W" 0.1us"W"]
 @named sys = SDESystem(eqs, noiseeqs, t, [P, E], [τ, Q])
 
 # Invalid noise matrix
-noiseeqs = [0.1u"W" 0.1u"W"
-            0.1u"W" 0.1u"s"]
+noiseeqs = [0.1us"W" 0.1us"W"
+            0.1us"W" 0.1us"s"]
 @test !MT.validate(eqs, noiseeqs)
 
 # Non-trivial simplifications
@@ -210,3 +217,10 @@ end
     @test prob = ODEProblem(
         pend, u0, (0.0, 1.0), p; guesses = guess, check_units = false) isa Any
 end
+
+@parameters p [unit = u"L/s"] d [unit = u"s^(-1)"]
+@parameters tt [unit = u"s"]
+@variables X(tt) [unit = u"L"]
+DD = Differential(tt)
+eqs = [DD(X) ~ p - d * X + d * X]
+@test ModelingToolkit.validate(eqs)