Merge pull request #200 from MakieOrg/hw/bumpnd

bump NetworkDynamics to 0.9

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"

docs/examples/truss.jl

@@ -6,38 +6,50 @@ 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!()
 
 #=
 ## 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