diff --git a/.github/workflows/FormatCheck.yml b/.github/workflows/FormatCheck.yml new file mode 100644 index 0000000..6762c6f --- /dev/null +++ b/.github/workflows/FormatCheck.yml @@ -0,0 +1,19 @@ +name: format-check + +on: + push: + branches: + - 'master' + - 'main' + - 'release-' + tags: '*' + pull_request: + +jobs: + runic: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - uses: fredrikekre/runic-action@v1 + with: + version: '1' diff --git a/docs/make.jl b/docs/make.jl index 8c9c307..24964f4 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,35 +1,43 @@ -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 +# 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 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.jl b/docs/src/lectures/lecture1.jl index 5189668..4af3ab2 100644 --- a/docs/src/lectures/lecture1.jl +++ b/docs/src/lectures/lecture1.jl @@ -1,515 +1,465 @@ -using ForwardDiff -using Plots - -d=1 -k=1000 -Δt=1e-3 -F = 100 - -function f(xᵢ, xᵢ₋₁) - - ẋᵢ = (xᵢ - xᵢ₋₁)/Δt - lhs = d*ẋᵢ + k*xᵢ^1.5 - rhs = F - - return lhs - rhs -end - -# Newton's Method -# first time step (i=2) -xᵢ₋₁ = 0.0 -xᵢ = xᵢ₋₁ #<-- guess -g(xᵢ) = f(xᵢ, xᵢ₋₁) -xᵢ -= g(xᵢ)/ForwardDiff.derivative(g, xᵢ) -xᵢ -= g(xᵢ)/ForwardDiff.derivative(g, xᵢ) -xᵢ -= g(xᵢ)/ForwardDiff.derivative(g, xᵢ) - - -# --------------------------------- - -tol = 1e-3 -x = zeros(10) -for i=2:10 - g(xᵢ) = f(xᵢ, x[i-1]) - Δx = Inf - while abs(Δx) > tol - Δx = g(x[i])/ForwardDiff.derivative(g, x[i]) - x[i] -= Δx - end -end - -plot(x; ylabel="x [m]", xlabel="time step") - -# --------------------------------- -using DifferentialEquations - -prob = NonlinearProblem(f, 0.0, 0.0) -sol=solve(prob, NewtonRaphson(); abstol=tol) - -# --------------------------------- - -x = zeros(10) -for i=2:10 - prob′ = remake(prob; u0=x[i], p=x[i-1]) - sol = solve(prob′, NewtonRaphson(); abstol=tol) - x[i] = sol[] -end -plot(x; ylabel="x [m]", xlabel="time step") - -# --------------------------------- - -function du_dt(u,p,t) - F, k, d = p - x = u - return (F - k*x^1.5)/d -end - -prob = ODEProblem(du_dt, 0.0, (0.0, 0.01), [F, k, d]) -sol = solve(prob) -plot(sol; xlabel="time [s]", ylabel="x [m]") - -# --------------------------------- -tol = 1e-3 -d=1 -k=1000 -F = 100 - -function du_dt1(u,p,t) - F, k, d = p - x, ẋ = u - - eqs = [ - ẋ # D(x) = ẋ - (d*ẋ + k*x^1.5) - (F) # 0 = ( lhs ) - ( rhs ) - ] - - return eqs -end - -fmm = ODEFunction(du_dt1; mass_matrix=[1 0;0 0]) -prob = ODEProblem(fmm, [0.0, F/d], (0.0, 0.01), [F, k, d]) -sol = solve(prob; abstol=tol) -plot(sol; layout=2) - - -function du_dt1(du,u,p,t) - F, k, d = p - x, ẋ, ẍ = u - - du[1] = ẋ - du[2] = ẍ - du[3] = (d*ẋ + k*(x^1.5)) - (F) - -end - -fmm = ODEFunction(du_dt1; mass_matrix=[1 0 0;0 1 0;0 0 0]) -prob = ODEProblem(fmm, [0.0, F/d, 0.0], (0.0, 0.01), [F, k, d]) - -sol = solve(prob, ImplicitEuler(autodiff=false); abstol=tol) -sol = solve(prob, ImplicitEuler(); abstol=tol, dt=0.001, adaptive=false) -sol = solve(prob, Rosenbrock23(); abstol=tol) - -plot(sol; layout=3) - -using ModelingToolkitComponents: SimpleImplicitEuler -sol = solve(prob, SimpleImplicitEuler(); abstol=tol) -plot(sol; layout=3) - -function du_dt2(u,p,t) - F, k, d = p - x, ẋ, ẍ = u - - eqs = [ - ẋ # D(x) = ẋ - ẍ # D(ẋ) = ẍ - (d*ẍ + 1.5*k*(x^0.5)*ẋ) - (0) # 0 = ( lhs ) - ( rhs ) - ] - - return eqs -end - -fmm = ODEFunction(du_dt2; mass_matrix=[1 0 0;0 1 0;0 0 0]) -prob = ODEProblem(fmm, [0.0, F/d, 0.0], (0.0, 0.01), [F, k, d]) -sol = solve(prob, ImplicitEuler(); abstol=tol) - - - -sol′ = solve(prob, RadauIIA5(); abstol=tol, reltol=10.0) - -plot(sol′ ; layout=3) - -plot(sol; idxs=1, xlabel="time [s]", ylabel="x [m]") -plot!(sol′; idxs=1, xlabel="time [s]", ylabel="x [m]") - -plot(sol; idxs=2, xlabel="time [s]", ylabel="ẋ [m/s]") -plot!(sol′; idxs=2, xlabel="time [s]", ylabel="ẋ [m/s]") - -plot(0:1e-5:0.01, x->ForwardDiff.derivative(t->sol(t), x)[2]; xlabel="time [s]", ylabel="ẍ [m/s^2]") -plot!(sol′; idxs=3, xlabel="time [s]", ylabel="ẍ [m/s^2]") -# --------------------------------- - -using ModelingToolkit -using ModelingToolkit: t_nounits as t, D_nounits as D - -vars = @variables x(t)=0.0 ẋ(t)=F/d ẍ(t)=0.0 -eqs = [ - D(x) ~ ẋ - D(ẋ) ~ ẍ - d*ẋ + k*x^1.5 ~ F -] -@mtkbuild odesys = ODESystem(eqs, t, vars, []) -sys = structural_simplify(odesys) -prob = ODEProblem(sys, [], (0.0, 0.01)) -sol = solve(prob; abstol=tol) -plot(sol; idxs=ẍ, xlabel="time [s]", ylabel="ẍ [m/s^2]") - -ẍ_sol = sol[ẍ] -t_sol = sol.t - - -# --------------------------------- -@connector MechanicalPort begin - v(t) - f(t), [connect = Flow] -end - -@mtkmodel Mass begin - @parameters begin - m = 10 - end - @variables begin - v(t) = 0 - f(t) = 0 - end - @components begin - port = MechanicalPort(;v=v, f=f) - end - @equations begin - # connectors - port.v ~ v - port.f ~ f - - # physics - f ~ m*D(v) - end -end - -@mtkmodel Damper begin - @parameters begin - d = 1 - end - @variables begin - v(t) = 0 - f(t) = d*v - end - @components begin - port = MechanicalPort(;v=v, f=f) - end - @equations begin - # connectors - port.v ~ v - port.f ~ f - - # physics - f ~ d*v - end -end - -@mtkmodel System begin - @parameters begin - v - m - d - end - @components begin - mass = Mass(;v,m) - damper = Damper(;v,d) - end - @equations begin - connect(mass.port, damper.port) - end -end - -@mtkbuild sys = System(;v=100, m=5, d=3) -full_equations(sys) - -defs = ModelingToolkit.defaults(sys) -defs[sys.damper.d] - - - - - -prob1 = ODEProblem(sys, [], (0,10)) -sol1 = solve(prob1) -prob2 = remake(prob1; p=[sys.m=>3, sys.d=>4]) -sol2 = solve(prob2) -plot(sol1; idxs=sys.mass.v) -plot!(sol2; idxs=sys.mass.v) - - - - - - - - - - - - - - - - - - - - - -using DAE2AE -aesys = dae_to_ae(odesys, Δt, 2) -sys = no_simplify(aesys) -prob = ODEProblem(sys, [], (0.0, 0.01)) -sol = solve(prob, ImplicitEuler(nlsolve=NLNewton(check_div=false)); abstol=tol, dt=Δt, adaptive=false, initializealg=NoInit()) -plot(sol; idxs=ẍ, xlabel="time [s]", ylabel="ẍ [m/s^2]") - -ẍ_sol = sol[ẍ] -t_sol = sol.t - - - - -# ------------------- -# ------------------- -# ------------------- -# --- PROJECT ----- -# ------------------- -# ------------------- -# ------------------- - -Δt = 1e-3 - - -function f(Xᵢ, P) - - Xᵢ₋₁, Xᵢ₋₂ = P - - xᵢ, ẋᵢ, ẍᵢ = Xᵢ - xᵢ₋₁, ẋᵢ₋₁, ẍᵢ₋₁ = Xᵢ₋₁ - xᵢ₋₂, ẋᵢ₋₂, ẍᵢ₋₂ = Xᵢ₋₂ - - eqs = [ - ( ẋᵢ ) - ( (3*xᵢ - 4*xᵢ₋₁ + xᵢ₋₂)/(2*Δt) ) - ( ẍᵢ ) - ( (3*ẋᵢ - 4*ẋᵢ₋₁ + ẋᵢ₋₂)/(2*Δt) ) - ( d*ẋᵢ + k*xᵢ^1.5 ) - ( F ) - ] - - return eqs -end - -T = (0:Δt:Δt*99) -X = zeros(100,3) -X[1,2] = F/d - -prob = NonlinearProblem(f, X[1,:], (X[1,:],X[1,:])) -sol = solve(prob, NewtonRaphson()) - -for i=2:100 - prob′ = remake(prob; u0=X[i-1,:], p=(X[i-1,:],X[i>2 ? i-2 : i-1, :])) - sol = solve(prob′, NewtonRaphson(); abstol=tol) - X[i,:] = sol[:] -end -plot(X[:,1]; ylabel="x [m]", xlabel="time step") -plot(X[:,2]; ylabel="ẋ [m/s]", xlabel="time step") - -plot(T, X[:,3]; ylabel="ẍ [m/s²]", xlabel="time [s]") -plot!(t_sol, ẍ_sol) - - - - - - - - - - - - - - - - - - - - - - - -xo = zeros(11) -ẋo = zeros(11) -ẋo[1] = F/d -function du_dt3(u,p,t) - F, k, d, Δt = p - x, ẋ, ẍ = u - - i = round(Int, t/Δt)+1 - xo[i] = ForwardDiff.value(x) - ẋo[i] = ForwardDiff.value(ẋ) - - eqs = [ - (ẋ) - (x-xo[i>1 ? i-1 : i])/Δt # D(x) = ẋ - (ẍ) - (ẋ-ẋo[i>1 ? i-1 : i])/Δt # D(ẋ) = ẍ - (d*ẋ + k*x^1.5) - (F) # 0 = ( lhs ) - ( rhs ) - ] - - - return eqs -end - -fmm = ODEFunction(du_dt3; mass_matrix=[0 0 0;0 0 0;0 0 0]) -prob = ODEProblem(fmm, [0.0, F/d, 0.0], (0.0, 0.01), [F, k, d, Δt]) -sol′′ = solve(prob, ImplicitEuler(nlsolve=NLNewton(always_new=true, check_div=false, relax=4//10)); abstol=tol, adaptive=false, dt=Δt, initializealg=NoInit()) -# -plot(sol′′ ; idxs=3, xlabel="time [s]", ylabel="ẍ [m/s^2]") - - - -fe(x,p,t) = error_tracker(f, x, p, t) - -function step_julia(xₒ, t) - - x = [xₒ] - p = [Δt, d, k, xₒ] - t += Δt - g(x) = fe(x,p,t) - - while norm(g(x)) > tol - x = x .- ForwardDiff.jacobian(g, x)\g(x) - end - - return x[1], t -end - -function solver_julia() - xs = zeros(10) - t = 0.0 - for i=2:10 - xs[i], t = step_julia(xs[i-1], t) - end - return xs -end - -@time x_julia=solver_julia(); -errors_julia = copy(errors) - -plot(x_julia) - -plot(iteration_error(errors_julia); yscale=:log10) -hline!([tol]) - -plot(iteration_number(errors_julia)) - - - - -# DifferentialEquations.jl -using DifferentialEquations - -xₒ=0.0 -p = [Δt, d, k, xₒ] -prob = SteadyStateProblem(fe, [xₒ], p) -sol = solve(prob, NewtonRaphson(); abstol=tol) - -function step_diffeq(xₒ, t) - - prob′ = remake(prob, p=[Δt, d, k, xₒ], u0=[xₒ]) - sol = solve(prob′, NewtonRaphson(); abstol=tol) - - for (i,er) in enumerate(errors) - if isinf(er[1]) - errors[i] = (t, er[2], er[3]) - end - end - - return sol.u[] -end - -empty!(errors) -function solver_diffeq() - xs = zeros(10) - t = 0.0 - for i=2:10 - t += Δt - xs[i] = step_diffeq(xs[i-1], t) - end - return xs -end - -@time x_diffeq=solver_diffeq(); -plot(x_diffeq) -plot!(x_julia) - -errors_diffeq = copy(errors) - -plot(iteration_error(errors_julia); yscale=:log10) -plot!(iteration_error(errors_diffeq); yscale=:log10) -hline!([tol]) - -plot(iteration_number(errors_julia)) -plot!(iteration_number(errors_diffeq)) - -fode = ODEFunction(fe, mass_matrix=zeros(1,1)) -prob = ODEProblem(fode, [0.0], (0, Inf), [Δt, d, k, xₒ]) -empty!(errors) -integrator = init(prob, ImplicitEuler(nlsolve=NLNewton(always_new=false, max_iter=100)); adaptive=false, dt=Δt, abstol=tol, initializealg=NoInit()) - -function solver_ie() - xs = zeros(10) - - for i=2:10 - integrator.p[end] = xs[i-1] - step!(integrator) - xs[i] = integrator.u[1] - end - return xs -end - -x_ie = solver_ie() - -errors_ie = copy(errors) - -plot(x_diffeq) -plot!(x_julia) -plot!(x_ie) - -plot(iteration_error(errors_julia); yscale=:log10) -plot!(iteration_error(errors_diffeq); yscale=:log10) -plot!(iteration_error(errors_ie); yscale=:log10) -hline!([tol]) - -plot(iteration_number(errors_julia)) -plot!(iteration_number(errors_diffeq)) -plot!(iteration_number(errors_ie)) - - - -using ModelingToolkit -using DAE2AE -using ModelingToolkit: t_nounits as t, D_nounits as D - -vars = @variables x(t)=0 - -eqs =[ - D(x) ~ (F - k*x^1.5)/d -] - -@mtkbuild odesys = ODESystem(eqs, t, vars, []) - -aesys = DAE2AE.dae_to_ae(odesys, Δt) -sys = structural_simplify(aesys) -prob = ODEProblem(sys,[], (0, Inf), []) -integrator = init(prob, ImplicitEuler(); adaptive=false, dt=Δt, abstol=tol, initializealg=NoInit()) -step!(integrator) \ No newline at end of file +using ForwardDiff +using Plots + +d = 1 +k = 1000 +Δt = 1.0e-3 +F = 100 + +function f(xᵢ, xᵢ₋₁) + + ẋᵢ = (xᵢ - xᵢ₋₁) / Δt + lhs = d * ẋᵢ + k * xᵢ^1.5 + rhs = F + + return lhs - rhs +end + +# Newton's Method +# first time step (i=2) +xᵢ₋₁ = 0.0 +xᵢ = xᵢ₋₁ #<-- guess +g(xᵢ) = f(xᵢ, xᵢ₋₁) +xᵢ -= g(xᵢ) / ForwardDiff.derivative(g, xᵢ) +xᵢ -= g(xᵢ) / ForwardDiff.derivative(g, xᵢ) +xᵢ -= g(xᵢ) / ForwardDiff.derivative(g, xᵢ) + + +# --------------------------------- + +tol = 1.0e-3 +x = zeros(10) +for i in 2:10 + g(xᵢ) = f(xᵢ, x[i - 1]) + Δx = Inf + while abs(Δx) > tol + Δx = g(x[i]) / ForwardDiff.derivative(g, x[i]) + x[i] -= Δx + end +end + +plot(x; ylabel = "x [m]", xlabel = "time step") + +# --------------------------------- +using DifferentialEquations + +prob = NonlinearProblem(f, 0.0, 0.0) +sol = solve(prob, NewtonRaphson(); abstol = tol) + +# --------------------------------- + +x = zeros(10) +for i in 2:10 + prob′ = remake(prob; u0 = x[i], p = x[i - 1]) + sol = solve(prob′, NewtonRaphson(); abstol = tol) + x[i] = sol[] +end +plot(x; ylabel = "x [m]", xlabel = "time step") + +# --------------------------------- + +function du_dt(u, p, t) + F, k, d = p + x = u + return (F - k * x^1.5) / d +end + +prob = ODEProblem(du_dt, 0.0, (0.0, 0.01), [F, k, d]) +sol = solve(prob) +plot(sol; xlabel = "time [s]", ylabel = "x [m]") + +# --------------------------------- +tol = 1.0e-3 +d = 1 +k = 1000 +F = 100 + +function du_dt1(u, p, t) + F, k, d = p + x, ẋ = u + + eqs = [ + ẋ # D(x) = ẋ + (d * ẋ + k * x^1.5) - (F) # 0 = ( lhs ) - ( rhs ) + ] + + return eqs +end + +fmm = ODEFunction(du_dt1; mass_matrix = [1 0;0 0]) +prob = ODEProblem(fmm, [0.0, F / d], (0.0, 0.01), [F, k, d]) +sol = solve(prob; abstol = tol) +plot(sol; layout = 2) + + +function du_dt1(du, u, p, t) + F, k, d = p + x, ẋ, ẍ = u + + du[1] = ẋ + du[2] = ẍ + return du[3] = (d * ẋ + k * (x^1.5)) - (F) + +end + +fmm = ODEFunction(du_dt1; mass_matrix = [1 0 0;0 1 0;0 0 0]) +prob = ODEProblem(fmm, [0.0, F / d, 0.0], (0.0, 0.01), [F, k, d]) + +sol = solve(prob, ImplicitEuler(autodiff = false); abstol = tol) +sol = solve(prob, ImplicitEuler(); abstol = tol, dt = 0.001, adaptive = false) +sol = solve(prob, Rosenbrock23(); abstol = tol) + +plot(sol; layout = 3) + +using ModelingToolkitComponents: SimpleImplicitEuler +sol = solve(prob, SimpleImplicitEuler(); abstol = tol) +plot(sol; layout = 3) + +function du_dt2(u, p, t) + F, k, d = p + x, ẋ, ẍ = u + + eqs = [ + ẋ # D(x) = ẋ + ẍ # D(ẋ) = ẍ + (d * ẍ + 1.5 * k * (x^0.5) * ẋ) - (0) # 0 = ( lhs ) - ( rhs ) + ] + + return eqs +end + +fmm = ODEFunction(du_dt2; mass_matrix = [1 0 0;0 1 0;0 0 0]) +prob = ODEProblem(fmm, [0.0, F / d, 0.0], (0.0, 0.01), [F, k, d]) +sol = solve(prob, ImplicitEuler(); abstol = tol) + + +sol′ = solve(prob, RadauIIA5(); abstol = tol, reltol = 10.0) + +plot(sol′; layout = 3) + +plot(sol; idxs = 1, xlabel = "time [s]", ylabel = "x [m]") +plot!(sol′; idxs = 1, xlabel = "time [s]", ylabel = "x [m]") + +plot(sol; idxs = 2, xlabel = "time [s]", ylabel = "ẋ [m/s]") +plot!(sol′; idxs = 2, xlabel = "time [s]", ylabel = "ẋ [m/s]") + +plot(0:1.0e-5:0.01, x -> ForwardDiff.derivative(t -> sol(t), x)[2]; xlabel = "time [s]", ylabel = "ẍ [m/s^2]") +plot!(sol′; idxs = 3, xlabel = "time [s]", ylabel = "ẍ [m/s^2]") +# --------------------------------- + +using ModelingToolkit +using ModelingToolkit: t_nounits as t, D_nounits as D + +vars = @variables x(t) = 0.0 ẋ(t) = F / d ẍ(t) = 0.0 +eqs = [ + D(x) ~ ẋ + D(ẋ) ~ ẍ + d * ẋ + k * x^1.5 ~ F +] +@mtkbuild odesys = ODESystem(eqs, t, vars, []) +sys = structural_simplify(odesys) +prob = ODEProblem(sys, [], (0.0, 0.01)) +sol = solve(prob; abstol = tol) +plot(sol; idxs = ẍ, xlabel = "time [s]", ylabel = "ẍ [m/s^2]") + +ẍ_sol = sol[ẍ] +t_sol = sol.t + + +# --------------------------------- +@connector MechanicalPort begin + v(t) + f(t), [connect = Flow] +end + +@mtkmodel Mass begin + @parameters begin + m = 10 + end + @variables begin + v(t) = 0 + f(t) = 0 + end + @components begin + port = MechanicalPort(; v = v, f = f) + end + @equations begin + # connectors + port.v ~ v + port.f ~ f + + # physics + f ~ m * D(v) + end +end + +@mtkmodel Damper begin + @parameters begin + d = 1 + end + @variables begin + v(t) = 0 + f(t) = d * v + end + @components begin + port = MechanicalPort(; v = v, f = f) + end + @equations begin + # connectors + port.v ~ v + port.f ~ f + + # physics + f ~ d * v + end +end + +@mtkmodel System begin + @parameters begin + v + m + d + end + @components begin + mass = Mass(; v, m) + damper = Damper(; v, d) + end + @equations begin + connect(mass.port, damper.port) + end +end + +@mtkbuild sys = System(; v = 100, m = 5, d = 3) +full_equations(sys) + +defs = ModelingToolkit.defaults(sys) +defs[sys.damper.d] + + +prob1 = ODEProblem(sys, [], (0, 10)) +sol1 = solve(prob1) +prob2 = remake(prob1; p = [sys.m => 3, sys.d => 4]) +sol2 = solve(prob2) +plot(sol1; idxs = sys.mass.v) +plot!(sol2; idxs = sys.mass.v) + + +using DAE2AE +aesys = dae_to_ae(odesys, Δt, 2) +sys = no_simplify(aesys) +prob = ODEProblem(sys, [], (0.0, 0.01)) +sol = solve(prob, ImplicitEuler(nlsolve = NLNewton(check_div = false)); abstol = tol, dt = Δt, adaptive = false, initializealg = NoInit()) +plot(sol; idxs = ẍ, xlabel = "time [s]", ylabel = "ẍ [m/s^2]") + +ẍ_sol = sol[ẍ] +t_sol = sol.t + + +# ------------------- +# ------------------- +# ------------------- +# --- PROJECT ----- +# ------------------- +# ------------------- +# ------------------- + +Δt = 1.0e-3 + + +function f(Xᵢ, P) + + Xᵢ₋₁, Xᵢ₋₂ = P + + xᵢ, ẋᵢ, ẍᵢ = Xᵢ + xᵢ₋₁, ẋᵢ₋₁, ẍᵢ₋₁ = Xᵢ₋₁ + xᵢ₋₂, ẋᵢ₋₂, ẍᵢ₋₂ = Xᵢ₋₂ + + eqs = [ + (ẋᵢ) - ((3 * xᵢ - 4 * xᵢ₋₁ + xᵢ₋₂) / (2 * Δt)) + (ẍᵢ) - ((3 * ẋᵢ - 4 * ẋᵢ₋₁ + ẋᵢ₋₂) / (2 * Δt)) + (d * ẋᵢ + k * xᵢ^1.5) - (F) + ] + + return eqs +end + +T = (0:Δt:(Δt * 99)) +X = zeros(100, 3) +X[1, 2] = F / d + +prob = NonlinearProblem(f, X[1, :], (X[1, :], X[1, :])) +sol = solve(prob, NewtonRaphson()) + +for i in 2:100 + prob′ = remake(prob; u0 = X[i - 1, :], p = (X[i - 1, :], X[i > 2 ? i - 2 : i - 1, :])) + sol = solve(prob′, NewtonRaphson(); abstol = tol) + X[i, :] = sol[:] +end +plot(X[:, 1]; ylabel = "x [m]", xlabel = "time step") +plot(X[:, 2]; ylabel = "ẋ [m/s]", xlabel = "time step") + +plot(T, X[:, 3]; ylabel = "ẍ [m/s²]", xlabel = "time [s]") +plot!(t_sol, ẍ_sol) + + +xo = zeros(11) +ẋo = zeros(11) +ẋo[1] = F / d +function du_dt3(u, p, t) + F, k, d, Δt = p + x, ẋ, ẍ = u + + i = round(Int, t / Δt) + 1 + xo[i] = ForwardDiff.value(x) + ẋo[i] = ForwardDiff.value(ẋ) + + eqs = [ + (ẋ) - (x - xo[i > 1 ? i - 1 : i]) / Δt # D(x) = ẋ + (ẍ) - (ẋ - ẋo[i > 1 ? i - 1 : i]) / Δt # D(ẋ) = ẍ + (d * ẋ + k * x^1.5) - (F) # 0 = ( lhs ) - ( rhs ) + ] + + + return eqs +end + +fmm = ODEFunction(du_dt3; mass_matrix = [0 0 0;0 0 0;0 0 0]) +prob = ODEProblem(fmm, [0.0, F / d, 0.0], (0.0, 0.01), [F, k, d, Δt]) +sol′′ = solve(prob, ImplicitEuler(nlsolve = NLNewton(always_new = true, check_div = false, relax = 4 // 10)); abstol = tol, adaptive = false, dt = Δt, initializealg = NoInit()) +# +plot(sol′′; idxs = 3, xlabel = "time [s]", ylabel = "ẍ [m/s^2]") + + +fe(x, p, t) = error_tracker(f, x, p, t) + +function step_julia(xₒ, t) + + x = [xₒ] + p = [Δt, d, k, xₒ] + t += Δt + g(x) = fe(x, p, t) + + while norm(g(x)) > tol + x = x .- ForwardDiff.jacobian(g, x) \ g(x) + end + + return x[1], t +end + +function solver_julia() + xs = zeros(10) + t = 0.0 + for i in 2:10 + xs[i], t = step_julia(xs[i - 1], t) + end + return xs +end + +@time x_julia = solver_julia(); +errors_julia = copy(errors) + +plot(x_julia) + +plot(iteration_error(errors_julia); yscale = :log10) +hline!([tol]) + +plot(iteration_number(errors_julia)) + + +# DifferentialEquations.jl +using DifferentialEquations + +xₒ = 0.0 +p = [Δt, d, k, xₒ] +prob = SteadyStateProblem(fe, [xₒ], p) +sol = solve(prob, NewtonRaphson(); abstol = tol) + +function step_diffeq(xₒ, t) + + prob′ = remake(prob, p = [Δt, d, k, xₒ], u0 = [xₒ]) + sol = solve(prob′, NewtonRaphson(); abstol = tol) + + for (i, er) in enumerate(errors) + if isinf(er[1]) + errors[i] = (t, er[2], er[3]) + end + end + + return sol.u[] +end + +empty!(errors) +function solver_diffeq() + xs = zeros(10) + t = 0.0 + for i in 2:10 + t += Δt + xs[i] = step_diffeq(xs[i - 1], t) + end + return xs +end + +@time x_diffeq = solver_diffeq(); +plot(x_diffeq) +plot!(x_julia) + +errors_diffeq = copy(errors) + +plot(iteration_error(errors_julia); yscale = :log10) +plot!(iteration_error(errors_diffeq); yscale = :log10) +hline!([tol]) + +plot(iteration_number(errors_julia)) +plot!(iteration_number(errors_diffeq)) + +fode = ODEFunction(fe, mass_matrix = zeros(1, 1)) +prob = ODEProblem(fode, [0.0], (0, Inf), [Δt, d, k, xₒ]) +empty!(errors) +integrator = init(prob, ImplicitEuler(nlsolve = NLNewton(always_new = false, max_iter = 100)); adaptive = false, dt = Δt, abstol = tol, initializealg = NoInit()) + +function solver_ie() + xs = zeros(10) + + for i in 2:10 + integrator.p[end] = xs[i - 1] + step!(integrator) + xs[i] = integrator.u[1] + end + return xs +end + +x_ie = solver_ie() + +errors_ie = copy(errors) + +plot(x_diffeq) +plot!(x_julia) +plot!(x_ie) + +plot(iteration_error(errors_julia); yscale = :log10) +plot!(iteration_error(errors_diffeq); yscale = :log10) +plot!(iteration_error(errors_ie); yscale = :log10) +hline!([tol]) + +plot(iteration_number(errors_julia)) +plot!(iteration_number(errors_diffeq)) +plot!(iteration_number(errors_ie)) + + +using ModelingToolkit +using DAE2AE +using ModelingToolkit: t_nounits as t, D_nounits as D + +vars = @variables x(t) = 0 + +eqs = [ + D(x) ~ (F - k * x^1.5) / d, +] + +@mtkbuild odesys = ODESystem(eqs, t, vars, []) + +aesys = DAE2AE.dae_to_ae(odesys, Δt) +sys = structural_simplify(aesys) +prob = ODEProblem(sys, [], (0, Inf), []) +integrator = init(prob, ImplicitEuler(); adaptive = false, dt = Δt, abstol = tol, initializealg = NoInit()) +step!(integrator) diff --git a/docs/src/lectures/lecture2.jl b/docs/src/lectures/lecture2.jl index 8e0cb9e..84040e2 100644 --- a/docs/src/lectures/lecture2.jl +++ b/docs/src/lectures/lecture2.jl @@ -1,183 +1,182 @@ -using ModelingToolkit -using DifferentialEquations -using Symbolics -using Plots -using ModelingToolkit: t_nounits as t, D_nounits as D - -# parameters ------- -pars = @parameters begin - r₀ = 876 #kg/s - β = 1.2e9 #Pa - A = 0.01 #m² - x₀ = 1.0 #m - M = 10_000 #kg - g = 9.807 #m/s² - amp = 5e-2 #m - f = 15 #Hz -end - -dt = 1e-4 #s -t_end = 0.2 #s -time = 0:dt:t_end - -x_fun(t,amp,f) = amp*sin(2π*t*f) + x₀ -ẋ_fun = build_function(expand_derivatives(D(x_fun(t,amp,f))), t, amp, f; expression=false) - - -vars = @variables begin - x(t) # = x₀ - ẋ(t) - ẍ(t) - p(t) # = m*g/A #Pa - ṁ(t) - r(t) - ṙ(t) -end - -function get_base_equations(density_type) - - eqs = [ - D(x) ~ ẋ - D(ẋ) ~ ẍ - D(r) ~ ṙ - - r ~ r₀*(1 + p/β) - - ṁ ~ ṙ*x*A + (density_type)*ẋ*A - M*ẍ ~ p*A - M*g - ] - - return eqs -end - -eqs_ṁ1 = [ - get_base_equations(r₀)... - ṁ ~ ẋ_fun(t,amp,f)*A*r -] - -eqs_ṁ2 = [ - get_base_equations(r)... - ṁ ~ ẋ_fun(t,amp,f)*A*r -] - -eqs_x = [ - get_base_equations(r)... - x ~ x_fun(t,amp,f) -] - -@mtkbuild odesys_x = ODESystem(eqs_x, t, vars, pars) -prob_x = ODEProblem(sys_x, [], (0, t_end)) -sol_x = solve(prob_x; saveat=time) - -@mtkbuild odesys_ṁ1 = ODESystem(eqs_ṁ1, t, vars, pars) -u0 = [sol_x[s][1] for s in unknowns(sys_ṁ1)] -prob_ṁ1 = ODEProblem(sys_ṁ1, u0, (0, t_end)) -sol_ṁ1 = solve(prob_ṁ1) - -plot(sol_ṁ1; idxs=ṁ, label="guess", ylabel="ṁ") -plot!(sol_x; idxs=ṁ, label="solution") - -@mtkbuild odesys_ṁ2 = ODESystem(eqs_ṁ2, t, vars, pars) -prob_ṁ2 = ODEProblem(sys_ṁ2, u0, (0, t_end)) -sol_ṁ2 = solve(prob_ṁ2) - -plot(sol_x; idxs=x, label="solution", ylabel="x") -plot!(sol_ṁ1; idxs=x, label="case 1") -plot!(sol_ṁ2; idxs=x, label="case 2") - -plot(time, (sol_ṁ1(time)[x] .- sol_ṁ2(time)[x])/1e-3, ylabel="x error [mm]", xlabel="t [s]") - -# ----------------------------------------------------- - -using DataInterpolations -mass_flow_fun = LinearInterpolation(sol_x[ṁ], sol_x.t) - -open("dm.jl","w") do io - print(io, "u = [") - join(io, sol_x[ṁ], ',') - print(io, "]") -end - - -plot(sol_x; idxs=ṁ) -plot!(time, -263.9489641054512*cos.(2π*time*15)) - -import ModelingToolkitStandardLibrary.Mechanical.Translational as T -import ModelingToolkitStandardLibrary.Hydraulic.IsothermalCompressible as IC -import ModelingToolkitStandardLibrary.Blocks as B - -function MassVolume(; name, dx, drho, dm) - - pars = @parameters begin - A = 0.01 #m² - x₀ = 1.0 #m - M = 10_000 #kg - g = 9.807 #m/s² - amp = 5e-2 #m - f = 15 #Hz - p_int=M*g/A - dx=dx - drho=drho - dm=dm - end - 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_flow_input = B.TimeVaryingFunction(;f = mass_flow_fun) - end - - eqs = [ - connect(mass.flange, vol.flange) - connect(vol.port, mass_flow.port) - connect(mass_flow.dm, mass_flow_input.output) - connect(mass_flow.port, fluid) - ] - - return ODESystem(eqs, t, vars, pars; systems, name) -end - -dx = sol_x[ẋ][1] -drho = sol_x[ṙ][1] -dm = sol_x[ṁ][1] - -@named odesys = MassVolume(; dx, drho, dm) -# using DAE2AE -sys = structural_simplify(odesys) -# sys = no_simplify(expand_connections(odesys)) |> complete -prob = ODEProblem(sys, [], (0, t_end)) -sol=solve(prob) -plot(sol; idxs=sys.vol.x, linewidth=2) -plot!(sol_x; idxs=x) - -plot(sol; idxs=sys.vol.dx, linewidth=2) -plot!(sol_x; idxs=ẋ) - - - -du0 = [ - 0 # mass₊v(t) - sol_x[ẋ][1] # vol₊x(t) - sol_x[ẋ][1] # vol₊moving_volume₊x(t) - sol_x[ṙ][1] # vol₊moving_volume₊rho(t) - sol_x[ṙ*β/r₀][1] # vol₊moving_volume₊port₊p(t) - sol_x[ṙ*β/r₀][1] # vol₊damper₊port_b₊p(t) - 0 # vol₊moving_volume₊drho(t) -] -prob = DAEProblem(sys, du0, [], (0, t_end)) -sol = solve(prob, DImplicitEuler(nlsolve=NLNewton(check_div=false))) - -prob = ODEProblem(sys, [], (0, t_end)) -sol = solve(prob, ImplicitEuler(nlsolve=NLNewton(check_div=false, max_iter=100, relax=4//10)); dt, adaptive=false) - -using SimpleEuler -sol = solve(prob, BackwardEuler(); dt, adaptive=false) - -plot(sol; idxs=sys.vol.x, linewidth=2) -plot!(sol_x; idxs=x) - -plot(sol; idxs=sys.mass_flow.dm.u) - -plot(sol; idxs=sys.vol.moving_volume.port.p) +using ModelingToolkit +using DifferentialEquations +using Symbolics +using Plots +using ModelingToolkit: t_nounits as t, D_nounits as D + +# parameters ------- +pars = @parameters begin + r₀ = 876 #kg/s + β = 1.2e9 #Pa + A = 0.01 #m² + x₀ = 1.0 #m + M = 10_000 #kg + g = 9.807 #m/s² + amp = 5.0e-2 #m + f = 15 #Hz +end + +dt = 1.0e-4 #s +t_end = 0.2 #s +time = 0:dt:t_end + +x_fun(t, amp, f) = amp * sin(2π * t * f) + x₀ +ẋ_fun = build_function(expand_derivatives(D(x_fun(t, amp, f))), t, amp, f; expression = false) + + +vars = @variables begin + x(t) # = x₀ + ẋ(t) + ẍ(t) + p(t) # = m*g/A #Pa + ṁ(t) + r(t) + ṙ(t) +end + +function get_base_equations(density_type) + + eqs = [ + D(x) ~ ẋ + D(ẋ) ~ ẍ + D(r) ~ ṙ + + r ~ r₀ * (1 + p / β) + + ṁ ~ ṙ * x * A + (density_type) * ẋ * A + M * ẍ ~ p * A - M * g + ] + + return eqs +end + +eqs_ṁ1 = [ + get_base_equations(r₀)... + ṁ ~ ẋ_fun(t, amp, f) * A * r +] + +eqs_ṁ2 = [ + get_base_equations(r)... + ṁ ~ ẋ_fun(t, amp, f) * A * r +] + +eqs_x = [ + get_base_equations(r)... + x ~ x_fun(t, amp, f) +] + +@mtkbuild odesys_x = ODESystem(eqs_x, t, vars, pars) +prob_x = ODEProblem(sys_x, [], (0, t_end)) +sol_x = solve(prob_x; saveat = time) + +@mtkbuild odesys_ṁ1 = ODESystem(eqs_ṁ1, t, vars, pars) +u0 = [sol_x[s][1] for s in unknowns(sys_ṁ1)] +prob_ṁ1 = ODEProblem(sys_ṁ1, u0, (0, t_end)) +sol_ṁ1 = solve(prob_ṁ1) + +plot(sol_ṁ1; idxs = ṁ, label = "guess", ylabel = "ṁ") +plot!(sol_x; idxs = ṁ, label = "solution") + +@mtkbuild odesys_ṁ2 = ODESystem(eqs_ṁ2, t, vars, pars) +prob_ṁ2 = ODEProblem(sys_ṁ2, u0, (0, t_end)) +sol_ṁ2 = solve(prob_ṁ2) + +plot(sol_x; idxs = x, label = "solution", ylabel = "x") +plot!(sol_ṁ1; idxs = x, label = "case 1") +plot!(sol_ṁ2; idxs = x, label = "case 2") + +plot(time, (sol_ṁ1(time)[x] .- sol_ṁ2(time)[x]) / 1.0e-3, ylabel = "x error [mm]", xlabel = "t [s]") + +# ----------------------------------------------------- + +using DataInterpolations +mass_flow_fun = LinearInterpolation(sol_x[ṁ], sol_x.t) + +open("dm.jl", "w") do io + print(io, "u = [") + join(io, sol_x[ṁ], ',') + print(io, "]") +end + + +plot(sol_x; idxs = ṁ) +plot!(time, -263.9489641054512 * cos.(2π * time * 15)) + +import ModelingToolkitStandardLibrary.Mechanical.Translational as T +import ModelingToolkitStandardLibrary.Hydraulic.IsothermalCompressible as IC +import ModelingToolkitStandardLibrary.Blocks as B + +function MassVolume(; name, dx, drho, dm) + + pars = @parameters begin + A = 0.01 #m² + x₀ = 1.0 #m + M = 10_000 #kg + g = 9.807 #m/s² + amp = 5.0e-2 #m + f = 15 #Hz + p_int = M * g / A + dx = dx + drho = drho + dm = dm + end + 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_flow_input = B.TimeVaryingFunction(; f = mass_flow_fun) + end + + eqs = [ + connect(mass.flange, vol.flange) + connect(vol.port, mass_flow.port) + connect(mass_flow.dm, mass_flow_input.output) + connect(mass_flow.port, fluid) + ] + + return ODESystem(eqs, t, vars, pars; systems, name) +end + +dx = sol_x[ẋ][1] +drho = sol_x[ṙ][1] +dm = sol_x[ṁ][1] + +@named odesys = MassVolume(; dx, drho, dm) +# using DAE2AE +sys = structural_simplify(odesys) +# sys = no_simplify(expand_connections(odesys)) |> complete +prob = ODEProblem(sys, [], (0, t_end)) +sol = solve(prob) +plot(sol; idxs = sys.vol.x, linewidth = 2) +plot!(sol_x; idxs = x) + +plot(sol; idxs = sys.vol.dx, linewidth = 2) +plot!(sol_x; idxs = ẋ) + + +du0 = [ + 0 # mass₊v(t) + sol_x[ẋ][1] # vol₊x(t) + sol_x[ẋ][1] # vol₊moving_volume₊x(t) + sol_x[ṙ][1] # vol₊moving_volume₊rho(t) + sol_x[ṙ * β / r₀][1] # vol₊moving_volume₊port₊p(t) + sol_x[ṙ * β / r₀][1] # vol₊damper₊port_b₊p(t) + 0 # vol₊moving_volume₊drho(t) +] +prob = DAEProblem(sys, du0, [], (0, t_end)) +sol = solve(prob, DImplicitEuler(nlsolve = NLNewton(check_div = false))) + +prob = ODEProblem(sys, [], (0, t_end)) +sol = solve(prob, ImplicitEuler(nlsolve = NLNewton(check_div = false, max_iter = 100, relax = 4 // 10)); dt, adaptive = false) + +using SimpleEuler +sol = solve(prob, BackwardEuler(); dt, adaptive = false) + +plot(sol; idxs = sys.vol.x, linewidth = 2) +plot!(sol_x; idxs = x) + +plot(sol; idxs = sys.mass_flow.dm.u) + +plot(sol; idxs = sys.vol.moving_volume.port.p) diff --git a/docs/src/lectures/run_away_train_hydraulic.jl b/docs/src/lectures/run_away_train_hydraulic.jl index 9f6ef89..ef1a7ce 100644 --- a/docs/src/lectures/run_away_train_hydraulic.jl +++ b/docs/src/lectures/run_away_train_hydraulic.jl @@ -1,311 +1,322 @@ -using ModelingToolkit -using DifferentialEquations -using ModelingToolkitStandardLibrary.Mechanical.Translational: MechanicalPort -using ModelingToolkitStandardLibrary.Hydraulic.IsothermalCompressible: Valve, DynamicVolume, HydraulicFluid, FixedPressure -using ModelingToolkitStandardLibrary.Blocks: Constant -using Plots -using ModelingToolkit: t_nounits as t, D_nounits as D - -# NOTE: How Decouple works to provide discontinuous behavior -@mtkmodel Decouple begin #Decouple(;k, d, v_a, v_b, x_a, x_b, f=0) | port_a, port_b - @parameters begin - k - d - end - @variables begin - v_a(t) - v_b(t) - x_a(t) - x_b(t) - f(t)=0 - end - @components begin - port_a = MechanicalPort(;v=v_a,f=+f) - port_b = MechanicalPort(;v=v_b,f=-f) - end - @equations begin - # differentials - D(x_a) ~ v_a - D(x_b) ~ v_b - - # connectors - port_a.v ~ v_a - port_b.v ~ v_b - port_a.f ~ +f - port_b.f ~ -f - - # physics - f ~ ifelse(x_a >= x_b, (v_a - v_b)*d + k*(x_a - x_b), 0) - end -end - -@mtkmodel Mass begin #Mass(;m, v, f, a=f/m) | port - @parameters begin - m = 10 - end - @variables begin - v(t) - f(t) = 0 - a(t) = f/m - end - @components begin - port = MechanicalPort(; v, f) - end - @equations begin - # derivatives - D(v) ~ a - - # connectors - port.v ~ v - port.f ~ f - - # physics - f ~ m*D(v) - end -end - -@mtkmodel Damper begin #Damper(;d=1,v,f) | port_a, port_b - @parameters begin - d = 1 - end - @variables begin - v(t) - f(t) - end - @components begin - port_a = MechanicalPort() - port_b = MechanicalPort() - end - @equations begin - # connectors - (port_a.v - port_b.v) ~ v - port_a.f ~ +f - port_b.f ~ -f - - # physics - f ~ d*v - end -end - -@mtkmodel Spring begin #Spring(;k=100,x,v,f) | port_a, port_b - @parameters begin - k = 100 - end - @variables begin - x(t) - v(t) - f(t) - end - @components begin - port_a = MechanicalPort() - port_b = MechanicalPort() - end - @equations begin - # derivatives - D(x) ~ v - - # connectors - (port_a.v - port_b.v) ~ v - port_a.f ~ +f - port_b.f ~ -f - - # physics - f ~ k*x - end -end - -@mtkmodel Reference begin #Reference(;f) | port - @parameters begin - - end - @variables begin - f(t) - end - @components begin - port = MechanicalPort() - end - @equations begin - # connectors - port.v ~ 0 - port.f ~ -f - end -end - -# NOTE: How TrainCar works to track absolute position -@mtkmodel TrainCar begin #TrainCar(;m, v, x_a, x_b) | port_a, port_b - @parameters begin - m - v - end - @variables begin - x_a(t) - x_b(t) - end - @components begin - port_a = MechanicalPort(;v) - port_b = MechanicalPort(;v) - mass = Mass(;m,v) - end - @equations begin - D(x_a) ~ port_a.v - D(x_b) ~ port_b.v - - connect(mass.port, port_a, port_b) - end -end - -@mtkmodel Coupler begin #Coupler(;k,d,v) | port_a, port_b - @parameters begin - k - d - v - end - @variables begin - - end - @components begin - port_a = MechanicalPort(;v) - port_b = MechanicalPort(;v) - spring = Spring(;k,x=0,v=0,f=0) - damper = Damper(;d,v=0,f=0) - end - @equations begin - connect(port_a, spring.port_a, damper.port_a) - connect(spring.port_b, damper.port_b, port_b) - end -end - -@mtkmodel Stopper begin #Stopper(;k,d) | port - @parameters begin - k - d - end - @variables begin - - end - @components begin - port = MechanicalPort() - spring = Spring(;k,x=0,v=0,f=0) - damper = Damper(;d,v=0,f=0) - ref = Reference() - end - @equations begin - connect(port, spring.port_a, damper.port_a) - connect(spring.port_b, damper.port_b, ref.port) - end -end - -# NOTE: How HydraulicStopper works -@component function HydraulicStopper(;name, area) #HydraulicStopper(;area) | port - pars = @parameters begin - area = area - end - vars = @variables begin - - end - systems = @named begin - port = MechanicalPort() - volume = DynamicVolume(5, true, false; p_int=0, area=0.1, x_int = 5, x_max = 5, x_min = 1, x_damp = 2, direction = -1) - valve = Valve(;p_a_int=0, p_b_int=0, area_int=area, Cd=0.7) - constarea = Constant(;k=area) - open = FixedPressure(;p=0) - fluid = HydraulicFluid(; density=876, bulk_modulus=1.2e9) - end - eqs = [ - connect(port, volume.flange) - connect(volume.port, valve.port_a) - connect(valve.port_b, open.port) - connect(constarea.output, valve.area) - - connect(fluid, volume.port) - ] - - return ODESystem(eqs, t, vars, pars; systems, name) -end - - -@mtkmodel System begin - @parameters begin - v=10 - end - @variables begin - - end - @components begin - car1 = TrainCar(;m=1000, v, x_a=0, x_b=0.9) - cx1 = Coupler(;k=1e5,d=1e5,v) - - car2 = TrainCar(;m=1000, v, x_a=1.1, x_b=1.9) - cx2 = Coupler(;k=1e5,d=1e5,v) - - car3 = TrainCar(;m=1000, v, x_a=2.1, x_b=2.9) - cx3 = Coupler(;k=1e5,d=1e5,v) - - engine = TrainCar(;m=2000, v=10, x_a=3.1, x_b=3.9) - decouple = Decouple(;k=1e6, d=1e6, v_a=v, v_b=0, x_a=3.9, x_b=10, f=0) - stopper = HydraulicStopper(; area=1e-2) - end - @equations begin - connect(car1.port_b, cx1.port_a) - connect(car2.port_a, cx1.port_b) - connect(car2.port_b, cx2.port_a) - connect(car3.port_a, cx2.port_b) - connect(car3.port_b, cx3.port_a) - connect(engine.port_a, cx3.port_b) - connect(engine.port_b, decouple.port_a) - connect(decouple.port_b, stopper.port) - end -end - -@named odesys = System() -using JuliaSimCompiler -sys = structural_simplify(IRSystem(odesys)) -@mtkbuild sys = System() - -# NOTE: missing_variable_defaults() -prob = ODEProblem(sys, ModelingToolkit.missing_variable_defaults(sys), (0, 5)) -prob = ODEProblem(sys, [], (0, 5)) -# NOTE: strategies for challenging solve -@time "solve Rodas" solve(prob, Rodas4()); -@time "solve" sol = solve(prob, ImplicitEuler(nlsolve=NLNewton(check_div=false, always_new=true, relax=4/10, max_iter=100)); adaptive=false, dt=1e-5); - - - -# NOTE: pressure wave and no negative pressure -plot(sol; - idxs=[sys.stopper.volume.v1.port.p/1e5, - sys.stopper.volume.v2.port.p/1e5, - sys.stopper.volume.v3.port.p/1e5, - sys.stopper.volume.v4.port.p/1e5, - sys.stopper.volume.v5.port.p/1e5], - xlims=(0.72, 0.74), - ylabel="pressure [bar]", - xlabel="time [s]") - -# NOTE: pressure bump at the end when the valve closes (1m x_min) -plot(sol; - idxs=[sys.stopper.volume.v1.port.p/1e5, - sys.stopper.volume.v2.port.p/1e5, - sys.stopper.volume.v3.port.p/1e5, - sys.stopper.volume.v4.port.p/1e5, - sys.stopper.volume.v5.port.p/1e5], - ylims=(0, 20), - ylabel="pressure [bar]", - xlabel="time [s]") - -plot(sol; - idxs=[sys.engine.mass.v]) - -# NOTE: How to design safe acceleration limit? -g = 9.807 -plot(sol; - idxs=[ - sys.car1.mass.a/g, - sys.car2.mass.a/g, - sys.car3.mass.a/g, - sys.engine.mass.a/g - ], - ylims=(-25,5), - xlims=(0.7,0.85), - ylabel="acceleration [g]", - xlabel="time [s]") +using ModelingToolkit +using DifferentialEquations +using ModelingToolkitStandardLibrary.Mechanical.Translational: MechanicalPort +using ModelingToolkitStandardLibrary.Hydraulic.IsothermalCompressible: Valve, DynamicVolume, HydraulicFluid, FixedPressure +using ModelingToolkitStandardLibrary.Blocks: Constant +using Plots +using ModelingToolkit: t_nounits as t, D_nounits as D + +# NOTE: How Decouple works to provide discontinuous behavior +@mtkmodel Decouple begin #Decouple(;k, d, v_a, v_b, x_a, x_b, f=0) | port_a, port_b + @parameters begin + k + d + end + @variables begin + v_a(t) + v_b(t) + x_a(t) + x_b(t) + f(t) = 0 + end + @components begin + port_a = MechanicalPort(; v = v_a, f = +f) + port_b = MechanicalPort(; v = v_b, f = -f) + end + @equations begin + # differentials + D(x_a) ~ v_a + D(x_b) ~ v_b + + # connectors + port_a.v ~ v_a + port_b.v ~ v_b + port_a.f ~ +f + port_b.f ~ -f + + # physics + f ~ ifelse(x_a >= x_b, (v_a - v_b) * d + k * (x_a - x_b), 0) + end +end + +@mtkmodel Mass begin #Mass(;m, v, f, a=f/m) | port + @parameters begin + m = 10 + end + @variables begin + v(t) + f(t) = 0 + a(t) = f / m + end + @components begin + port = MechanicalPort(; v, f) + end + @equations begin + # derivatives + D(v) ~ a + + # connectors + port.v ~ v + port.f ~ f + + # physics + f ~ m * D(v) + end +end + +@mtkmodel Damper begin #Damper(;d=1,v,f) | port_a, port_b + @parameters begin + d = 1 + end + @variables begin + v(t) + f(t) + end + @components begin + port_a = MechanicalPort() + port_b = MechanicalPort() + end + @equations begin + # connectors + (port_a.v - port_b.v) ~ v + port_a.f ~ +f + port_b.f ~ -f + + # physics + f ~ d * v + end +end + +@mtkmodel Spring begin #Spring(;k=100,x,v,f) | port_a, port_b + @parameters begin + k = 100 + end + @variables begin + x(t) + v(t) + f(t) + end + @components begin + port_a = MechanicalPort() + port_b = MechanicalPort() + end + @equations begin + # derivatives + D(x) ~ v + + # connectors + (port_a.v - port_b.v) ~ v + port_a.f ~ +f + port_b.f ~ -f + + # physics + f ~ k * x + end +end + +@mtkmodel Reference begin #Reference(;f) | port + @parameters begin + + end + @variables begin + f(t) + end + @components begin + port = MechanicalPort() + end + @equations begin + # connectors + port.v ~ 0 + port.f ~ -f + end +end + +# NOTE: How TrainCar works to track absolute position +@mtkmodel TrainCar begin #TrainCar(;m, v, x_a, x_b) | port_a, port_b + @parameters begin + m + v + end + @variables begin + x_a(t) + x_b(t) + end + @components begin + port_a = MechanicalPort(; v) + port_b = MechanicalPort(; v) + mass = Mass(; m, v) + end + @equations begin + D(x_a) ~ port_a.v + D(x_b) ~ port_b.v + + connect(mass.port, port_a, port_b) + end +end + +@mtkmodel Coupler begin #Coupler(;k,d,v) | port_a, port_b + @parameters begin + k + d + v + end + @variables begin + + end + @components begin + port_a = MechanicalPort(; v) + port_b = MechanicalPort(; v) + spring = Spring(; k, x = 0, v = 0, f = 0) + damper = Damper(; d, v = 0, f = 0) + end + @equations begin + connect(port_a, spring.port_a, damper.port_a) + connect(spring.port_b, damper.port_b, port_b) + end +end + +@mtkmodel Stopper begin #Stopper(;k,d) | port + @parameters begin + k + d + end + @variables begin + + end + @components begin + port = MechanicalPort() + spring = Spring(; k, x = 0, v = 0, f = 0) + damper = Damper(; d, v = 0, f = 0) + ref = Reference() + end + @equations begin + connect(port, spring.port_a, damper.port_a) + connect(spring.port_b, damper.port_b, ref.port) + end +end + +# NOTE: How HydraulicStopper works +@component function HydraulicStopper(; name, area) #HydraulicStopper(;area) | port + pars = @parameters begin + area = area + end + vars = @variables begin + + end + systems = @named begin + port = MechanicalPort() + volume = DynamicVolume(5, true, false; p_int = 0, area = 0.1, x_int = 5, x_max = 5, x_min = 1, x_damp = 2, direction = -1) + valve = Valve(; p_a_int = 0, p_b_int = 0, area_int = area, Cd = 0.7) + constarea = Constant(; k = area) + open = FixedPressure(; p = 0) + fluid = HydraulicFluid(; density = 876, bulk_modulus = 1.2e9) + end + eqs = [ + connect(port, volume.flange) + connect(volume.port, valve.port_a) + connect(valve.port_b, open.port) + connect(constarea.output, valve.area) + + connect(fluid, volume.port) + ] + + return ODESystem(eqs, t, vars, pars; systems, name) +end + + +@mtkmodel System begin + @parameters begin + v = 10 + end + @variables begin + + end + @components begin + car1 = TrainCar(; m = 1000, v, x_a = 0, x_b = 0.9) + cx1 = Coupler(; k = 1.0e5, d = 1.0e5, v) + + car2 = TrainCar(; m = 1000, v, x_a = 1.1, x_b = 1.9) + cx2 = Coupler(; k = 1.0e5, d = 1.0e5, v) + + car3 = TrainCar(; m = 1000, v, x_a = 2.1, x_b = 2.9) + cx3 = Coupler(; k = 1.0e5, d = 1.0e5, v) + + engine = TrainCar(; m = 2000, v = 10, x_a = 3.1, x_b = 3.9) + decouple = Decouple(; k = 1.0e6, d = 1.0e6, v_a = v, v_b = 0, x_a = 3.9, x_b = 10, f = 0) + stopper = HydraulicStopper(; area = 1.0e-2) + end + @equations begin + connect(car1.port_b, cx1.port_a) + connect(car2.port_a, cx1.port_b) + connect(car2.port_b, cx2.port_a) + connect(car3.port_a, cx2.port_b) + connect(car3.port_b, cx3.port_a) + connect(engine.port_a, cx3.port_b) + connect(engine.port_b, decouple.port_a) + connect(decouple.port_b, stopper.port) + end +end + +@named odesys = System() +using JuliaSimCompiler +sys = structural_simplify(IRSystem(odesys)) +@mtkbuild sys = System() + +# NOTE: missing_variable_defaults() +prob = ODEProblem(sys, ModelingToolkit.missing_variable_defaults(sys), (0, 5)) +prob = ODEProblem(sys, [], (0, 5)) +# NOTE: strategies for challenging solve +@time "solve Rodas" solve(prob, Rodas4()); +@time "solve" sol = solve(prob, ImplicitEuler(nlsolve = NLNewton(check_div = false, always_new = true, relax = 4 / 10, max_iter = 100)); adaptive = false, dt = 1.0e-5); + + +# NOTE: pressure wave and no negative pressure +plot( + sol; + idxs = [ + sys.stopper.volume.v1.port.p / 1.0e5, + sys.stopper.volume.v2.port.p / 1.0e5, + sys.stopper.volume.v3.port.p / 1.0e5, + sys.stopper.volume.v4.port.p / 1.0e5, + sys.stopper.volume.v5.port.p / 1.0e5, + ], + xlims = (0.72, 0.74), + ylabel = "pressure [bar]", + xlabel = "time [s]" +) + +# NOTE: pressure bump at the end when the valve closes (1m x_min) +plot( + sol; + idxs = [ + sys.stopper.volume.v1.port.p / 1.0e5, + sys.stopper.volume.v2.port.p / 1.0e5, + sys.stopper.volume.v3.port.p / 1.0e5, + sys.stopper.volume.v4.port.p / 1.0e5, + sys.stopper.volume.v5.port.p / 1.0e5, + ], + ylims = (0, 20), + ylabel = "pressure [bar]", + xlabel = "time [s]" +) + +plot( + sol; + idxs = [sys.engine.mass.v] +) + +# NOTE: How to design safe acceleration limit? +g = 9.807 +plot( + sol; + idxs = [ + sys.car1.mass.a / g, + sys.car2.mass.a / g, + sys.car3.mass.a / g, + sys.engine.mass.a / g, + ], + ylims = (-25, 5), + xlims = (0.7, 0.85), + ylabel = "acceleration [g]", + xlabel = "time [s]" +) diff --git a/docs/src/lectures/run_away_train_utils.jl b/docs/src/lectures/run_away_train_utils.jl index ca88278..1933246 100644 --- a/docs/src/lectures/run_away_train_utils.jl +++ b/docs/src/lectures/run_away_train_utils.jl @@ -1,167 +1,180 @@ -using CairoMakie -using Makie.GeometryBasics - -function plot_train(sol, sys, tx=0.0) - tm = Observable(tx) - idx = Dict(reverse.(enumerate(unknowns(sys)))) - - fig = Figure() - a = Axis(fig[1,1], aspect=DataAspect(), ) - hidedecorations!(a) - slt(t,idxs) = sol(t; idxs) # =[sys.car1.x_a, sys.car1.x_b, sys.car2.x_a, sys.car2.x_b, sys.engine.x_a, sys.engine.x_b, sys.stopper.spring.x+sys.engine.barrier]) - - c1xa(t) = slt(t, sys.car1.x_a) - c1xb(t) = slt(t, sys.car1.x_b) - - flb1(t) = (c1xa(t), 0) - frb1(t) = (c1xb(t), 0) - frt1(t) = (c1xb(t), 0.25) - delta1(t) = c1xb(t)-c1xa(t) - flt1a(t) = (c1xb(t)-delta1(t)*0.75, 0.25) - flt1b(t) = (c1xa(t), 0.1) - - ply1(t) = Polygon(Point2f[flb1(t), frb1(t), frt1(t), flt1a(t), flt1b(t)]) - poly!(a, lift(ply1, tm), color=:blue) - - c2xa(t) = slt(t, sys.car2.x_a) - c2xb(t) = slt(t, sys.car2.x_b) - - flb2(t) = (c2xa(t), 0) - frb2(t) = (c2xb(t), 0) - frt2(t) = (c2xb(t), 0.25) - flt2(t) = (c2xa(t), 0.25) - - ply2(t) = Polygon(Point2f[flb2(t), frb2(t), frt2(t), flt2(t)]) - poly!(a, lift(ply2, tm), color=:blue) - - - c3xa(t) = slt(t, sys.car3.x_a) - c3xb(t) = slt(t, sys.car3.x_b) - - flb3(t) = (c3xa(t), 0) - frb3(t) = (c3xb(t), 0) - frt3(t) = (c3xb(t), 0.25) - flt3(t) = (c3xa(t), 0.25) - - ply3(t) = Polygon(Point2f[flb3(t), frb3(t), frt3(t), flt3(t)]) - poly!(a, lift(ply3, tm), color=:blue) - - - - exa(t) = slt(t, sys.engine.x_a) - exb(t) = slt(t, sys.engine.x_b) - - flbe(t) = (exa(t), 0) - frbe(t) = (exb(t), 0) - deltae(t) = exb(t)-exa(t) - frtea(t) = (exb(t), 0.1) - frteb(t) = (exa(t)+deltae(t)*0.75, 0.25) - flte(t) = (exa(t), 0.25) - - plye(t) = Polygon(Point2f[flbe(t), frbe(t), frtea(t), frteb(t), flte(t)]) - poly!(a, lift(plye, tm), color=:blue) - - - br(t) = slt(t, (5-sys.stopper.volume.x)+10) - brp(t) = Polygon(Point2f[ - (br(t),0) - (br(t)+0.1,0) - (br(t)+0.1,0.1) - (br(t)+0.5,0.1) - (br(t)+0.5,0.2) - (br(t)+0.1,0.2) - (br(t)+0.1,0.3) - (br(t),0.3) - ]) - - poly!(a, lift(brp, tm), color=(:red, 0.5)) - - scaling(t,p) = (10+log((slt(t,p)/1e5+0.01)/100))/10 - - v5(t) = slt(t, 10-sys.stopper.volume.v5.x) - v5c(t) = (:red, scaling(t, sys.stopper.volume.v5.port.p) ) - v5p(t) = Polygon(Point2f[ - (v5(t),0) - (11,0) - (11,0.3) - (v5(t),0.3) - ]) - poly!(a, lift(v5p, tm); color=lift(v5c, tm)) #, colormap = :jet, colorrange = (1, 20)) - - - v4(t) = slt(t, 10+1-sys.stopper.volume.v4.x) - v4c(t) = (:red, scaling(t, sys.stopper.volume.v4.port.p) ) - v4p(t) = Polygon(Point2f[ - (v4(t),0) - (12,0) - (12,0.3) - (v4(t),0.3) - ]) - poly!(a, lift(v4p, tm); color=lift(v4c, tm)) #, colormap = :jet, colorrange = (1, 20)) - - - v3(t) = slt(t, 10+2-sys.stopper.volume.v3.x) - v3c(t) = (:red, scaling(t, sys.stopper.volume.v3.port.p) ) - v3p(t) = Polygon(Point2f[ - (v3(t),0) - (13,0) - (13,0.3) - (v3(t),0.3) - ]) - poly!(a, lift(v3p, tm); color=lift(v3c, tm)) #, colormap = :jet, colorrange = (1, 20)) - - v2(t) = slt(t, 10+3-sys.stopper.volume.v2.x) - v2c(t) = (:red, scaling(t, sys.stopper.volume.v2.port.p) ) - v2p(t) = Polygon(Point2f[ - (v2(t),0) - (14,0) - (14,0.3) - (v2(t),0.3) - ]) - poly!(a, lift(v2p, tm); color=lift(v2c, tm)) #, colormap = :jet, colorrange = (1, 20)) - - v1(t) = slt(t, 10+4-sys.stopper.volume.v1.x) - v1c(t) = (:red, scaling(t, sys.stopper.volume.v1.port.p) ) - v1p(t) = Polygon(Point2f[ - (v1(t),0) - (15,0) - (15,0.3) - (v1(t),0.3) - ]) - poly!(a, lift(v1p, tm); color=lift(v1c, tm)) #, colormap = :jet, colorrange = (1, 20)) - - # track - lines!(a, [0,15], [-0.05,-0.05]; linewidth=1, color=:black) - - - CairoMakie.ylims!(a, -5, 5) - CairoMakie.xlims!(a, 0, 15) - fig, tm -end - -# fig,tm = plot_train(sol, sys); fig -# fig, = plot_train(sol, sys, 0.76); fig - -function record_train(fig, tm; start, stop, step, framerate, filename) - - timestamps = range(start, stop; step) - - record(fig, filename, timestamps; - framerate) do t - tm[] = t - end - - nothing -end - -fig, tm = plot_train(sol, sys) -record_train(fig, tm; start=0.72, stop=0.729, step=1e-4, framerate=20, filename="run_away_train_slow.mp4") -record_train(fig, tm; start=0.0, stop=5, step=1/30, framerate=30, filename="run_away_train.mp4") - -# plot(sol; idxs=sys.stopper.volume.moving_volume.port.p/1e5, xlims=(0.7, 0.9)) -# plot!(sol; idxs=sys.stopper.volume.v1.port.p/1e5, xlims=(0.7, 0.9)) -# plot!(sol; idxs=sys.stopper.volume.v5.port.p/1e5, xlims=(0.7, 0.9)) - -# plot(sol; idxs=sys.stopper.volume.moving_volume.dx, xlims=(0.7, 0.9)) -# plot!(sol; idxs=sys.stopper.volume.v1.dx, xlims=(0.7, 0.9)) -# plot!(sol; idxs=sys.stopper.volume.v5.dx, xlims=(0.7, 0.9)) \ No newline at end of file +using CairoMakie +using Makie.GeometryBasics + +function plot_train(sol, sys, tx = 0.0) + tm = Observable(tx) + idx = Dict(reverse.(enumerate(unknowns(sys)))) + + fig = Figure() + a = Axis(fig[1, 1], aspect = DataAspect()) + hidedecorations!(a) + slt(t, idxs) = sol(t; idxs) # =[sys.car1.x_a, sys.car1.x_b, sys.car2.x_a, sys.car2.x_b, sys.engine.x_a, sys.engine.x_b, sys.stopper.spring.x+sys.engine.barrier]) + + c1xa(t) = slt(t, sys.car1.x_a) + c1xb(t) = slt(t, sys.car1.x_b) + + flb1(t) = (c1xa(t), 0) + frb1(t) = (c1xb(t), 0) + frt1(t) = (c1xb(t), 0.25) + delta1(t) = c1xb(t) - c1xa(t) + flt1a(t) = (c1xb(t) - delta1(t) * 0.75, 0.25) + flt1b(t) = (c1xa(t), 0.1) + + ply1(t) = Polygon(Point2f[flb1(t), frb1(t), frt1(t), flt1a(t), flt1b(t)]) + poly!(a, lift(ply1, tm), color = :blue) + + c2xa(t) = slt(t, sys.car2.x_a) + c2xb(t) = slt(t, sys.car2.x_b) + + flb2(t) = (c2xa(t), 0) + frb2(t) = (c2xb(t), 0) + frt2(t) = (c2xb(t), 0.25) + flt2(t) = (c2xa(t), 0.25) + + ply2(t) = Polygon(Point2f[flb2(t), frb2(t), frt2(t), flt2(t)]) + poly!(a, lift(ply2, tm), color = :blue) + + + c3xa(t) = slt(t, sys.car3.x_a) + c3xb(t) = slt(t, sys.car3.x_b) + + flb3(t) = (c3xa(t), 0) + frb3(t) = (c3xb(t), 0) + frt3(t) = (c3xb(t), 0.25) + flt3(t) = (c3xa(t), 0.25) + + ply3(t) = Polygon(Point2f[flb3(t), frb3(t), frt3(t), flt3(t)]) + poly!(a, lift(ply3, tm), color = :blue) + + + exa(t) = slt(t, sys.engine.x_a) + exb(t) = slt(t, sys.engine.x_b) + + flbe(t) = (exa(t), 0) + frbe(t) = (exb(t), 0) + deltae(t) = exb(t) - exa(t) + frtea(t) = (exb(t), 0.1) + frteb(t) = (exa(t) + deltae(t) * 0.75, 0.25) + flte(t) = (exa(t), 0.25) + + plye(t) = Polygon(Point2f[flbe(t), frbe(t), frtea(t), frteb(t), flte(t)]) + poly!(a, lift(plye, tm), color = :blue) + + + br(t) = slt(t, (5 - sys.stopper.volume.x) + 10) + brp(t) = Polygon( + Point2f[ + (br(t), 0) + (br(t) + 0.1, 0) + (br(t) + 0.1, 0.1) + (br(t) + 0.5, 0.1) + (br(t) + 0.5, 0.2) + (br(t) + 0.1, 0.2) + (br(t) + 0.1, 0.3) + (br(t), 0.3) + ] + ) + + poly!(a, lift(brp, tm), color = (:red, 0.5)) + + scaling(t, p) = (10 + log((slt(t, p) / 1.0e5 + 0.01) / 100)) / 10 + + v5(t) = slt(t, 10 - sys.stopper.volume.v5.x) + v5c(t) = (:red, scaling(t, sys.stopper.volume.v5.port.p)) + v5p(t) = Polygon( + Point2f[ + (v5(t), 0) + (11, 0) + (11, 0.3) + (v5(t), 0.3) + ] + ) + poly!(a, lift(v5p, tm); color = lift(v5c, tm)) #, colormap = :jet, colorrange = (1, 20)) + + + v4(t) = slt(t, 10 + 1 - sys.stopper.volume.v4.x) + v4c(t) = (:red, scaling(t, sys.stopper.volume.v4.port.p)) + v4p(t) = Polygon( + Point2f[ + (v4(t), 0) + (12, 0) + (12, 0.3) + (v4(t), 0.3) + ] + ) + poly!(a, lift(v4p, tm); color = lift(v4c, tm)) #, colormap = :jet, colorrange = (1, 20)) + + + v3(t) = slt(t, 10 + 2 - sys.stopper.volume.v3.x) + v3c(t) = (:red, scaling(t, sys.stopper.volume.v3.port.p)) + v3p(t) = Polygon( + Point2f[ + (v3(t), 0) + (13, 0) + (13, 0.3) + (v3(t), 0.3) + ] + ) + poly!(a, lift(v3p, tm); color = lift(v3c, tm)) #, colormap = :jet, colorrange = (1, 20)) + + v2(t) = slt(t, 10 + 3 - sys.stopper.volume.v2.x) + v2c(t) = (:red, scaling(t, sys.stopper.volume.v2.port.p)) + v2p(t) = Polygon( + Point2f[ + (v2(t), 0) + (14, 0) + (14, 0.3) + (v2(t), 0.3) + ] + ) + poly!(a, lift(v2p, tm); color = lift(v2c, tm)) #, colormap = :jet, colorrange = (1, 20)) + + v1(t) = slt(t, 10 + 4 - sys.stopper.volume.v1.x) + v1c(t) = (:red, scaling(t, sys.stopper.volume.v1.port.p)) + v1p(t) = Polygon( + Point2f[ + (v1(t), 0) + (15, 0) + (15, 0.3) + (v1(t), 0.3) + ] + ) + poly!(a, lift(v1p, tm); color = lift(v1c, tm)) #, colormap = :jet, colorrange = (1, 20)) + + # track + lines!(a, [0, 15], [-0.05, -0.05]; linewidth = 1, color = :black) + + + CairoMakie.ylims!(a, -5, 5) + CairoMakie.xlims!(a, 0, 15) + return fig, tm +end + +# fig,tm = plot_train(sol, sys); fig +# fig, = plot_train(sol, sys, 0.76); fig + +function record_train(fig, tm; start, stop, step, framerate, filename) + + timestamps = range(start, stop; step) + + record( + fig, filename, timestamps; + framerate + ) do t + tm[] = t + end + + return nothing +end + +fig, tm = plot_train(sol, sys) +record_train(fig, tm; start = 0.72, stop = 0.729, step = 1.0e-4, framerate = 20, filename = "run_away_train_slow.mp4") +record_train(fig, tm; start = 0.0, stop = 5, step = 1 / 30, framerate = 30, filename = "run_away_train.mp4") + +# plot(sol; idxs=sys.stopper.volume.moving_volume.port.p/1e5, xlims=(0.7, 0.9)) +# plot!(sol; idxs=sys.stopper.volume.v1.port.p/1e5, xlims=(0.7, 0.9)) +# plot!(sol; idxs=sys.stopper.volume.v5.port.p/1e5, xlims=(0.7, 0.9)) + +# plot(sol; idxs=sys.stopper.volume.moving_volume.dx, xlims=(0.7, 0.9)) +# plot!(sol; idxs=sys.stopper.volume.v1.dx, xlims=(0.7, 0.9)) +# plot!(sol; idxs=sys.stopper.volume.v5.dx, xlims=(0.7, 0.9)) diff --git a/docs/src/lectures/volume.jl b/docs/src/lectures/volume.jl index 0195acd..e5d13e8 100644 --- a/docs/src/lectures/volume.jl +++ b/docs/src/lectures/volume.jl @@ -1,56 +1,57 @@ -using ModelingToolkit -using ModelingToolkitStandardLibrary.Mechanical.Translational -using ModelingToolkitStandardLibrary.Hydraulic.IsothermalCompressible - -using ModelingToolkitStandardLibrary.Hydraulic.IsothermalCompressible: liquid_density - - - -@component function Volume(; - #initial conditions - x, - dx=0, - p, - drho=0, - dm=0, - - #parameters - area, - direction=+1, name) - pars = @parameters begin - area = area - end - - vars = @variables begin - x(t) = x - dx(t) = dx - p(t) = p - f(t) = p * area - rho(t) - drho(t) = drho - dm(t) = dm - end - - systems = @named begin - port = HydraulicPort(; p_int=p) - flange = MechanicalPort(; f, v=dx) - end - - eqs = [ - # connectors - port.p ~ p - port.dm ~ dm - flange.v * direction ~ dx - flange.f * direction ~ -f - - # differentials - D(x) ~ dx - D(rho) ~ drho - - # physics - rho ~ liquid_density(port, p) - f ~ p * area - dm ~ drho * x * area + rho * dx * area] - - ODESystem(eqs, t, vars, pars; name, systems, defaults=[rho => liquid_density(port)]) -end \ No newline at end of file +using ModelingToolkit +using ModelingToolkitStandardLibrary.Mechanical.Translational +using ModelingToolkitStandardLibrary.Hydraulic.IsothermalCompressible + +using ModelingToolkitStandardLibrary.Hydraulic.IsothermalCompressible: liquid_density + + +@component function Volume(; + #initial conditions + x, + dx = 0, + p, + drho = 0, + dm = 0, + + #parameters + area, + direction = +1, name + ) + pars = @parameters begin + area = area + end + + vars = @variables begin + x(t) = x + dx(t) = dx + p(t) = p + f(t) = p * area + rho(t) + drho(t) = drho + dm(t) = dm + end + + systems = @named begin + port = HydraulicPort(; p_int = p) + flange = MechanicalPort(; f, v = dx) + end + + eqs = [ + # connectors + port.p ~ p + port.dm ~ dm + flange.v * direction ~ dx + flange.f * direction ~ -f + + # differentials + D(x) ~ dx + D(rho) ~ drho + + # physics + rho ~ liquid_density(port, p) + f ~ p * area + dm ~ drho * x * area + rho * dx * area + ] + + ODESystem(eqs, t, vars, pars; name, systems, defaults = [rho => liquid_density(port)]) +end diff --git a/test/explicit_imports.jl b/test/explicit_imports.jl index 324fa4e..d495bc1 100644 --- a/test/explicit_imports.jl +++ b/test/explicit_imports.jl @@ -1,8 +1,8 @@ -using ExplicitImports -using ModelingToolkitCourse -using Test - -@testset "ExplicitImports" begin - @test check_no_implicit_imports(ModelingToolkitCourse) === nothing - @test check_no_stale_explicit_imports(ModelingToolkitCourse) === nothing -end +using ExplicitImports +using ModelingToolkitCourse +using Test + +@testset "ExplicitImports" begin + @test check_no_implicit_imports(ModelingToolkitCourse) === nothing + @test check_no_stale_explicit_imports(ModelingToolkitCourse) === nothing +end