WIP

Author: YingboMa

Project.toml

@@ -17,6 +17,7 @@ Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 DomainSets = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
+DynamicQuantities = "06fc5a27-2a28-4c7c-a15d-362465fb6821"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 FunctionWrappersWrappers = "77dc65aa-8811-40c2-897b-53d922fa7daf"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
@@ -71,6 +72,7 @@ DiffRules = "0.1, 1.0"
 Distributions = "0.23, 0.24, 0.25"
 DocStringExtensions = "0.7, 0.8, 0.9"
 DomainSets = "0.6"
+DynamicQuantities = "0.8"
 ForwardDiff = "0.10.3"
 FunctionWrappersWrappers = "0.1"
 Graphs = "1.5.2"

src/ModelingToolkit.jl

@@ -152,6 +152,7 @@ include("systems/pde/pdesystem.jl")
 
 include("systems/sparsematrixclil.jl")
 include("systems/discrete_system/discrete_system.jl")
+include("systems/unit_check.jl")
 include("systems/validation.jl")
 include("systems/dependency_graphs.jl")
 include("clock.jl")

src/systems/diffeqs/odesystem.jl

@@ -163,7 +163,8 @@ struct ODESystem <: AbstractODESystem
             check_equations(equations(cevents), iv)
         end
         if checks == true || (checks & CheckUnits) > 0
-            all_dimensionless([dvs; ps; iv]) || check_units(deqs)
+            u = __get_unit_type(dvs, ps, iv)
+            check_units(u, deqs)
         end
         new(tag, deqs, iv, dvs, ps, tspan, var_to_name, ctrls, observed, tgrad, jac,
             ctrl_jac, Wfact, Wfact_t, name, systems, defaults, torn_matching,

src/systems/unit_check.jl

@@ -0,0 +1,42 @@
+import DynamicQuantities
+const DQ = DynamicQuantities
+
+struct ValidationError <: Exception
+    message::String
+end
+
+check_units(::Nothing, _...) = true
+
+__get_literal_unit(x) = getmetadata(x, VariableUnit, nothing)
+function __get_unit_type(vs′...)
+    vs = Iterators.flatten(vs′)
+    for v in vs
+        u = __get_literal_unit(v)
+        if u isa DQ.AbstractQuantity
+            return Val(:Unitful)
+        else
+            return Val(:DynamicQuantities)
+        end
+    end
+    return nothing
+end
+
+function check_units(::Val{:DynamicQuantities}, eqs...)
+    validate(eqs...) ||
+        throw(ValidationError("Some equations had invalid units. See warnings for details."))
+end
+
+function screen_units(result)
+    if result isa DQ.AbstractQuantity
+        d = DQ.dimension(result)
+        if d isa DQ.Dimensions
+            return result
+        elseif d isa DQ.SymbolicDimensions
+            throw(ValidationError("$result uses SymbolicDimensions, please use `u\"m\"` to instantiate SI unit only."))
+        else
+            throw(ValidationError("$result doesn't use SI unit, please use `u\"m\"` to instantiate SI unit only."))
+        end
+    end
+end
+
+

src/systems/validation.jl

@@ -1,10 +1,12 @@
+module UnitfulUnitCheck
+
+using .ModelingToolkit, Symbolics, SciMLBase
+using .ModelingToolkit: ValidationError
+const MT = ModelingToolkit
+
 Base.:*(x::Union{Num, Symbolic}, y::Unitful.AbstractQuantity) = x * y
 Base.:/(x::Union{Num, Symbolic}, y::Unitful.AbstractQuantity) = x / y
 
-struct ValidationError <: Exception
-    message::String
-end
-
 """
 Throw exception on invalid unit types, otherwise return argument.
 """
@@ -60,7 +62,11 @@ get_literal_unit(x) = screen_unit(getmetadata(x, VariableUnit, unitless))
 function get_unit(op, args) # Fallback
     result = op(1 .* get_unit.(args)...)
     try
-        unit(result)
+        if result isa DQ.AbstractQuantity
+            oneunit(result)
+        else
+            unit(result)
+        end
     catch
         throw(ValidationError("Unable to get unit for operation $op with arguments $args."))
     end
@@ -211,15 +217,15 @@ function _validate(conn::Connection; info::String = "")
     valid
 end
 
-function validate(jump::Union{ModelingToolkit.VariableRateJump,
-            ModelingToolkit.ConstantRateJump}, t::Symbolic;
+function validate(jump::Union{MT.VariableRateJump,
+            MT.ConstantRateJump}, t::Symbolic;
         info::String = "")
     newinfo = replace(info, "eq." => "jump")
     _validate([jump.rate, 1 / t], ["rate", "1/t"], info = newinfo) && # Assuming the rate is per time units
         validate(jump.affect!, info = newinfo)
 end
 
-function validate(jump::ModelingToolkit.MassActionJump, t::Symbolic; info::String = "")
+function validate(jump::MT.MassActionJump, t::Symbolic; info::String = "")
     left_symbols = [x[1] for x in jump.reactant_stoch] #vector of pairs of symbol,int -> vector symbols
     net_symbols = [x[1] for x in jump.net_stoch]
     all_symbols = vcat(left_symbols, net_symbols)
@@ -235,18 +241,18 @@ function validate(jumps::ArrayPartition{<:Union{Any, Vector{<:JumpType}}}, t::Sy
     all([validate(jumps.x[idx], t, info = labels[idx]) for idx in 1:3])
 end
 
-function validate(eq::ModelingToolkit.Equation; info::String = "")
+function validate(eq::MT.Equation; info::String = "")
     if typeof(eq.lhs) == Connection
         _validate(eq.rhs; info)
     else
         _validate([eq.lhs, eq.rhs], ["left", "right"]; info)
     end
 end
-function validate(eq::ModelingToolkit.Equation,
+function validate(eq::MT.Equation,
         term::Union{Symbolic, Unitful.Quantity, Num}; info::String = "")
     _validate([eq.lhs, eq.rhs, term], ["left", "right", "noise"]; info)
 end
-function validate(eq::ModelingToolkit.Equation, terms::Vector; info::String = "")
+function validate(eq::MT.Equation, terms::Vector; info::String = "")
     _validate(vcat([eq.lhs, eq.rhs], terms),
         vcat(["left", "right"], "noise  #" .* string.(1:length(terms))); info)
 end
@@ -273,8 +279,10 @@ validate(term::Symbolics.SymbolicUtils.Symbolic) = safe_get_unit(term, "") !== n
 """
 Throws error if units of equations are invalid.
 """
-function check_units(eqs...)
+function MT.check_units(::Val{:Unitful}, eqs...)
     validate(eqs...) ||
         throw(ValidationError("Some equations had invalid units. See warnings for details."))
 end
 all_dimensionless(states) = all(x -> safe_get_unit(x, "") in (unitless, nothing), states)
+
+end # module

test/runtests.jl

@@ -13,7 +13,14 @@ using SafeTestsets, Test
 @safetestset "Clock Test" include("clock.jl")
 @safetestset "DiscreteSystem Test" include("discretesystem.jl")
 @safetestset "ODESystem Test" include("odesystem.jl")
-@safetestset "Unitful Quantities Test" include("units.jl")
+@safetestset "Dynamic Quantities Test" begin
+    using DynamicQuantities
+    include("units.jl")
+end
+@safetestset "Unitful Quantities Test" begin
+    using Unitful
+    include("units.jl")
+end
 @safetestset "LabelledArrays Test" include("labelledarrays.jl")
 @safetestset "Mass Matrix Test" include("mass_matrix.jl")
 @safetestset "SteadyStateSystem Test" include("steadystatesystems.jl")

test/units.jl

@@ -1,4 +1,4 @@
-using ModelingToolkit, Unitful, OrdinaryDiffEq, JumpProcesses, IfElse
+using ModelingToolkit, OrdinaryDiffEq, JumpProcesses, IfElse
 using Test
 MT = ModelingToolkit
 @parameters τ [unit = u"ms"] γ