Merge branch 'master' into kf/state_selection

Author: YingboMa

NEWS.md

@@ -0,0 +1,10 @@
+# ModelingToolkit v8 Release Notes
+
+
+### Upgrade guide
+
+- `connect` should not be overloaded by users anymore. `[connect = Flow]`
+  informs ModelingToolkit that particular variable in a connector ought to sum
+  to zero, and by default, variables are equal in a connection. Please check out
+  [acausal components tutorial](https://mtk.sciml.ai/dev/tutorials/acausal_components/)
+  for examples.

Project.toml

@@ -1,7 +1,7 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
 authors = ["Chris Rackauckas <[email protected]>"]
-version = "8.0.0"
+version = "8.3.1"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -45,10 +45,10 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
 AbstractTrees = "0.3"
-ArrayInterface = "3.1.39"
+ArrayInterface = "3.1.39, 4"
 ConstructionBase = "1"
 DataStructures = "0.17, 0.18"
-DiffEqBase = "6.54.0"
+DiffEqBase = "6.81.0"
 DiffEqCallbacks = "2.16"
 DiffEqJump = "7.0, 8"
 DiffRules = "0.1, 1.0"
@@ -57,7 +57,7 @@ DocStringExtensions = "0.7, 0.8"
 DomainSets = "0.5"
 Graphs = "1.4"
 IfElse = "0.1"
-JuliaFormatter = "0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.20"
+JuliaFormatter = "0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.20, 0.21"
 LabelledArrays = "1.3"
 Latexify = "0.11, 0.12, 0.13, 0.14, 0.15"
 MacroTools = "0.5"

docs/src/basics/Composition.md

@@ -167,13 +167,17 @@ before numerically solving. The `structural_simplify` function removes
 these trivial equality relationships and trivial singularity equations,
 i.e. equations which result in `0~0` expressions, in over-specified systems.
 
-## Inheritance and Combine (TODO)
+## Inheritance and Combine
 
-Model inheritance can be done in two ways: explicitly or implicitly.
-The explicit way is to shadow variables with equality expressions.
-For example, let's assume we have three separate systems which we
-want to compose to a single one. This is how one could explicitly
-forward all states and parameters to the higher level system:
+Model inheritance can be done in two ways: implicitly or explicitly. First, one
+can use the `extend` function to extend a base model with another set of
+equations, states, and parameters. An example can be found in the
+[acausal components tutorial](@ref acausal).
+
+The explicit way is to shadow variables with equality expressions. For example,
+let's assume we have three separate systems which we want to compose to a single
+one. This is how one could explicitly forward all states and parameters to the
+higher level system:
 
 ```julia
 using ModelingToolkit, OrdinaryDiffEq, Plots
@@ -259,7 +263,7 @@ ODESystem(eqs, ...; continuous_events::Pair{Vector{Equation}, Vector{Equation}})
 ```
 where equations can be added that evaluate to 0 at discontinuities.
 
-To model events that have an affect on the state, provide `events::Pair{Vector{Equation}, Vector{Equation}}` where the first entry in the pair is a vector of equations describing event conditions, and the second vector of equations describe the affect on the state. The affect equations must be on the form
+To model events that have an effect on the state, provide `events::Pair{Vector{Equation}, Vector{Equation}}` where the first entry in the pair is a vector of equations describing event conditions, and the second vector of equations describe the effect on the state. The effect equations must be of the form
 ```
 single_state_variable ~ expression_involving_any_variables
 ```
@@ -275,7 +279,7 @@ function UnitMassWithFriction(k; name)
     D(x) ~ v
     D(v) ~ sin(t) - k*sign(v) # f = ma, sinusoidal force acting on the mass, and Coulomb friction opposing the movement
   ]
-  ODESystem(eqs, t, continuous_events=[v ~ 0], name=name) # when v = 0 there is a discontinuity
+  ODESystem(eqs, t; continuous_events=[v ~ 0], name) # when v = 0 there is a discontinuity
 end
 @named m = UnitMassWithFriction(0.7)
 prob = ODEProblem(m, Pair[], (0, 10pi))
@@ -296,7 +300,7 @@ affect   = [v ~ -v] # the effect is that the velocity changes sign
 @named ball = ODESystem([
     D(x) ~ v
     D(v) ~ -9.8
-], t, continuous_events = root_eqs => affect) # equation => affect
+], t; continuous_events = root_eqs => affect) # equation => affect
 
 ball = structural_simplify(ball)
 
@@ -323,7 +327,7 @@ continuous_events = [ # This time we have a vector of pairs
     D(y)  ~ vy,
     D(vx) ~ -9.8-0.1vx, # gravity + some small air resistance
     D(vy) ~ -0.1vy,
-], t, continuous_events = continuous_events)
+], t; continuous_events)
 
 
 ball = structural_simplify(ball)
@@ -340,4 +344,4 @@ tv = sort([LinRange(0, 10, 200); sol.t])
 plot(sol(tv)[y], sol(tv)[x], line_z=tv)
 vline!([-1.5, 1.5], l=(:black, 5), primary=false)
 hline!([0], l=(:black, 5), primary=false)
-```
+```

docs/src/tutorials/acausal_components.md

@@ -16,22 +16,10 @@ using ModelingToolkit, Plots, DifferentialEquations
 
 @variables t
 @connector function Pin(;name)
-    sts = @variables v(t)=1.0 i(t)=1.0
+    sts = @variables v(t)=1.0 i(t)=1.0 [connect = Flow]
     ODESystem(Equation[], t, sts, []; name=name)
 end
 
-function ModelingToolkit.connect(::Type{Pin}, ps...)
-    eqs = [
-           0 ~ sum(p->p.i, ps) # KCL
-          ]
-    # KVL
-    for i in 1:length(ps)-1
-        push!(eqs, ps[i].v ~ ps[i+1].v)
-    end
-
-    return eqs
-end
-
 function Ground(;name)
     @named g = Pin()
     eqs = [g.v ~ 0]
@@ -92,7 +80,8 @@ V = 1.0
 rc_eqs = [
           connect(source.p, resistor.p)
           connect(resistor.n, capacitor.p)
-          connect(capacitor.n, source.n, ground.g)
+          connect(capacitor.n, source.n)
+          connect(capacitor.n, ground.g)
          ]
 
 @named _rc_model = ODESystem(rc_eqs, t)
@@ -117,11 +106,15 @@ For each of our components we use a Julia function which emits an `ODESystem`.
 At the top we start with defining the fundamental qualities of an electrical
 circuit component. At every input and output pin a circuit component has
 two values: the current at the pin and the voltage. Thus we define the `Pin`
-component (connector) to simply be the values there:
+component (connector) to simply be the values there. Whenever two `Pin`s in a
+circuit are connected together, the system satisfies [Kirchoff's laws](https: //en.wikipedia.org/wiki/Kirchhoff%27s_circuit_laws),
+i.e. that currents sum to zero and voltages across the pins are equal.
+`[connect = Flow]` informs MTK that currents ought to sum to zero, and by
+default, variables are equal in a connection.
 
 ```julia
 @connector function Pin(;name)
-    sts = @variables v(t)=1.0 i(t)=1.0
+    sts = @variables v(t)=1.0 i(t)=1.0 [connect = Flow]
     ODESystem(Equation[], t, sts, []; name=name)
 end
 ```
@@ -253,26 +246,6 @@ V = 1.0
 @named ground = Ground()
 ```
 
-Next we have to define how we connect the circuit. Whenever two `Pin`s in a
-circuit are connected together, the system satisfies
-[Kirchoff's laws](https://en.wikipedia.org/wiki/Kirchhoff%27s_circuit_laws),
-i.e. that currents sum to zero and voltages across the pins are equal. Thus
-we will build a helper function `connect_pins` which implements these rules:
-
-```julia
-function ModelingToolkit.connect(::Type{Pin}, ps...)
-    eqs = [
-           0 ~ sum(p->p.i, ps) # KCL
-          ]
-    # KVL
-    for i in 1:length(ps)-1
-        push!(eqs, ps[i].v ~ ps[i+1].v)
-    end
-
-    return eqs
-end
-```
-
 Finally we will connect the pieces of our circuit together. Let's connect the
 positive pin of the resistor to the source, the negative pin of the resistor
 to the capacitor, and the negative pin of the capacitor to a junction between
@@ -282,7 +255,8 @@ the source and the ground. This would mean our connection equations are:
 rc_eqs = [
           connect(source.p, resistor.p)
           connect(resistor.n, capacitor.p)
-          connect(capacitor.n, source.n, ground.g)
+          connect(capacitor.n, source.n)
+          connect(capacitor.n, ground.g)
          ]
 ```
 

docs/src/tutorials/spring_mass.md

@@ -6,14 +6,15 @@ In this tutorial we will build a simple component-based model of a spring-mass s
 
 ```julia
 using ModelingToolkit, Plots, DifferentialEquations, LinearAlgebra
+using Symbolics: scalarize
 
 @variables t
 D = Differential(t)
 
 function Mass(; name, m = 1.0, xy = [0., 0.], u = [0., 0.])
     ps = @parameters m=m
     sts = @variables pos[1:2](t)=xy v[1:2](t)=u
-    eqs = collect(D.(pos) .~ v)
+    eqs = scalarize(D.(pos) .~ v)
     ODESystem(eqs, t, [pos..., v...], ps; name)
 end
 
@@ -25,12 +26,12 @@ end
 
 function connect_spring(spring, a, b)
     [
-        spring.x ~ norm(collect(a .- b))
-        collect(spring.dir .~ collect(a .- b))
+        spring.x ~ norm(scalarize(a .- b))
+        scalarize(spring.dir .~ scalarize(a .- b))
     ]
 end
 
-spring_force(spring) = -spring.k .* collect(spring.dir) .* (spring.x - spring.l)  ./ spring.x
+spring_force(spring) = -spring.k .* scalarize(spring.dir) .* (spring.x - spring.l)  ./ spring.x
 
 m = 1.0
 xy = [1., -1.]
@@ -43,7 +44,7 @@ g = [0., -9.81]
 
 eqs = [
     connect_spring(spring, mass.pos, center)
-    collect(D.(mass.v) .~ spring_force(spring) / mass.m .+ g)
+    scalarize(D.(mass.v) .~ spring_force(spring) / mass.m .+ g)
 ]
 
 @named _model = ODESystem(eqs, t)
@@ -65,7 +66,7 @@ For each component we use a Julia function that returns an `ODESystem`. At the t
 function Mass(; name, m = 1.0, xy = [0., 0.], u = [0., 0.])
     ps = @parameters m=m
     sts = @variables pos[1:2](t)=xy v[1:2](t)=u
-    eqs = collect(D.(pos) .~ v)
+    eqs = scalarize(D.(pos) .~ v)
     ODESystem(eqs, t, [pos..., v...], ps; name)
 end
 ```
@@ -97,16 +98,16 @@ We now define functions that help construct the equations for a mass-spring syst
 ```julia
 function connect_spring(spring, a, b)
     [
-        spring.x ~ norm(collect(a .- b))
-        collect(spring.dir .~ collect(a .- b))
+        spring.x ~ norm(scalarize(a .- b))
+        scalarize(spring.dir .~ scalarize(a .- b))
     ]
 end
 ```
 
 Lastly, we define the `spring_force` function that takes a `spring` and returns the force exerted by this spring.
 
 ```julia
-spring_force(spring) = -spring.k .* collect(spring.dir) .* (spring.x - spring.l)  ./ spring.x
+spring_force(spring) = -spring.k .* scalarize(spring.dir) .* (spring.x - spring.l)  ./ spring.x
 ```
 
 To create our system, we will first create the components: a mass and a spring. This is done as follows:
@@ -127,7 +128,7 @@ We can now create the equations describing this system, by connecting `spring` t
 ```julia
 eqs = [
     connect_spring(spring, mass.pos, center)
-    collect(D.(mass.v) .~ spring_force(spring) / mass.m .+ g)
+    scalarize(D.(mass.v) .~ spring_force(spring) / mass.m .+ g)
 ]
 ```
 
@@ -230,4 +231,4 @@ We can also plot the path of the mass:
 plot(sol, vars = (mass.pos[1], mass.pos[2]))
 ```
 
-![plotpos](https://user-images.githubusercontent.com/23384717/130322197-cff35eb7-0739-471d-a3d9-af83d87f1cc7.png)
+![plotpos](https://user-images.githubusercontent.com/23384717/130322197-cff35eb7-0739-471d-a3d9-af83d87f1cc7.png)

docs/src/tutorials/tearing_parallelism.md

@@ -12,34 +12,11 @@ electrical circuits:
 ```julia
 using ModelingToolkit, OrdinaryDiffEq
 
-function connect_pin(ps...)
-    eqs = [
-           0 ~ sum(p->p.i, ps) # KCL
-          ]
-    # KVL
-    for i in 1:length(ps)-1
-        push!(eqs, ps[i].v ~ ps[i+1].v)
-    end
-
-    return eqs
-end
-
-function connect_heat(ps...)
-    eqs = [
-           0 ~ sum(p->p.Q_flow, ps)
-          ]
-    for i in 1:length(ps)-1
-        push!(eqs, ps[i].T ~ ps[i+1].T)
-    end
-
-    return eqs
-end
-
 # Basic electric components
 @variables t
 const D = Differential(t)
-function Pin(;name)
-    @variables v(t)=1.0 i(t)=1.0
+@connector function Pin(;name)
+    @variables v(t)=1.0 i(t)=1.0 [connect = Flow]
     ODESystem(Equation[], t, [v, i], [], name=name)
 end
 
@@ -61,8 +38,8 @@ function ConstantVoltage(;name, V = 1.0)
     compose(ODESystem(eqs, t, [], [V], name=name), p, n)
 end
 
-function HeatPort(;name)
-    @variables T(t)=293.15 Q_flow(t)=0.0
+@connector function HeatPort(;name)
+    @variables T(t)=293.15 Q_flow(t)=0.0 [connect = Flow]
     ODESystem(Equation[], t, [T, Q_flow], [], name=name)
 end
 
@@ -120,10 +97,10 @@ function parallel_rc_model(i; name, source, ground, R, C)
     heat_capacitor = HeatCapacitor(name=Symbol(:heat_capacitor, i))
 
     rc_eqs = [
-              connect_pin(source.p, resistor.p)
-              connect_pin(resistor.n, capacitor.p)
-              connect_pin(capacitor.n, source.n, ground.g)
-              connect_heat(resistor.h, heat_capacitor.h)
+              connect(source.p, resistor.p)
+              connect(resistor.n, capacitor.p)
+              connect(capacitor.n, source.n, ground.g)
+              connect(resistor.h, heat_capacitor.h)
              ]
 
     compose(ODESystem(rc_eqs, t, name=Symbol(name, i)),
@@ -139,8 +116,8 @@ we can connect a bunch of RC components as follows:
 
 ```julia
 V = 2.0
-source = ConstantVoltage(name=:source, V=V)
-ground = Ground(name=:ground)
+@named source = ConstantVoltage(V=V)
+@named ground = Ground()
 N = 50
 Rs = 10 .^range(0, stop=-4, length=N)
 Cs = 10 .^range(-3, stop=0, length=N)
@@ -170,19 +147,9 @@ Yes, that's a good question! Let's investigate a little bit more what had happen
 If you look at the system we defined:
 
 ```julia
-equations(big_rc)
-
-1051-element Vector{Equation}:
- Differential(t)(E(t)) ~ rc10₊resistor10₊h₊Q_flow(t) + rc11₊resistor11₊h₊Q_flow(t) + rc12₊resistor12₊h₊Q_flow(t) + rc13₊resistor13₊h₊Q_flow(t) + rc14₊resistor14₊h₊Q_flow(t) + rc15₊resistor15₊h₊Q_flow(t) + rc16₊resistor16₊h₊Q_flow(t) + rc17₊resistor17₊h₊Q_flow(t) + rc18₊resistor18₊h₊Q_flow(t) + rc19₊resistor19₊h₊Q_flow(t) + rc1₊resistor1₊h₊Q_flow(t) + rc20₊resistor20₊h₊Q_flow(t) + rc21₊resistor21₊h₊Q_flow(t) + rc22₊resistor22₊h₊Q_flow(t) + rc23₊resistor23₊h₊Q_flow(t) + rc24₊resistor24₊h₊Q_flow(t) + rc25₊resistor25₊h₊Q_flow(t) + rc26₊resistor26₊h₊Q_flow(t) + rc27₊resistor27₊h₊Q_flow(t) + rc28₊resistor28₊h₊Q_flow(t) + rc29₊resistor29₊h₊Q_flow(t) + rc2₊resistor2₊h₊Q_flow(t) + rc30₊resistor30₊h₊Q_flow(t) + rc31₊resistor31₊h₊Q_flow(t) + rc32₊resistor32₊h₊Q_flow(t) + rc33₊resistor33₊h₊Q_flow(t) + rc34₊resistor34₊h₊Q_flow(t) + rc35₊resistor35₊h₊Q_flow(t) + rc36₊resistor36₊h₊Q_flow(t) + rc37₊resistor37₊h₊Q_flow(t) + rc38₊resistor38₊h₊Q_flow(t) + rc39₊resistor39₊h₊Q_flow(t) + rc3₊resistor3₊h₊Q_flow(t) + rc40₊resistor40₊h₊Q_flow(t) + rc41₊resistor41₊h₊Q_flow(t) + rc42₊resistor42₊h₊Q_flow(t) + rc43₊resistor43₊h₊Q_flow(t) + rc44₊resistor44₊h₊Q_flow(t) + rc45₊resistor45₊h₊Q_flow(t) + rc46₊resistor46₊h₊Q_flow(t) + rc47₊resistor47₊h₊Q_flow(t) + rc48₊resistor48₊h₊Q_flow(t) + rc49₊resistor49₊h₊Q_flow(t) + rc4₊resistor4₊h₊Q_flow(t) + rc50₊resistor50₊h₊Q_flow(t) + rc5₊resistor5₊h₊Q_flow(t) + rc6₊resistor6₊h₊Q_flow(t) + rc7₊resistor7₊h₊Q_flow(t) + rc8₊resistor8₊h₊Q_flow(t) + rc9₊resistor9₊h₊Q_flow(t)
- 0 ~ rc1₊resistor1₊p₊i(t) + rc1₊source₊p₊i(t)
- rc1₊source₊p₊v(t) ~ rc1₊resistor1₊p₊v(t)
- 0 ~ rc1₊capacitor1₊p₊i(t) + rc1₊resistor1₊n₊i(t)
- rc1₊resistor1₊n₊v(t) ~ rc1₊capacitor1₊p₊v(t)
- ⋮
- rc50₊source₊V ~ rc50₊source₊p₊v(t) - rc50₊source₊n₊v(t)
- 0 ~ rc50₊source₊n₊i(t) + rc50₊source₊p₊i(t)
- rc50₊ground₊g₊v(t) ~ 0
- Differential(t)(rc50₊heat_capacitor50₊h₊T(t)) ~ rc50₊heat_capacitor50₊h₊Q_flow(t)*(rc50₊heat_capacitor50₊V^-1)*(rc50₊heat_capacitor50₊cp^-1)*(rc50₊heat_capacitor50₊rho^-1)
+length(equations(big_rc))
+
+801
 ```
 
 You see it started as a massive 1051 set of equations. However, after eliminating
@@ -211,7 +178,7 @@ investigate what this means:
 
 ```julia
 using ModelingToolkit.BipartiteGraphs
-ts = TearingState(big_rc)
+ts = TearingState(expand_connections(big_rc))
 inc_org = BipartiteGraphs.incidence_matrix(ts.graph)
 blt_org = StructuralTransformations.sorted_incidence_matrix(big_rc, only_algeqs=true, only_algvars=true)
 blt_reduced = StructuralTransformations.sorted_incidence_matrix(sys, only_algeqs=true, only_algvars=true)
@@ -223,7 +190,7 @@ The figure on the left is the original incidence matrix of the algebraic equatio
 Notice that the original formulation of the model has dependencies between different
 equations, and so the full set of equations must be solved together. That exposes
 no parallelism. However, the Block Lower Triangular (BLT) transformation exposes
-independent blocks. This is then further impoved by the tearing process, which
+independent blocks. This is then further improved by the tearing process, which
 removes 90% of the equations and transforms the nonlinear equations into 50
 independent blocks *which can now all be solved in parallel*. The conclusion
 is that, your attempts to parallelize are neigh: performing parallelism after