From f359235c28571ff86bad4457a0822605ae66cdca Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Hans=20W=C3=BCrfel?= Date: Thu, 17 Oct 2024 10:25:26 +0200 Subject: [PATCH 1/7] add hydrogen example --- docs/Project.toml | 5 + docs/examples/gas_network.jl | 379 +++++++++++++++++++++++++++++++++++ docs/make.jl | 3 +- 3 files changed, 386 insertions(+), 1 deletion(-) create mode 100644 docs/examples/gas_network.jl diff --git a/docs/Project.toml b/docs/Project.toml index fc4c563b..18ccf228 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,11 +1,15 @@ +name = "ND.jl docs" + [deps] CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" GraphMakie = "1ecd5474-83a3-4783-bb4f-06765db800d2" +DynamicQuantities = "06fc5a27-2a28-4c7c-a15d-362465fb6821" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" +LinearInterpolations = "b20c7882-2490-4592-a606-fbbfe9e745e8" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" NetworkDynamics = "22e9dc34-2a0d-11e9-0de0-8588d035468b" @@ -20,6 +24,7 @@ SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" SimpleWeightedGraphs = "47aef6b3-ad0c-573a-a1e2-d07658019622" StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0" StyledStrings = "f489334b-da3d-4c2e-b8f0-e476e12c162b" SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" diff --git a/docs/examples/gas_network.jl b/docs/examples/gas_network.jl new file mode 100644 index 00000000..350bba90 --- /dev/null +++ b/docs/examples/gas_network.jl @@ -0,0 +1,379 @@ +#= +# Dynamic Flow in simple Gas Network + +This Example is based on the paper + + Albertus J. Malan, Lukas Rausche, Felix Strehle, Sören Hohmann, + Port-Hamiltonian Modelling for Analysis and Control of Gas Networks, + IFAC-PapersOnLine, Volume 56, Issue 2, 2023, https://doi.org/10.1016/j.ifacol.2023.10.193. + +and tries replicate a simple simulation of flow in a 3-node gas network. + +This example can be dowloaded as a normal Julia script [here](@__NAME__.jl). #md + +We start by importing the necessary packages: +=# +using NetworkDynamics +using ModelingToolkit +using DynamicQuantities +using ModelingToolkit: D as Dt, t as t +using Test +using StaticArrays +using LinearInterpolations +using OrdinaryDiffEqTsit5 +using CairoMakie +CairoMakie.activate!(type="svg") #hide + +#= +## Node Models + +There are 2 node models used in the paper. The first node type has a constant +pressure. We model this as an first order ODE with zero rhs. + +In this example, we use the equation based modeling using `ModelingToolkit.jl`. +To verify the equations on a basic level we also provide units to eveything to +perform dimensionality checks. +=# +@mtkmodel ConstantPressureNode begin + @parameters begin + p_set, [description="Constant pressure setpoint", unit=u"Pa"] + end + @variables begin + p(t) = p_set, [description="Pressure", unit=u"Pa", output=true] + q̃_nw(t), [description="aggregated flow from pipes into node", unit=u"m^3/s", input=true] + q̃_inj(t), [description="internal help variable", unit=u"m^3/s"] + end + @equations begin + Dt(p) ~ 0 + q̃_inj ~ -q̃_nw + end +end +nothing #hide + +#= +The second node model is a variable pressure node. It has one output state, the pressure and one input state, +the aggregated flows from the connected pipes. +As an internal state we have the injected flow from our source/load. +The source/load behaviour itself is provided via a time dependent function. +=# +@mtkmodel VariablePressureNode begin + @structural_parameters begin + load_profile # time dependent load profile + end + @parameters begin + C, [description="Lumped capacitance of connected pipes", unit=u"m^4 * s^2 / kg"] + end + @variables begin + p(t)=5e6, [description="Pressure", unit=u"Pa", output=true] + q̃_inj(t), [description="external injection into node", unit=u"m^3/s"] + q̃_nw(t), [description="aggregated flow from pipes into node", unit=u"m^3/s", input=true] + end + @equations begin + q̃_inj ~ load_profile(t) + C * Dt(p) ~ q̃_inj + q̃_nw # (30) + end +end +nothing #hide + +#= +## Pipe Model + +The pipe is modeld as a first order ODE for the volumetric flow at the `dst` end. It has two inputs: +the pressure at the source and and the pressure at the destination end. +Later on, we'll specify the model to be antisymmetric, thus the flow is calculated explicitly for the +destination end, but the source end will just recive just that times (-1). +=# +@mtkmodel Pipe begin + @parameters begin + L, [description="Length of pipe", unit=u"m"] + A, [description="Cross-sectional area of pipe", unit=u"m^2"] + ρ̃, [description="constant standard density", unit=u"kg/m^3"] + D, [description="Diameter of pipe", unit=u"m"] + g, [description="Gravitational acceleration", unit=u"m/s^2"] + sinθ, [description="Angle of inclination" ] + c, [description="Speed of sound in fluid", unit=u"m/s"] + γ, [description="Friction efficiency factor"] + η, [description="Dynamic viscosity", unit=u"kg/(m*s)"] + r, [description="Pipe roughness", unit=u"m"] + end + @variables begin + p_src(t), [description="Pressure at source end", unit=us"Pa", input=true] + p_dst(t), [description="Pressure at destination end", unit=us"Pa", input=true] + q̃(t)=1, [description="Flow through pipe", unit=u"m^3/s", output=true] + Re(t), [description="Reynolds number"] + λ(t), [description="Friction factor"] + λe(t), [description="Effective friction factor"] + pM(t), [description="mean pressure", unit=us"Pa"] + end + @equations begin + Re ~ (ρ̃ * abs(q̃) * D) / (η * A) # (6) + λ ~ ifelse(Re < 2300, + 64/Re, # laminar (7) + (2*log10(4.518/Re * log10(Re/7) + r/(3.71*D)))^(-2) # turbulent (8) + ) + λe ~ λ/γ^2 # (10) + pM ~ 2/3*(p_src + p_dst - (p_src*p_dst)/(p_src + p_dst)) # (20) + + ρ̃ * L / A * Dt(q̃) ~ -(λe * ρ̃^2 * c^2 * L * abs(q̃))/(2 * D * A^2 * pM) * q̃ - (g * L * sinθ)/(c^2) * pM + (p_src - p_dst) # (31) + end +end +nothing #hide + +#= +## Parametrization + +The parameterization turned out to be a bit tricky. There might be errors in there. + +Some of them are quite cleare and explicitly given. +=# +g = 9.81u"m/s^2" # that we just know +Rs = 518.28u"J/(kg*K)" # Specific gas constant for hydrogen +η = 1e-5u"kg/(m*s)" # Dynamic viscosity +pc = 46.5u"bar" # Critical pressure +p̃ = 1.01325u"bar" # standard pressure +Tc = 190.55u"K" # critical temperature +T̃ = 273.15u"K" # standard temperature +γ = 0.98 # friction efficiency factor +r = 0.012u"mm" # pipe roughness +D = 0.6u"m" # pipe diameter + +L₁₂ = 80u"km" +L₁₃ = 90u"km" +L₂₃ = 100u"km" +Δh₁ = 0u"km" # this value is different for different sims in the paper +p₁_set = 50u"bar" +nothing # hide + +#= +Some neede parameters for the equations are easily derived. +=# +A = π/4 * D^2 +sinθ₁₂ = ustrip(Δh₁ / L₁₂) +sinθ₁₃ = ustrip(Δh₁ / L₁₃) +sinθ₂₃ = 0.0 +nothing # hide + +#= +Other parameters are more complicated. The compressibility factor Z according to (5) +depends on temperature and pressure. While the temperature is not considered in the model +at all, I am not sure whether we sould include the pressure dependency. +Here, i just calculated it using the "standard" pressure and "standard" temperature. +=# +Z̃ = 1 - 3.52 * p̃/pc * exp(-2.26*(T̃/Tc)) + 0.274 * (p̃/pc)^2 * exp(-1.878*(T̃/Tc)) # (5) +nothing # hide + +#= +Similarily, we can use (4) to derive the speed of sound. This again is +temperature and compressibility dependent. Here I use the "standard" compressibility again +to calculate a constant "standard" sound speed. +Equivalently, we use the "standard" compressibility to calculate the standard density. +=# +c = sqrt(T̃ * Rs * Z̃) # (4) but in theory this is also time dependend +ρ̃ = p̃ / (Rs * T̃ * Z̃) # (4) +nothing #hide + +#= +The equivalent "pressure capacity" at the nodes is calculated as a sum of the connected +pipe parameters according ot (28). +=# +C₂ = L₁₂*A/(2*ρ̃*c^2) + L₂₃*A/(2*ρ̃*c^2) # (28) +C₃ = L₁₃*A/(2*ρ̃*c^2) + L₂₃*A/(2*ρ̃*c^2) # (28) +nothing #hide + +#= +## Load Profile + +The paper specifies the load profile at two nodes. We use the package `LinearInterpolations` +to get a callable object which represents this picewise linear interpolation. + +However, this function is not Symbolics.jl compatible, so we need to stop Symbolics.jl/ModelingToolkit.jl +from tracing it. To do so, we use `@register_symbolic` to declare it as a symbolic function which is treated +as a blackbox. + +Additionally, we need to tell ModelingToolkit about the units of this object. This is just used +for the static unit check during construction of the model. Later one, when we generate the +Julia code from the symbolic reepresentation all units will be stripped. +=# +load2(t) = -Interpolate(SA[0, 4, 12, 20, 24]*3600, SA[20, 30, 10, 30, 20], extrapolate=LinearInterpolations.Constant(20))(t) +load3(t) = -Interpolate(SA[0, 4, 12, 20, 24]*3600, SA[40, 50, 30, 50, 40], extrapolate=LinearInterpolations.Constant(40))(t) +@register_symbolic load2(t) +@register_symbolic load3(t) +ModelingToolkit.get_unit(op::typeof(load2), _) = u"m^3/s" +ModelingToolkit.get_unit(op::typeof(load3), _) = u"m^3/s" +nothing #hide + +#= +## Building the Network + +To bild the Network we first need to define the components. This is a two step process: + +- first create the symbolic `ODESystem` using ModelingToolkit +- secondly build a NetworkDynamics component function (`ODEVertex`/`ODEEdge`) based on the symbolic system. + +In the first step we can use the keyword arguments to pass "default" values for our parameters and states. +Those values will be automaticially transfered to the metadata of the component function the second step. + +The second step requires to define the interface variables, i.e. what are the "input" states of your +component function and what are the "output" states. +For `ODEVertex` the input state is the aggregated flow of all connected pipes. The output state is the pressure of the node. +=# +@named v1_mtk = ConstantPressureNode(p_set=p₁_set) +v1 = ODEVertex(v1_mtk, [:q̃_nw], [:p]; name=:v1, vidx=1) +# + +@named v2_mtk = VariablePressureNode(C=C₂, load_profile=load2) +v2 = ODEVertex(v2_mtk, [:q̃_nw], [:p]; name=:v2, vidx=2) +# + +@named v3_mtk = VariablePressureNode(C=C₃, load_profile=load3) +v3 = ODEVertex(v3_mtk, [:q̃_nw], [:p]; name=:v3, vidx=3) + +#= +For `ODEEdge` we have two inputs: the pressure on both source and destination end. +There is a single output state: the volumetric flow. However we also need to tell +NetworkDynamics about the coupling type. In this case we use `AntiSymmetric`, which +meas that the source end will recieve the same flow, just inverted sign. +=# +@named e12_mtk = Pipe(; L=L₁₂, sinθ=sinθ₁₂, A, ρ̃, D, g, c, γ, η, r) +@named e13_mtk = Pipe(; L=L₁₃, sinθ=sinθ₁₃, A, ρ̃, D, g, c, γ, η, r) +@named e23_mtk = Pipe(; L=L₂₃, sinθ=sinθ₂₃, A, ρ̃, D, g, c, γ, η, r) + +e12 = ODEEdge(e12_mtk, [:p_src], [:p_dst], [:q̃], AntiSymmetric(); name=:e12, src=1, dst=2) +e13 = ODEEdge(e13_mtk, [:p_src], [:p_dst], [:q̃], AntiSymmetric(); name=:e13, src=1, dst=3) +e23 = ODEEdge(e23_mtk, [:p_src], [:p_dst], [:q̃], AntiSymmetric(); name=:e23, src=2, dst=3) + +#= +To build the network object we just need to pass the vertices and edges to the constructor. + +Note that we've used the `vidx` and `src`/`dst` keywords in the constructors to define for +each component to which "part" of the network it belongs. + +This means, the constructor can automaticially construct a graph based on those informations and +we don't need to pass it explicitly. +=# + +nw = Network([v1, v2, v3], [e12, e13, e23]) + +#= +As a result, we recive a network with 3 unique types (v2 and v3 are similar but +structurally different, because both functions capure a unique loadprofile function). + +## Finding a Steady State + +To simulate the systme, we first need to find a steadystate. As a "guess" for that +we create a `NWState` object from the network. +This will allocate flat arrays for states `u` and parameters `p` and fill them with the +default values. +=# +uguess = NWState(nw) + +#= +This is not a steadystate of the system however. To find a true steadystate we want to ensure +that the lhs of the system is zero. +We can solve for a steady state numerically by defining a Nonlinear Rootfind problem. + +To do so, we need to wrap the Network object in a closure. This closure is important to +overwrite all changes of the solver to `u[1]`, to ensure that it is not altered while +solving. +=# + +nwwrap = (du, u, p) -> begin + u[1] = uflat(uguess)[1] + nw(du, u, p, 0) + nothing +end +initprob = NonlinearProblem(nwwrap, uflat(uguess), pflat(uguess)) +initsol = solve(initprob) + +#= +The solver complains a a bit about rank deficient matrices but is able to find a solution. + +We can create a new `NWState` object by wrapping the solution from the nonlinear problem and the +original prameters in a new `NWState` object. +=# +u0 = NWState(nw, initsol.u, uguess.p) + +#= +## Solving the ODE + +Using this as our initial state we can create the actual `ODEProblem`. +Since the ode allways operates on flat state and aprameter arrays we use `uflat` and `pflat` to extract them. +=# +prob = ODEProblem(nw, uflat(u0), (0.0,24*3600), copy(pflat(u0))) +sol = solve(prob, Tsit5(), tstops=[0,4,12,20,24]*3600) +nothing #hide + +#= +## Inspect the Solution + +Inspecting the solution is all which is left to do. +=# +xticks = ((0:4:24)*3600, string.(0:4:24)) # its nice to display hours +fig = begin + _fig = Figure() + row = 1 + ax = Axis(_fig[row, 1]; xlabel="time [h]", ylabel="pressure [Pa]", title="Pressure at nodes", xticks) + xlims!(ax, sol.t[begin], sol.t[end]) + for i in 1:3 + lines!(ax, sol, idxs=vidxs(nw, i, :p); label="v$i", color=Cycled(i)) + end + axislegend(ax) + row += 1 + + ax = Axis(_fig[row, 1]; xlabel="time [h]", ylabel="flow [m³/s]", title="Flow through pipes", xticks) + xlims!(ax, sol.t[begin], sol.t[end]) + for i in 1:2 + lines!(ax, sol, idxs=eidxs(nw, i, :q̃); label="e$i flow", color=Cycled(i)) + end + axislegend(ax) + row += 1 + _fig +end + + +#= +Notably, the "internal" states defined in the symbolic models are not "states" in the sense of the ODE. +For example, we captured the load profile in the `q̃_inj` state of the `VariablePressureNode`. +The only dynamic state of the model however is `p`. +Using the "observables" mechanism from SciML, which is implemented by NetworkDynamics, we can reconstruct +those "optimized" states which have been removed symbolicially. +Here we plot the reconstructed load profile of nodes 2 and 3. Also, we know that node 1 is infinetly stiff, +acting as an infinite source of volumetric flow. We can reconstruct this flow too. +=# +fig = begin + _fig = Figure() + row = 1 + ax = Axis(_fig[row, 1]; xlabel="time [h]", ylabel="flow [m³/s]", title="Flow at nodes", xticks) + xlims!(ax, sol.t[begin], sol.t[end]) + lines!(ax, sol, idxs=vidxs(nw, 1, :q̃_inj); label="v1 compensation", color=Cycled(1)) + for i in 2:3 + lines!(ax, sol, idxs=vidxs(nw, i, :q̃_inj); label="v$i load profile", color=Cycled(i)) + end + axislegend(ax, position=:rc) + _fig +end + +#= +Lastly we want to observe two internal states of the pipes: the Reynolds number and the mean pressure. +We see, that we're purely in the turbulent flow regime. +=# +fig = begin + _fig = Figure() + row = 1 + ax = Axis(_fig[row, 1]; xlabel="time [h]", ylabel="Reynolds number", title="Reynolds number", xticks) + xlims!(ax, sol.t[begin], sol.t[end]) + for i in 1:3 + lines!(ax, sol, idxs=eidxs(nw, i, :Re); label="e $i", color=Cycled(i)) + end + hlines!(ax, 2300, color=:black, linestyle=:dash, label="L/T transition") + axislegend(ax, position=:rb) + row += 1 + + ax = Axis(_fig[row, 1]; xlabel="time [h]", ylabel="Mean pressure [Pa]", title="Mean pressure in pipes", xticks) + xlims!(ax, sol.t[begin], sol.t[end]) + for i in 1:3 + lines!(ax, sol, idxs=eidxs(nw, i, :pM); label="e $i", color=Cycled(i)) + end + axislegend(ax, position=:rb) + _fig +end diff --git a/docs/make.jl b/docs/make.jl index d935f357..a83a7bae 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -42,7 +42,8 @@ makedocs(; # "Stochastic differential equations" => "generated/StochasticSystem.md", # "Delay differential equations" => "generated/kuramoto_delay.md", "Cascading Failure" => "generated/cascading_failure.md", - "Stress on Truss" => "generated/stress_on_truss.md" + "Stress on Truss" => "generated/stress_on_truss.md", + "Gas Network" => "generated/gas_network.md", ] ], draft=false, From d076663a68d41ce88b7f353d4b0910f370313a88 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Hans=20W=C3=BCrfel?= Date: Thu, 17 Oct 2024 10:45:15 +0200 Subject: [PATCH 2/7] add CairoMakie dep and add note on discontinuities --- docs/examples/gas_network.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/docs/examples/gas_network.jl b/docs/examples/gas_network.jl index 350bba90..11dabeaa 100644 --- a/docs/examples/gas_network.jl +++ b/docs/examples/gas_network.jl @@ -193,6 +193,11 @@ as a blackbox. Additionally, we need to tell ModelingToolkit about the units of this object. This is just used for the static unit check during construction of the model. Later one, when we generate the Julia code from the symbolic reepresentation all units will be stripped. + +!!! note "Discontinuities in RHS" + The picewise linear interpolated function creates discontinuities in the RHS of the system. + However since we know the times exactly, we can handle this by simply giving a list of explicit + tstops to the solve command, to make sure those are hit exactly. =# load2(t) = -Interpolate(SA[0, 4, 12, 20, 24]*3600, SA[20, 30, 10, 30, 20], extrapolate=LinearInterpolations.Constant(20))(t) load3(t) = -Interpolate(SA[0, 4, 12, 20, 24]*3600, SA[40, 50, 30, 50, 40], extrapolate=LinearInterpolations.Constant(40))(t) From 76c933d7dad340258f14c2e0ea08b7493a373175 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Hans=20W=C3=BCrfel?= Date: Tue, 22 Oct 2024 07:13:54 +0200 Subject: [PATCH 3/7] update p --- docs/examples/gas_network.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/examples/gas_network.jl b/docs/examples/gas_network.jl index 11dabeaa..a3a3e38c 100644 --- a/docs/examples/gas_network.jl +++ b/docs/examples/gas_network.jl @@ -159,7 +159,7 @@ depends on temperature and pressure. While the temperature is not considered in at all, I am not sure whether we sould include the pressure dependency. Here, i just calculated it using the "standard" pressure and "standard" temperature. =# -Z̃ = 1 - 3.52 * p̃/pc * exp(-2.26*(T̃/Tc)) + 0.274 * (p̃/pc)^2 * exp(-1.878*(T̃/Tc)) # (5) +Z̃ = 1 - 3.52 * p₁_set/pc * exp(-2.26*(T̃/Tc)) + 0.274 * (p₁_set/pc)^2 * exp(-1.878*(T̃/Tc)) # (5) nothing # hide #= From 4dd36e394dc7abaad10843b76154be76b06dc999 Mon Sep 17 00:00:00 2001 From: Benchmark Bot Date: Thu, 14 Nov 2024 13:28:16 +0100 Subject: [PATCH 4/7] single sided MTK edges with Annotation(symvec) syntax --- ext/MTKExt.jl | 41 ++++++++++++++++++++++------- src/component_functions.jl | 54 +++++++++++++++++++++++++++++++++++--- 2 files changed, 82 insertions(+), 13 deletions(-) diff --git a/ext/MTKExt.jl b/ext/MTKExt.jl index 337d9164..2e1a036c 100644 --- a/ext/MTKExt.jl +++ b/ext/MTKExt.jl @@ -11,7 +11,7 @@ using LinearAlgebra: Diagonal, I using NetworkDynamics: NetworkDynamics, set_metadata!, PureFeedForward, FeedForward, NoFeedForward, PureStateMap -import NetworkDynamics: VertexFunction, EdgeFunction +import NetworkDynamics: VertexFunction, EdgeFunction, AnnotatedSym include("MTKUtils.jl") @@ -49,17 +49,33 @@ function VertexFunction(sys::ODESystem, inputs, outputs; verbose=false, name=get c end +EdgeFunction(sys::ODESystem, srcin, dstin, dstout; kwargs...) = EdgeFunction(sys, srcin, dstin, nothing, dstout; kwargs...) function EdgeFunction(sys::ODESystem, srcin, dstin, srcout, dstout; verbose=false, name=getname(sys), kwargs...) warn_events(sys) srcin = srcin isa AbstractVector ? srcin : [srcin] dstin = dstin isa AbstractVector ? dstin : [dstin] - srcout = srcout isa AbstractVector ? srcout : [srcout] - dstout = dstout isa AbstractVector ? dstout : [dstin] - gen = generate_io_function(sys, (srcin, dstin), (srcout, dstout); verbose) + singlesided = isnothing(srcout) + if singlesided && !(dstout isa AnnotatedSym) + throw(ArgumentError("If you only provide one output (single sided \ + model), it musst be wrapped either in `AntiSymmetric`, `Symmetric` or \ + `Directed`!")) + end + + if singlesided + gwrap = NetworkDynamics.wrapper(dstout) + dstout = NetworkDynamics.sym(dstout) + outs = (dstout, ) + else + srcout = srcout isa AbstractVector ? srcout : [srcout] + dstout = dstout isa AbstractVector ? dstout : [dstin] + outs = (srcout, dstout) + end + + gen = generate_io_function(sys, (srcin, dstin), outs; verbose) f = gen.f - g = gen.g + g = singlesided ? gwrap(gen.g; ff=gen.fftype) : gen.g obsf = gen.obsf _sym = getname.(gen.states) @@ -77,11 +93,16 @@ function EdgeFunction(sys::ODESystem, srcin, dstin, srcout, dstout; verbose=fals insym_dst = [s => _get_metadata(sys, s) for s in _insym_dst] insym = (;src=insym_src, dst=insym_dst) - _outsym_src = getname.(srcout) - outsym_src = [s => _get_metadata(sys, s) for s in _outsym_src] - _outsym_dst = getname.(dstout) - outsym_dst = [s => _get_metadata(sys, s) for s in _outsym_dst] - outsym = (;src=outsym_src, dst=outsym_dst) + if singlesided + _outsym_dst = getname.(dstout) + outsym = [s => _get_metadata(sys, s) for s in _outsym_dst] + else + _outsym_src = getname.(srcout) + outsym_src = [s => _get_metadata(sys, s) for s in _outsym_src] + _outsym_dst = getname.(dstout) + outsym_dst = [s => _get_metadata(sys, s) for s in _outsym_dst] + outsym = (;src=outsym_src, dst=outsym_dst) + end mass_matrix = gen.mass_matrix c = EdgeFunction(;f, g, sym, insym, outsym, psym, obssym, diff --git a/src/component_functions.jl b/src/component_functions.jl index 0899f4c9..2bd3da3c 100644 --- a/src/component_functions.jl +++ b/src/component_functions.jl @@ -55,11 +55,14 @@ Retrieve the feed forward trait of `x`. """ function fftype end +""" + hasfftype(x) + +Defaults to false. Musst be overloaded for objects for which `fftype(x)` is defined. +""" hasfftype(::Any) = false hasff(x) = fftype(x) isa Union{PureFeedForward, FeedForward} -_has_sym_to_outsym_mapping(::Any) = false - """ StateMask(i::AbstractArray) StateMaks(i::Number) @@ -224,6 +227,49 @@ Fiducial(;src, dst, ff=nothing) = Fiducial(src, dst; ff) nothing end +""" + AnnotatedSym{F} + +Wrapper to annotate a vector of symbols as AntiSymmetric, Symmetric or Directed. +""" +struct AnnotatedSym{F} + s::Vector{Symbol} +end +function AnnotatedSym(wrapper, s::AbstractVector{<:Symbol}) + if wrapper ∉ (AntiSymmetric, Symmetric, Directed) + throw(ArgumentError("Only AntiSymmetric, Symmetric and Directed are allowed.")) + end + AnnotatedSym{wrapper}(s) +end +sym(s::AnnotatedSym) = s.s +wrapper(s::AnnotatedSym{W}) where {W} = W +Base.show(io::IO, ::MIME"text/plain", s::AnnotatedSym) = print(io, wrapper(s), "(", sym(s), ")") + +""" + AntiSymmetric(s::AbstractVector{<:Symbol}) + +Annotate a vector of output-symbols as `AntiSymmetric`, used when creating `EdgeFunctions` from +single-sided MTK models. +""" +AntiSymmetric(s::Symbol) = AntiSymmetric([s]) +AntiSymmetric(s::AbstractVector{<:Symbol}) = AnnotatedSym(AntiSymmetric, s) +""" + Symmetric(s::AbstractVector{<:Symbol}) + +Annotate a vector of output-symbols as `Symmetric`, used when creating `EdgeFunctions` from +single-sided MTK models. +""" +Symmetric(s::Symbol) = Symmetric([s]) +Symmetric(s::AbstractVector{<:Symbol}) = AnnotatedSym(Symmetric, s) +""" + Directed(s::AbstractVector{<:Symbol}) + +Annotate a vector of output-symbols as `Directed`, used when creating `EdgeFunctions` from +single-sided MTK models. +""" +Directed(s::Symbol) = Directed([s]) +Directed(s::AbstractVector{<:Symbol}) = AnnotatedSym(Directed, s) + abstract type ComponentFunction end @@ -750,7 +796,8 @@ function _fill_defaults(T, @nospecialize(kwargs)) outsym = dict[:outsym] if T <: EdgeFunction && outsym isa AbstractVector if _has_metadata(outsym) - throw(ArgumentError("Metadata for outsyms can only be provided when using full (;src, dst) form")) + outsym, _metadata = _split_metadata(outsym) + mergewith!(merge!, symmetadata, _metadata) end outsym = dict[:outsym] = _symvec_to_sym_tup(g, outsym) end @@ -937,6 +984,7 @@ function _fill_defaults(T, @nospecialize(kwargs)) end # define the symbolmapping to infer output symbols from state symbols +_has_sym_to_outsym_mapping(::Any) = false _has_sym_to_outsym_mapping(::StateMask) = true _has_sym_to_outsym_mapping(::Directed{<:Any, <:StateMask}) = true _has_sym_to_outsym_mapping(::AntiSymmetric{<:Any, <:StateMask}) = true From 78c9cbde641f2892bcf21dea486eeac657afd339 Mon Sep 17 00:00:00 2001 From: Benchmark Bot Date: Thu, 14 Nov 2024 13:30:33 +0100 Subject: [PATCH 5/7] update gas network example for new f/g interface --- docs/examples/gas_network.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/docs/examples/gas_network.jl b/docs/examples/gas_network.jl index a3a3e38c..da4260fb 100644 --- a/docs/examples/gas_network.jl +++ b/docs/examples/gas_network.jl @@ -213,28 +213,28 @@ nothing #hide To bild the Network we first need to define the components. This is a two step process: - first create the symbolic `ODESystem` using ModelingToolkit -- secondly build a NetworkDynamics component function (`ODEVertex`/`ODEEdge`) based on the symbolic system. +- secondly build a NetworkDynamics component function ([`VertexFunction`](@ref)/[`EdgeFunsction`](@ref)) based on the symbolic system. In the first step we can use the keyword arguments to pass "default" values for our parameters and states. Those values will be automaticially transfered to the metadata of the component function the second step. The second step requires to define the interface variables, i.e. what are the "input" states of your component function and what are the "output" states. -For `ODEVertex` the input state is the aggregated flow of all connected pipes. The output state is the pressure of the node. +For `VertexFunction` the input state is the aggregated flow of all connected pipes. The output state is the pressure of the node. =# @named v1_mtk = ConstantPressureNode(p_set=p₁_set) -v1 = ODEVertex(v1_mtk, [:q̃_nw], [:p]; name=:v1, vidx=1) +v1 = VertexFunction(v1_mtk, [:q̃_nw], [:p]; name=:v1, vidx=1) # @named v2_mtk = VariablePressureNode(C=C₂, load_profile=load2) -v2 = ODEVertex(v2_mtk, [:q̃_nw], [:p]; name=:v2, vidx=2) +v2 = VertexFunction(v2_mtk, [:q̃_nw], [:p]; name=:v2, vidx=2) # @named v3_mtk = VariablePressureNode(C=C₃, load_profile=load3) -v3 = ODEVertex(v3_mtk, [:q̃_nw], [:p]; name=:v3, vidx=3) +v3 = VertexFunction(v3_mtk, [:q̃_nw], [:p]; name=:v3, vidx=3) #= -For `ODEEdge` we have two inputs: the pressure on both source and destination end. +For the edge Model we have two inputs: the pressure on both source and destination end. There is a single output state: the volumetric flow. However we also need to tell NetworkDynamics about the coupling type. In this case we use `AntiSymmetric`, which meas that the source end will recieve the same flow, just inverted sign. @@ -243,9 +243,9 @@ meas that the source end will recieve the same flow, just inverted sign. @named e13_mtk = Pipe(; L=L₁₃, sinθ=sinθ₁₃, A, ρ̃, D, g, c, γ, η, r) @named e23_mtk = Pipe(; L=L₂₃, sinθ=sinθ₂₃, A, ρ̃, D, g, c, γ, η, r) -e12 = ODEEdge(e12_mtk, [:p_src], [:p_dst], [:q̃], AntiSymmetric(); name=:e12, src=1, dst=2) -e13 = ODEEdge(e13_mtk, [:p_src], [:p_dst], [:q̃], AntiSymmetric(); name=:e13, src=1, dst=3) -e23 = ODEEdge(e23_mtk, [:p_src], [:p_dst], [:q̃], AntiSymmetric(); name=:e23, src=2, dst=3) +e12 = EdgeFunction(e12_mtk, [:p_src], [:p_dst], AntiSymmetric([:q̃]); name=:e12, src=1, dst=2) +e13 = EdgeFunction(e13_mtk, [:p_src], [:p_dst], AntiSymmetric([:q̃]); name=:e13, src=1, dst=3) +e23 = EdgeFunction(e23_mtk, [:p_src], [:p_dst], AntiSymmetric([:q̃]); name=:e23, src=2, dst=3) #= To build the network object we just need to pass the vertices and edges to the constructor. From f8e82c7b1a2f05ea430ec23e14d4719262337b52 Mon Sep 17 00:00:00 2001 From: Benchmark Bot Date: Thu, 14 Nov 2024 13:39:48 +0100 Subject: [PATCH 6/7] add docstrings on MTK constructors --- ext/MTKExt.jl | 31 +++++++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) diff --git a/ext/MTKExt.jl b/ext/MTKExt.jl index 2e1a036c..975fe2a1 100644 --- a/ext/MTKExt.jl +++ b/ext/MTKExt.jl @@ -15,6 +15,14 @@ import NetworkDynamics: VertexFunction, EdgeFunction, AnnotatedSym include("MTKUtils.jl") +""" + VertexFunction(sys::ODESystem, inputs, outputs; kwargs...) + +Create a vertex function object from a given `ODESystem` created with ModelingToolkit. +You need to provide 2 lists of symbolic names (`Symbol` or `Vector{Symbols}`): +- `inputs`: names of variables in you equation representing the aggregated edge states +- `outputs`: names of variables in you equation representing the node output +""" function VertexFunction(sys::ODESystem, inputs, outputs; verbose=false, name=getname(sys), kwargs...) warn_events(sys) inputs = inputs isa AbstractVector ? inputs : [inputs] @@ -49,7 +57,30 @@ function VertexFunction(sys::ODESystem, inputs, outputs; verbose=false, name=get c end +""" + EdgeFunction(sys::ODESystem, srcin, dstin, AntiSymmetric(dstout); kwargs...) + +Create a edge function object from a given `ODESystem` created with ModelingToolkit. + +Here you only need to provide one list of output symbols: `dstout`. +To make it clear how to handle the single-sided output definiton, you musst wrap +the symbol vector in +- `AntiSymmetric(dstout)`, +- `Symmetric(dstout)`, or +- `Directed(dstout)`. +""" EdgeFunction(sys::ODESystem, srcin, dstin, dstout; kwargs...) = EdgeFunction(sys, srcin, dstin, nothing, dstout; kwargs...) + +""" + EdgeFunction(sys::ODESystem, srcin, srcout, dstin, dstout; kwargs...) + +Create a edge function object from a given `ODESystem` created with ModelingToolkit. +You need to provide 4 lists of symbolic names (`Symbol` or `Vector{Symbols}`): +- `srcin`: names of variables in you equation representing the node state at the source +- `dstin`: names of variables in you equation representing the node state at the destination +- `srcout`: names of variables in you equation representing the output at the source +- `dstout`: names of variables in you equation representing the output at the destination +""" function EdgeFunction(sys::ODESystem, srcin, dstin, srcout, dstout; verbose=false, name=getname(sys), kwargs...) warn_events(sys) srcin = srcin isa AbstractVector ? srcin : [srcin] From 1e2e89e7202c9477ac70d503603078df1b69d186 Mon Sep 17 00:00:00 2001 From: Benchmark Bot Date: Thu, 14 Nov 2024 13:57:51 +0100 Subject: [PATCH 7/7] remove outdated test --- test/construction_test.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/test/construction_test.jl b/test/construction_test.jl index bf0808c0..4d576865 100644 --- a/test/construction_test.jl +++ b/test/construction_test.jl @@ -203,8 +203,6 @@ end cf = EdgeFunction(f=f, g=Directed(1:2), dim=2, outdim=2) @test isempty(cf.outsym.src) @test cf.outsym.dst == cf.sym - # no metadata without src/dst - @test_throws ArgumentError EdgeFunction(f=f, g=Symmetric(g_single_ff), dim=1, outsym=[:x=>1]) cf = EdgeFunction(g=AntiSymmetric(g_single_pff), outdim=2) @test cf.outsym == (; src=[:₋o₁, :₋o₂], dst=[:o₁, :o₂])