From 547cce815369448acff3da42f2b5a8220f708bda Mon Sep 17 00:00:00 2001
From: Christopher Tessum <ctessum@gmail.com>
Date: Wed, 21 Aug 2024 19:55:25 +0900
Subject: [PATCH 1/3] Update to MTK v9

---
 Project.toml                     | 13 +++++------
 src/EarthSciMLBase.jl            |  2 +-
 src/advection.jl                 | 11 ++++-----
 src/coupled_system.jl            |  4 ++--
 src/domaininfo.jl                | 39 ++++---------------------------
 src/simulator.jl                 |  4 ++--
 src/simulator_strategies.jl      |  2 +-
 src/simulator_strategy_strang.jl |  2 +-
 src/simulator_utils.jl           | 14 +++++------
 test/advection_test.jl           |  7 ++----
 test/coord_trans_test.jl         |  9 +++----
 test/coupled_system_test.jl      | 26 +++++++++------------
 test/domaininfo_test.jl          | 32 ++++++++++++-------------
 test/operator_compose_test.jl    | 40 +++++++++++---------------------
 test/param_to_var_test.jl        | 17 +++++++-------
 test/simulator_test.jl           |  2 +-
 test/simulator_utils_test.jl     | 17 +++++++-------
 17 files changed, 93 insertions(+), 148 deletions(-)

diff --git a/Project.toml b/Project.toml
index 533044e1..169a06d3 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "EarthSciMLBase"
 uuid = "e53f1632-a13c-4728-9402-0c66d48804b0"
 authors = ["EarthSciML Authors and Contributors"]
-version = "0.15.3"
+version = "0.16.0"
 
 [deps]
 BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
@@ -9,6 +9,7 @@ Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 DomainSets = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
+DynamicQuantities = "06fc5a27-2a28-4c7c-a15d-362465fb6821"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MetaGraphsNext = "fa8bd995-216d-47f1-8a91-f3b68fbeb377"
@@ -18,25 +19,23 @@ ProgressLogging = "33c8b6b6-d38a-422a-b730-caa89a2f386c"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
-Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
 BlockBandedMatrices = "0.13"
-Catalyst = "13"
+Catalyst = "14"
 DiffEqCallbacks = "2"
 DifferentialEquations = "7"
 DocStringExtensions = "0.9"
-DomainSets = "0.5, 0.6"
+DomainSets = "0.7"
 Graphs = "1"
 MetaGraphsNext = "0.5, 0.6, 0.7"
-MethodOfLines = "0.10"
-ModelingToolkit = "8"
+MethodOfLines = "0.11"
+ModelingToolkit = "9"
 OrdinaryDiffEq = "6"
 ProgressLogging = "0.1"
 SciMLBase = "2"
 SciMLOperators = "0.3"
 Symbolics = "5"
-Unitful = "1"
 julia = "1.6"
 
 [extras]
diff --git a/src/EarthSciMLBase.jl b/src/EarthSciMLBase.jl
index c35da70e..307ec62e 100644
--- a/src/EarthSciMLBase.jl
+++ b/src/EarthSciMLBase.jl
@@ -2,7 +2,7 @@ module EarthSciMLBase
 using ModelingToolkit, Symbolics, Catalyst
 using Graphs, MetaGraphsNext
 using DocStringExtensions
-using Unitful
+using DynamicQuantities
 using OrdinaryDiffEq, DomainSets
 using SciMLBase: DECallback, CallbackSet
 using DiffEqCallbacks
diff --git a/src/advection.jl b/src/advection.jl
index 870f057d..72b8fa59 100644
--- a/src/advection.jl
+++ b/src/advection.jl
@@ -95,18 +95,15 @@ function ConstantWind(t, vals...; name=:ConstantWind)
     uvars = Num[]
     for (i, val) ∈ enumerate(vals)
         sym = Symbol("v_$i")
-        uv = (@variables $sym(t))[1]
-        uv = Symbolics.setmetadata(uv, ModelingToolkit.VariableUnit, unit(val))
-        uv = Symbolics.setmetadata(uv, ModelingToolkit.VariableDescription,
-            "Constant wind speed in the $(i)$(counts[i]) direction.")
+        desc = "Constant wind speed in the $(i)$(counts[i]) direction."
+        uv = only(@variables $sym(t), [unit=val/ustrip(val), description=desc])
         push!(uvars, uv)
     end
     uvals = []
     for i in eachindex(vals)
-        u = unit(vals[i])
         v = ustrip(vals[i])
         sym = Symbol("c_v$i")
-        c = (@constants $sym = v [unit = u])[1]
+        c = only(@constants $sym = v [unit = vals[i]/v])
         push!(uvals, c)
     end
     eqs = convert(Vector{Equation}, Symbolics.scalarize(uvars .~ uvals))
@@ -120,7 +117,7 @@ function couple2(mw::MeanWindCoupler, w::ConstantWindCoupler)
     @named a = ODESystem(Equation[], ModelingToolkit.get_iv(mw), [], [], systems=[mw])
     @named b = ODESystem(Equation[], ModelingToolkit.get_iv(w), [], [], systems=[w])
     ConnectorSystem(
-        Symbolics.scalarize(states(a) .~ states(b)),
+        Symbolics.scalarize(unknowns(a) .~ unknowns(b)),
         mw, w,
     )
 end
\ No newline at end of file
diff --git a/src/coupled_system.jl b/src/coupled_system.jl
index 4f8184af..89313f9b 100644
--- a/src/coupled_system.jl
+++ b/src/coupled_system.jl
@@ -93,14 +93,14 @@ function couple(systems...)::CoupledSystem
             end
         elseif sys isa DECallback
             push!(o.callbacks, sys)
+        elseif (sys isa Tuple) || (sys isa AbstractVector)
+            o = couple(o, sys...)
         elseif hasmethod(init_callback, (typeof(sys), Simulator))
             push!(o.init_callbacks, sys)
         elseif applicable(couple, o, sys)
             o = couple(o, sys)
         elseif applicable(couple, sys, o)
             o = couple(sys, o)
-        elseif (sys isa Tuple) || (sys isa AbstractVector)
-            o = couple(o, sys...)
         else
             error("Cannot couple a $(typeof(sys)).")
         end
diff --git a/src/domaininfo.jl b/src/domaininfo.jl
index 0b18ef89..ace5f2a0 100644
--- a/src/domaininfo.jl
+++ b/src/domaininfo.jl
@@ -310,22 +310,14 @@ domains(di::DomainInfo) = unique(vcat(domains.(di.icbc)...))
 
 function Base.:(+)(sys::ModelingToolkit.ODESystem, di::DomainInfo)::ModelingToolkit.PDESystem
     dimensions = dims(di)
-    allvars = states(sys)
-    statevars = states(structural_simplify(sys))
-    # TODO(CT): Update once the MTK get_defaults function can get defaults for composed system.
-    # defaults = ModelingToolkit.get_defaults(sys)
-    toreplace, replacements = replacement_params(parameters(sys), pvars(di))
-    defaults = get_defaults_all(sys)
-    ps = [k => v for (k,v) in defaults] # Add parameters and their default values
-    parameterstokeep = setdiff(parameters(sys), toreplace)
-    if !all([p ∈ keys(defaults) for p in parameterstokeep])
-        error("All parameters in the system of equations must have default values, but these ones don't: $(setdiff(parameterstokeep, keys(defaults))).")
-    end
-    ivs = dims(di) # New dimensions are the independent variables.
+    allvars = unknowns(sys)
+    statevars = unknowns(structural_simplify(sys))
+    ps = parameters(sys)
+    toreplace, replacements = replacement_params(ps, pvars(di))
     dvs = add_dims(allvars, dimensions) # Add new dimensions to dependent variables.
     eqs = substitute(equations(sys), Dict(zip(toreplace, replacements))) # Substitute local coordinate parameters for global ones.
     eqs = Vector{Equation}([add_dims(eq, allvars, dimensions) for eq in eqs]) # Add new dimensions to equations.
-    PDESystem(eqs, icbc(di, statevars), domains(di), ivs, dvs, ps, name=nameof(sys), defaults=defaults)
+    PDESystem(eqs, icbc(di, statevars), domains(di), dimensions, dvs, ps, name=nameof(sys))
 end
 
 Base.:(+)(di::DomainInfo, sys::ModelingToolkit.ODESystem)::ModelingToolkit.PDESystem = sys + di
@@ -336,27 +328,6 @@ end
 
 Base.:(+)(di::DomainInfo, sys::Catalyst.ReactionSystem)::ModelingToolkit.PDESystem = sys + di
 
-# TODO(CT): Delete once the MTK get_defaults function can get defaults for composed system.
-get_defaults_all(sys) = get_defaults_all(sys, "", 0)
-function get_defaults_all(sys, prefix, depth)
-    dmap = Dict()
-    if depth == 1
-        prefix = Symbol("$(nameof(sys))₊")
-    elseif depth > 1
-        prefix = Symbol("$(prefix)₊$(nameof(sys))₊")
-    end
-    for (p, d) ∈ ModelingToolkit.get_defaults(sys)
-        n = Symbol("$(prefix)$(p)")
-        pp = (@parameters $(n))[1]
-        pp = add_metadata(pp, p)
-        dmap[pp] = d
-    end
-    for child ∈ ModelingToolkit.get_systems(sys)
-        dmap = merge(dmap, get_defaults_all(child, prefix, depth+1))
-    end
-    dmap
-end
-
 # Match local parameters with the global parameters of the same name.
 function replacement_params(localcoords::AbstractVector, globalcoords::AbstractVector)
     gcstr = string.(globalcoords)
diff --git a/src/simulator.jl b/src/simulator.jl
index a05adc61..269e1673 100644
--- a/src/simulator.jl
+++ b/src/simulator.jl
@@ -42,7 +42,7 @@ struct Simulator{T,FT1,FT2,TG}
 
         mtk_sys, obs_eqs = prune_observed(mtk_sys) # Remove unused variables to speed up computation.
 
-        vars = states(mtk_sys)
+        vars = unknowns(mtk_sys)
         ps = parameters(mtk_sys)
 
         dflts = ModelingToolkit.get_defaults(mtk_sys)
@@ -100,5 +100,5 @@ function get_callbacks(s::Simulator)
     [s.sys.callbacks; extra_cb]
 end
 
-Base.size(s::Simulator) = (length(states(s.sys_mtk)), [length(g) for g ∈ s.grid]...)
+Base.size(s::Simulator) = (length(unknowns(s.sys_mtk)), [length(g) for g ∈ s.grid]...)
 Base.length(s::Simulator) = *(size(s)...)
\ No newline at end of file
diff --git a/src/simulator_strategies.jl b/src/simulator_strategies.jl
index 2ff1b46a..c2dfe968 100644
--- a/src/simulator_strategies.jl
+++ b/src/simulator_strategies.jl
@@ -59,7 +59,7 @@ function mtk_op(s::Simulator)
 end
 
 function mtk_func(s::Simulator)
-    b = repeat([length(states(s.sys_mtk))], length(s) ÷ size(s)[1])
+    b = repeat([length(unknowns(s.sys_mtk))], length(s) ÷ size(s)[1])
     j = BlockBandedMatrix{Float64}(undef, b, b, (0,0)) # Jacobian prototype
     ODEFunction(mtk_op(s); jac_prototype=j)
 end
diff --git a/src/simulator_strategy_strang.jl b/src/simulator_strategy_strang.jl
index 7fd47f6c..018506f2 100644
--- a/src/simulator_strategy_strang.jl
+++ b/src/simulator_strategy_strang.jl
@@ -95,7 +95,7 @@ end
 "A callback to periodically run the stiff solver."
 function stiff_callback(s::Simulator{T}, st::SimulatorStrang, IIchunks, integrators) where T
     function affect!(integrator)
-        u = reshape(integrator.u, length(states(s.sys_mtk)), [length(g) for g in s.grid]...)
+        u = reshape(integrator.u, length(unknowns(s.sys_mtk)), [length(g) for g in s.grid]...)
         ode_step!(s, st, u, IIchunks, integrators, T(integrator.t), T(st.timestep))
         @views integrator.u .= u[:]
     end
diff --git a/src/simulator_utils.jl b/src/simulator_utils.jl
index a67bbc6e..3706fd78 100644
--- a/src/simulator_utils.jl
+++ b/src/simulator_utils.jl
@@ -31,7 +31,7 @@ function observed_expression(eqs, x)
     for v ∈ Symbolics.get_variables(expr)
         v_expr = observed_expression(eqs, v)
         if !isnothing(v_expr)
-            expr = Symbolics.replace(expr, v => v_expr)
+            expr = Symbolics.substitute(expr, v => v_expr)
         end
     end
     # Do it again to catch extra variables TODO(CT): Theoretically this could recurse forever; when to stop?
@@ -145,14 +145,14 @@ function get_needed_vars(sys::ODESystem)
         ModelingToolkit.asgraph(sys), 
         ModelingToolkit.variable_dependencies(sys),
     )
-    g = SimpleDiGraph(length(states(sys)))
+    g = SimpleDiGraph(length(unknowns(sys)))
     for (i, es) in enumerate(varvardeps.badjlist)
         for e in es
             add_edge!(g, i, e)
         end
     end
-    allst = states(sys)
-    simpst = states(structural_simplify(sys))
+    allst = unknowns(sys)
+    simpst = unknowns(structural_simplify(sys))
     stidx = [only(findall(isequal(s), allst)) for s in simpst]
     collect(Graphs.DFSIterator(g, stidx))    
 end
@@ -166,8 +166,8 @@ remove computationally intensive equations that are not used in the final model.
 """
 function prune_observed(sys::ODESystem)
     needed_var_idxs = get_needed_vars(sys)
-    needed_vars = Symbolics.tosymbol.(states(sys)[needed_var_idxs]; escape=true)
-    sys = structural_simplify(sys)
+    needed_vars = Symbolics.tosymbol.(unknowns(sys)[needed_var_idxs]; escape=true)
+    sys = structural_simplify(sys, split=false)
     deleteindex = []
     for (i, eq) ∈ enumerate(observed(sys))
         lhsvars = Symbolics.tosymbol.(Symbolics.get_variables(eq.lhs); escape=true)
@@ -180,6 +180,6 @@ function prune_observed(sys::ODESystem)
     obs = observed(sys)
     deleteat!(obs, deleteindex)
     sys2 = structural_simplify(ODESystem([equations(sys); obs], 
-        ModelingToolkit.get_iv(sys), name=nameof(sys)))
+        ModelingToolkit.get_iv(sys), name=nameof(sys)), split=false)
     return sys2, observed(sys)
 end
\ No newline at end of file
diff --git a/test/advection_test.jl b/test/advection_test.jl
index d0a37f7a..20a3ea88 100644
--- a/test/advection_test.jl
+++ b/test/advection_test.jl
@@ -1,19 +1,18 @@
 using Main.EarthSciMLBase
 using Test
 using DomainSets, MethodOfLines, ModelingToolkit, DifferentialEquations
+using ModelingToolkit: t, D
 import SciMLBase
-using Unitful
+using DynamicQuantities
 using Dates, DomainSets
 
 @testset "Composed System" begin
-    @parameters t [unit = u"s"]
     @parameters x [unit = u"m"]
 
     struct ExampleSysCoupler2 sys end
     function ExampleSys()
         @variables y(t) [unit = u"kg"]
         @parameters p = 1.0 [unit = u"kg/s"]
-        D = Differential(t)
         ODESystem([D(y) ~ p], t; name=:examplesys,
             metadata=Dict(:coupletype=>ExampleSysCoupler2))
     end
@@ -22,7 +21,6 @@ using Dates, DomainSets
     function ExampleSysCopy()
         @variables y(t) [unit = u"kg"]
         @parameters p = 1.0 [unit = u"kg/s"]
-        D = Differential(t)
         ODESystem([D(y) ~ p], t; name=:examplesyscopy,
             metadata=Dict(:coupletype=>ExampleSysCopyCoupler2))
     end
@@ -60,7 +58,6 @@ using Dates, DomainSets
 end
 
 @testset "Coordinate transform" begin
-    @parameters t [unit = u"s"]
     @parameters lon [unit = u"rad"]
     @parameters lat [unit = u"rad"]
     @constants c_unit = 180 / π / 6 [unit = u"rad" description = "constant to make units cancel out"]
diff --git a/test/coord_trans_test.jl b/test/coord_trans_test.jl
index 163122f0..10487c42 100644
--- a/test/coord_trans_test.jl
+++ b/test/coord_trans_test.jl
@@ -1,9 +1,10 @@
 using EarthSciMLBase
 using ModelingToolkit, DomainSets
-using Unitful
+using ModelingToolkit: t, D
+using DynamicQuantities
 
 @testset "varindex" begin
-    @parameters lon lat x y t
+    @parameters lon lat x y
     @test EarthSciMLBase.varindex([lon, lat, x, y, t], :lat) == 2
 end
 
@@ -17,7 +18,6 @@ end
     @parameters lat [unit = u"rad"]
     @parameters x [unit = u"m"]
     @parameters y [unit = u"m"]
-    @parameters t [unit = u"s"]
     pd = partialderivatives_δxyδlonlat([lon, x, lat, y, t])
     @test isequal(pd, Dict(3 => 1.0 / EarthSciMLBase.lat2meters, 1 => 1.0 / (EarthSciMLBase.lon2m * cos(lat))))
 end
@@ -26,13 +26,11 @@ end
     @parameters lon [unit = u"rad"]
     @parameters lat [unit = u"rad"]
     @parameters lev [unit = u"m"]
-    @parameters t [unit = u"s"]
 
     function Example()
         @variables c(t) = 5.0 [unit = u"kg"]
         @constants t_c = 1.0 [unit = u"s"] # constant to make `sin` unitless
         @constants c_c = 1.0 [unit = u"kg/s"] # constant to make equation units work out
-        D = Differential(t)
         ODESystem([D(c) ~ sin(t / t_c) * c_c], t, name=:examplesys)
     end
     examplesys = Example()
@@ -65,7 +63,6 @@ end
     @parameters lon [unit = u"rad"]
     @parameters lat [unit = u"rad"]
     @parameters lev [unit = u"m"]
-    @parameters t [unit = u"s"]
 
     deg2rad(x) = x * π / 180.0
     domain = DomainInfo(
diff --git a/test/coupled_system_test.jl b/test/coupled_system_test.jl
index 8990a366..119479d1 100644
--- a/test/coupled_system_test.jl
+++ b/test/coupled_system_test.jl
@@ -1,10 +1,9 @@
 using Main.EarthSciMLBase
 using ModelingToolkit
+using ModelingToolkit: t, D
 using Test
 using Catalyst
-using Unitful
-
-@parameters t [unit = u"s"]
+using DynamicQuantities
 
 @testset "Composed System" begin
     struct SEqnCoupler
@@ -12,10 +11,9 @@ using Unitful
     end
     function SEqn()
         @variables S(t), I(t), R(t)
-        D = Differential(t)
         N = S + I + R
         @parameters β [unit = u"s^-1"]
-        @named seqn = ODESystem([D(S) ~ -β * S * I / N],
+        @named seqn = ODESystem([D(S) ~ -β * S * I / N], t,
             metadata=Dict(:coupletype => SEqnCoupler))
     end
 
@@ -24,11 +22,10 @@ using Unitful
     end
     function IEqn()
         @variables S(t), I(t), R(t)
-        D = Differential(t)
         N = S + I + R
         @parameters β [unit = u"s^-1"]
         @parameters γ [unit = u"s^-1"]
-        @named ieqn = ODESystem([D(I) ~ β * S * I / N - γ * I],
+        @named ieqn = ODESystem([D(I) ~ β * S * I / N - γ * I], t,
             metadata=Dict(:coupletype => IEqnCoupler))
     end
 
@@ -37,9 +34,8 @@ using Unitful
     end
     function REqn()
         @variables I(t), R(t)
-        D = Differential(t)
         @parameters γ [unit = u"s^-1"]
-        @named reqn = ODESystem([D(R) ~ γ * I],
+        @named reqn = ODESystem([D(R) ~ γ * I], t,
             metadata=Dict(:coupletype => REqnCoupler))
     end
 
@@ -71,9 +67,9 @@ using Unitful
     sir_simple = structural_simplify(sirfinal)
 
     want_eqs = [
-        Differential(t)(reqn.R) ~ reqn.γ * reqn.I,
-        Differential(t)(seqn.S) ~ (-seqn.β * seqn.I * seqn.S) / (seqn.I + seqn.R + seqn.S),
-        Differential(t)(ieqn.I) ~ (ieqn.β * ieqn.I * ieqn.S) / (ieqn.I + ieqn.R + ieqn.S) - ieqn.γ * ieqn.I,
+        D(reqn.R) ~ reqn.γ * reqn.I,
+        D(seqn.S) ~ (-seqn.β * seqn.I * seqn.S) / (seqn.I + seqn.R + seqn.S),
+        D(ieqn.I) ~ (ieqn.β * ieqn.I * ieqn.S) / (ieqn.I + ieqn.R + ieqn.S) - ieqn.γ * ieqn.I,
     ]
 
     have_eqs = equations(sir_simple)
@@ -120,9 +116,9 @@ end
         @parameters jNO2 = 0.0149 [unit = u"s^-1"]
         @species NO2(t) = 10.0
         rxs = [
-            Reaction(jNO2, [NO2], [], [1], [1])
+            Reaction(jNO2, [NO2], [], [1], [])
         ]
-        rs = ReactionSystem(rxs, t; combinatoric_ratelaws=false, name=:b)
+        rs = complete(ReactionSystem(rxs, t; combinatoric_ratelaws=false, name=:b))
         convert(ODESystem, rs; metadata=Dict(:coupletype=>BCoupler))
     end
 
@@ -165,7 +161,7 @@ end
             eqstr = string(equations(m))
             @test occursin("b₊c_NO2(t)", eqstr)
             @test occursin("b₊jNO2(t)", eqstr)
-            @test occursin("b₊NO2(t)", string(states(m)))
+            @test occursin("b₊NO2(t)", string(unknowns(m)))
             obstr = string(observed(m))
             @test occursin("a₊j_NO2(t) ~ a₊j_unit", obstr)
             @test occursin("c₊NO2(t) ~ c₊emis", obstr)
diff --git a/test/domaininfo_test.jl b/test/domaininfo_test.jl
index 31c52181..e3b50e92 100644
--- a/test/domaininfo_test.jl
+++ b/test/domaininfo_test.jl
@@ -1,21 +1,23 @@
-using EarthSciMLBase
+using Test
+using Main.EarthSciMLBase
 using ModelingToolkit, Catalyst
 using MethodOfLines, DifferentialEquations, DomainSets
+using ModelingToolkit: t_nounits; t=t_nounits
+using ModelingToolkit: D_nounits; D=D_nounits
 import SciMLBase
 
-@parameters x y t α = 10.0
+@parameters x y α = 10.0
 @variables u(t) v(t)
-Dt = Differential(t)
 
 x_min = y_min = t_min = 0.0
 x_max = y_max = 1.0
 t_max = 11.5
 
-eqs = [Dt(u) ~ -α * √abs(v),
-    Dt(v) ~ -α * √abs(u),
+eqs = [D(u) ~ -α * √abs(v),
+    D(v) ~ -α * √abs(u),
 ]
 
-@named sys = ODESystem(eqs)
+@named sys = ODESystem(eqs, t)
 
 indepdomain = t ∈ Interval(t_min, t_max)
 
@@ -46,16 +48,15 @@ end
 
 @testset "pde" begin
     pde_want = let
-        @parameters x y t α = 10.0
+        @parameters x y α = 10.0
         @variables u(..) v(..)
-        Dt = Differential(t)
 
         x_min = y_min = t_min = 0.0
         x_max = y_max = 1.0
         t_max = 11.5
 
-        eqs = [Dt(u(t, x, y)) ~ -α * √abs(v(t, x, y)),
-            Dt(v(t, x, y)) ~ -α * √abs(u(t, x, y)),
+        eqs = [D(u(t, x, y)) ~ -α * √abs(v(t, x, y)),
+            D(v(t, x, y)) ~ -α * √abs(u(t, x, y)),
         ]
 
         domains = [
@@ -78,7 +79,7 @@ end
             v(t, x, y_max) ~ 16.0,
         ]
 
-        @named pdesys = PDESystem(eqs, bcs, domains, [t, x, y], [u(t, x, y), v(t, x, y)], [α => 10.0])
+        @named pdesys = PDESystem(eqs, bcs, domains, [t, x, y], [u(t, x, y), v(t, x, y)], [α])
     end
 
     pde_result = sys + domain
@@ -93,12 +94,11 @@ end
 
 @testset "ReactionSystem" begin
     pde_want = let
-        @parameters x y t
+        @parameters x y
         @variables m₁(..) m₂(..)
-        Dt = Differential(t)
         eqs = [
-            Dt(m₁(t, x, y)) ~ -10.0 * m₁(t, x, y),
-            Dt(m₂(t, x, y)) ~ 10.0 * m₁(t, x, y),
+            D(m₁(t, x, y)) ~ -10.0 * m₁(t, x, y),
+            D(m₂(t, x, y)) ~ 10.0 * m₁(t, x, y),
         ]
 
         bcs = [
@@ -173,7 +173,7 @@ end
 end
 
 @testset "Simplify" begin
-    @parameters x, t
+    @parameters x
     domain = DomainInfo(
         constIC(16.0, t ∈ Interval(0, 1)),
         periodicBC(x ∈ Interval(0, 1)),
diff --git a/test/operator_compose_test.jl b/test/operator_compose_test.jl
index e3e54438..8fc559db 100644
--- a/test/operator_compose_test.jl
+++ b/test/operator_compose_test.jl
@@ -1,19 +1,17 @@
 using Main.EarthSciMLBase
 using ModelingToolkit
+using ModelingToolkit: t, D, t_nounits, D_nounits
 using Catalyst
 using Test
-using Unitful
-
-@parameters t
+using DynamicQuantities
 
 struct ExampleSysCoupler
     sys
 end
 function ExampleSys()
-    @variables x(t)
+    @variables x(t_nounits)
     @parameters p
-    D = Differential(t)
-    ODESystem([D(x) ~ p], t; name=:sys1,
+    ODESystem([D_nounits(x) ~ p], t; name=:sys1,
         metadata=Dict(:coupletype => ExampleSysCoupler))
 end
 
@@ -21,10 +19,9 @@ struct ExampleSysCopyCoupler
     sys
 end
 function ExampleSysCopy()
-    @variables x(t)
+    @variables x(t_nounits)
     @parameters p
-    D = Differential(t)
-    ODESystem([D(x) ~ p], t; name=:syscopy,
+    ODESystem([D_nounits(x) ~ p], t; name=:syscopy,
         metadata=Dict(:coupletype => ExampleSysCopyCoupler))
 end
 
@@ -32,10 +29,9 @@ struct ExampleSys2Coupler
     sys
 end
 function ExampleSys2(; name=:sys2)
-    @variables y(t)
+    @variables y(t_nounits)
     @parameters p
-    D = Differential(t)
-    ODESystem([D(y) ~ p], t; name=name,
+    ODESystem([D_nounits(y) ~ p], t; name=name,
         metadata=Dict(:coupletype => ExampleSys2Coupler))
 end
 
@@ -82,7 +78,7 @@ end
         sys
     end
     function ExampleSysNonODE()
-        @variables y(t)
+        @variables y(t_nounits)
         @parameters p
         ODESystem([y ~ p], t; name=:sysnonode,
             metadata=Dict(:coupletype => ExampleSysNonODECoupler))
@@ -124,14 +120,12 @@ end
 end
 
 @testset "Units" begin
-    @parameters t [unit = u"s"]
     struct U1Coupler
         sys
     end
     function U1()
         @variables x(t) [unit = u"kg"]
         @parameters p [unit = u"kg/s"]
-        D = Differential(t)
         ODESystem([D(x) ~ p], t; name=:sys1,
             metadata=Dict(:coupletype => U1Coupler))
     end
@@ -146,7 +140,6 @@ end
             metadata=Dict(:coupletype => U2Coupler))
     end
     
-    
     sys1 = U1()
     sys2 = U2()
 
@@ -166,14 +159,12 @@ end
 end
 
 @testset "Units Non-ODE" begin
-    @parameters t [unit = u"s"]
     struct U1Coupler
         sys
     end
     function U1()
         @variables x(t) [unit = u"kg"]
         @parameters p [unit = u"kg/s"]
-        D = Differential(t)
         ODESystem([D(x) ~ p], t; name=:sys1,
             metadata=Dict(:coupletype => U1Coupler))
     end
@@ -207,13 +198,11 @@ end
 end
 
 @testset "Units 2" begin
-    @parameters t [unit = u"s"]
     struct U1Coupler
         sys
     end
     function U1()
         @variables x(t) [unit = u"kg*m^-3"]
-        D = Differential(t)
         ODESystem([D(x) ~ 0], t; name=:sys1,
             metadata=Dict(:coupletype => U1Coupler))
     end
@@ -247,12 +236,12 @@ end
         sys
     end
     function Chem()
-        @species SO2(t) O2(t) SO4(t)
+        @species SO2(t_nounits) O2(t_nounits) SO4(t_nounits)
         @parameters α β
         rxns = [
             Reaction(α, [SO2, O2], [SO4], [1, 1], [1])
         ]
-        rs = ReactionSystem(rxns, t; name=:chem)
+        rs = complete(ReactionSystem(rxns, t_nounits; name=:chem))
         convert(ODESystem, rs; metadata=Dict(:coupletype => ChemCoupler))
     end
 
@@ -260,14 +249,13 @@ end
         sys
     end
     function Deposition()
-        @variables SO2(t)
+        @variables SO2(t_nounits)
         @parameters k = 2
-        D = Differential(t)
 
         eqs = [
-            D(SO2) ~ -k * SO2
+            D_nounits(SO2) ~ -k * SO2
         ]
-        ODESystem(eqs, t, [SO2], [k]; name=:deposition,
+        ODESystem(eqs, t_nounits, [SO2], [k]; name=:deposition,
             metadata=Dict(:coupletype => DepositionCoupler))
     end
 
diff --git a/test/param_to_var_test.jl b/test/param_to_var_test.jl
index 0af1b6ca..566d68e9 100644
--- a/test/param_to_var_test.jl
+++ b/test/param_to_var_test.jl
@@ -1,13 +1,12 @@
 using Test
-using Main.EarthSciMLBase, ModelingToolkit, Unitful, Symbolics
-
+using Main.EarthSciMLBase, ModelingToolkit, DynamicQuantities, Symbolics
+using ModelingToolkit: t, D
 
 @parameters α=1 [unit = u"kg", description="α description"] 
 @parameters β=2 [unit = u"kg*s", description="β description"]
-@variables t [unit=u"s", description="time"]
 @variables x(t) [unit=u"m", description="x description"]
-eq = Differential(t)(x) ~ α * x / β
-@named sys = ODESystem([eq]; metadata=:metatest)
+eq = D(x) ~ α * x / β
+@named sys = ODESystem([eq], t; metadata=:metatest)
 
 ii(x, y) = findfirst(isequal(x), y)
 isin(x, y) = ii(x, y) !== nothing
@@ -16,10 +15,10 @@ isin(x, y) = ii(x, y) !== nothing
 
 @testset "Single substitution" begin
     sys2 = param_to_var(sys, :β)
-    @test isin(β, states(sys2)) == true
+    @test isin(β, unknowns(sys2)) == true
     @test isin(β, parameters(sys2)) == false
     @test isin(β, Symbolics.get_variables(equations(sys2)[1])) == true
-    var = states(sys2)[ii(β, states(sys2))]
+    var = unknowns(sys2)[ii(β, unknowns(sys2))]
     @test Symbolics.getmetadata(var, ModelingToolkit.VariableUnit) == u"kg*s"
     @test Symbolics.getmetadata(var, ModelingToolkit.VariableDescription) == "β description"
     @test ModelingToolkit.get_metadata(sys2) == :metatest
@@ -29,10 +28,10 @@ end
 
 @testset "Multiple substitutions" begin
     sys3 = param_to_var(sys, :β, :α)
-    @test isin(β, states(sys3)) == true
+    @test isin(β, unknowns(sys3)) == true
     @test isin(β, parameters(sys3)) == false
     @test isin(β, Symbolics.get_variables(equations(sys3)[1])) == true
-    @test isin(α, states(sys3)) == true
+    @test isin(α, unknowns(sys3)) == true
     @test isin(α, parameters(sys3)) == false
     @test isin(α, Symbolics.get_variables(equations(sys3)[1])) == true
 end
\ No newline at end of file
diff --git a/test/simulator_test.jl b/test/simulator_test.jl
index 6e02eb64..a9f07f95 100644
--- a/test/simulator_test.jl
+++ b/test/simulator_test.jl
@@ -1,4 +1,4 @@
-using EarthSciMLBase
+using Main.EarthSciMLBase
 using Test
 using ModelingToolkit, DomainSets, OrdinaryDiffEq
 using SciMLOperators
diff --git a/test/simulator_utils_test.jl b/test/simulator_utils_test.jl
index 4199d6fe..a544c572 100644
--- a/test/simulator_utils_test.jl
+++ b/test/simulator_utils_test.jl
@@ -1,24 +1,25 @@
-using EarthSciMLBase: steplength, observed_expression, observed_function, utype, grid, timesteps, icbc
-using EarthSciMLBase
+using Main.EarthSciMLBase: steplength, observed_expression, observed_function, utype, grid, timesteps, icbc
+using Main.EarthSciMLBase
 using Test
 using ModelingToolkit, DomainSets
+using ModelingToolkit: t_nounits; t=t_nounits
+using ModelingToolkit: D_nounits; D=D_nounits
 
 @test steplength([0, 0.1, 0.2]) == 0.1
 
-@parameters x y lon = 0.0 lat = 0.0 lev = 1.0 t α = 10.0
+@parameters x y lon = 0.0 lat = 0.0 lev = 1.0 α = 10.0
 @constants p = 1.0
 @variables(
     u(t) = 1.0, v(t) = 1.0, x(t) = 1.0, y(t) = 1.0
 )
-Dt = Differential(t)
 
-eqs = [Dt(u) ~ -α * √abs(v) + lon,
-    Dt(v) ~ -α * √abs(u) + lat,
+eqs = [D(u) ~ -α * √abs(v) + lon,
+    D(v) ~ -α * √abs(u) + lat,
     x ~ 2α + p + y,
     y ~ 4α - 2p
 ]
 
-@named sys = ODESystem(eqs)
+@named sys = ODESystem(eqs, t)
 sys = structural_simplify(sys)
 
 xx = observed_expression(observed(sys), x)
@@ -55,7 +56,7 @@ domain = DomainInfo(
     partialderivatives_δxyδlonlat,
     constIC(16.0, indepdomain), constBC(16.0, partialdomains...))
 
-vars = states(sys)
+vars = unknowns(sys)
 
 bcs = icbc(domain, vars)
 

From b5ac27c1230d5f4076564b728774a4500e7dea3a Mon Sep 17 00:00:00 2001
From: Christopher Tessum <ctessum@gmail.com>
Date: Wed, 21 Aug 2024 20:44:42 +0900
Subject: [PATCH 2/3] Change `get_mtk` to `convert`

---
 src/coupled_system.jl         | 17 ++++++++---------
 src/simulator.jl              |  2 +-
 test/advection_test.jl        |  4 ++--
 test/coord_trans_test.jl      |  2 +-
 test/coupled_system_test.jl   |  8 ++++----
 test/domaininfo_test.jl       |  2 +-
 test/operator_compose_test.jl | 16 ++++++++--------
 7 files changed, 25 insertions(+), 26 deletions(-)

diff --git a/src/coupled_system.jl b/src/coupled_system.jl
index 89313f9b..435248c1 100644
--- a/src/coupled_system.jl
+++ b/src/coupled_system.jl
@@ -1,5 +1,4 @@
-export CoupledSystem, ConnectorSystem, get_mtk, couple
-
+export CoupledSystem, ConnectorSystem, couple
 
 """
 Types that implement an:
@@ -143,11 +142,11 @@ struct ACoupler sys end
 couple2() = nothing
 
 """
-    $(SIGNATURES)
+$(SIGNATURES)
 
-Get the ODE ModelingToolkit representation of a [`CoupledSystem`](@ref).
+Get the ODE ModelingToolkit ODESystem representation of a [`CoupledSystem`](@ref).
 """
-function get_mtk_ode(sys::CoupledSystem; name=:model)::ModelingToolkit.AbstractSystem
+function Base.convert(::Type{<:ODESystem}, sys::CoupledSystem; name=:model, kwargs...)::ModelingToolkit.AbstractSystem
     connector_eqs = []
     systems = copy(sys.systems)
     for (i, a) ∈ enumerate(systems)
@@ -167,7 +166,7 @@ function get_mtk_ode(sys::CoupledSystem; name=:model)::ModelingToolkit.AbstractS
         end
     end
     iv = ModelingToolkit.get_iv(first(systems))
-    connectors = ODESystem(connector_eqs, iv; name=name)
+    connectors = ODESystem(connector_eqs, iv; name=name, kwargs...)
 
     # Compose everything together.
     compose(connectors, systems...)
@@ -176,10 +175,10 @@ end
 """
     $(SIGNATURES)
 
-Get the ModelingToolkit representation of a [`CoupledSystem`](@ref).
+Get the ModelingToolkit PDESystem representation of a [`CoupledSystem`](@ref).
 """
-function get_mtk(sys::CoupledSystem; name=:model)::ModelingToolkit.AbstractSystem
-    o = get_mtk_ode(sys; name)
+function Base.convert(::Type{<:PDESystem}, sys::CoupledSystem; name=:model, kwargs...)::ModelingToolkit.AbstractSystem
+    o = convert(ODESystem, sys; name, kwargs...)
 
     if sys.domaininfo !== nothing
         o += sys.domaininfo
diff --git a/src/simulator.jl b/src/simulator.jl
index 269e1673..2c4f25f0 100644
--- a/src/simulator.jl
+++ b/src/simulator.jl
@@ -38,7 +38,7 @@ struct Simulator{T,FT1,FT2,TG}
 
     function Simulator(sys::CoupledSystem, Δs::AbstractVector{T2}) where {T2<:AbstractFloat}
         @assert !isnothing(sys.domaininfo) "The system must have a domain specified; see documentation for EarthSciMLBase.DomainInfo."
-        mtk_sys = get_mtk_ode(sys; name=:model)
+        mtk_sys = convert(ODESystem, sys; name=:model)
 
         mtk_sys, obs_eqs = prune_observed(mtk_sys) # Remove unused variables to speed up computation.
 
diff --git a/test/advection_test.jl b/test/advection_test.jl
index 20a3ea88..c76b9a98 100644
--- a/test/advection_test.jl
+++ b/test/advection_test.jl
@@ -36,7 +36,7 @@ using Dates, DomainSets
 
     combined = couple(sys1, sys2)
     combined_pde = couple(combined, domain, ConstantWind(t, 1.0u"m/s"), Advection())
-    combined_mtk = get_mtk(combined_pde)
+    combined_mtk = convert(PDESystem, combined_pde)
 
     @test length(equations(combined_mtk)) == 6
     @test length(combined_mtk.ivs) == 2
@@ -81,7 +81,7 @@ end
     )
 
     composed_sys = couple(examplesys, domain, Advection(), wind)
-    pde_sys = get_mtk(composed_sys)
+    pde_sys = convert(PDESystem, composed_sys)
 
     eqs = equations(pde_sys)
 
diff --git a/test/coord_trans_test.jl b/test/coord_trans_test.jl
index 10487c42..f30ff353 100644
--- a/test/coord_trans_test.jl
+++ b/test/coord_trans_test.jl
@@ -46,7 +46,7 @@ end
 
     composed_sys = couple(examplesys, domain, Advection())
 
-    sys_mtk = get_mtk(composed_sys)
+    sys_mtk = convert(PDESystem, composed_sys)
 
     have_eq = equations(sys_mtk)
     @assert length(have_eq) == 1
diff --git a/test/coupled_system_test.jl b/test/coupled_system_test.jl
index 119479d1..8e5a3a52 100644
--- a/test/coupled_system_test.jl
+++ b/test/coupled_system_test.jl
@@ -62,7 +62,7 @@ using DynamicQuantities
 
     sir = couple(seqn, ieqn, reqn)
 
-    sirfinal = get_mtk(sir)
+    sirfinal = convert(PDESystem, sir)
 
     sir_simple = structural_simplify(sirfinal)
 
@@ -156,7 +156,7 @@ end
     ]
     for (i, model) ∈ enumerate(models)
         @testset "permutation $i" begin
-            model_mtk = get_mtk(model)
+            model_mtk = convert(ODESystem, model)
             m = structural_simplify(model_mtk)
             eqstr = string(equations(m))
             @test occursin("b₊c_NO2(t)", eqstr)
@@ -172,9 +172,9 @@ end
 
     @testset "Stable evaluation" begin
         sys = couple(A(), B(), C())
-        eqs1 = string(equations(get_mtk(sys)))
+        eqs1 = string(equations(convert(ODESystem, sys)))
         @test occursin("b₊c_NO2(t)", eqs1)
-        eqs2 = string(equations(get_mtk(sys)))
+        eqs2 = string(equations(convert(ODESystem, sys)))
         @test occursin("b₊c_NO2(t)", eqs2)
         @test eqs1 == eqs2
     end
diff --git a/test/domaininfo_test.jl b/test/domaininfo_test.jl
index e3b50e92..03d9a4b9 100644
--- a/test/domaininfo_test.jl
+++ b/test/domaininfo_test.jl
@@ -189,7 +189,7 @@ end
     end
 
     sys_domain = couple(ExSys(), domain)
-    sys_mtk = get_mtk(sys_domain)
+    sys_mtk = convert(PDESystem, sys_domain)
 
     discretization = MOLFiniteDifference([x => 10], t, approx_order=2)
     prob = discretize(sys_mtk, discretization)
diff --git a/test/operator_compose_test.jl b/test/operator_compose_test.jl
index 8fc559db..17115bc9 100644
--- a/test/operator_compose_test.jl
+++ b/test/operator_compose_test.jl
@@ -46,7 +46,7 @@ end
 
     combined = couple(sys1, sys2)
 
-    ox = get_mtk(combined)
+    ox = convert(ODESystem, combined)
     op = structural_simplify(ox)
     eq = equations(op)
 
@@ -66,7 +66,7 @@ end
 
     combined = couple(sys1, sys2)
 
-    ox = get_mtk(combined)
+    ox = convert(ODESystem, combined)
     op = structural_simplify(ox)
     eq = equations(op)
     eqstr = replace(string(eq), "Symbolics." => "")
@@ -93,7 +93,7 @@ end
     end
 
     combined = couple(sys1, sys2)
-    combined_mtk = get_mtk(combined)
+    combined_mtk = convert(ODESystem, combined)
     sys_combined = structural_simplify(combined_mtk)
 
     streq = string(equations(sys_combined))
@@ -112,7 +112,7 @@ end
 
     combined = couple(sys1, sys2)
 
-    ox = get_mtk(combined)
+    ox = convert(ODESystem, combined)
     op = structural_simplify(ox)
     streq = string(equations(op))
     @test occursin("sys1₊p", streq)
@@ -151,7 +151,7 @@ end
 
     combined = couple(sys1, sys2)
 
-    ox = get_mtk(combined)
+    ox = convert(ODESystem, combined)
     op = structural_simplify(ox)
     streq = string(equations(op))
     @test occursin("sys1₊p", streq)
@@ -190,7 +190,7 @@ end
 
     combined = couple(sys1, sys2)
 
-    ox = get_mtk(combined)
+    ox = convert(ODESystem, combined)
     op = structural_simplify(ox)
     streq = string(equations(op))
     @test occursin("sys1₊p", streq)
@@ -227,7 +227,7 @@ end
 
     combined = couple(sys1, sys2)
 
-    sys = get_mtk(combined)
+    sys = convert(ODESystem, combined)
     @test occursin("sys1₊sys2_x(t)", string(equations(sys)))
 end
 
@@ -268,7 +268,7 @@ end
     end
 
     combined = couple(rn, dep)
-    cs = structural_simplify(get_mtk(combined))
+    cs = structural_simplify(convert(ODESystem, combined))
     eq = equations(cs)
 
     eqstr = replace(string(eq), "Symbolics." => "")

From a71359bc96a96ce64c654aa43215b70bf9887e43 Mon Sep 17 00:00:00 2001
From: Christopher Tessum <ctessum@gmail.com>
Date: Wed, 21 Aug 2024 22:34:34 +0900
Subject: [PATCH 3/3] Update docs

---
 Project.toml                     |  1 +
 docs/Project.toml                | 10 +++----
 docs/src/advection.md            | 16 ++++++------
 docs/src/composition.md          | 25 +++++++++---------
 docs/src/coord_transforms.md     |  6 ++---
 docs/src/example_all_together.md | 26 +++++++++---------
 docs/src/icbc.md                 | 18 +++++++------
 docs/src/operator_compose.md     | 45 ++++++++++++++++----------------
 docs/src/param_to_var.md         | 29 ++++++++++----------
 docs/src/simulator.md            | 10 ++++---
 src/advection.jl                 | 12 +++++++--
 11 files changed, 106 insertions(+), 92 deletions(-)

diff --git a/Project.toml b/Project.toml
index 169a06d3..9fda3e0d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -27,6 +27,7 @@ DiffEqCallbacks = "2"
 DifferentialEquations = "7"
 DocStringExtensions = "0.9"
 DomainSets = "0.7"
+DynamicQuantities = "0.13"
 Graphs = "1"
 MetaGraphsNext = "0.5, 0.6, 0.7"
 MethodOfLines = "0.11"
diff --git a/docs/Project.toml b/docs/Project.toml
index dab095ea..8fa00294 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -4,6 +4,7 @@ Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DomainSets = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
+DynamicQuantities = "06fc5a27-2a28-4c7c-a15d-362465fb6821"
 EarthSciMLBase = "e53f1632-a13c-4728-9402-0c66d48804b0"
 GraphMakie = "1ecd5474-83a3-4783-bb4f-06765db800d2"
 MetaGraphsNext = "fa8bd995-216d-47f1-8a91-f3b68fbeb377"
@@ -12,19 +13,18 @@ ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
-Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
 CairoMakie = "0.12"
-Catalyst = "13"
+Catalyst = "14"
 DifferentialEquations = "7"
 Documenter = "1"
 DomainSets = "0.6,0.7"
+DynamicQuantities = "0.13"
 GraphMakie = "0.5"
 MetaGraphsNext = "0.7"
-MethodOfLines = "0.8,0.9,0.10"
-ModelingToolkit = "8"
+MethodOfLines = "0.11"
+ModelingToolkit = "9"
 Plots = "1"
 SciMLOperators = "0.3"
 Symbolics = "5"
-Unitful = "1"
diff --git a/docs/src/advection.md b/docs/src/advection.md
index 6b93dba7..07041cff 100644
--- a/docs/src/advection.md
+++ b/docs/src/advection.md
@@ -9,18 +9,18 @@ Advection is implemented with the [`Advection`](@ref) type.
 To demonstrate how this can work, we will start with a simple system of equations:
 
 ```@example advection
-using EarthSciMLBase, ModelingToolkit, Unitful
+using EarthSciMLBase, ModelingToolkit
+using ModelingToolkit: t_nounits, D_nounits
+t = t_nounits
+D = D_nounits
 
-@parameters t
-
-function ExampleSys(t)
+function ExampleSys()
     @variables y(t)
     @parameters p=2.0
-    D = Differential(t)
     ODESystem([D(y) ~ p], t; name=:ExampleSys)
 end
 
-ExampleSys(t)
+ExampleSys()
 ```
 
 We also need to create our initial and boundary conditions.
@@ -35,8 +35,8 @@ Now we convert add advection to each of the state variables.
 We're also adding a constant wind ([`ConstantWind`](@ref)) in the x-direction, with a speed of 1.0.
 
 ```@example advection
-sys_advection = couple(ExampleSys(t), domain, ConstantWind(t, 1.0), Advection())
-sys_mtk = get_mtk(sys_advection)
+sys_advection = couple(ExampleSys(), domain, ConstantWind(t, 1.0), Advection())
+sys_mtk = convert(PDESystem, sys_advection)
 ```
 
 Finally, we can discretize the system and solve it:
diff --git a/docs/src/composition.md b/docs/src/composition.md
index 3bb1100f..be96c728 100644
--- a/docs/src/composition.md
+++ b/docs/src/composition.md
@@ -7,10 +7,11 @@ To demonstrate how this works, below we define three model components, `Photolys
 
 ```@example composition
 using ModelingToolkit, Catalyst, EarthSciMLBase
-@parameters t
+using ModelingToolkit: t_nounits
+t = t_nounits
 
 struct PhotolysisCoupler sys end
-function Photolysis(t; name=:Photolysis)
+function Photolysis(; name=:Photolysis)
     @variables j_NO2(t)
     eqs = [
         j_NO2 ~ max(sin(t/86400),0)
@@ -19,7 +20,7 @@ function Photolysis(t; name=:Photolysis)
         metadata=Dict(:coupletype=>PhotolysisCoupler))
 end
 
-Photolysis(t)
+Photolysis()
 ```
 
 You can see that the system defined above is mostly a standard ModelingToolkit ODE system,
@@ -44,18 +45,18 @@ Let's follow the same process for some additional components:
 
 ```@example composition
 struct ChemistryCoupler sys end
-function Chemistry(t; name=:Chemistry)
+function Chemistry(; name=:Chemistry)
     @parameters jNO2
     @species NO2(t)
     rxs = [
-        Reaction(jNO2, [NO2], [], [1], [1])
+        Reaction(jNO2, [NO2], [], [1], [])
     ]
     rsys = ReactionSystem(rxs, t, [NO2], [jNO2]; 
         combinatoric_ratelaws=false, name=name)
-    convert(ODESystem, rsys, metadata=Dict(:coupletype=>ChemistryCoupler))
+    convert(ODESystem, complete(rsys), metadata=Dict(:coupletype=>ChemistryCoupler))
 end
 
-Chemistry(t)
+Chemistry()
 ```
 
 For our chemistry component above, because it's is originally a ReactionSystem instead of an
@@ -63,7 +64,7 @@ ODESystem, we convert it to an ODESystem before adding the metadata.
 
 ```@example composition
 struct EmissionsCoupler sys end
-function Emissions(t; name=:Emissions)
+function Emissions(; name=:Emissions)
     @parameters emis = 1
     @variables NO2(t)
     eqs = [NO2 ~ emis]
@@ -71,7 +72,7 @@ function Emissions(t; name=:Emissions)
         metadata=Dict(:coupletype=>EmissionsCoupler))
 end
 
-Emissions(t)
+Emissions()
 ```
 
 Now, we need to define ways to couple the model components together.
@@ -104,11 +105,11 @@ end
 nothing # hide
 ```
 
-Finally, we can compose the model components together to create our complete model. To do so, we just initialize each model component and add the components together using the [`couple`](@ref) function. We can use the [`get_mtk`](@ref) function to convert the composed model to a [ModelingToolkit](https://mtk.sciml.ai/dev/) model so we can see what the combined equations look like.
+Finally, we can compose the model components together to create our complete model. To do so, we just initialize each model component and add the components together using the [`couple`](@ref) function. We can use the `convert` function to convert the composed model to a [ModelingToolkit](https://mtk.sciml.ai/dev/) `ODESystem` so we can see what the combined equations look like.
 
 ```@example composition
-model = couple(Photolysis(t), Chemistry(t), Emissions(t))
-get_mtk(model)
+model = couple(Photolysis(), Chemistry(), Emissions())
+convert(ODESystem, model)
 ```
 
 Finally, we can use the [`graph`](@ref) function to visualize the model components and their connections.
diff --git a/docs/src/coord_transforms.md b/docs/src/coord_transforms.md
index 3f2ef32f..877f8b56 100644
--- a/docs/src/coord_transforms.md
+++ b/docs/src/coord_transforms.md
@@ -8,12 +8,12 @@ To handle this, the [`DomainInfo`](@ref) type can be used to define coordinate s
 ```@example trans
 using EarthSciMLBase
 using ModelingToolkit
+using ModelingToolkit: t, D
 using DomainSets
-using Unitful
+using DynamicQuantities
 
 @parameters lon [unit = u"rad"]
 @parameters lat [unit = u"rad"]
-@parameters t [unit = u"s"]
 @parameters lev
 
 partialderivatives_δxyδlonlat([lev, lon, lat])
@@ -57,7 +57,7 @@ This returns a list of functions, one corresponding to each coordinate in our do
 Then we can calculate the symbolic partial derivative of a variable by just calling each function:
 
 ```@example trans
-@variables u
+@variables u(t)
 
 [δs[i](u) for i ∈ eachindex(δs)]
 ```
diff --git a/docs/src/example_all_together.md b/docs/src/example_all_together.md
index 17ad5a2f..a4c546bb 100644
--- a/docs/src/example_all_together.md
+++ b/docs/src/example_all_together.md
@@ -11,37 +11,39 @@ Our first example system is a simple reaction system:
 ```@example ex1
 using EarthSciMLBase
 using ModelingToolkit, Catalyst, DomainSets, MethodOfLines, DifferentialEquations
+using ModelingToolkit: t_nounits, D_nounits
+t = t_nounits
+D = D_nounits
 using Plots
 
 # Create our independent variable `t` and our partially-independent variables `x` and `y`.
-@parameters t x y
+@parameters x y
 
 struct ExampleSys1Coupler sys end
-function ExampleSys1(t)
+function ExampleSys1()
     @species c₁(t)=5.0 c₂(t)=5.0
     rs = ReactionSystem(
         [Reaction(2.0, [c₁], [c₂])],
         t; name=:Sys1, combinatoric_ratelaws=false)
-    convert(ODESystem, rs, metadata=Dict(:coupletype=>ExampleSys1Coupler))
+    convert(ODESystem, complete(rs), metadata=Dict(:coupletype=>ExampleSys1Coupler))
 end
 
-ExampleSys1(t)
+ExampleSys1()
 ```
 
 Our second example system is a simple ODE system, with the same two variables.
 
 ```@example ex1
 struct ExampleSys2Coupler sys end
-function ExampleSys2(t)
+function ExampleSys2()
     @variables c₁(t)=5.0 c₂(t)=5.0
     @parameters p₁=1.0 p₂=0.5
-    D = Differential(t)
     ODESystem(
         [D(c₁) ~ -p₁, D(c₂) ~ p₂],
         t; name=:Sys2, metadata=Dict(:coupletype=>ExampleSys2Coupler))
 end
 
-ExampleSys2(t)
+ExampleSys2()
 ```
 
 Now, we specify what should happen when we couple the two systems together.
@@ -62,18 +64,18 @@ nothing # hide
 Once we specify all of the above, it is simple to create our two individual systems and then couple them together. 
 
 ```@example ex1
-sys1 = ExampleSys1(t)
-sys2 = ExampleSys2(t)
+sys1 = ExampleSys1()
+sys2 = ExampleSys2()
 sys = couple(sys1, sys2)
 
-get_mtk(sys)
+convert(ODESystem, sys)
 ```
 
 At this point we have an ODE system that is composed of two other ODE systems.
 We can inspect its observed variables using the `observed` function.
 
 ```@example ex1
-simplified_sys = structural_simplify(get_mtk(sys))
+simplified_sys = structural_simplify(convert(ODESystem, sys))
 ```
 
 ```@example ex1
@@ -105,7 +107,7 @@ domain = DomainInfo(
 
 sys_pde = couple(sys, domain, ConstantWind(t, 1.0, 1.0), Advection())
 
-sys_pde_mtk = get_mtk(sys_pde)
+sys_pde_mtk = convert(PDESystem, sys_pde)
 ```
 
 Now we can inspect this new system that we've created:
diff --git a/docs/src/icbc.md b/docs/src/icbc.md
index ad688547..20e6b324 100644
--- a/docs/src/icbc.md
+++ b/docs/src/icbc.md
@@ -8,20 +8,22 @@ To demonstrate how to do this, we will use the following simple system of ordina
 ```@example icbc
 using EarthSciMLBase
 using ModelingToolkit
+using ModelingToolkit: t_nounits, D_nounits
+t = t_nounits
+D = D_nounits
 
-@parameters x y t
+@parameters x y
 
-function ExampleSys(t)
+function ExampleSys()
     @variables u(t) v(t)
-    Dt = Differential(t)
     eqs = [
-        Dt(u) ~ √abs(v),
-        Dt(v) ~ √abs(u),
+        D(u) ~ √abs(v),
+        D(v) ~ √abs(u),
     ]
     ODESystem(eqs, t; name=:Docs₊Example)
 end
 
-ExampleSys(t)
+ExampleSys()
 ```
 Next, we specify our initial and boundary conditions using the [`DomainInfo`](@ref) type.
 We initialize [`DomainInfo`](@ref) with sets of initial and boundary conditions.
@@ -49,9 +51,9 @@ It is also possible to use periodic boundary conditions with [`periodicBC`](@ref
 Finally, we combine our initial and boundary conditions with our system of equations using the [`couple`](@ref) function.
 
 ```@example icbc
-model = couple(ExampleSys(t), icbc)
+model = couple(ExampleSys(), icbc)
 
-eq_sys = get_mtk(model)
+eq_sys = convert(PDESystem, model)
 ```
 
 We can also look at the expanded boundary conditions of the resulting equation system:
diff --git a/docs/src/operator_compose.md b/docs/src/operator_compose.md
index 59961608..bffedd6c 100644
--- a/docs/src/operator_compose.md
+++ b/docs/src/operator_compose.md
@@ -25,37 +25,37 @@ The example below shows that when we `operator_compose` two systems together tha
 ```@example operator_compose
 using EarthSciMLBase
 using ModelingToolkit
+using ModelingToolkit: t_nounits, D_nounits
+t = t_nounits
+D = D_nounits
 
-@parameters t
 
 struct ExampleSysCoupler sys end
-function ExampleSys(t)
+function ExampleSys()
     @variables x(t)
     @parameters p
-    D = Differential(t)
     ODESystem([D(x) ~ p], t; name=:ExampleSys,
         metadata=Dict(:coupletype=>ExampleSysCoupler))
 end
 
-ExampleSys(t)
+ExampleSys()
 ```
 
 ```@example operator_compose
 struct ExampleSys2Coupler sys end
-function ExampleSys2(t)
+function ExampleSys2()
     @variables x(t)
     @parameters p
-    D = Differential(t)
     ODESystem([D(x) ~ 2p], t; name=:ExampleSys2,
         metadata=Dict(:coupletype=>ExampleSys2Coupler))
 end
 
-ExampleSys2(t)
+ExampleSys2()
 ```
 
 ```@example operator_compose
-sys1 = ExampleSys(t)
-sys2 = ExampleSys2(t)
+sys1 = ExampleSys()
+sys2 = ExampleSys2()
 
 function EarthSciMLBase.couple2(sys1::ExampleSysCoupler, sys2::ExampleSys2Coupler)
     sys1, sys2 = sys1.sys, sys2.sys
@@ -64,7 +64,7 @@ end
 
 combined = couple(sys1, sys2)
 
-combined_mtk = get_mtk(combined)
+combined_mtk = convert(ODESystem, combined)
 ```
 
 The simplified equation should be D(x) = p + sys2_xˍt:
@@ -84,16 +84,15 @@ This example demonstrates a case where one variable in the first system is equal
 
 ```@example operator_compose
 struct ExampleSys3Coupler sys end
-function ExampleSys3(t)
+function ExampleSys3()
     @variables y(t)
     @parameters p
-    D = Differential(t)
     ODESystem([D(y) ~ p], t; name=:ExampleSys3,
         metadata=Dict(:coupletype=>ExampleSys3Coupler))
 end
 
-sys1 = ExampleSys(t)
-sys2 = ExampleSys3(t)
+sys1 = ExampleSys()
+sys2 = ExampleSys3()
 
 function EarthSciMLBase.couple2(sys1::ExampleSysCoupler, sys2::ExampleSys3Coupler)
     sys1, sys2 = sys1.sys, sys2.sys
@@ -101,7 +100,7 @@ function EarthSciMLBase.couple2(sys1::ExampleSysCoupler, sys2::ExampleSys3Couple
 end
 
 combined = couple(sys1, sys2)
-combined_simplified = structural_simplify(get_mtk(combined))
+combined_simplified = structural_simplify(convert(ODESystem, combined))
 ```
 
 ```@example operator_compose
@@ -116,15 +115,15 @@ This could happen if we are extracting emissions from a file, and those emission
 
 ```@example operator_compose
 struct ExampleSysNonODECoupler sys end
-function ExampleSysNonODE(t)
+function ExampleSysNonODE()
     @variables y(t)
     @parameters p
     ODESystem([y ~ p], t; name=:ExampleSysNonODE,
         metadata=Dict(:coupletype=>ExampleSysNonODECoupler))
 end
 
-sys1 = ExampleSys(t)
-sys2 = ExampleSysNonODE(t)
+sys1 = ExampleSys()
+sys2 = ExampleSysNonODE()
 
 function EarthSciMLBase.couple2(sys1::ExampleSysCoupler, sys2::ExampleSysNonODECoupler)
     sys1, sys2 = sys1.sys, sys2.sys
@@ -132,7 +131,7 @@ function EarthSciMLBase.couple2(sys1::ExampleSysCoupler, sys2::ExampleSysNonODEC
 end
 
 combined = couple(sys1, sys2)
-sys_combined = structural_simplify(get_mtk(combined))
+sys_combined = structural_simplify(convert(ODESystem, combined))
 ```
 
 ```@example operator_compose
@@ -145,15 +144,15 @@ Finally, this last example shows the fourth case, where a conversion factor is i
 
 ```@example operator_compose
 struct ExampleSysNonODE2Coupler sys end
-function ExampleSysNonODE2(t)
+function ExampleSysNonODE2()
     @variables y(t)
     @parameters p
     ODESystem([y ~ p], t; name=:Docs₊ExampleSysNonODE2,
         metadata=Dict(:coupletype=>ExampleSysNonODE2Coupler))
 end
 
-sys1 = ExampleSys(t)
-sys2 = ExampleSysNonODE2(t)
+sys1 = ExampleSys()
+sys2 = ExampleSysNonODE2()
 
 function EarthSciMLBase.couple2(sys1::ExampleSysCoupler, sys2::ExampleSysNonODE2Coupler)
     sys1, sys2 = sys1.sys, sys2.sys
@@ -161,7 +160,7 @@ function EarthSciMLBase.couple2(sys1::ExampleSysCoupler, sys2::ExampleSysNonODE2
 end
 
 combined = couple(sys1, sys2)
-combined_simplified = structural_simplify(get_mtk(combined))
+combined_simplified = structural_simplify(convert(ODESystem, combined))
 ```
 
 ```@example operator_compose
diff --git a/docs/src/param_to_var.md b/docs/src/param_to_var.md
index 3d24b5ec..c1b92739 100644
--- a/docs/src/param_to_var.md
+++ b/docs/src/param_to_var.md
@@ -10,36 +10,35 @@ As an example, we will create a loss equation that depends on the temperature, s
 So first, let's specify the original system with constant temperature.
 
 ```@example param_to_var
-using ModelingToolkit, EarthSciMLBase, Unitful
-
-@variables t [unit=u"s", description="time"]
+using ModelingToolkit, EarthSciMLBase, DynamicQuantities
+using ModelingToolkit: t, D
 
 struct LossCoupler sys end
-function Loss(t)
+function Loss()
     @variables A(t)=1 [unit=u"kg"]
     @parameters k=1 [unit=u"s^-1"]
     @parameters T=300 [unit=u"K"]
     @constants T₀=300 [unit=u"K"]
-    eq = Differential(t)(A) ~ -k*exp(T/T₀) * A
-    ODESystem([eq]; name=:Loss, metadata=Dict(:coupletype=>LossCoupler))
+    eq = D(A) ~ -k*exp(T/T₀) * A
+    ODESystem([eq], t; name=:Loss, metadata=Dict(:coupletype=>LossCoupler))
 end
 
-Loss(t)
+Loss()
 ```
 
 Next, we specify the temperature that varies in time.
 
 ```@example param_to_var
 struct TemperatureCoupler sys end
-function Temperature(t)
+function Temperature()
     @variables T(t)=300 [unit=u"K"]
     @constants Tc=1.0 [unit=u"K/s"]
     @constants tc=1.0 [unit=u"s"]
-    eq = Differential(t)(T) ~ sin(t/tc)*Tc
-    ODESystem([eq]; name=:Temperature, metadata=Dict(:coupletype=>TemperatureCoupler))
+    eq = D(T) ~ sin(t/tc)*Tc
+    ODESystem([eq], t; name=:Temperature, metadata=Dict(:coupletype=>TemperatureCoupler))
 end
 
-Temperature(t)
+Temperature()
 ```
 
 Now, we specify how to compose the two systems using `param_to_var`.
@@ -54,11 +53,11 @@ end
 
 Finally, we create the system components and the composed system.
 ```@example param_to_var
-l = Loss(t)
-t = Temperature(t)
-variable_loss = couple(l, t)
+l = Loss()
+temp = Temperature()
+variable_loss = couple(l, temp)
 
-get_mtk(variable_loss)
+convert(ODESystem, variable_loss)
 ```
 
 If we wanted to, we could then run a simulation with the composed system.
\ No newline at end of file
diff --git a/docs/src/simulator.md b/docs/src/simulator.md
index dd2508c7..80d7bd49 100644
--- a/docs/src/simulator.md
+++ b/docs/src/simulator.md
@@ -15,16 +15,18 @@ As an example, let's first define a system of ODEs:
 using EarthSciMLBase
 using ModelingToolkit, DomainSets, DifferentialEquations
 using SciMLOperators, Plots
+using ModelingToolkit: t_nounits, D_nounits
+t = t_nounits
+D = D_nounits
 
-@parameters y lon = 0.0 lat = 0.0 lev = 1.0 t α = 10.0
+@parameters y lon = 0.0 lat = 0.0 lev = 1.0 α = 10.0
 @constants p = 1.0
 @variables(
     u(t) = 1.0, v(t) = 1.0, x(t) = 1.0, y(t) = 1.0, windspeed(t) = 1.0
 )
-Dt = Differential(t)
 
-eqs = [Dt(u) ~ -α * √abs(v) + lon,
-    Dt(v) ~ -α * √abs(u) + lat + 1e-14 * lev,
+eqs = [D(u) ~ -α * √abs(v) + lon,
+    D(v) ~ -α * √abs(u) + lat + 1e-14 * lev,
     windspeed ~ lat + lon + lev,
 ]
 sys = ODESystem(eqs, t; name=:Docs₊sys)
diff --git a/src/advection.jl b/src/advection.jl
index 72b8fa59..217dc3fe 100644
--- a/src/advection.jl
+++ b/src/advection.jl
@@ -96,14 +96,22 @@ function ConstantWind(t, vals...; name=:ConstantWind)
     for (i, val) ∈ enumerate(vals)
         sym = Symbol("v_$i")
         desc = "Constant wind speed in the $(i)$(counts[i]) direction."
-        uv = only(@variables $sym(t), [unit=val/ustrip(val), description=desc])
+        if val isa DynamicQuantities.AbstractQuantity
+            uv = only(@variables $sym(t), [unit=val/ustrip(val), description=desc])
+        else
+            uv = only(@variables $sym(t), [description=desc])
+        end
         push!(uvars, uv)
     end
     uvals = []
     for i in eachindex(vals)
         v = ustrip(vals[i])
         sym = Symbol("c_v$i")
-        c = only(@constants $sym = v [unit = vals[i]/v])
+        if vals[i] isa DynamicQuantities.AbstractQuantity
+            c = only(@constants $sym = v [unit = vals[i]/v])
+        else
+            c = only(@constants $sym = v)
+        end
         push!(uvals, c)
     end
     eqs = convert(Vector{Equation}, Symbolics.scalarize(uvars .~ uvals))