From 20f4e1d1f1dc74d02b1c69f5c7add88f93382e65 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Wed, 26 Aug 2026 08:41:29 -0400 Subject: [PATCH 1/3] Migrate the lecture docs to ModelingToolkit v11 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The Documentation build has been failing since MTK v11 became resolvable in the docs environment: the component DSL moved to SciCompDSL.jl, solvers are no longer re-exported by DifferentialEquations v8, `missing_variable_defaults` and `defaults` were removed, several standard-library components changed their constructor signatures, and symbolic defaults are now bindings rather than initial conditions. This updates the lecture examples to the current APIs so `makedocs` runs every `@example` block again, instead of capping MTK to v9. The hydraulic debugging model in lecture 6 is rewritten in the same functional `ODESystem` style used by the rest of that lecture, and its flattened variable names use `_` instead of the `₊` namespace separator, which MTK now strips defaults from. Co-authored-by: Cursor --- docs/Project.toml | 18 +- docs/make.jl | 79 ++++---- docs/src/lectures/lecture1.md | 19 +- docs/src/lectures/lecture2.md | 25 +-- docs/src/lectures/lecture6.md | 332 ++++++++++++++++++---------------- docs/src/lectures/lecture7.md | 1 + docs/src/lectures/lecture8.md | 10 +- docs/src/lectures/volume.jl | 7 +- 8 files changed, 259 insertions(+), 232 deletions(-) diff --git a/docs/Project.toml b/docs/Project.toml index f83ccb0..8cae03c 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -6,21 +6,31 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" ModelingToolkitCourse = "d0ecabdd-fe99-481e-9814-e16fa9c541df" ModelingToolkitStandardLibrary = "16a59e39-deab-5bd0-87e4-056b12336739" +NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +OrdinaryDiffEqBDF = "6ad6398a-0878-4a85-9266-38940aa047c8" +OrdinaryDiffEqNonlinearSolve = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8" +OrdinaryDiffEqSDIRK = "2d112036-d095-4a1e-ab9a-08536f3ecdbf" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +SciCompDSL = "91a8cdf1-4ca6-467b-a780-87fda3fff15e" Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46" Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" [compat] DataInterpolations = "5, 6, 8.10, 9.0" -DifferentialEquations = "7, 8.0" +DifferentialEquations = "8" Documenter = "1" ForwardDiff = "0.10, 1.4" -ModelingToolkit = "9, 11.28" -ModelingToolkitStandardLibrary = "2" -OrdinaryDiffEq = "=6.74.1, 7.1" +ModelingToolkit = "11.28" +ModelingToolkitStandardLibrary = "2.29" +NonlinearSolve = "4" +OrdinaryDiffEq = "7.1" +OrdinaryDiffEqBDF = "2" +OrdinaryDiffEqNonlinearSolve = "2" +OrdinaryDiffEqSDIRK = "2" Plots = "1" +SciCompDSL = "1" Setfield = "1" Sundials = "4, 6.2" Symbolics = "5, 6, 7.28" diff --git a/docs/make.jl b/docs/make.jl index 24964f4..d31a9a0 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,43 +1,42 @@ -using Documenter, ModelingToolkitCourse -# NOTE: OrdinaryDiffEq limited to v6.74.1 because of bug https://github.com/SciML/OrdinaryDiffEq.jl/issues/2250 - -pages = [ - "Home" => "index.md", - "lectures/lecture1.md", - "lectures/lecture2.md", - "lectures/lecture3.md", - "lectures/lecture4.md", - "lectures/lecture6.md", - "lectures/lecture7.md", - "lectures/lecture8.md", -] - -ENV["GKSwstype"] = "100" -using Plots - -makedocs( - sitename = "ModelingToolkit Course", - authors = "Chris Rackauckas", - modules = [ModelingToolkitCourse], - clean = true, doctest = false, linkcheck = true, - linkcheck_ignore = [ - "https://epubs.siam.org/doi/10.1137/0903023", - "https://link.springer.com/book/10.1007/978-3-642-05221-7", - "http://www.siam.org/journals/auth-info.php", - ], - format = Documenter.HTML( - assets = ["assets/favicon.ico"], - canonical = "https://docs.sciml.ai/ModelingToolkitCourse/stable/" - ), - pages = pages -) - +using Documenter, ModelingToolkitCourse + +pages = [ + "Home" => "index.md", + "lectures/lecture1.md", + "lectures/lecture2.md", + "lectures/lecture3.md", + "lectures/lecture4.md", + "lectures/lecture6.md", + "lectures/lecture7.md", + "lectures/lecture8.md", +] + +ENV["GKSwstype"] = "100" +using Plots + +makedocs( + sitename = "ModelingToolkit Course", + authors = "Chris Rackauckas", + modules = [ModelingToolkitCourse], + clean = true, doctest = false, linkcheck = true, + linkcheck_ignore = [ + "https://epubs.siam.org/doi/10.1137/0903023", + "https://link.springer.com/book/10.1007/978-3-642-05221-7", + "http://www.siam.org/journals/auth-info.php", + ], + format = Documenter.HTML( + assets = ["assets/favicon.ico"], + canonical = "https://docs.sciml.ai/ModelingToolkitCourse/stable/" + ), + pages = pages +) + #= using LiveServer serve(dir="build") -=# - -deploydocs( - repo = "github.com/SciML/ModelingToolkitCourse.git"; - push_preview = true -) +=# + +deploydocs( + repo = "github.com/SciML/ModelingToolkitCourse.git"; + push_preview = true +) diff --git a/docs/src/lectures/lecture1.md b/docs/src/lectures/lecture1.md index dcbc439..70f7bcf 100644 --- a/docs/src/lectures/lecture1.md +++ b/docs/src/lectures/lecture1.md @@ -86,6 +86,7 @@ Then we can solve by specifying the method, in this case we specify `NewtonRaphs ```@example l1 using DifferentialEquations +using NonlinearSolve p = xᵢ₋₁ = 0.0 # initial condition if i=2, x[1]=0 u0 = xᵢ = xᵢ₋₁ # guess value for x[i] @@ -186,6 +187,7 @@ ModelingToolkit.jl uses symbolic math from Symbolics.jl to provide automatic ind ```@example l1 using ModelingToolkit using ModelingToolkit: t_nounits as t, D_nounits as D +using SciCompDSL # provides the `@connector`/`@mtkmodel` component DSL nothing # hide ``` @@ -238,9 +240,8 @@ Notice how the 2nd derivative term `ẍ(t)` has been automatically determined fr We can now assemble a problem and solve it. The initial conditions do not need to be supplied here because the `sys` contains the variable defaults from `vars`. The solution object `sol` can now be indexed symbolically from any symbol of the system regardless if it's a solved variable, observable, or even a parameter. This way, if for example doing a batch of simulations, each respective solution object can easily retrieve all respective information about the simulation. ```@example l1 -u0 = [] # <-- used to override defaults of ODESystem variables -p = [] # <-- used to override defaults of ODESystem parameters -prob = ODEProblem(odesys, u0, tspan, p) +op = [] # <-- symbolic map used to override defaults of ODESystem variables and parameters +prob = ODEProblem(odesys, op, tspan) sol = solve(prob; abstol=tol) plot(sol; idxs=ẍ, xlabel="time [s]", ylabel="ẍ [m/s^2]") ``` @@ -368,7 +369,7 @@ nothing # hide Now the `Mass` and `Damper` components can be assembled in a system and connected together (note: the `connect` equation). Also note the parameters `v`, `m`, and `d` are defined to expose the properties which can be set as keyword arguments of the same name. ```@example l1 -@mtkmodel System begin +@mtkmodel MassDamperSystem begin @parameters begin v m @@ -383,7 +384,7 @@ Now the `Mass` and `Damper` components can be assembled in a system and connecte end end -@mtkbuild sys = System(;v=100, m=5, d=3) +@mtkbuild sys = MassDamperSystem(;v=100, m=5, d=3) nothing # hide ``` @@ -533,7 +534,7 @@ As with the `Reference` component, the force is a boundary condition and is leav Now let's assemble a *mass-spring-damper* system with the full collection of components. ```@example l1 -@mtkmodel System begin +@mtkmodel MassSpringDamperSystem begin @parameters begin v=0 x=0 @@ -555,7 +556,7 @@ Now let's assemble a *mass-spring-damper* system with the full collection of com end end -@mtkbuild sys = System() +@mtkbuild sys = MassSpringDamperSystem() ``` ![mass-spring-damper](../img/System2.png) @@ -638,7 +639,7 @@ nothing # hide As an example, the `MassSpringDamper` component can be connected in series to make a complex system. One can imagine then how this enables easy construction of complex models that can be quickly modified, extremely useful for the application of model based design. ```@example l1 -@mtkmodel System begin +@mtkmodel SeriesSystem begin @parameters begin v = 0 x = 0 @@ -658,7 +659,7 @@ As an example, the `MassSpringDamper` component can be connected in series to ma end end -@mtkbuild sys = System() +@mtkbuild sys = SeriesSystem() nothing #hide ``` diff --git a/docs/src/lectures/lecture2.md b/docs/src/lectures/lecture2.md index 0778ab8..ef1a94f 100644 --- a/docs/src/lectures/lecture2.md +++ b/docs/src/lectures/lecture2.md @@ -290,9 +290,9 @@ odesys_ṁ1 Notice that now, with a simple change of the system input variable, `structural_simplify()` outputs a system with 4 states to be solved. We can find the initial conditions needed for these states from `sol_x` and solve. ```@example l2 -u0 = [sol_x[s][1] for s in unknowns(odesys_ṁ1)] -prob_ṁ1 = ODEProblem(odesys_ṁ1, u0, (0, t_end)) -@time sol_ṁ1 = solve(prob_ṁ1; initializealg=NoInit()); +u0 = [s => sol_x[s][1] for s in unknowns(odesys_ṁ1)] +prob_ṁ1 = ODEProblem(odesys_ṁ1, u0, (0, t_end); build_initializeprob = false) +@time sol_ṁ1 = solve(prob_ṁ1; initializealg=SciMLBase.NoInit()); nothing # hide ``` @@ -307,8 +307,8 @@ If we now solve for case 2, we can study the impact the compressibility derivati ```@example l2 @mtkbuild odesys_ṁ2 = ODESystem(eqs_ṁ2, t, vars, pars) -prob_ṁ2 = ODEProblem(odesys_ṁ2, u0, (0, t_end)) -@time sol_ṁ2 = solve(prob_ṁ2; initializealg=NoInit()); +prob_ṁ2 = ODEProblem(odesys_ṁ2, u0, (0, t_end); build_initializeprob = false) +@time sol_ṁ2 = solve(prob_ṁ2; initializealg=SciMLBase.NoInit()); nothing # hide ``` @@ -333,13 +333,13 @@ plot(time, (sol_ṁ1(time)[x] .- sol_ṁ2(time)[x])/1e-3, Also note the difference in computation. ```@repl l2 -sol_ṁ1.destats +sol_ṁ1.stats ``` As can be seen, including the detail of full compressibility resulted in more computation: more function evaluations, Jacobians, solves, and steps. ```@repl l2 -sol_ṁ2.destats +sol_ṁ2.stats ``` ### ModelingToolkitStandardLibrary.jl @@ -350,6 +350,8 @@ import ModelingToolkitStandardLibrary.Mechanical.Translational as T import ModelingToolkitStandardLibrary.Hydraulic.IsothermalCompressible as IC import ModelingToolkitStandardLibrary.Blocks as B +include("volume.jl") # <-- moving-wall `Volume` component + using DataInterpolations mass_flow_fun = LinearInterpolation(sol_x[ṁ], sol_x.t) @@ -370,9 +372,9 @@ function MassVolume(; name, dx, drho, dm) vars = [] systems = @named begin fluid = IC.HydraulicFluid(; density = 876, bulk_modulus = 1.2e9) - mass = T.Mass(;v=dx,m=M,g=-g) - vol = IC.Volume(;area=A, x=x₀, p=p_int, dx, drho, dm) - mass_flow = IC.MassFlow(;p_int) + mass = T.Mass(;m=M,g=-g) + vol = Volume(;area=A, x=x₀, p=p_int, dx, drho, dm) + mass_flow = IC.MassFlow() mass_flow_input = B.TimeVaryingFunction(;f = mass_flow_fun) end @@ -383,7 +385,8 @@ function MassVolume(; name, dx, drho, dm) connect(mass_flow.port, fluid) ] - return ODESystem(eqs, t, vars, pars; systems, name) + return ODESystem(eqs, t, vars, pars; systems, name, + initial_conditions = [mass.v => dx, mass.s => x₀]) end dx = sol_x[ẋ][1] diff --git a/docs/src/lectures/lecture6.md b/docs/src/lectures/lecture6.md index 2a2baa8..77499da 100644 --- a/docs/src/lectures/lecture6.md +++ b/docs/src/lectures/lecture6.md @@ -17,6 +17,8 @@ It is not always the case, but for most models, the unsimplified system should g ```@example l6 using ModelingToolkit, DifferentialEquations, Plots using ModelingToolkit: t_nounits as t, D_nounits as D +using OrdinaryDiffEqNonlinearSolve # ShampineCollocationInit, NLNewton +using OrdinaryDiffEqSDIRK # ImplicitEuler pars = @parameters m = 1 g = 1 L = 1 Φ=0 @@ -62,7 +64,7 @@ If we attempt to solve this system we can see that it only solves up to the poin ```@example l6 sys = complete(structural_simplify(pendulum)) -prob = ODEProblem(sys, ModelingToolkit.missing_variable_defaults(sys), (0, 10)) +prob = ODEProblem(sys, [], (0, 10)) sol = solve(prob)# gives retcode: DtLessThanMin plot(sol; idxs=[x,y]) ``` @@ -95,7 +97,7 @@ eqs = [ @named stiffness_pendulum = ODESystem(eqs, t, vars, pars) sys = structural_simplify(stiffness_pendulum) -prob = ODEProblem(sys, ModelingToolkit.missing_variable_defaults(sys), (0, 10)) +prob = ODEProblem(sys, [], (0, 10)) sol = solve(prob)# Success plot(sol; idxs=[x,y]) ``` @@ -107,7 +109,7 @@ In some cases we can apply `dae_index_lowering()` to further simplify the proble ```@example l6 sys = structural_simplify(dae_index_lowering(pendulum)) -prob = ODEProblem(sys, ModelingToolkit.missing_variable_defaults(sys), (0, 10)) +prob = ODEProblem(sys, [], (0, 10)) ref = solve(prob) plot(ref; idxs=x, label="dae_index_lowering") plot!(sol; idxs=x, label="compliance") @@ -157,9 +159,9 @@ function MassVolume(solves_force = true; name) vars = [] systems = @named begin fluid = IC.HydraulicFluid(; density = 876, bulk_modulus = 1.2e9) - mass = T.Mass(;v=dx,m=M,g=-g) + mass = T.Mass(;m=M,g=-g) vol = Volume(;area=A, x=x₀, p=p_int, dx, drho, dm) # <-- missing Volume component from MTKSL (will be released in new version) - mass_flow = IC.Open(;p_int) + mass_flow = IC.Open(;p=p_int) position = T.Position(solves_force) position_input = B.TimeVaryingFunction(;f = t -> amp*sin(2π*t*f) + x₀) end @@ -170,7 +172,8 @@ function MassVolume(solves_force = true; name) connect(position.s, position_input.output) ] - return ODESystem(eqs, t, vars, pars; systems, name) + return ODESystem(eqs, t, vars, pars; systems, name, + initial_conditions = [mass.v => dx]) end @named odesys = MassVolume() @@ -194,151 +197,160 @@ The reason for the mismatch is that the input boundary condition `Position()` ne It's very difficult to identify what is wrong with a model if it's not outputting any data. This section discusses ways to force a model solution. It's still possible that something with the model is wrong, but the best way to know that is to see what the equations are outputting. For example if the model is simulating negative pressure, but negative pressure is impossible, then this is a good clue of what is wrong with the model! The strategies for forcing a model solve will come from a simple hydraulic system that is attempting to start a hydraulic cylinder at a high pressure differential. See [ModelingToolkit Industrial Example](https://github.com/bradcarman/ModelingToolkitWebinar) for more information about the model. ```@example l6 -@mtkmodel System begin - @parameters begin - res₁₊Cₒ = 2.7 - res₁₊Aₒ = 0.00094 - res₁₊ρ₀ = 1000 - res₁₊p′ = 3.0e7 - res₂₊Cₒ = 2.7 - res₂₊Aₒ = 0.00094 - res₂₊ρ₀ = 1000 - res₂₊p′ = 0 - act₊p₁′ = 3.0e7 - act₊p₂′ = 0 - act₊vol₁₊A = 0.1 - act₊vol₁₊ρ₀ = 1000 - act₊vol₁₊β = 2.0e9 - act₊vol₁₊direction = -1 - act₊vol₁₊p′ = act₊p₁′ - act₊vol₁₊x′ = 0.5 - act₊vol₂₊A = 0.1 - act₊vol₂₊ρ₀ = 1000 - act₊vol₂₊β = 2.0e9 - act₊vol₂₊direction = 1 - act₊vol₂₊p′ = act₊p₂′ - act₊vol₂₊x′ = 0.5 - act₊mass₊m = 100 - act₊mass₊f′ = 0.1(-act₊p₁′ + act₊p₂′) - src₊p′ = 3.0e7 - snk₊p′ = 0 - dmp₊c = 1000 - end - @variables begin - res₁₊ṁ(t) = 0 - res₁₊p₁(t) = res₁₊p′ - res₁₊p₂(t) = res₁₊p′ - res₁₊port₁₊p(t) = res₁₊p′ - res₁₊port₁₊ṁ(t) = 0 - res₁₊port₂₊p(t) = res₁₊p′ - res₁₊port₂₊ṁ(t) = 0 - res₂₊ṁ(t) = 0 - res₂₊p₁(t) = res₂₊p′ - res₂₊p₂(t) = res₂₊p′ - res₂₊port₁₊p(t) = res₂₊p′ - res₂₊port₁₊ṁ(t) = 0 - res₂₊port₂₊p(t) = res₂₊p′ - res₂₊port₂₊ṁ(t) = 0 - act₊port₁₊p(t) = act₊p₁′ - act₊port₁₊ṁ(t) = 0 - act₊port₂₊p(t) = act₊p₂′ - act₊port₂₊ṁ(t) = 0 - act₊vol₁₊p(t) = act₊vol₁₊p′ - act₊vol₁₊x(t) = act₊vol₁₊x′ - act₊vol₁₊ṁ(t) = 0 - act₊vol₁₊f(t) = act₊vol₁₊A * act₊vol₁₊p′ - act₊vol₁₊ẋ(t) = 0 - act₊vol₁₊r(t) = act₊vol₁₊ρ₀ * (1 + act₊vol₁₊p′ / act₊vol₁₊β) - act₊vol₁₊ṙ(t) = 0 - act₊vol₁₊port₊p(t) = act₊vol₁₊p′ - act₊vol₁₊port₊ṁ(t) = 0 - act₊vol₁₊flange₊ẋ(t) = 0 - act₊vol₁₊flange₊f(t) = -act₊vol₁₊A * act₊vol₁₊direction * act₊vol₁₊p′ - act₊vol₂₊p(t) = act₊vol₂₊p′ - act₊vol₂₊x(t) = act₊vol₂₊x′ - act₊vol₂₊ṁ(t) = 0 - act₊vol₂₊f(t) = act₊vol₂₊A * act₊vol₂₊p′ - act₊vol₂₊ẋ(t) = 0 - act₊vol₂₊r(t) = act₊vol₂₊ρ₀ * (1 + act₊vol₂₊p′ / act₊vol₂₊β) - act₊vol₂₊ṙ(t) = 0 - act₊vol₂₊port₊p(t) = act₊vol₂₊p′ - act₊vol₂₊port₊ṁ(t) = 0 - act₊vol₂₊flange₊ẋ(t) = 0 - act₊vol₂₊flange₊f(t) = -act₊vol₂₊A * act₊vol₂₊direction * act₊vol₂₊p′ - act₊mass₊f(t) = act₊mass₊f′ - act₊mass₊x(t) = 0 - act₊mass₊ẋ(t) = 0 - act₊mass₊ẍ(t) = act₊mass₊f′ / act₊mass₊m - act₊mass₊flange₊ẋ(t) = 0 - act₊mass₊flange₊f(t) = act₊mass₊f′ - act₊flange₊ẋ(t) = 0 - act₊flange₊f(t) = 0 - src₊port₊p(t) = src₊p′ - src₊port₊ṁ(t) = 0 - snk₊port₊p(t) = snk₊p′ - snk₊port₊ṁ(t) = 0 - dmp₊flange₊ẋ(t) = 0 - dmp₊flange₊f(t) = 0 - end - @equations begin - res₁₊ṁ ~ res₁₊port₁₊ṁ - res₁₊ṁ ~ -res₁₊port₂₊ṁ - res₁₊p₁ ~ res₁₊port₁₊p - res₁₊p₂ ~ res₁₊port₂₊p - -res₁₊p₂ + res₁₊p₁ ~ 0.5res₁₊Cₒ * res₁₊ρ₀ * ((res₁₊ṁ / (res₁₊Aₒ * res₁₊ρ₀))^2) - res₂₊ṁ ~ res₂₊port₁₊ṁ - res₂₊ṁ ~ -res₂₊port₂₊ṁ - res₂₊p₁ ~ res₂₊port₁₊p - res₂₊p₂ ~ res₂₊port₂₊p - -res₂₊p₂ + res₂₊p₁ ~ 0.5res₂₊Cₒ * res₂₊ρ₀ * ((res₂₊ṁ / (res₂₊Aₒ * res₂₊ρ₀))^2) - D(act₊vol₁₊x) ~ act₊vol₁₊ẋ - D(act₊vol₁₊r) ~ act₊vol₁₊ṙ - act₊vol₁₊p ~ act₊vol₁₊port₊p - act₊vol₁₊ṁ ~ act₊vol₁₊port₊ṁ - act₊vol₁₊f ~ -act₊vol₁₊direction * act₊vol₁₊flange₊f - act₊vol₁₊ẋ ~ act₊vol₁₊direction * act₊vol₁₊flange₊ẋ - act₊vol₁₊r ~ act₊vol₁₊ρ₀ * (1 + act₊vol₁₊p / act₊vol₁₊β) - act₊vol₁₊ṁ ~ act₊vol₁₊A * act₊vol₁₊ẋ * act₊vol₁₊r + act₊vol₁₊A * act₊vol₁₊x * act₊vol₁₊ṙ - act₊vol₁₊f ~ act₊vol₁₊A * act₊vol₁₊p - D(act₊vol₂₊x) ~ act₊vol₂₊ẋ - D(act₊vol₂₊r) ~ act₊vol₂₊ṙ - act₊vol₂₊p ~ act₊vol₂₊port₊p - act₊vol₂₊ṁ ~ act₊vol₂₊port₊ṁ - act₊vol₂₊f ~ -act₊vol₂₊direction * act₊vol₂₊flange₊f - act₊vol₂₊ẋ ~ act₊vol₂₊direction * act₊vol₂₊flange₊ẋ - act₊vol₂₊r ~ act₊vol₂₊ρ₀ * (1 + act₊vol₂₊p / act₊vol₂₊β) - act₊vol₂₊ṁ ~ act₊vol₂₊A * act₊vol₂₊r * act₊vol₂₊ẋ + act₊vol₂₊A * act₊vol₂₊ṙ * act₊vol₂₊x - act₊vol₂₊f ~ act₊vol₂₊A * act₊vol₂₊p - D(act₊mass₊x) ~ act₊mass₊ẋ - D(act₊mass₊ẋ) ~ act₊mass₊ẍ - act₊mass₊f ~ act₊mass₊flange₊f - act₊mass₊ẋ ~ act₊mass₊flange₊ẋ - act₊mass₊m * act₊mass₊ẍ ~ act₊mass₊f - src₊port₊p ~ src₊p′ - snk₊port₊p ~ snk₊p′ - dmp₊flange₊f ~ dmp₊c * dmp₊flange₊ẋ - src₊port₊p ~ res₁₊port₁₊p - 0 ~ res₁₊port₁₊ṁ + src₊port₊ṁ - res₁₊port₂₊p ~ act₊port₁₊p - 0 ~ act₊port₁₊ṁ + res₁₊port₂₊ṁ - act₊port₂₊p ~ res₂₊port₁₊p - 0 ~ act₊port₂₊ṁ + res₂₊port₁₊ṁ - res₂₊port₂₊p ~ snk₊port₊p - 0 ~ res₂₊port₂₊ṁ + snk₊port₊ṁ - dmp₊flange₊ẋ ~ act₊flange₊ẋ - 0 ~ act₊flange₊f + dmp₊flange₊f - act₊port₁₊p ~ act₊vol₁₊port₊p - 0 ~ act₊vol₁₊port₊ṁ - act₊port₁₊ṁ - act₊port₂₊p ~ act₊vol₂₊port₊p - 0 ~ -act₊port₂₊ṁ + act₊vol₂₊port₊ṁ - act₊vol₁₊flange₊ẋ ~ act₊vol₂₊flange₊ẋ - act₊vol₁₊flange₊ẋ ~ act₊mass₊flange₊ẋ - act₊vol₁₊flange₊ẋ ~ act₊flange₊ẋ - 0 ~ act₊vol₁₊flange₊f - act₊flange₊f + act₊vol₂₊flange₊f + act₊mass₊flange₊f - end +pars = @parameters begin + res₁_Cₒ = 2.7 + res₁_Aₒ = 0.00094 + res₁_ρ₀ = 1000 + res₁_p′ = 3.0e7 + res₂_Cₒ = 2.7 + res₂_Aₒ = 0.00094 + res₂_ρ₀ = 1000 + res₂_p′ = 0 + act_p₁′ = 3.0e7 + act_p₂′ = 0 + act_vol₁_A = 0.1 + act_vol₁_ρ₀ = 1000 + act_vol₁_β = 2.0e9 + act_vol₁_direction = -1 + act_vol₁_p′ = act_p₁′ + act_vol₁_x′ = 0.5 + act_vol₂_A = 0.1 + act_vol₂_ρ₀ = 1000 + act_vol₂_β = 2.0e9 + act_vol₂_direction = 1 + act_vol₂_p′ = act_p₂′ + act_vol₂_x′ = 0.5 + act_mass_m = 100 + act_mass_f′ = 0.1(-act_p₁′ + act_p₂′) + src_p′ = 3.0e7 + snk_p′ = 0 + dmp_c = 1000 +end + +vars = @variables begin + res₁_ṁ(t) = 0 + res₁_p₁(t) = res₁_p′ + res₁_p₂(t) = res₁_p′ + res₁_port₁_p(t) = res₁_p′ + res₁_port₁_ṁ(t) = 0 + res₁_port₂_p(t) = res₁_p′ + res₁_port₂_ṁ(t) = 0 + res₂_ṁ(t) = 0 + res₂_p₁(t) = res₂_p′ + res₂_p₂(t) = res₂_p′ + res₂_port₁_p(t) = res₂_p′ + res₂_port₁_ṁ(t) = 0 + res₂_port₂_p(t) = res₂_p′ + res₂_port₂_ṁ(t) = 0 + act_port₁_p(t) = act_p₁′ + act_port₁_ṁ(t) = 0 + act_port₂_p(t) = act_p₂′ + act_port₂_ṁ(t) = 0 + act_vol₁_p(t) = act_vol₁_p′ + act_vol₁_x(t) + act_vol₁_ṁ(t) = 0 + act_vol₁_f(t) = act_vol₁_A * act_vol₁_p′ + act_vol₁_ẋ(t) = 0 + act_vol₁_r(t) + act_vol₁_ṙ(t) = 0 + act_vol₁_port_p(t) = act_vol₁_p′ + act_vol₁_port_ṁ(t) = 0 + act_vol₁_flange_ẋ(t) = 0 + act_vol₁_flange_f(t) = -act_vol₁_A * act_vol₁_direction * act_vol₁_p′ + act_vol₂_p(t) = act_vol₂_p′ + act_vol₂_x(t) + act_vol₂_ṁ(t) = 0 + act_vol₂_f(t) = act_vol₂_A * act_vol₂_p′ + act_vol₂_ẋ(t) = 0 + act_vol₂_r(t) + act_vol₂_ṙ(t) = 0 + act_vol₂_port_p(t) = act_vol₂_p′ + act_vol₂_port_ṁ(t) = 0 + act_vol₂_flange_ẋ(t) = 0 + act_vol₂_flange_f(t) = -act_vol₂_A * act_vol₂_direction * act_vol₂_p′ + act_mass_f(t) = act_mass_f′ + act_mass_x(t) = 0 + act_mass_ẋ(t) = 0 + act_mass_ẍ(t) = act_mass_f′ / act_mass_m + act_mass_flange_ẋ(t) = 0 + act_mass_flange_f(t) = act_mass_f′ + act_flange_ẋ(t) = 0 + act_flange_f(t) = 0 + src_port_p(t) = src_p′ + src_port_ṁ(t) = 0 + snk_port_p(t) = snk_p′ + snk_port_ṁ(t) = 0 + dmp_flange_ẋ(t) = 0 + dmp_flange_f(t) = 0 end -@mtkbuild sys = System() +eqs = [ + res₁_ṁ ~ res₁_port₁_ṁ + res₁_ṁ ~ -res₁_port₂_ṁ + res₁_p₁ ~ res₁_port₁_p + res₁_p₂ ~ res₁_port₂_p + -res₁_p₂ + res₁_p₁ ~ 0.5res₁_Cₒ * res₁_ρ₀ * ((res₁_ṁ / (res₁_Aₒ * res₁_ρ₀))^2) + res₂_ṁ ~ res₂_port₁_ṁ + res₂_ṁ ~ -res₂_port₂_ṁ + res₂_p₁ ~ res₂_port₁_p + res₂_p₂ ~ res₂_port₂_p + -res₂_p₂ + res₂_p₁ ~ 0.5res₂_Cₒ * res₂_ρ₀ * ((res₂_ṁ / (res₂_Aₒ * res₂_ρ₀))^2) + D(act_vol₁_x) ~ act_vol₁_ẋ + D(act_vol₁_r) ~ act_vol₁_ṙ + act_vol₁_p ~ act_vol₁_port_p + act_vol₁_ṁ ~ act_vol₁_port_ṁ + act_vol₁_f ~ -act_vol₁_direction * act_vol₁_flange_f + act_vol₁_ẋ ~ act_vol₁_direction * act_vol₁_flange_ẋ + act_vol₁_r ~ act_vol₁_ρ₀ * (1 + act_vol₁_p / act_vol₁_β) + act_vol₁_ṁ ~ act_vol₁_A * act_vol₁_ẋ * act_vol₁_r + act_vol₁_A * act_vol₁_x * act_vol₁_ṙ + act_vol₁_f ~ act_vol₁_A * act_vol₁_p + D(act_vol₂_x) ~ act_vol₂_ẋ + D(act_vol₂_r) ~ act_vol₂_ṙ + act_vol₂_p ~ act_vol₂_port_p + act_vol₂_ṁ ~ act_vol₂_port_ṁ + act_vol₂_f ~ -act_vol₂_direction * act_vol₂_flange_f + act_vol₂_ẋ ~ act_vol₂_direction * act_vol₂_flange_ẋ + act_vol₂_r ~ act_vol₂_ρ₀ * (1 + act_vol₂_p / act_vol₂_β) + act_vol₂_ṁ ~ act_vol₂_A * act_vol₂_r * act_vol₂_ẋ + act_vol₂_A * act_vol₂_ṙ * act_vol₂_x + act_vol₂_f ~ act_vol₂_A * act_vol₂_p + D(act_mass_x) ~ act_mass_ẋ + D(act_mass_ẋ) ~ act_mass_ẍ + act_mass_f ~ act_mass_flange_f + act_mass_ẋ ~ act_mass_flange_ẋ + act_mass_m * act_mass_ẍ ~ act_mass_f + src_port_p ~ src_p′ + snk_port_p ~ snk_p′ + dmp_flange_f ~ dmp_c * dmp_flange_ẋ + src_port_p ~ res₁_port₁_p + 0 ~ res₁_port₁_ṁ + src_port_ṁ + res₁_port₂_p ~ act_port₁_p + 0 ~ act_port₁_ṁ + res₁_port₂_ṁ + act_port₂_p ~ res₂_port₁_p + 0 ~ act_port₂_ṁ + res₂_port₁_ṁ + res₂_port₂_p ~ snk_port_p + 0 ~ res₂_port₂_ṁ + snk_port_ṁ + dmp_flange_ẋ ~ act_flange_ẋ + 0 ~ act_flange_f + dmp_flange_f + act_port₁_p ~ act_vol₁_port_p + 0 ~ act_vol₁_port_ṁ - act_port₁_ṁ + act_port₂_p ~ act_vol₂_port_p + 0 ~ -act_port₂_ṁ + act_vol₂_port_ṁ + act_vol₁_flange_ẋ ~ act_vol₂_flange_ẋ + act_vol₁_flange_ẋ ~ act_mass_flange_ẋ + act_vol₁_flange_ẋ ~ act_flange_ẋ + 0 ~ act_vol₁_flange_f - act_flange_f + act_vol₂_flange_f + act_mass_flange_f +] + +# Initial conditions given as expressions of parameters are supplied separately, since a +# symbolic default in `@variables` is a binding that cannot be overridden later on. +ics = [ + act_vol₁_x => act_vol₁_x′ + act_vol₁_r => act_vol₁_ρ₀ * (1 + act_vol₁_p′ / act_vol₁_β) + act_vol₂_x => act_vol₂_x′ + act_vol₂_r => act_vol₂_ρ₀ * (1 + act_vol₂_p′ / act_vol₂_β) +] + +@mtkbuild sys = ODESystem(eqs, t, vars, pars; initial_conditions = ics) prob = ODEProblem(sys, [], (0, 0.1)) sol = solve(prob) ``` @@ -350,11 +362,11 @@ First, let's check the initial conditions to see if at time 0 we are starting wi ```@example l6 eqs = full_equations(sys) -defs = ModelingToolkit.defaults(sys) +defs = merge(Dict(ModelingToolkit.bindings(sys)), Dict(ModelingToolkit.initial_conditions(sys))) residuals = Float64[] for eq in eqs if !ModelingToolkit.isdifferential(eq.lhs) - push!(residuals, ModelingToolkit.fixpoint_sub(eq.rhs, defs)) + push!(residuals, ModelingToolkit.value(ModelingToolkit.fixpoint_sub(eq.rhs, defs))) end end residuals @@ -371,12 +383,12 @@ The `ShampineCollocationInit` solves the initial conditions by essentially takin ```@example l6 prob = ODEProblem(sys, [], (0, dt)) -sol = solve(prob, ImplicitEuler(nlsolve=NLNewton(check_div=false, always_new=true, relax=4/10, max_iter=100)); dt, adaptive=false) +sol = solve(prob, ImplicitEuler(nlsolve=OrdinaryDiffEqNonlinearSolve.NLNewton(check_div=false, always_new=true, relax=4/10, max_iter=100)); dt, adaptive=false) # update u0 with the ImplicitEuler non-adaptive step -prob′ = ODEProblem(sys, sol[2], (0, 0.1)) -sol = solve(prob′); -plot(sol; idxs=sys.act₊mass₊ẋ) +prob′ = ODEProblem(sys, unknowns(sys) .=> sol.u[2], (0, 0.1)) +sol = solve(prob′; initializealg=SciMLBase.NoInit()); +plot(sol; idxs=sys.act_mass_ẋ) ``` As can be seen, now we have a successful solve. We can see the change to the initial conditions is very minimal. As can be seen, the solver needs the derivative terms to be offset by a small amount. @@ -388,12 +400,12 @@ println(join(["$s : $(round(x; digits=3)) -> $(round(y; digits=3))" for (s,x,y) Another strategy that can help issues with initial conditions is to offset or perturb any initial conditions from 0 by a small value. ### Adjust tolerance -Here we get a solve by increasing the `abstol` and `reltol` to very large values. This is therefore understood to give us a very low resolution solution that is far from the true solution, but we can now at least see if the model is calculating generally correct values, at least with the correct sign. Here we expect the `act₊mass₊ẋ` to be around -1 and that's exactly what we get. However, as can be seen the tolerance is too open to resolve the dynamics. +Here we get a solve by increasing the `abstol` and `reltol` to very large values. This is therefore understood to give us a very low resolution solution that is far from the true solution, but we can now at least see if the model is calculating generally correct values, at least with the correct sign. Here we expect the `act_mass_ẋ` to be around -1 and that's exactly what we get. However, as can be seen the tolerance is too open to resolve the dynamics. ```@example l6 prob = ODEProblem(sys, [], (0, 0.1)) sol = solve(prob, ImplicitEuler(); abstol=10000.0, reltol=100.0) -plot(sol; idxs=sys.act₊mass₊ẋ) +plot(sol; idxs=sys.act_mass_ẋ) ``` @@ -402,8 +414,8 @@ Another strategy similar to adjusting tolerance is to turn off adaptivity. This ```@example l6 prob = ODEProblem(sys, [], (0, 0.1)) -sol = solve(prob, ImplicitEuler(nlsolve=NLNewton(check_div=false, always_new=true, relax=4/10, max_iter=100)); initializealg=NoInit(), adaptive=false, dt=1e-6) -plot(sol; idxs=sys.act₊mass₊ẋ) +sol = solve(prob, ImplicitEuler(nlsolve=OrdinaryDiffEqNonlinearSolve.NLNewton(check_div=false, always_new=true, relax=4/10, max_iter=100)); initializealg=SciMLBase.NoInit(), adaptive=false, dt=1e-6) +plot(sol; idxs=sys.act_mass_ẋ) ``` Note the use of keywords: @@ -494,7 +506,7 @@ Implementing this for the hydraulic system works well, giving an adaptive time s odesys = dae_to_ode(sys) prob = ODEProblem(odesys, [], (0,0.1)) sol = solve(prob) -plot(sol; idxs=sys.act₊mass₊ẋ) +plot(sol; idxs=sys.act_mass_ẋ) ``` Note this problem, as we've seen, has a lot of trouble with initialization. Note how the first 200 steps are taken with a very small time step. The `Tsit5` solver is able to successfully push through the model initialization and then solve the remaining steps at a reasonable time step. diff --git a/docs/src/lectures/lecture7.md b/docs/src/lectures/lecture7.md index 3df58c4..0bc1b12 100644 --- a/docs/src/lectures/lecture7.md +++ b/docs/src/lectures/lecture7.md @@ -14,6 +14,7 @@ numerical solvers and ModelingToolkit behave. ```@example l7 using DifferentialEquations, Sundials, ModelingToolkit, Plots, LinearAlgebra using ModelingToolkit: t_nounits as t, D_nounits as D +using OrdinaryDiffEqBDF # DFBDF function f!(out, du, u, p, t) # u[1]: x, du[1]: x' diff --git a/docs/src/lectures/lecture8.md b/docs/src/lectures/lecture8.md index 080b410..eeb6bbf 100644 --- a/docs/src/lectures/lecture8.md +++ b/docs/src/lectures/lecture8.md @@ -256,7 +256,7 @@ eqs = [ @named pend = ODESystem(eqs,t) pend = complete(pend) ss = structural_simplify(pend) -prob_ir = ODEProblem(ss, [ModelingToolkit.missing_variable_defaults(ss); x => 1], (0.0, 25.0), [g => 1]) +prob_ir = ODEProblem(ss, [x => 1, g => 1], (0.0, 25.0)) sol = solve(prob_ir, Rodas5P()) plot(sol, idxs = (x, y), lab = "pendulum", aspect_ratio = 1) ``` @@ -288,7 +288,7 @@ eqs = [ @named pend = ODESystem(eqs,t) pend = complete(pend) ss = structural_simplify(pend) -prob_ir = ODEProblem(ss, ModelingToolkit.missing_variable_defaults(ss), (0.0, 25.0), [g => 1]) +prob_ir = ODEProblem(ss, [g => 1], (0.0, 25.0)) sol = solve(prob_ir, Rodas5P()) plot(sol, idxs = (x, y), lab = "pendulum", aspect_ratio = 1) ``` @@ -333,16 +333,14 @@ eqs = [ @named pend = ODESystem(eqs,t) pend = complete(pend) ss = structural_simplify(pend) -prob_ir = ODEProblem(ss, - [ModelingToolkit.missing_variable_defaults(ss); θ2=>1.4], - (0.0, 25.0), [g => 1]) +prob_ir = ODEProblem(ss, [θ2 => 1.4, g => 1], (0.0, 25.0)) sol = solve(prob_ir, Rodas5P(), reltol=1e-7, abstol=1e-9) plot(sol, idxs = (x1, y1)) plot!(sol, idxs = (x2, y2), xlab = "x", ylab = "y", aspect_ratio=1, dpi=400) ``` ```@example l8 -plot(sol, idxs = [T+V-sol[T+V, 1]]) +plot(sol, idxs = [T+V-(sol[T, 1]+sol[V, 1])]) ``` ![](../img/double_pendulum.mp4) diff --git a/docs/src/lectures/volume.jl b/docs/src/lectures/volume.jl index e5d13e8..a30326b 100644 --- a/docs/src/lectures/volume.jl +++ b/docs/src/lectures/volume.jl @@ -32,7 +32,7 @@ using ModelingToolkitStandardLibrary.Hydraulic.IsothermalCompressible: liquid_de end systems = @named begin - port = HydraulicPort(; p_int = p) + port = HydraulicPort() flange = MechanicalPort(; f, v = dx) end @@ -53,5 +53,8 @@ using ModelingToolkitStandardLibrary.Hydraulic.IsothermalCompressible: liquid_de dm ~ drho * x * area + rho * dx * area ] - ODESystem(eqs, t, vars, pars; name, systems, defaults = [rho => liquid_density(port)]) + ODESystem( + eqs, t, vars, pars; name, systems, + initial_conditions = [rho => liquid_density(port)] + ) end From 1e95a79a8a4cab4c9b0b0d10689c7eead4e34e0f Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Wed, 26 Aug 2026 08:54:17 -0400 Subject: [PATCH 2/3] Keep docs/make.jl line endings unchanged Co-authored-by: Cursor --- docs/make.jl | 78 ++++++++++++++++++++++++++-------------------------- 1 file changed, 39 insertions(+), 39 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index d31a9a0..0b9314c 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,42 +1,42 @@ -using Documenter, ModelingToolkitCourse - -pages = [ - "Home" => "index.md", - "lectures/lecture1.md", - "lectures/lecture2.md", - "lectures/lecture3.md", - "lectures/lecture4.md", - "lectures/lecture6.md", - "lectures/lecture7.md", - "lectures/lecture8.md", -] - -ENV["GKSwstype"] = "100" -using Plots - -makedocs( - sitename = "ModelingToolkit Course", - authors = "Chris Rackauckas", - modules = [ModelingToolkitCourse], - clean = true, doctest = false, linkcheck = true, - linkcheck_ignore = [ - "https://epubs.siam.org/doi/10.1137/0903023", - "https://link.springer.com/book/10.1007/978-3-642-05221-7", - "http://www.siam.org/journals/auth-info.php", - ], - format = Documenter.HTML( - assets = ["assets/favicon.ico"], - canonical = "https://docs.sciml.ai/ModelingToolkitCourse/stable/" - ), - pages = pages -) - +using Documenter, ModelingToolkitCourse + +pages = [ + "Home" => "index.md", + "lectures/lecture1.md", + "lectures/lecture2.md", + "lectures/lecture3.md", + "lectures/lecture4.md", + "lectures/lecture6.md", + "lectures/lecture7.md", + "lectures/lecture8.md", +] + +ENV["GKSwstype"] = "100" +using Plots + +makedocs( + sitename = "ModelingToolkit Course", + authors = "Chris Rackauckas", + modules = [ModelingToolkitCourse], + clean = true, doctest = false, linkcheck = true, + linkcheck_ignore = [ + "https://epubs.siam.org/doi/10.1137/0903023", + "https://link.springer.com/book/10.1007/978-3-642-05221-7", + "http://www.siam.org/journals/auth-info.php", + ], + format = Documenter.HTML( + assets = ["assets/favicon.ico"], + canonical = "https://docs.sciml.ai/ModelingToolkitCourse/stable/" + ), + pages = pages +) + #= using LiveServer serve(dir="build") -=# - -deploydocs( - repo = "github.com/SciML/ModelingToolkitCourse.git"; - push_preview = true -) +=# + +deploydocs( + repo = "github.com/SciML/ModelingToolkitCourse.git"; + push_preview = true +) From 8dda68d27566214145f4447365a6717349cc68a6 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Fri, 28 Aug 2026 06:01:02 -0400 Subject: [PATCH 3/3] Ignore bot-blocked MathWorks documentation link MathWorks returns HTTP 403 to Documenter's link checker for the referenced hydraulic chamber page. Exclude only that URL while retaining link checks for all other external references. Co-Authored-By: Chris Rackauckas Co-Authored-By: Codex Agent-Harness: Codex (version unknown) Agent-Model: gpt-5 Agent-Session: 01a04793-6175-7cf2-b514-9bcca858906d (local session ID) --- docs/make.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/make.jl b/docs/make.jl index 0b9314c..b685d23 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -22,6 +22,7 @@ makedocs( linkcheck_ignore = [ "https://epubs.siam.org/doi/10.1137/0903023", "https://link.springer.com/book/10.1007/978-3-642-05221-7", + "https://www.mathworks.com/help/simscape/ref/variablehydraulicchamber.html", "http://www.siam.org/journals/auth-info.php", ], format = Documenter.HTML(