Basic functionality is implemented.

Author: lamorton

src/ModelingToolkit.jl

@@ -115,6 +115,7 @@ using .BipartiteGraphs
 
 include("variables.jl")
 include("parameters.jl")
+include("constants.jl")
 
 include("utils.jl")
 include("domains.jl")
@@ -206,7 +207,7 @@ export toexpr, get_variables
 export simplify, substitute
 export build_function
 export modelingtoolkitize
-export @variables, @parameters
+export @variables, @parameters, @constants
 export @named, @nonamespace, @namespace, extend, compose
 
 end # module

src/constants.jl

@@ -5,23 +5,23 @@ function isconstant(x)
     x = unwrap(x)
     x isa Symbolic &&  getmetadata(x, MTKConstantCtx, false)
 end
-
+isconstant(x::Num) = isconstant(unwrap(x))
 """
     toconst(s::Sym)
 
 Maps the parameter to a constant. The parameter must have a default.
 """
-function toconst(s)
+function toconstant(s)
     if s isa Symbolics.Arr
-        Symbolics.wrap(toconst(Symbolics.unwrap(s)))
+        Symbolics.wrap(toconstant(Symbolics.unwrap(s)))
     elseif s isa AbstractArray
-        map(toconst, s)
+        map(toconstant, s)
     else
-        assert(hasmetadata(s,VariableDefaultValue))
-        setmetadata(s, MTKConstCtx, true)
+        hasmetadata(s, Symbolics.VariableDefaultValue) || throw(ArgumentError("Constant `$(s)` must be assigned a default value."))
+        setmetadata(s, MTKConstantCtx, true)
     end
 end
-toconst(s::Num) = wrap(toconst(value(s)))
+toconstant(s::Num) = wrap(toconstant(value(s)))
 
 """
 $(SIGNATURES)
@@ -32,5 +32,5 @@ macro constants(xs...)
     Symbolics._parse_vars(:constants,
                           Real,
                           xs,
-                          toconst) |> esc
+                          toconstant) |> esc
 end

src/systems/abstractsystem.jl

@@ -167,6 +167,7 @@ for prop in [:eqs
              :iv
              :states
              :ps
+             :cs
              :var_to_name
              :ctrls
              :defaults
@@ -376,6 +377,7 @@ end
 namespace_variables(sys::AbstractSystem) = states(sys, states(sys))
 namespace_parameters(sys::AbstractSystem) = parameters(sys, parameters(sys))
 namespace_controls(sys::AbstractSystem) = controls(sys, controls(sys))
+namespace_constants(sys::AbstractSystem) = constants(sys, constants(sys))
 
 function namespace_defaults(sys)
     defs = defaults(sys)
@@ -437,6 +439,12 @@ function parameters(sys::AbstractSystem)
     unique(isempty(systems) ? ps : [ps; reduce(vcat, namespace_parameters.(systems))])
 end
 
+function constants(sys::AbstractSystem)
+    cs = get_cs(sys)
+    systems = get_systems(sys)
+    unique(isempty(systems) ? cs : [cs; reduce(vcat, namespace_constants.(systems))])
+end
+
 function controls(sys::AbstractSystem)
     ctrls = get_ctrls(sys)
     systems = get_systems(sys)

src/systems/diffeqs/abstractodesystem.jl

@@ -127,6 +127,13 @@ function generate_function(sys::AbstractODESystem, dvs = states(sys), ps = param
     rhss = implicit_dae ? [_iszero(eq.lhs) ? eq.rhs : eq.rhs - eq.lhs for eq in eqs] :
            [eq.rhs for eq in eqs]
 
+    # Swap constants for their values
+    cs = constants(sys)
+    if !isempty(cs) > 0
+        cmap = map(x -> x => getdefault(x), cs)
+        rhss = map(x -> substitute(x, cmap), rhss)
+    end
+
     # TODO: add an optional check on the ordering of observed equations
     u = map(x -> time_varying_as_func(value(x), sys), dvs)
     p = map(x -> time_varying_as_func(value(x), sys), ps)

src/systems/diffeqs/odesystem.jl

@@ -36,6 +36,8 @@ struct ODESystem <: AbstractODESystem
     states::Vector
     """Parameter variables. Must not contain the independent variable."""
     ps::Vector
+    """Symbolic constants."""
+    cs::Vector
     """Array variables."""
     var_to_name::Any
     """Control parameters (some subset of `ps`)."""
@@ -120,7 +122,7 @@ struct ODESystem <: AbstractODESystem
     """
     metadata::Any
 
-    function ODESystem(deqs, iv, dvs, ps, var_to_name, ctrls, observed, tgrad,
+    function ODESystem(deqs, iv, dvs, ps, cs, var_to_name, ctrls, observed, tgrad,
                        jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults,
                        torn_matching, connector_type, connections, preface, cevents,
                        devents, tearing_state = nothing, substitutions = nothing,
@@ -133,16 +135,16 @@ struct ODESystem <: AbstractODESystem
             check_equations(equations(cevents), iv)
         end
         if checks == true || (checks & CheckUnits) > 0
-            all_dimensionless([dvs; ps; iv]) || check_units(deqs)
+            all_dimensionless([dvs; ps; iv; cs]) || check_units(deqs)
         end
-        new(deqs, iv, dvs, ps, var_to_name, ctrls, observed, tgrad, jac,
+        new(deqs, iv, dvs, ps, cs, var_to_name, ctrls, observed, tgrad, jac,
             ctrl_jac, Wfact, Wfact_t, name, systems, defaults, torn_matching,
             connector_type, connections, preface, cevents, devents, tearing_state,
             substitutions, metadata)
     end
 end
 
-function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps;
+function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps, cs;
                    controls = Num[],
                    observed = Equation[],
                    systems = ODESystem[],
@@ -164,6 +166,7 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps;
     iv′ = value(scalarize(iv))
     dvs′ = value.(scalarize(dvs))
     ps′ = value.(scalarize(ps))
+    cs′ = value.(scalarize(cs))
     ctrl′ = value.(scalarize(controls))
 
     if !(isempty(default_u0) && isempty(default_p))
@@ -176,6 +179,7 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps;
     var_to_name = Dict()
     process_variables!(var_to_name, defaults, dvs′)
     process_variables!(var_to_name, defaults, ps′)
+    process_variables!(var_to_name, defaults, cs′)
     isempty(observed) || collect_var_to_name!(var_to_name, (eq.lhs for eq in observed))
 
     tgrad = RefValue(EMPTY_TGRAD)
@@ -189,7 +193,7 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps;
     end
     cont_callbacks = SymbolicContinuousCallbacks(continuous_events)
     disc_callbacks = SymbolicDiscreteCallbacks(discrete_events)
-    ODESystem(deqs, iv′, dvs′, ps′, var_to_name, ctrl′, observed, tgrad, jac,
+    ODESystem(deqs, iv′, dvs′, ps′, cs′, var_to_name, ctrl′, observed, tgrad, jac,
               ctrl_jac, Wfact, Wfact_t, name, systems, defaults, nothing,
               connector_type, nothing, preface, cont_callbacks, disc_callbacks,
               metadata, checks = checks)
@@ -201,6 +205,7 @@ function ODESystem(eqs, iv = nothing; kwargs...)
     diffvars = OrderedSet()
     allstates = OrderedSet()
     ps = OrderedSet()
+    cs = OrderedSet() #Constants
     # reorder equations such that it is in the form of `diffeq, algeeq`
     diffeq = Equation[]
     algeeq = Equation[]
@@ -218,8 +223,8 @@ function ODESystem(eqs, iv = nothing; kwargs...)
     compressed_eqs = Equation[] # equations that need to be expanded later, like `connect(a, b)`
     for eq in eqs
         eq.lhs isa Union{Symbolic, Number} || (push!(compressed_eqs, eq); continue)
-        collect_vars!(allstates, ps, eq.lhs, iv)
-        collect_vars!(allstates, ps, eq.rhs, iv)
+        collect_vars!(allstates, ps, cs, eq.lhs, iv)
+        collect_vars!(allstates, ps, cs, eq.rhs, iv)
         if isdiffeq(eq)
             diffvar, _ = var_from_nested_derivative(eq.lhs)
             isequal(iv, iv_from_nested_derivative(eq.lhs)) ||
@@ -235,7 +240,7 @@ function ODESystem(eqs, iv = nothing; kwargs...)
     algevars = setdiff(allstates, diffvars)
     # the orders here are very important!
     return ODESystem(Equation[diffeq; algeeq; compressed_eqs], iv,
-                     collect(Iterators.flatten((diffvars, algevars))), ps; kwargs...)
+                     collect(Iterators.flatten((diffvars, algevars))), ps, cs; kwargs...)
 end
 
 # NOTE: equality does not check cached Jacobian
@@ -260,6 +265,7 @@ function flatten(sys::ODESystem, noeqs = false)
                          get_iv(sys),
                          states(sys),
                          parameters(sys),
+                         constants(sys),
                          observed = observed(sys),
                          continuous_events = continuous_events(sys),
                          discrete_events = discrete_events(sys),

src/utils.jl

@@ -466,38 +466,40 @@ function find_derivatives!(vars, expr, f)
     return vars
 end
 
-function collect_vars!(states, parameters, expr, iv)
+function collect_vars!(states, parameters, constants, expr, iv)
     if expr isa Sym
-        collect_var!(states, parameters, expr, iv)
+        collect_var!(states, parameters, constants, expr, iv)
     else
         for var in vars(expr)
             if istree(var) && operation(var) isa Differential
                 var, _ = var_from_nested_derivative(var)
             end
-            collect_var!(states, parameters, var, iv)
+            collect_var!(states, parameters, constants, var, iv)
         end
     end
     return nothing
 end
 
-function collect_vars_difference!(states, parameters, expr, iv)
+function collect_vars_difference!(states, parameters, constants, expr, iv)
     if expr isa Sym
-        collect_var!(states, parameters, expr, iv)
+        collect_var!(states, parameters, constants, expr, iv)
     else
         for var in vars(expr)
             if istree(var) && operation(var) isa Difference
                 var, _ = var_from_nested_difference(var)
             end
-            collect_var!(states, parameters, var, iv)
+            collect_var!(states, parameters, constants, var, iv)
         end
     end
     return nothing
 end
 
-function collect_var!(states, parameters, var, iv)
+function collect_var!(states, parameters, constants, var, iv)
     isequal(var, iv) && return nothing
     if isparameter(var) || (istree(var) && isparameter(operation(var)))
         push!(parameters, var)
+    elseif isconstant(var)
+        push!(constants,var)
     else
         push!(states, var)
     end

test/constants.jl

@@ -0,0 +1,14 @@
+using ModelingToolkit, OrdinaryDiffEq
+using Test
+MT = ModelingToolkit
+
+@constants a = 1
+@test_throws MT.MissingDefaultError @constants b
+
+@variables t x(t) w(t)
+D = Differential(t)
+eqs = [D(x) ~ a]
+@named sys = ODESystem(eqs)
+prob = ODEProblem(sys, [0, ], [0.0, 1.0],[])
+sol = solve(prob,Tsit5())
+

test/runtests.jl

@@ -48,4 +48,4 @@ println("Last test requires gcc available in the path!")
 @safetestset "FuncAffect Test" begin include("funcaffect.jl") end
 
 # Reference tests go Last
-@safetestset "Latexify recipes Test" begin include("latexify.jl") end
+#@safetestset "Latexify recipes Test" begin include("latexify.jl") end