Updates

Author: YingboMa

docs/src/basics/AbstractSystem.md

@@ -22,7 +22,7 @@ same keyword arguments, which are:
 
 - `system`: This is used for specifying subsystems for hierarchical modeling with
   reusable components. For more information, see the [components page](@ref components)
-- Defaults: Keyword arguments like `default_u0` are used for specifying default
+- Defaults: Keyword arguments like `defaults` are used for specifying default
   values which are used. If a value is not given at the `SciMLProblem` construction
   time, its numerical value will be the default.
 
@@ -49,9 +49,7 @@ Optionally, a system could have:
 
 - `observed(sys)`: All observed equations of the system and its subsystems.
 - `get_observed(sys)`: Observed equations of the current-level system.
-- `get_default_u0(sys)`: A `Dict` that maps states into their default initial
-  condition.
-- `get_default_p(sys)`: A `Dict` that maps parameters into their default value.
+- `get_defaults(sys)`: A `Dict` that maps variables into their default values.
 - `independent_variable(sys)`: The independent variable of a system.
 - `get_noiseeqs(sys)`: Noise equations of the current-level system.
 
@@ -132,7 +130,7 @@ u0 = [
 ## Default Value Handling
 
 The `AbstractSystem` types allow for specifying default values, for example
-`default_p` inside of them. At problem construction time, these values are merged
+`defaults` inside of them. At problem construction time, these values are merged
 into the value maps, where for any repeats the value maps override the default.
 In addition, defaults of a higher level in the system override the defaults of
-a lower level in the system. 
+a lower level in the system.

docs/src/tutorials/acausal_components.md

@@ -19,7 +19,7 @@ using ModelingToolkit, Plots, DifferentialEquations
 # Basic electric components
 function Pin(;name)
     @variables v(t) i(t)
-    ODESystem(Equation[], t, [v, i], [], name=name, default_u0=[v=>1.0, i=>1.0])
+    ODESystem(Equation[], t, [v, i], [], name=name, defaults=[v=>1.0, i=>1.0])
 end
 
 function Ground(;name)
@@ -39,7 +39,7 @@ function Resistor(;name, R = 1.0)
            0 ~ p.i + n.i
            v ~ p.i * R
           ]
-    ODESystem(eqs, t, [v], [R], systems=[p, n], default_p=Dict(R => val), name=name)
+    ODESystem(eqs, t, [v], [R], systems=[p, n], defaults=Dict(R => val), name=name)
 end
 
 function Capacitor(; name, C = 1.0)
@@ -54,7 +54,7 @@ function Capacitor(; name, C = 1.0)
            0 ~ p.i + n.i
            D(v) ~ p.i / C
           ]
-    ODESystem(eqs, t, [v], [C], systems=[p, n], default_p=Dict(C => val), name=name)
+    ODESystem(eqs, t, [v], [C], systems=[p, n], defaults=Dict(C => val), name=name)
 end
 
 function ConstantVoltage(;name, V = 1.0)
@@ -66,7 +66,7 @@ function ConstantVoltage(;name, V = 1.0)
            V ~ p.v - n.v
            0 ~ p.i + n.i
           ]
-    ODESystem(eqs, t, [], [V], systems=[p, n], default_p=Dict(V => val), name=name)
+    ODESystem(eqs, t, [], [V], systems=[p, n], defaults=Dict(V => val), name=name)
 end
 
 R = 1.0
@@ -122,7 +122,7 @@ component to simply be the values there:
 ```julia
 function Pin(;name)
     @variables v(t) i(t)
-    ODESystem(Equation[], t, [v, i], [], name=name, default_u0=[v=>1.0, i=>1.0])
+    ODESystem(Equation[], t, [v, i], [], name=name, defaults=[v=>1.0, i=>1.0])
 end
 ```
 
@@ -176,11 +176,11 @@ function Resistor(;name, R = 1.0)
            0 ~ p.i + n.i
            v ~ p.i * R
           ]
-    ODESystem(eqs, t, [v], [R], systems=[p, n], default_p=Dict(R => val), name=name)
+    ODESystem(eqs, t, [v], [R], systems=[p, n], defaults=Dict(R => val), name=name)
 end
 ```
 
-Notice that we have created this system with a `default_p` for the resistor's
+Notice that we have created this system with a `defaults` for the resistor's
 resistance. By doing so, if the resistance of this resistor is not overridden
 by a higher level default or overridden at `ODEProblem` construction time, this
 will be the value of the resistance.
@@ -200,7 +200,7 @@ function Capacitor(; name, C = 1.0)
            0 ~ p.i + n.i
            D(v) ~ p.i / C
           ]
-    ODESystem(eqs, t, [v], [C], systems=[p, n], default_p=Dict(C => val), name=name)
+    ODESystem(eqs, t, [v], [C], systems=[p, n], defaults=Dict(C => val), name=name)
 end
 ```
 
@@ -219,7 +219,7 @@ function ConstantVoltage(;name, V = 1.0)
            V ~ p.v - n.v
            0 ~ p.i + n.i
           ]
-    ODESystem(eqs, t, [], [V], systems=[p, n], default_p=Dict(V => val), name=name)
+    ODESystem(eqs, t, [], [V], systems=[p, n], defaults=Dict(V => val), name=name)
 end
 ```
 

docs/src/tutorials/tearing_parallelism.md

@@ -41,7 +41,7 @@ const t = Sym{ModelingToolkit.Parameter{Real}}(:t)
 const D = Differential(t)
 function Pin(;name)
     @variables v(t) i(t)
-    ODESystem(Equation[], t, [v, i], [], name=name, default_u0=Dict([v=>1.0, i=>1.0]))
+    ODESystem(Equation[], t, [v, i], [], name=name, defaults=Dict([v=>1.0, i=>1.0]))
 end
 
 function Ground(;name)
@@ -59,12 +59,12 @@ function ConstantVoltage(;name, V = 1.0)
            V ~ p.v - n.v
            0 ~ p.i + n.i
           ]
-    ODESystem(eqs, t, [], [V], systems=[p, n], default_p=Dict(V => val), name=name)
+    ODESystem(eqs, t, [], [V], systems=[p, n], default=Dict(V => val), name=name)
 end
 
 function HeatPort(;name)
     @variables T(t) Q_flow(t)
-    return ODESystem(Equation[], t, [T, Q_flow], [], default_u0=Dict(T=>293.15, Q_flow=>0.0), name=name)
+    return ODESystem(Equation[], t, [T, Q_flow], [], default=Dict(T=>293.15, Q_flow=>0.0), name=name)
 end
 
 function HeatingResistor(;name, R=1.0, TAmbient=293.15, alpha=1.0)
@@ -83,8 +83,10 @@ function HeatingResistor(;name, R=1.0, TAmbient=293.15, alpha=1.0)
           ]
     ODESystem(
         eqs, t, [v, RTherm], [R, TAmbient, alpha], systems=[p, n, h],
-        default_p=Dict(R=>R_val, TAmbient=>TAmbient_val, alpha=>alpha_val),
-        default_u0=Dict(v=>0.0, RTherm=>R_val),
+        defaults=Dict(
+            R=>R_val, TAmbient=>TAmbient_val, alpha=>alpha_val,
+            v=>0.0, RTherm=>R_val
+        ),
         name=name,
     )
 end
@@ -99,7 +101,7 @@ function HeatCapacitor(;name, rho=8050, V=1, cp=460, TAmbient=293.15)
           ]
     ODESystem(
         eqs, t, [], [rho, V, cp], systems=[h],
-        default_p=Dict(rho=>rho_val, V=>V_val, cp=>cp_val),
+        default=Dict(rho=>rho_val, V=>V_val, cp=>cp_val),
         name=name,
     )
 end
@@ -117,8 +119,7 @@ function Capacitor(;name, C = 1.0)
           ]
     ODESystem(
         eqs, t, [v], [C], systems=[p, n],
-        default_u0=Dict(v => 0.0),
-        default_p=Dict(C => val),
+        defaults=Dict(v => 0.0, C => val),
         name=name
     )
 end
@@ -159,7 +160,7 @@ end
 eqs = [
        D(E) ~ sum(((i, sys),)->getproperty(sys, Symbol(:resistor, i)).h.Q_flow, enumerate(rc_systems))
       ]
-big_rc = ODESystem(eqs, t, [], [], systems=rc_systems, default_u0=Dict(E=>0.0))
+big_rc = ODESystem(eqs, t, [], [], systems=rc_systems, default=Dict(E=>0.0))
 ```
 
 Now let's say we want to expose a bit more parallelism via running tearing.

src/structural_transformation/StructuralTransformations.jl

@@ -16,7 +16,7 @@ using ModelingToolkit: ODESystem, var_from_nested_derivative, Differential,
                        states, equations, vars, Symbolic, diff2term, value,
                        operation, arguments, Sym, Term, simplify, solve_for,
                        isdiffeq, isdifferential,
-                       get_structure, get_reduced_states, default_u0, default_p
+                       get_structure, get_reduced_states, defaults
 
 using ModelingToolkit.BipartiteGraphs
 using ModelingToolkit.SystemStructures

src/structural_transformation/codegen.jl

@@ -131,7 +131,7 @@ function gen_nlsolve(sys, eqs, vars)
     allvars = unique(collect(Iterators.flatten(map(ModelingToolkit.vars, rhss))))
     params = setdiff(allvars, vars)
 
-    u0map = default_u0(sys)
+    u0map = defaults(sys)
     # splatting to tighten the type
     u0 = [map(var->get(u0map, var, 1e-3), vars)...]
     # specialize on the scalar case
@@ -338,22 +338,11 @@ function ODAEProblem{iip}(
     s = structure(sys)
     @unpack fullvars = s
     dvs = fullvars[diffvars_range(s)]
-    defaults = merge(default_p(sys), default_u0(sys))
-    u0map′ = ModelingToolkit.lower_mapnames(u0map, independent_variable(sys))
-    u0 = ModelingToolkit.varmap_to_vars(u0map′, dvs; defaults=defaults)
-
     ps = parameters(sys)
-    if parammap isa DiffEqBase.NullParameters && isempty(default_p(sys))
-        isempty(ps) || throw(ArgumentError("The model has non-empty parameters but no parameters are specified in the problem."))
-        p = parammap
-    else
-        if parammap isa DiffEqBase.NullParameters
-            pp = Pair[]
-        else
-            pp = ModelingToolkit.lower_mapnames(parammap)
-        end
-        p = ModelingToolkit.varmap_to_vars(pp, ps; defaults=defaults)
-    end
+    defs = defaults(sys)
+
+    u0 = ModelingToolkit.varmap_to_vars(u0map, dvs; defaults=defs)
+    p = ModelingToolkit.varmap_to_vars(parammap, ps; defaults=defs)
 
     ODEProblem{iip}(build_torn_function(sys; kw...), u0, tspan, p; kw...)
 end

src/utils.jl

@@ -88,21 +88,3 @@ Maps the variable to a state.
 """
 tovar(s::Sym{<:Parameter}) = Sym{symtype(s)}(s.name)
 tovar(s::Sym) = s
-
-function lower_mapnames(umap::AbstractArray{T}) where {T<:Pair}
-    T[value(k) => value(v) for (k, v) in umap]
-end
-function lower_mapnames(umap::AbstractArray{T},name) where {T<:Pair}
-    T[lower_varname(value(k), name) => value(v) for (k, v) in umap]
-end
-function lower_mapnames(umap::NTuple{N,T}) where {N,T<:Pair}
-    ntuple(i->value(umap[i][1]) => value(umap[i][2]),N)
-end
-function lower_mapnames(umap::NTuple{N,T},name) where {N,T<:Pair}
-    ntuple(i->lower_varname(value(umap[i][1]), name) => value(umap[i][2]),N)
-end
-
-lower_mapnames(umap::AbstractArray) = umap # Ambiguity
-lower_mapnames(umap::AbstractArray,name) = umap
-lower_mapnames(umap::Tuple) = umap
-lower_mapnames(umap::Tuple, name) = umap

test/components.jl

@@ -5,7 +5,7 @@ using ModelingToolkit, OrdinaryDiffEq
 const t = Sym{ModelingToolkit.Parameter{Real}}(:t)
 function Pin(;name)
     @variables v(t) i(t)
-    ODESystem(Equation[], t, [v, i], [], name=name, default_u0=[v=>1.0, i=>1.0])
+    ODESystem(Equation[], t, [v, i], [], name=name, default=[v=>1.0, i=>1.0])
 end
 
 function Ground(name)
@@ -23,7 +23,7 @@ function ConstantVoltage(name; V = 1.0)
            V ~ p.v - n.v
            0 ~ p.i + n.i
           ]
-    ODESystem(eqs, t, [], [V], systems=[p, n], default_p=Dict(V => val), name=name)
+    ODESystem(eqs, t, [], [V], systems=[p, n], default=Dict(V => val), name=name)
 end
 
 function Resistor(name; R = 1.0)
@@ -37,7 +37,7 @@ function Resistor(name; R = 1.0)
            0 ~ p.i + n.i
            v ~ p.i * R
           ]
-    ODESystem(eqs, t, [v], [R], systems=[p, n], default_p=Dict(R => val), name=name)
+    ODESystem(eqs, t, [v], [R], systems=[p, n], default=Dict(R => val), name=name)
 end
 
 function Capacitor(name; C = 1.0)
@@ -52,7 +52,7 @@ function Capacitor(name; C = 1.0)
            0 ~ p.i + n.i
            D(v) ~ p.i / C
           ]
-    ODESystem(eqs, t, [v], [C], systems=[p, n], default_p=Dict(C => val), name=name)
+    ODESystem(eqs, t, [v], [C], systems=[p, n], default=Dict(C => val), name=name)
 end
 
 R = 1.0
@@ -82,7 +82,7 @@ rc_eqs = [
 
 rc_model = ODESystem(rc_eqs, t, systems=[resistor, capacitor, source, ground], name=:rc)
 sys = structural_simplify(rc_model)
-@test ModelingToolkit.default_p(sys) == Dict(
+@test ModelingToolkit.default(sys) == Dict(
                                              capacitor.C => 1.0,
                                              source.V => 1.0,
                                              resistor.R => 1.0,

test/odesystem.jl

@@ -205,7 +205,7 @@ D = Differential(t)
 eqs = [D(y₁) ~ -k₁*y₁+k₃*y₂*y₃,
        0     ~  y₁ + y₂ + y₃ - 1,
        D(y₂) ~  k₁*y₁-k₂*y₂^2-k₃*y₂*y₃]
-sys = ODESystem(eqs, default_p=[k₁ => 100, k₂ => 3e7], default_u0=[y₁ => 1.0])
+sys = ODESystem(eqs, defaults=[k₁ => 100, k₂ => 3e7, y₁ => 1.0])
 u0 = Pair[]
 push!(u0, y₂ => 0.0)
 push!(u0, y₃ => 0.0)
@@ -311,8 +311,8 @@ eqs = [D(D(x)) ~ -b/M*D(x) - k/M*x]
 ps = [M, b, k]
 default_u0 = [D(x) => 0.0, x => 10.0]
 default_p = [M => 1.0, b => 1.0, k => 1.0]
-@named sys = ODESystem(eqs, t, [x], ps, default_u0=default_u0, default_p=default_p)
+@named sys = ODESystem(eqs, t, [x], ps, defaults=[default_u0; default_p])
 sys = ode_order_lowering(sys)
-prob = ODEProblem(sys, Pair[], tspan)
+prob = ODEProblem(sys, [], tspan)
 sol = solve(prob, Tsit5())
 @test sum(abs, sol[end]) < 1

test/structural_transformation/tearing.jl

@@ -156,7 +156,7 @@ du = [0.0]; u = [1.0]; pr = 0.2; tt = 0.1
 @test du ≈ [-asin(u[1] - pr * tt)] atol=1e-5
 
 # test the initial guess is respected
-infprob = ODAEProblem(tearing(ODESystem(eqs, t, default_u0=Dict(z=>Inf))), [x=>1.0], (0, 1.0), [p=>0.2])
+infprob = ODAEProblem(tearing(ODESystem(eqs, t, defaults=Dict(z=>Inf))), [x=>1.0], (0, 1.0), [p=>0.2])
 @test_throws DomainError infprob.f(du, u, pr, tt)
 
 sol1 = solve(prob, Tsit5())

test/symbolic_parameters.jl

@@ -19,10 +19,9 @@ u0 = Pair{Num, Any}[
     y => σ, # default u0 from default p
     z => u-0.1,
 ]
-ns = NonlinearSystem(eqs, [x,y,z],[σ,ρ,β], name=:ns, default_p=par, default_u0=u0)
+ns = NonlinearSystem(eqs, [x,y,z],[σ,ρ,β], name=:ns, defaults=merge(par, u0))
 ns.y = u*1.1
-ModelingToolkit.default_p(ns)
-resolved = ModelingToolkit.varmap_to_vars(Dict(), parameters(ns), defaults=ModelingToolkit.default_p(ns))
+resolved = ModelingToolkit.varmap_to_vars(Dict(), parameters(ns), defaults=ModelingToolkit.defaults(ns))
 @test resolved == [1, 0.1+1, (0.1+1)*1.1]
 
 prob = NonlinearProblem(ns, [u=>1.0], Pair[])
@@ -35,7 +34,7 @@ top = NonlinearSystem([0 ~ -a + ns.x+b], [a], [b], systems=[ns], name=:top)
 top.b = ns.σ*0.5
 top.ns.x = u*0.5
 
-res = ModelingToolkit.varmap_to_vars(Dict(), parameters(top), defaults=ModelingToolkit.default_p(top))
+res = ModelingToolkit.varmap_to_vars(Dict(), parameters(top), defaults=ModelingToolkit.defaults(top))
 @test res == [0.5, 1, 0.1+1, (0.1+1)*1.1]
 
 prob = NonlinearProblem(top, [states(ns, u)=>1.0, a=>1.0], Pair[])
@@ -45,4 +44,4 @@ prob = NonlinearProblem(top, [states(ns, u)=>1.0, a=>1.0], Pair[])
 # test NullParameters+defaults
 prob = NonlinearProblem(top, [states(ns, u)=>1.0, a=>1.0])
 @test prob.u0 == [1.0, 0.5, 1.1, 0.9]
-@show sol = solve(prob,NewtonRaphson())
+@show sol = solve(prob,NewtonRaphson())