Merge pull request #2442 from AayushSabharwal/as/common-t-D

feat!: add common definitions of t and D

examples/electrical_components.jl

@@ -1,7 +1,7 @@
 using Test
 using ModelingToolkit, OrdinaryDiffEq
+using ModelingToolkit: t_nounits as t, D_nounits as D
 
-@isdefined(t) || @parameters t
 @connector function Pin(; name)
     sts = @variables v(t)=1.0 i(t)=1.0 [connect = Flow]
     ODESystem(Equation[], t, sts, []; name = name)
@@ -37,7 +37,6 @@ end
     @named oneport = OnePort()
     @unpack v, i = oneport
     ps = @parameters C = C
-    D = Differential(t)
     eqs = [
         D(v) ~ i / C,
     ]
@@ -58,7 +57,6 @@ end
     @named oneport = OnePort()
     @unpack v, i = oneport
     ps = @parameters L = L
-    D = Differential(t)
     eqs = [
         D(i) ~ v / L,
     ]
@@ -89,7 +87,6 @@ end
     @parameters rho=rho V=V cp=cp
     C = rho * V * cp
     @named h = HeatPort()
-    D = Differential(t)
     eqs = [
         D(h.T) ~ h.Q_flow / C,
     ]

src/ModelingToolkit.jl

@@ -74,6 +74,9 @@ end
 @reexport using UnPack
 RuntimeGeneratedFunctions.init(@__MODULE__)
 
+import DynamicQuantities, Unitful
+const DQ = DynamicQuantities
+
 export @derivatives
 
 for fun in [:toexpr]
@@ -172,11 +175,24 @@ for S in subtypes(ModelingToolkit.AbstractSystem)
     @eval convert_system(::Type{<:$S}, sys::$S) = sys
 end
 
+const t_nounits = let
+    only(@parameters t)
+end
+const t_unitful = let
+    only(@parameters t [unit = Unitful.u"s"])
+end
+const t = let
+    only(@parameters t [unit = DQ.u"s"])
+end
+
+const D_nounits = Differential(t_nounits)
+const D_unitful = Differential(t_unitful)
+const D = Differential(t)
+
 PrecompileTools.@compile_workload begin
     using ModelingToolkit
-    @variables t x(t)
-    D = Differential(t)
-    @named sys = ODESystem([D(x) ~ -x])
+    @variables x(ModelingToolkit.t_nounits)
+    @named sys = ODESystem([ModelingToolkit.D_nounits(x) ~ -x])
     prob = ODEProblem(structural_simplify(sys), [x => 30.0], (0, 100), [], jac = true)
 end
 

src/structural_transformation/codegen.jl

@@ -299,7 +299,6 @@ function build_torn_function(sys;
 
     unknown_vars = Any[fullvars[i] for i in unknowns_idxs]
     @set! sys.solved_unknowns = unknown_vars
-    syms = map(Symbol, unknown_vars)
 
     pre = get_postprocess_fbody(sys)
     cpre = get_preprocess_constants(rhss)
@@ -354,9 +353,6 @@ function build_torn_function(sys;
                 eqs_idxs,
                 unknowns_idxs) :
                        nothing,
-            syms = syms,
-            paramsyms = Symbol.(parameters(sys)),
-            indepsym = Symbol(get_iv(sys)),
             observed = observedfun,
             mass_matrix = mass_matrix,
             sys = sys),

src/systems/abstractsystem.jl

@@ -921,6 +921,7 @@ function toexpr(sys::AbstractSystem)
                 name = $name, checks = false)))
     end
 
+    expr = :(let; $expr; end)
     Base.remove_linenums!(expr) # keeping the line numbers is never helpful
 end
 

src/systems/diffeqs/abstractodesystem.jl

@@ -471,9 +471,6 @@ function DiffEqBase.ODEFunction{iip, specialize}(sys::AbstractODESystem, dvs = u
         tgrad = _tgrad === nothing ? nothing : _tgrad,
         mass_matrix = _M,
         jac_prototype = jac_prototype,
-        syms = collect(Symbol.(unknowns(sys))),
-        indepsym = Symbol(get_iv(sys)),
-        paramsyms = collect(Symbol.(ps)),
         observed = observedfun,
         sparsity = sparsity ? jacobian_sparsity(sys) : nothing,
         analytic = analytic)
@@ -569,9 +566,6 @@ function DiffEqBase.DAEFunction{iip}(sys::AbstractODESystem, dvs = unknowns(sys)
     DAEFunction{iip}(f,
         sys = sys,
         jac = _jac === nothing ? nothing : _jac,
-        syms = Symbol.(dvs),
-        indepsym = Symbol(get_iv(sys)),
-        paramsyms = Symbol.(ps),
         jac_prototype = jac_prototype,
         observed = observedfun)
 end
@@ -596,11 +590,7 @@ function DiffEqBase.DDEFunction{iip}(sys::AbstractODESystem, dvs = unknowns(sys)
     f(u, h, p, t) = f_oop(u, h, p, t)
     f(du, u, h, p, t) = f_iip(du, u, h, p, t)
 
-    DDEFunction{iip}(f,
-        sys = sys,
-        syms = Symbol.(dvs),
-        indepsym = Symbol(get_iv(sys)),
-        paramsyms = Symbol.(ps))
+    DDEFunction{iip}(f, sys = sys)
 end
 
 function DiffEqBase.SDDEFunction(sys::AbstractODESystem, args...; kwargs...)
@@ -628,11 +618,7 @@ function DiffEqBase.SDDEFunction{iip}(sys::AbstractODESystem, dvs = unknowns(sys
     g(u, h, p, t) = g_oop(u, h, p, t)
     g(du, u, h, p, t) = g_iip(du, u, h, p, t)
 
-    SDDEFunction{iip}(f, g,
-        sys = sys,
-        syms = Symbol.(dvs),
-        indepsym = Symbol(get_iv(sys)),
-        paramsyms = Symbol.(ps))
+    SDDEFunction{iip}(f, g, sys = sys)
 end
 
 """
@@ -718,9 +704,6 @@ function ODEFunctionExpr{iip}(sys::AbstractODESystem, dvs = unknowns(sys),
             tgrad = $tgradsym,
             mass_matrix = M,
             jac_prototype = $jp_expr,
-            syms = $(Symbol.(unknowns(sys))),
-            indepsym = $(QuoteNode(Symbol(get_iv(sys)))),
-            paramsyms = $(Symbol.(parameters(sys))),
             sparsity = $(sparsity ? jacobian_sparsity(sys) : nothing),
             observed = $observedfun_exp)
     end

src/systems/diffeqs/sdesystem.jl

@@ -474,9 +474,6 @@ function DiffEqBase.SDEFunction{iip}(sys::SDESystem, dvs = unknowns(sys),
         Wfact = _Wfact === nothing ? nothing : _Wfact,
         Wfact_t = _Wfact_t === nothing ? nothing : _Wfact_t,
         mass_matrix = _M,
-        syms = Symbol.(unknowns(sys)),
-        indepsym = Symbol(get_iv(sys)),
-        paramsyms = Symbol.(ps),
         observed = observedfun)
 end
 
@@ -563,10 +560,7 @@ function SDEFunctionExpr{iip}(sys::SDESystem, dvs = unknowns(sys),
             tgrad = tgrad,
             Wfact = Wfact,
             Wfact_t = Wfact_t,
-            mass_matrix = M,
-            syms = $(Symbol.(unknowns(sys))),
-            indepsym = $(Symbol(get_iv(sys))),
-            paramsyms = $(Symbol.(parameters(sys))))
+            mass_matrix = M)
     end
     !linenumbers ? Base.remove_linenums!(ex) : ex
 end

src/systems/jumps/jumpsystem.jl

@@ -324,10 +324,7 @@ function DiffEqBase.DiscreteProblem(sys::JumpSystem, u0map, tspan::Union{Tuple,
         end
     end
 
-    df = DiscreteFunction{true, true}(f; syms = Symbol.(unknowns(sys)),
-        indepsym = Symbol(get_iv(sys)),
-        paramsyms = Symbol.(ps), sys = sys,
-        observed = observedfun)
+    df = DiscreteFunction{true, true}(f; sys = sys, observed = observedfun)
     DiscreteProblem(df, u0, tspan, p; kwargs...)
 end
 
@@ -371,10 +368,9 @@ function DiscreteProblemExpr{iip}(sys::JumpSystem, u0map, tspan::Union{Tuple, No
         f = DiffEqBase.DISCRETE_INPLACE_DEFAULT
         u0 = $u0
         p = $p
+        sys = $sys
         tspan = $tspan
-        df = DiscreteFunction{true, true}(f; syms = $(Symbol.(unknowns(sys))),
-            indepsym = $(Symbol(get_iv(sys))),
-            paramsyms = $(Symbol.(parameters(sys))))
+        df = DiscreteFunction{true, true}(f; sys = sys)
         DiscreteProblem(df, u0, tspan, p)
     end
 end

src/systems/nonlinear/nonlinearsystem.jl

@@ -271,8 +271,6 @@ function SciMLBase.NonlinearFunction{iip}(sys::NonlinearSystem, dvs = unknowns(s
         jac_prototype = sparse ?
                         similar(calculate_jacobian(sys, sparse = sparse),
             Float64) : nothing,
-        syms = Symbol.(unknowns(sys)),
-        paramsyms = Symbol.(parameters(sys)),
         observed = observedfun)
 end
 
@@ -320,9 +318,7 @@ function NonlinearFunctionExpr{iip}(sys::NonlinearSystem, dvs = unknowns(sys),
         jac = $_jac
         NonlinearFunction{$iip}(f,
             jac = jac,
-            jac_prototype = $jp_expr,
-            syms = $(Symbol.(unknowns(sys))),
-            paramsyms = $(Symbol.(parameters(sys))))
+            jac_prototype = $jp_expr)
     end
     !linenumbers ? Base.remove_linenums!(ex) : ex
 end
@@ -347,7 +343,6 @@ function process_NonlinearProblem(constructor, sys::NonlinearSystem, u0map, para
 
     f = constructor(sys, dvs, ps, u0; jac = jac, checkbounds = checkbounds,
         linenumbers = linenumbers, parallel = parallel, simplify = simplify,
-        syms = Symbol.(dvs), paramsyms = Symbol.(ps),
         sparse = sparse, eval_expression = eval_expression, kwargs...)
     return f, u0, p
 end

src/systems/optimization/optimizationsystem.jl

@@ -362,8 +362,6 @@ function DiffEqBase.OptimizationProblem{iip}(sys::OptimizationSystem, u0map,
             grad = _grad,
             hess = _hess,
             hess_prototype = hess_prototype,
-            syms = Symbol.(unknowns(sys)),
-            paramsyms = Symbol.(parameters(sys)),
             cons = cons[2],
             cons_j = _cons_j,
             cons_h = _cons_h,
@@ -380,8 +378,6 @@ function DiffEqBase.OptimizationProblem{iip}(sys::OptimizationSystem, u0map,
             SciMLBase.NoAD();
             grad = _grad,
             hess = _hess,
-            syms = Symbol.(unknowns(sys)),
-            paramsyms = Symbol.(parameters(sys)),
             hess_prototype = hess_prototype,
             expr = obj_expr,
             observed = observedfun)
@@ -552,13 +548,9 @@ function OptimizationProblemExpr{iip}(sys::OptimizationSystem, u0map,
             ucons = $ucons
             cons_j = $_cons_j
             cons_h = $_cons_h
-            syms = $(Symbol.(unknowns(sys)))
-            paramsyms = $(Symbol.(parameters(sys)))
             _f = OptimizationFunction{iip}(f, SciMLBase.NoAD();
                 grad = grad,
                 hess = hess,
-                syms = syms,
-                paramsyms = paramsyms,
                 hess_prototype = hess_prototype,
                 cons = cons,
                 cons_j = cons_j,
@@ -580,13 +572,9 @@ function OptimizationProblemExpr{iip}(sys::OptimizationSystem, u0map,
             lb = $lb
             ub = $ub
             int = $int
-            syms = $(Symbol.(unknowns(sys)))
-            paramsyms = $(Symbol.(parameters(sys)))
             _f = OptimizationFunction{iip}(f, SciMLBase.NoAD();
                 grad = grad,
                 hess = hess,
-                syms = syms,
-                paramsyms = paramsyms,
                 hess_prototype = hess_prototype,
                 expr = obj_expr)
             OptimizationProblem{$iip}(_f, u0, p; lb = lb, ub = ub, int = int, kwargs...)

src/systems/unit_check.jl

@@ -1,6 +1,3 @@
-import DynamicQuantities, Unitful
-const DQ = DynamicQuantities
-
 #For dispatching get_unit
 const Conditional = Union{typeof(ifelse), typeof(IfElse.ifelse)}
 const Comparison = Union{typeof.([==, !=, ≠, <, <=, ≤, >, >=, ≥])...}

test/ccompile.jl

@@ -1,7 +1,8 @@
 using ModelingToolkit, Test
-@parameters t a
+using ModelingToolkit: t_nounits as t, D_nounits as D
+
+@parameters a
 @variables x y
-D = Differential(t)
 eqs = [D(x) ~ a * x - x * y,
     D(y) ~ -3y + x * y]
 f = build_function([x.rhs for x in eqs], [x, y], [a], t, expression = Val{false},

test/clock.jl

@@ -1,5 +1,6 @@
 using ModelingToolkit, Test, Setfield, OrdinaryDiffEq, DiffEqCallbacks
 using ModelingToolkit: Continuous
+using ModelingToolkit: t_nounits as t, D_nounits as D
 
 function infer_clocks(sys)
     ts = TearingState(sys)
@@ -9,9 +10,8 @@ end
 
 @info "Testing hybrid system"
 dt = 0.1
-@variables t x(t) y(t) u(t) yd(t) ud(t) r(t)
+@variables x(t) y(t) u(t) yd(t) ud(t) r(t)
 @parameters kp
-D = Differential(t)
 # u(n + 1) := f(u(n))
 
 eqs = [yd ~ Sample(t, dt)(y)
@@ -89,9 +89,8 @@ d = Clock(t, dt)
 
 @info "Testing shift normalization"
 dt = 0.1
-@variables t x(t) y(t) u(t) yd(t) ud(t) r(t) z(t)
+@variables x(t) y(t) u(t) yd(t) ud(t) r(t) z(t)
 @parameters kp
-D = Differential(t)
 d = Clock(t, dt)
 k = ShiftIndex(d)
 
@@ -161,9 +160,8 @@ sol2 = solve(prob, Tsit5())
 @info "Testing multi-rate hybrid system"
 dt = 0.1
 dt2 = 0.2
-@variables t x(t) y(t) u(t) r(t) yd1(t) ud1(t) yd2(t) ud2(t)
+@variables x(t) y(t) u(t) r(t) yd1(t) ud1(t) yd2(t) ud2(t)
 @parameters kp
-D = Differential(t)
 
 eqs = [
 # controller (time discrete part `dt=0.1`)
@@ -196,14 +194,12 @@ d2 = Clock(t, dt2)
 @info "test composed systems"
 
 dt = 0.5
-@variables t
 d = Clock(t, dt)
 k = ShiftIndex(d)
 timevec = 0:0.1:4
 
 function plant(; name)
     @variables x(t)=1 u(t)=0 y(t)=0
-    D = Differential(t)
     eqs = [D(x) ~ -x + u
         y ~ x]
     ODESystem(eqs, t; name = name)
@@ -253,9 +249,8 @@ ci, varmap = infer_clocks(cl)
 @info "Testing multi-rate hybrid system"
 dt = 0.1
 dt2 = 0.2
-@variables t x(t)=0 y(t)=0 u(t)=0 yd1(t)=0 ud1(t)=0 yd2(t)=0 ud2(t)=0
+@variables x(t)=0 y(t)=0 u(t)=0 yd1(t)=0 ud1(t)=0 yd2(t)=0 ud2(t)=0
 @parameters kp=1 r=1
-D = Differential(t)
 
 eqs = [
 # controller (time discrete part `dt=0.1`)
@@ -326,7 +321,6 @@ end
 using ModelingToolkitStandardLibrary.Blocks
 
 dt = 0.05
-@variables t
 d = Clock(t, dt)
 k = ShiftIndex(d)