Adding validation of Unitful quantities. Doesn't handle symbolic arrays yet.

Author: lamorton

src/systems/diffeqs/odesystem.jl

@@ -88,6 +88,7 @@ struct ODESystem <: AbstractODESystem
         check_variables(dvs,iv)
         check_parameters(ps,iv)
         check_equations(deqs,iv)
+        validate(deqs)
         new(deqs, iv, dvs, ps, var_to_name, ctrls, observed, tgrad, jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults, structure, connection_type)
     end
 end

src/systems/diffeqs/sdesystem.jl

@@ -90,6 +90,7 @@ struct SDESystem <: AbstractODESystem
         check_variables(dvs,iv)
         check_parameters(ps,iv)
         check_equations(deqs,iv)
+        validate(deqs)
         new(deqs, neqs, iv, dvs, ps, var_to_name, ctrls, observed, tgrad, jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults, connection_type)
     end
 end

src/systems/diffeqs/validation.jl

@@ -1,26 +1,72 @@
 Base.:*(x::Union{Num,Symbolic},y::Unitful.AbstractQuantity) = x * y
 
-instantiate(x::Sym{Real}) = 1.0
-instantiate(x::Symbolic) = oneunit(1*ModelingToolkit.vartype(x))
-function instantiate(x::Num)
-    x = value(x)
-    if operation(x) isa Sym
-        return instantiate(operation(x))
-    elseif operation(x) isa Differential
-        instantiate(arguments(x)[1])/instantiate(arguments(x)[1].args[1])
+
+function vartype(x::Symbolic)
+    if !(x.metadata isa Nothing)
+        return haskey(x.metadata,VariableUnit) ? x.metadata[VariableUnit] : 1.0
+    end
+    1.0
+end
+vartype(x::Num) = vartype(value(x))
+
+instantiate(x) = 1.0
+instantiate(x::Num) = instantiate(value(x))
+function instantiate(x::Symbolic)
+    vx = value(x)
+    if vx isa Sym || operation(vx) isa Sym
+        return oneunit(1 * ModelingToolkit.vartype(x))
+    elseif operation(vx) isa Differential
+        return instantiate(arguments(vx)[1]) / instantiate(arguments(arguments(vx)[1])[1])
+    elseif vx isa Pow
+        pargs = arguments(vx)
+        base,expon = instantiate.(pargs)
+        uconvert(NoUnits, expon) # This acts as an assertion
+        return base == 1.0 ? 1.0 : operation(vx)(base, pargs[2])
+    elseif vx isa Add # Cannot simply add the units b/c they may differ in magnitude (eg, kg vs g)
+        terms = instantiate.(arguments(vx))
+        firstunit = unit(terms[1])
+        @assert all(map(x -> ustrip(firstunit, x) == 1, terms[2:end]))
+        return 1.0 * firstunit
     else
-        operation(x)(instantiate.(arguments(x))...)
+        return oneunit(operation(vx)(instantiate.(arguments(vx))...))
     end
 end
 
-function validate(eq::ModelingToolkit.Equation)
+function validate(eq::ModelingToolkit.Equation; eqnum = 1)
+    lhs = rhs = nothing
     try
-        return typeof(instantiate(eq.lhs)) == typeof(instantiate(eq.rhs))
-    catch
-        return false
+        lhs = instantiate(eq.lhs)
+    catch err
+        if err isa Unitful.DimensionError 
+            @warn("In left-hand side of eq. #$eqnum: $(eq.lhs),  $(err.x) and $(err.y) are not dimensionally compatible.")
+        elseif err isa MethodError
+            @warn("In right-hand side of eq. #$eqnum: $(err.f) doesn't accept $(err.args).")
+        else
+            rethrow()
+        end
     end
+    try
+        rhs = instantiate(eq.rhs)
+    catch err
+        if err isa Unitful.DimensionError
+            @warn("In right-hand side of eq. #$eqnum: $(eq.rhs), $(err.x) and $(err.y) are not dimensionally compatible.")
+        elseif err isa MethodError
+            @warn("In right-hand side of eq. #$eqnum: $(err.f) doesn't accept $(err.args).")
+        else
+            rethrow()
+        end
+    end
+    if (rhs !== nothing) && (lhs !== nothing)
+        if !isequal(lhs, rhs)
+            @warn("In eq. #$eqnum, left-side units ($lhs) and right-side units ($rhs) don't match.")
+        end
+    end
+    (rhs !== nothing) && (lhs !== nothing) && isequal(lhs, rhs)
 end
 
-function validate(sys::AbstractODESystem)
-    all(validate.(equations(sys)))
+function validate(eqs::Vector{ModelingToolkit.Equation})
+    correct = [validate(eqs[idx],eqnum=idx) for idx in 1:length(eqs)]
+    all(correct) || throw(ArgumentError("Invalid equations, see warnings for details."))
 end
+
+validate(sys::AbstractODESystem) = validate(equations(sys))

test/runtests.jl

@@ -8,6 +8,7 @@ using SafeTestsets, Test
 @safetestset "System Linearity Test" begin include("linearity.jl") end
 @safetestset "DiscreteSystem Test" begin include("discretesystem.jl") end
 @safetestset "ODESystem Test" begin include("odesystem.jl") end
+@safetestset "Unitful Quantities Test" begin include("units.jl") end
 @safetestset "LabelledArrays Test" begin include("labelledarrays.jl") end
 @safetestset "Mass Matrix Test" begin include("mass_matrix.jl") end
 @safetestset "SteadyStateSystem Test" begin include("steadystatesystems.jl") end

test/units.jl

@@ -1,25 +1,58 @@
 using ModelingToolkit, Unitful
 using Test
-
-t = Variable{u"s"}(:t)()
-x = Variable{u"kg"}(:x)(t)
-y = Variable{u"kg"}(:y)(t)
+MT = ModelingToolkit
+@parameters τ [unit=u"ms"]
+@variables t [unit=u"ms"] E(t) [unit=u"kJ"] P(t) [unit=u"MW"]
 D = Differential(t)
 
-eq1 = x ~ y*t
-eq2 = x*10u"s" ~ y*t
+@test MT.vartype(t) == u"ms"
+@test MT.vartype(E) == u"kJ"
+@test MT.vartype(τ) == u"ms"
+
+eqs = [D(E) ~ P-E/τ ]
+sys = ODESystem(eqs)
+
+@test MT.instantiate(eqs[1].lhs) == 1.0u"MW"
+@test MT.instantiate(eqs[1].rhs) == 1.0u"MW"
+@test MT.validate(eqs[1])
+@test MT.validate(sys)
+
+@test MT.instantiate(0.5) == 1.0
+@test MT.instantiate(t) == 1.0u"ms"
+@test MT.instantiate(P) == 1.0u"MW"
+@test MT.instantiate(τ) == 1.0u"ms"
 
-@test ModelingToolkit.instantiate(t) == 1u"s"
-@test ModelingToolkit.instantiate(x) == 1u"kg"
-@test ModelingToolkit.instantiate(y) == 1u"kg"
+@test MT.instantiate(τ^-1) == 1/u"ms"
+@test MT.instantiate(D(E)) == 1.0u"MW"
+@test MT.instantiate(E/τ) == 1.0u"MW" 
+@test MT.instantiate(2*P) == 1.0u"MW"
+@test MT.instantiate(t/τ) == 1.0
+@test MT.instantiate(P-E/τ)/1.0u"MW" == 1.0
 
-@test !ModelingToolkit.validate(eq1)
-@test ModelingToolkit.validate(eq2)
+@test MT.instantiate(1.0^(t/τ)) == 1.0
+@test MT.instantiate(exp(t/τ)) == 1.0
+@test MT.instantiate(sin(t/τ)) == 1.0
+@test MT.instantiate(sin(1u"rad")) == 1.0
+@test MT.instantiate(t^2) == 1.0u"ms"^2
 
-eqs = [
-        D(x) ~ y/t
-        D(y) ~ (x*y)/(t*10u"kg")
-]
+@test !MT.validate(E^1.5~ E^(t/τ))
+@test MT.validate(E^(t/τ)~ E^(t/τ))
 
-sys = ODESystem(eqs,t,[x,y],[])
-@test ModelingToolkit.validate(sys)
+sys = ODESystem(eqs,t,[P,E],[τ])
+@test MT.validate(sys)
+
+@test !MT.validate(D(D(E))~P)
+@test !MT.validate(0~P+E*τ)
+@test_logs (:warn,) MT.validate(0 ~ P + E*τ)
+@test_logs (:warn,) MT.validate(P + E*τ ~ 0)
+@test_logs (:warn,) MT.validate(P ~ 0)
+
+@variables x y z u
+@parameters σ ρ β
+eqs = [0 ~ σ*(y-x)]
+@test MT.validate(eqs) #should cope with unit-free
+
+@variables t x[1:3,1:3](t) #should cope with arrays
+D = Differential(t)
+eqs = D.(x) .~ x
+ODESystem(eqs)