From 2ac51db69d40c6373dd3375b4b0e101fc8ae4f7f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Hans=20W=C3=BCrfel?= Date: Thu, 5 Dec 2024 11:27:56 +0100 Subject: [PATCH] bump NetworkDynamics to 0.9 --- docs/Project.toml | 6 +- docs/examples/truss.jl | 185 +++++++++++++++++++---------------------- 2 files changed, 87 insertions(+), 104 deletions(-) diff --git a/docs/Project.toml b/docs/Project.toml index ab3ba248..e5238f75 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -15,7 +15,7 @@ MLJ = "add582a8-e3ab-11e8-2d5e-e98b27df1bc7" Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" NetworkDynamics = "22e9dc34-2a0d-11e9-0de0-8588d035468b" NetworkLayout = "46757867-2c16-5918-afeb-47bfcb05e46a" -OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a" RegistryInstances = "2792f1a3-b283-48e8-9a74-f99dce5104f3" StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" WGLMakie = "276b4fcb-3e11-5398-bf8b-a0c2d153d008" @@ -33,9 +33,9 @@ LayeredLayouts = "0.2" Literate = "2" MLJ = "0.20" Makie = "0.21" -NetworkDynamics = "0.8.3" +NetworkDynamics = "0.9" NetworkLayout = "0.4.7" -OrdinaryDiffEq = "6" +OrdinaryDiffEqTsit5 = "1" RegistryInstances = "0.1" StableRNGs = "1" WGLMakie = "0.10" diff --git a/docs/examples/truss.jl b/docs/examples/truss.jl index 9abc8e4e..ca27683f 100644 --- a/docs/examples/truss.jl +++ b/docs/examples/truss.jl @@ -6,10 +6,10 @@ using `GaphMakie.jl` and [`NetworkDynamics.jl`](https://github.com/PIK-ICoN/Netw ![truss animation](truss.mp4) =# using NetworkDynamics -using OrdinaryDiffEq +using OrdinaryDiffEqTsit5 using Graphs using GraphMakie -using LinearAlgebra: norm, ⋅ +using LinearAlgebra: norm using Printf using CairoMakie CairoMakie.activate!() @@ -17,27 +17,39 @@ CairoMakie.activate!() #= ## Definition of the dynamical system Define dynamic model for NetworkDynamics +For more information on the models see the corresponding [example in the NetworkDynamics.jl docs](https://juliadynamics.github.io/NetworkDynamics.jl/stable/generated/stress_on_truss/). =# -function edge_f!(e, v_s, v_d, (K, L), _) - Δd = view(v_d, 1:2) - view(v_s, 1:2) - nd = norm(Δd) - F = K * ( L - nd ) - e[1:2] = F .* view(Δd, 1:2) ./ nd -end -function vertex_free!(du, u, edges, (M, γ, g), _) - F = sum(edges) - γ * view(u, 3:4) + [0, -M*g] - du[3:4] .= F ./ M - du[1:2] = view(u, 3:4) +function fixed_g(pos, x, p, t) + pos .= p end - -function vertex_fixed!(du, u, edges, _, _) - du[1:2] .= 0.0 +vertex_fix = VertexModel(g=fixed_g, psym=[:xfix, :yfix], outsym=[:x, :y], ff=NoFeedForward()) +function free_f(dx, x, Fsum, (M, γ, g), t) + v = view(x, 1:2) + dx[1:2] .= (Fsum .- γ .* v) ./ M + dx[2] -= g + dx[3:4] .= v + nothing end - -edge = StaticEdge(f = edge_f!, dim=2, coupling=:antisymmetric) -vertex_free = ODEVertex(f = vertex_free!, dim=4, sym=[:x, :y, :v, :w]) -vertex_fix = ODEVertex(f = vertex_fixed!, dim=2, sym=[:x, :y]) +vertex_free = VertexModel(f=free_f, g=3:4, sym=[:vx=>0, :vy=>0, :x, :y], + psym=[:M=>10, :γ=>200, :g=>9.81], insym=[:Fx, :Fy]) +function edge_g!(F, pos_src, pos_dst, (K, L), t) + dx = pos_dst[1] - pos_src[1] + dy = pos_dst[2] - pos_src[2] + d = sqrt(dx^2 + dy^2) + Fabs = K * (L - d) + F[1] = Fabs * dx / d + F[2] = Fabs * dy / d + nothing +end +function observedf(obsout, u, pos_src, pos_dst, (K, L), t) + dx = pos_dst[1] .- pos_src[1] + dy = pos_dst[2] .- pos_src[2] + d = sqrt(dx^2 + dy^2) + obsout[1] = K * (L - d) + nothing +end +beam = EdgeModel(g=AntiSymmetric(edge_g!), psym=[:K=>0.5e6, :L], outsym=[:Fx, :Fy], obsf=observedf, obssym=[:Fabs]) nothing #hide #= @@ -58,123 +70,94 @@ pos0 = zeros(Point2f, 2N + 1) pos0[1:N] = [Point((i-1)dx,0) for i in 1:N] pos0[N+1:2*N] = [Point(i*dx + shift, 1) for i in 1:N] pos0[2N+1] = Point(N*dx + 1, -1) - fixed = [1,4] # set fixed vertices nothing #hide -# create networkdynamics object -verts = ODEVertex[vertex_free for i in 1:nv(g)] -for i in fixed - verts[i] = vertex_fix # use the fixed vertex for the fixed -end -nd = network_dynamics(verts, edge, g); -nothing #hide - #= -write some auxiliary functions to map between state vector of DGL `u` and graph data +Now we can create: +- the `Network` object (rhs of the differential equation), +- the initial state `u0` and +- the parameters for the individual components. =# -x_idxs = [findfirst(isequal(Symbol("x_$i")), nd.syms) for i in 1:nv(g)] - -"Set positions `pos` inside dgl state `u`" -function set_pos!(u, pos) - for (i, idx) in enumerate(x_idxs) - u[idx] = pos[i][1] - u[idx+1] = pos[i][2] - end +verts = VertexModel[vertex_free for i in 1:nv(g)] +for i in fixed + verts[i] = vertex_fix # use the fixed vertex for the fixed points end - -"Extract vector of Points from dgl state `u`" -function to_pos(u) - pos = Vector{Point2f}(undef, nv(g)) - for (i, idx) in enumerate(x_idxs) - pos[i] = Point(u[idx], u[idx+1]) +nw = Network(g, verts, beam) +u0 = NWState(nw) +## set x/y initial conditions and xfix/yfix parameters +for i in eachindex(pos0, verts) + if i in fixed + u0.p.v[i, :xfix] = pos0[i][1] + u0.p.v[i, :yfix] = pos0[i][2] + else + u0.v[i, :x] = pos0[i][1] + u0.v[i, :y] = pos0[i][2] end - return pos end - -"Calculate load on edge for given state." -function get_load(u, p, t=0.0) - gd_nd = nd(u, p, t, GetGD) # exposes underlying graph data struct - force = Vector{Float64}(undef, ne(g)) - pos = to_pos(u) - for (i,e) in enumerate(edges(g)) - edgeval = get_edge(gd_nd, i) - fvec = Point(edgeval[1], edgeval[2]) - dir = pos[dst(e)] .- pos[src(e)] - force[i] = sign(fvec ⋅ dir) * norm(fvec) - end - return force +## set L for edges +for (i,e) in enumerate(edges(g)) + u0.p.e[i, :L] = norm(pos0[src(e)] - pos0[dst(e)]) end +## set damping and mass for "big mass" at the end +u0.p.v[11, :M] = 200 +u0.p.v[11, :γ] = 100 nothing #hide #= -Set parameters. +With rhs, parameters and initial conditions constructed we can integrate the system. =# -M = [10 for i in 1:nv(g)] # mass of the nodes -M[end] = 200 # heavy mass -gc = [9.81 for i in 1:nv(g)] # gravitational constant -γ = [200.0 for i in 1:nv(g)] # damping parameter -γ[end] = 100.0 - -L = [norm(pos0[src(e)] - pos0[dst(e)]) for e in edges(g)] # length of edges -K = [0.5e6 for i in 1:ne(g)] # spring constant of edges - -## bundle parameters for NetworkDynamics -para = collect.((zip(M, γ, gc), - zip(K, L))) -nothing #hide - -#= -Set initial state and solve the system -=# -u0 = zeros(length(nd.syms)) -set_pos!(u0, pos0) - tspan = (0.0, 12.0) -prob = ODEProblem(nd, u0, tspan, para) -sol = solve(prob, Tsit5()); +prob = ODEProblem(nw, uflat(u0), tspan, pflat(u0)) +sol = solve(prob, Tsit5()) nothing #hide #= ## Plot the solution =# -fig = Figure(size=(1000,550)) +fig = Figure(size=(1000,550)); fig[1,1] = title = Label(fig, "Stress on truss", fontsize=30) title.tellwidth = false fig[2,1] = ax = Axis(fig) +ax.aspect = DataAspect(); +hidespines!(ax); # no borders +hidedecorations!(ax); # no grid, axis, ... +limits!(ax, -0.1, pos0[end][1]+0.3, pos0[end][2]-0.5, 1.15) # axis limits to show full plot + +## get the maximum force during the simulation to get the color scale +## It is only possible to access `:Fabs` directly becaus we've define the observable function for it! +(fmin, fmax) = 0.3 .* extrema(Iterators.flatten(sol(sol.t, idxs=eidxs(nw, :, :Fabs)))) -## calculate some values for colorscaling -(fmin, fmax) = 0.3 .* extrema([(map(u->get_load(u, para), sol)...)...]) -nlabels = [" " for i in 1:nv(g)] -nlabels[fixed] .= "Δ" -elabels = ["edge $i" for i in 1:ne(g)] p = graphplot!(ax, g; - edge_width=4.0, - node_size=sqrt.(M)*3, - nlabels=nlabels, - nlabels_align=(:center,:top), - nlabels_fontsize=30, - elabels=elabels, - elabels_side=Dict(ne(g) => :right), - edge_color=[0.0 for i in 1:ne(g)], - edge_attr=(colorrange=(fmin,fmax), + edge_width = 4.0, + node_size = 3*sqrt.(try u0.p.v[i, :M] catch; 10.0 end for i in 1:nv(g)), + nlabels = [i in fixed ? "Δ" : "" for i in 1:nv(g)], + nlabels_align = (:center,:top), + nlabels_fontsize = 30, + elabels = ["edge $i" for i in 1:ne(g)], + elabels_side = Dict(ne(g) => :right), + edge_color = [0.0 for i in 1:ne(g)], + edge_attr = (colorrange=(fmin,fmax), colormap=:diverging_bkr_55_10_c35_n256)) -hidespines!(ax); hidedecorations!(ax); p[:node_pos][]=to_pos(u0); ax.aspect = DataAspect() -limits!(ax, -0.1, pos0[end][1]+0.3, pos0[end][2]-0.5, 1.15) ## draw colorbar -fig[3,1] = cb = Colorbar(fig, p.plots[1].plots[1], label = "Axial force", vertical=false) +fig[3,1] = cb = Colorbar(fig, get_edge_plot(p), label = "Axial force", vertical=false) T = tspan[2] fps = 30 trange = range(0.0, sol.t[end], length=Int(T * fps)) record(fig, "truss.mp4", trange; framerate=fps) do t title.text = @sprintf "Stress on truss (t = %.2f )" t - p[:node_pos][] = to_pos(sol(t)) - load = get_load(sol(t), para) - p.elabels = [@sprintf("%.0f", l) for l in load] + s_at_t = NWState(sol, t) + for i in eachindex(pos0) + p[:node_pos][][i] = (s_at_t.v[i, :x], s_at_t.v[i, :y]) + end + notify(p[:node_pos]) + load = s_at_t.e[:, :Fabs] p.edge_color[] = load + p.elabels = [@sprintf("%.0f", l) for l in load] + fig end nothing #hide