Merge pull request #1410 from SciML/myb/release

v8 Release

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.

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
@@ -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
-big_rc = initialize_system_structure(big_rc)
+big_rc = initialize_system_structure(expand_connections(big_rc))
 inc_org = BipartiteGraphs.incidence_matrix(structure(big_rc).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)