From 5b5e99f1cf07e978c0d1ff892e920a64323e57a7 Mon Sep 17 00:00:00 2001 From: Bart de Koning <74617371+SouthEndMusic@users.noreply.github.com> Date: Tue, 3 Dec 2024 16:05:23 +0100 Subject: [PATCH] Get rid of specifying sources in allocation by edges (#1956) Fixes https://github.com/Deltares/Ribasim/issues/1953. To do: come up with a unified way to handle deprecation, and incorporate this in the developer docs and release process. Maybe also use a 'deprecation' label instead of a 'breaking' label? --------- Co-authored-by: Martijn Visser --- core/src/allocation_init.jl | 71 ++++++++++-------- core/src/allocation_optim.jl | 60 +++++++--------- core/src/callback.jl | 18 +++-- core/src/graph.jl | 36 ++-------- core/src/model.jl | 1 + core/src/parameter.jl | 17 ++--- core/src/read.jl | 72 +++++++++++-------- core/src/util.jl | 35 +++++---- core/src/validation.jl | 55 +++----------- core/test/allocation_test.jl | 57 +++++++-------- core/test/validation_test.jl | 14 +--- docs/changelog.qmd | 3 + docs/concept/allocation.qmd | 4 +- docs/dev/callstacks.qmd | 4 +- docs/guide/examples.ipynb | 4 +- docs/reference/usage.qmd | 1 - python/ribasim/ribasim/geometry/edge.py | 8 +-- python/ribasim/ribasim/migrations.py | 3 + python/ribasim/tests/test_model.py | 1 - .../ribasim_testmodels/allocation.py | 27 ++++--- .../ribasim_testmodels/invalid.py | 1 - ribasim_qgis/core/nodes.py | 1 - 22 files changed, 229 insertions(+), 264 deletions(-) diff --git a/core/src/allocation_init.jl b/core/src/allocation_init.jl index 90e033229..d80660861 100644 --- a/core/src/allocation_init.jl +++ b/core/src/allocation_init.jl @@ -95,7 +95,7 @@ function get_subnetwork_capacity( return capacity end -const allocation_source_nodetypes = +const boundary_source_nodetypes = Set{NodeType.T}([NodeType.LevelBoundary, NodeType.FlowBoundary]) """ @@ -138,10 +138,6 @@ function get_capacity( capacity = get_subnetwork_capacity(p, subnetwork_id) add_subnetwork_connections!(capacity, p, subnetwork_id) - if !valid_sources(p, capacity, subnetwork_id) - error("Errors in sources in allocation network.") - end - return capacity end @@ -272,36 +268,55 @@ function add_constraints_user_source!( end """ -Add the source constraints to the allocation problem. +Add the boundary source constraints to the allocation problem. The actual threshold values will be set before each allocation solve. The constraint indices are (edge_source_id, edge_dst_id). Constraint: -flow over source edge <= source flow in subnetwork +flow over source edge <= source flow in physical layer """ -function add_constraints_source!( +function add_constraints_boundary_source!( problem::JuMP.Model, p::Parameters, subnetwork_id::Int32, )::Nothing - (; graph) = p - edges_source = Tuple{NodeID, NodeID}[] + # Source edges (without the basins) + edges_source = + [edge for edge in source_edges_subnetwork(p, subnetwork_id) if edge[1] != edge[2]] F = problem[:F] - # Find the edges in the whole model which are a source for - # this subnetwork - for edge_metadata in values(graph.edge_data) - (; edge) = edge_metadata - if graph[edge...].subnetwork_id_source == subnetwork_id - push!(edges_source, edge) - end - end + problem[:source_boundary] = JuMP.@constraint( + problem, + [edge_id = edges_source], + F[edge_id] <= 0.0, + base_name = "source_boundary" + ) + return nothing +end - problem[:source] = JuMP.@constraint( +""" +Add main network source constraints to the allocation problem. +The actual threshold values will be set before each allocation solve. +The constraint indices are (edge_source_id, edge_dst_id). + +Constraint: +flow over main network to subnetwork connection edge <= either 0 or allocated amount from the main network +""" +function add_constraints_main_network_source!( + problem::JuMP.Model, + p::Parameters, + subnetwork_id::Int32, +)::Nothing + F = problem[:F] + (; main_network_connections, subnetwork_ids) = p.allocation + subnetwork_id = searchsortedfirst(subnetwork_ids, subnetwork_id) + edges_source = main_network_connections[subnetwork_id] + + problem[:source_main_network] = JuMP.@constraint( problem, [edge_id = edges_source], F[edge_id] <= 0.0, - base_name = "source" + base_name = "source_main_network" ) return nothing end @@ -451,9 +466,9 @@ function allocation_problem( # Add constraints to problem add_constraints_conservation_node!(problem, p, subnetwork_id) - add_constraints_capacity!(problem, capacity, p, subnetwork_id) - add_constraints_source!(problem, p, subnetwork_id) + add_constraints_boundary_source!(problem, p, subnetwork_id) + add_constraints_main_network_source!(problem, p, subnetwork_id) add_constraints_user_source!(problem, p, subnetwork_id) add_constraints_basin_flow!(problem) add_constraints_buffer!(problem) @@ -484,12 +499,12 @@ function get_sources_in_order( sources[edge] = AllocationSource(; edge, type = AllocationSourceType.user_return) end - # Source edges (within subnetwork) - for edge in - sort(only(problem[:source].axes); by = edge -> (edge[1].value, edge[2].value)) - if graph[edge[1]].subnetwork_id == graph[edge[2]].subnetwork_id - sources[edge] = AllocationSource(; edge, type = AllocationSourceType.edge) - end + # Boundary node sources + for edge in sort( + only(problem[:source_boundary].axes); + by = edge -> (edge[1].value, edge[2].value), + ) + sources[edge] = AllocationSource(; edge, type = AllocationSourceType.boundary_node) end # Basins with level demand diff --git a/core/src/allocation_optim.jl b/core/src/allocation_optim.jl index 9a3f0d4e0..6e2404862 100644 --- a/core/src/allocation_optim.jl +++ b/core/src/allocation_optim.jl @@ -127,8 +127,7 @@ function assign_allocations!( # If this edge is a source edge from the main network to a subnetwork, # and demands are being collected, add its flow to the demand of this edge if optimization_type == OptimizationType.collect_demands - if graph[edge...].subnetwork_id_source == subnetwork_id && - edge ∈ main_network_source_edges + if edge in main_network_source_edges allocated = flow[edge] subnetwork_demands[edge][priority_idx] += allocated end @@ -203,22 +202,13 @@ function set_initial_capacities_source!( p::Parameters, )::Nothing (; sources) = allocation_model - (; graph, allocation) = p - (; mean_input_flows) = allocation (; subnetwork_id) = allocation_model - main_network_source_edges = get_main_network_connections(p, subnetwork_id) - for edge_metadata in values(graph.edge_data) - (; edge) = edge_metadata - if graph[edge...].subnetwork_id_source == subnetwork_id - # If it is a source edge for this allocation problem - if edge ∉ main_network_source_edges - # Reset the source to the averaged flow over the last allocation period - source = sources[edge] - @assert source.type == AllocationSourceType.edge - source.capacity = mean_input_flows[edge][] - end - end + mean_input_flows_subnetwork_ = mean_input_flows_subnetwork(p, subnetwork_id) + + for edge in keys(mean_input_flows_subnetwork_) + source = sources[edge] + source.capacity = mean_input_flows_subnetwork_[edge] end return nothing end @@ -231,7 +221,7 @@ function reduce_source_capacity!(problem::JuMP.Model, source::AllocationSource): used_capacity = if source.type in ( - AllocationSourceType.edge, + AllocationSourceType.boundary_node, AllocationSourceType.main_to_sub, AllocationSourceType.user_return, ) @@ -343,11 +333,10 @@ function get_basin_data( u::ComponentVector, node_id::NodeID, ) - (; graph, allocation, basin) = p - (; Δt_allocation) = allocation_model - (; mean_input_flows) = allocation + (; graph, basin) = p + (; Δt_allocation, subnetwork_id) = allocation_model @assert node_id.type == NodeType.Basin - influx = mean_input_flows[(node_id, node_id)][] + influx = mean_input_flows_subnetwork(p, subnetwork_id)[(node_id, node_id)] storage_basin = basin.current_properties.current_storage[parent(u)][node_id.idx] control_inneighbors = inneighbor_labels_type(graph, node_id, EdgeType.control) if isempty(control_inneighbors) @@ -841,9 +830,10 @@ function set_source_capacity!( optimization_type::OptimizationType.T, )::Nothing (; problem, sources) = allocation_model - constraints_source_edge = problem[:source] - constraints_source_basin = problem[:basin_outflow] + constraints_source_boundary = problem[:source_boundary] constraints_source_user_out = problem[:source_user] + constraints_source_main_network = problem[:source_main_network] + constraints_source_basin = problem[:basin_outflow] constraints_source_buffer = problem[:flow_buffer_outflow] for source in values(sources) @@ -860,16 +850,17 @@ function set_source_capacity!( 0.0 end - constraint = - if source.type in (AllocationSourceType.edge, AllocationSourceType.main_to_sub) - constraints_source_edge[edge] - elseif source.type == AllocationSourceType.basin - constraints_source_basin[edge[1]] - elseif source.type == AllocationSourceType.user_return - constraints_source_user_out[edge[1]] - elseif source.type == AllocationSourceType.buffer - constraints_source_buffer[edge[1]] - end + constraint = if source.type == AllocationSourceType.boundary_node + constraints_source_boundary[edge] + elseif source.type == AllocationSourceType.main_to_sub + constraints_source_main_network[edge] + elseif source.type == AllocationSourceType.basin + constraints_source_basin[edge[1]] + elseif source.type == AllocationSourceType.user_return + constraints_source_user_out[edge[1]] + elseif source.type == AllocationSourceType.buffer + constraints_source_buffer[edge[1]] + end JuMP.set_normalized_rhs(constraint, capacity_effective) end @@ -1062,7 +1053,8 @@ function empty_sources!(allocation_model::AllocationModel, allocation::Allocatio (; problem) = allocation_model (; subnetwork_demands) = allocation - for constraint_set_name in [:source, :source_user, :basin_outflow, :flow_buffer_outflow] + for constraint_set_name in + [:source_boundary, :source_user, :basin_outflow, :flow_buffer_outflow] constraint_set = problem[constraint_set_name] for key in only(constraint_set.axes) # Do not set the capacity to 0.0 if the edge diff --git a/core/src/callback.jl b/core/src/callback.jl index 213e97284..467e71060 100644 --- a/core/src/callback.jl +++ b/core/src/callback.jl @@ -153,8 +153,11 @@ function update_cumulative_flows!(u, t, integrator)::Nothing end # Update realized flows for allocation input - for edge in keys(allocation.mean_input_flows) - allocation.mean_input_flows[edge] += flow_update_on_edge(integrator, edge) + for subnetwork_id in allocation.subnetwork_ids + mean_input_flows_subnetwork_ = mean_input_flows_subnetwork(p, subnetwork_id) + for edge in keys(mean_input_flows_subnetwork_) + mean_input_flows_subnetwork_[edge] += flow_update_on_edge(integrator, edge) + end end # Update realized flows for allocation output @@ -773,8 +776,10 @@ function update_allocation!(integrator)::Nothing # Divide by the allocation Δt to get the mean input flows from the cumulative flows (; Δt_allocation) = allocation_models[1] - for edge in keys(mean_input_flows) - mean_input_flows[edge] /= Δt_allocation + for mean_input_flows_subnetwork in values(mean_input_flows) + for edge in keys(mean_input_flows_subnetwork) + mean_input_flows_subnetwork[edge] /= Δt_allocation + end end # Divide by the allocation Δt to get the mean realized flows from the cumulative flows @@ -797,11 +802,14 @@ function update_allocation!(integrator)::Nothing end # Reset the mean flows - for mean_flows in (mean_input_flows, mean_realized_flows) + for mean_flows in mean_input_flows for edge in keys(mean_flows) mean_flows[edge] = 0.0 end end + for edge in keys(mean_realized_flows) + mean_realized_flows[edge] = 0.0 + end end "Load updates from 'TabulatedRatingCurve / time' into the parameters" diff --git a/core/src/graph.jl b/core/src/graph.jl index d7e32d6a8..e2a3273e0 100644 --- a/core/src/graph.jl +++ b/core/src/graph.jl @@ -19,8 +19,7 @@ function create_graph(db::DB, config::Config)::MetaGraph FromNode.node_type AS from_node_type, ToNode.node_id AS to_node_id, ToNode.node_type AS to_node_type, - Edge.edge_type, - Edge.subnetwork_id + Edge.edge_type FROM Edge LEFT JOIN Node AS FromNode ON FromNode.node_id = Edge.from_node_id LEFT JOIN Node AS ToNode ON ToNode.node_id = Edge.to_node_id @@ -28,8 +27,7 @@ function create_graph(db::DB, config::Config)::MetaGraph ) # Node IDs per subnetwork node_ids = Dict{Int32, Set{NodeID}}() - # Source edges per subnetwork - edges_source = Dict{Int32, Set{EdgeMetadata}}() + # The metadata of the flow edges in the order in which they are in the input # and will be in the output flow_edges = EdgeMetadata[] @@ -57,15 +55,8 @@ function create_graph(db::DB, config::Config)::MetaGraph end errors = false - for (; - edge_id, - from_node_type, - from_node_id, - to_node_type, - to_node_id, - edge_type, - subnetwork_id, - ) in edge_rows + for (; edge_id, from_node_type, from_node_id, to_node_type, to_node_id, edge_type) in + edge_rows try # hasfield does not work edge_type = getfield(EdgeType, Symbol(edge_type)) @@ -74,15 +65,8 @@ function create_graph(db::DB, config::Config)::MetaGraph end id_src = NodeID(from_node_type, from_node_id, db) id_dst = NodeID(to_node_type, to_node_id, db) - if ismissing(subnetwork_id) - subnetwork_id = 0 - end - edge_metadata = EdgeMetadata(; - id = edge_id, - type = edge_type, - subnetwork_id_source = subnetwork_id, - edge = (id_src, id_dst), - ) + edge_metadata = + EdgeMetadata(; id = edge_id, type = edge_type, edge = (id_src, id_dst)) if edge_type == EdgeType.flow push!(flow_edges, edge_metadata) end @@ -91,12 +75,6 @@ function create_graph(db::DB, config::Config)::MetaGraph @error "Duplicate edge" id_src id_dst end graph[id_src, id_dst] = edge_metadata - if subnetwork_id != 0 - if !haskey(edges_source, subnetwork_id) - edges_source[subnetwork_id] = Set{EdgeMetadata}() - end - push!(edges_source[subnetwork_id], edge_metadata) - end end if errors error("Invalid edges found") @@ -106,7 +84,7 @@ function create_graph(db::DB, config::Config)::MetaGraph error("Incomplete connectivity in subnetwork") end - graph_data = (; node_ids, edges_source, flow_edges, config.solver.saveat) + graph_data = (; node_ids, flow_edges, config.solver.saveat) @reset graph.graph_data = graph_data return graph diff --git a/core/src/model.jl b/core/src/model.jl index e0047ca0d..2197d72b7 100644 --- a/core/src/model.jl +++ b/core/src/model.jl @@ -47,6 +47,7 @@ function Model(config::Config)::Model # so we can directly close it again. db = SQLite.DB(db_path) + database_warning(db) if !valid_nodes(db) error("Invalid nodes found.") end diff --git a/core/src/parameter.jl b/core/src/parameter.jl index 4c56611be..03b909799 100644 --- a/core/src/parameter.jl +++ b/core/src/parameter.jl @@ -128,7 +128,7 @@ The caches are always initialized with zeros """ cache(len::Int)::Cache = LazyBufferCache(Returns(len); initializer! = set_zero!) -@enumx AllocationSourceType edge basin main_to_sub user_return buffer +@enumx AllocationSourceType boundary_node basin main_to_sub user_return buffer """ Data structure for a single source within an allocation subnetwork. @@ -180,20 +180,21 @@ main_network_connections: (from_id, to_id) from the main network to the subnetwo priorities: All used priority values. subnetwork_demands: The demand of an edge from the main network to a subnetwork subnetwork_allocateds: The allocated flow of an edge from the main network to a subnetwork -mean_input_flows: Flows averaged over Δt_allocation over edges that are allocation sources +mean_input_flows: Per subnetwork, flows averaged over Δt_allocation over edges that are allocation sources mean_realized_flows: Flows averaged over Δt_allocation over edges that realize a demand record_demand: A record of demands and allocated flows for nodes that have these record_flow: A record of all flows computed by allocation optimization, eventually saved to output file """ @kwdef struct Allocation - subnetwork_ids::Vector{Int32} = [] - allocation_models::Vector{AllocationModel} = [] - main_network_connections::Vector{Vector{Tuple{NodeID, NodeID}}} = [] + subnetwork_ids::Vector{Int32} = Int32[] + allocation_models::Vector{AllocationModel} = AllocationModel[] + main_network_connections::Vector{Vector{Tuple{NodeID, NodeID}}} = + Vector{Tuple{NodeID, NodeID}}[] priorities::Vector{Int32} subnetwork_demands::Dict{Tuple{NodeID, NodeID}, Vector{Float64}} = Dict() subnetwork_allocateds::Dict{Tuple{NodeID, NodeID}, Vector{Float64}} = Dict() - mean_input_flows::Dict{Tuple{NodeID, NodeID}, Float64} + mean_input_flows::Vector{Dict{Tuple{NodeID, NodeID}, Float64}} mean_realized_flows::Dict{Tuple{NodeID, NodeID}, Float64} record_demand::@NamedTuple{ time::Vector{Float64}, @@ -255,14 +256,11 @@ end Type for storing metadata of edges in the graph: id: ID of the edge (only used for labeling flow output) type: type of the edge -subnetwork_id_source: ID of subnetwork where this edge is a source - (0 if not a source) edge: (from node ID, to node ID) """ @kwdef struct EdgeMetadata id::Int32 type::EdgeType.T - subnetwork_id_source::Int32 edge::Tuple{NodeID, NodeID} end @@ -902,7 +900,6 @@ const ModelGraph = MetaGraph{ EdgeMetadata, @NamedTuple{ node_ids::Dict{Int32, Set{NodeID}}, - edges_source::Dict{Int32, Set{EdgeMetadata}}, flow_edges::Vector{EdgeMetadata}, saveat::Float64, }, diff --git a/core/src/read.jl b/core/src/read.jl index edcce108c..e99d207e0 100644 --- a/core/src/read.jl +++ b/core/src/read.jl @@ -1268,40 +1268,56 @@ function Subgrid(db::DB, config::Config, basin::Basin)::Subgrid end function Allocation(db::DB, config::Config, graph::MetaGraph)::Allocation - mean_input_flows = Dict{Tuple{NodeID, NodeID}, Float64}() + mean_input_flows = Dict{Tuple{NodeID, NodeID}, Float64}[] + mean_realized_flows = Dict{Tuple{NodeID, NodeID}, Float64}() + subnetwork_ids = sort(collect(keys(graph[].node_ids))) - # Find edges which serve as sources in allocation - for edge_metadata in values(graph.edge_data) - (; subnetwork_id_source, edge) = edge_metadata - if subnetwork_id_source != 0 - mean_input_flows[edge] = 0.0 + if config.allocation.use_allocation + for subnetwork_id in subnetwork_ids + push!(mean_input_flows, Dict{Tuple{NodeID, NodeID}, Float64}()) end - end - # Find basins with a level demand - for node_id in values(graph.vertex_labels) - if has_external_demand(graph, node_id, :level_demand)[1] - mean_input_flows[(node_id, node_id)] = 0.0 + # Find edges which serve as sources in allocation + for edge_metadata in values(graph.edge_data) + (; edge) = edge_metadata + id_source, _ = edge + if id_source.type in boundary_source_nodetypes + (; subnetwork_id) = graph[id_source] + # Check whether the source node is part of a subnetwork + if subnetwork_id ≠ 0 + subnetwork_idx = searchsortedfirst(subnetwork_ids, subnetwork_id) + mean_input_flows[subnetwork_idx][edge] = 0.0 + end + end end - end - mean_realized_flows = Dict{Tuple{NodeID, NodeID}, Float64}() + # Find basins with a level demand + for node_id in values(graph.vertex_labels) + if has_external_demand(graph, node_id, :level_demand)[1] + subnetwork_id = graph[node_id].subnetwork_id + subnetwork_idx = searchsortedfirst(subnetwork_ids, subnetwork_id) + mean_input_flows[subnetwork_idx][(node_id, node_id)] = 0.0 + end + end - # Find edges that realize a demand - for edge_metadata in values(graph.edge_data) - (; type, edge) = edge_metadata - - src_id, dst_id = edge - user_demand_inflow = (type == EdgeType.flow) && (dst_id.type == NodeType.UserDemand) - level_demand_inflow = - (type == EdgeType.control) && (src_id.type == NodeType.LevelDemand) - flow_demand_inflow = - (type == EdgeType.flow) && has_external_demand(graph, dst_id, :flow_demand)[1] - - if user_demand_inflow || flow_demand_inflow - mean_realized_flows[edge] = 0.0 - elseif level_demand_inflow - mean_realized_flows[(dst_id, dst_id)] = 0.0 + # Find edges that realize a demand + for edge_metadata in values(graph.edge_data) + (; type, edge) = edge_metadata + + src_id, dst_id = edge + user_demand_inflow = + (type == EdgeType.flow) && (dst_id.type == NodeType.UserDemand) + level_demand_inflow = + (type == EdgeType.control) && (src_id.type == NodeType.LevelDemand) + flow_demand_inflow = + (type == EdgeType.flow) && + has_external_demand(graph, dst_id, :flow_demand)[1] + + if user_demand_inflow || flow_demand_inflow + mean_realized_flows[edge] = 0.0 + elseif level_demand_inflow + mean_realized_flows[(dst_id, dst_id)] = 0.0 + end end end diff --git a/core/src/util.jl b/core/src/util.jl index 379f2abd5..2eac3f267 100644 --- a/core/src/util.jl +++ b/core/src/util.jl @@ -587,15 +587,17 @@ function set_initial_allocation_mean_flows!(integrator)::Nothing du = get_du(integrator) water_balance!(du, u, p, t) - for edge in keys(mean_input_flows) - if edge[1] == edge[2] - q = get_influx(du, edge[1], p) - else - q = get_flow(du, p, t, edge) + for mean_input_flows_subnetwork in values(mean_input_flows) + for edge in keys(mean_input_flows_subnetwork) + if edge[1] == edge[2] + q = get_influx(du, edge[1], p) + else + q = get_flow(du, p, t, edge) + end + # Multiply by Δt_allocation as averaging divides by this factor + # in update_allocation! + mean_input_flows_subnetwork[edge] = q * Δt_allocation end - # Multiply by Δt_allocation as averaging divides by this factor - # in update_allocation! - mean_input_flows[edge] = q * Δt_allocation end # Mean realized demands for basins are calculated as Δstorage/Δt @@ -1004,12 +1006,8 @@ function set_state_flow_edges(p::Parameters, u0::ComponentVector)::Parameters state_inflow_edges = Vector{EdgeMetadata}[] state_outflow_edges = Vector{EdgeMetadata}[] - placeholder_edge = EdgeMetadata( - 0, - EdgeType.flow, - 0, - (NodeID(:Terminal, 0, 0), NodeID(:Terminal, 0, 0)), - ) + placeholder_edge = + EdgeMetadata(0, EdgeType.flow, (NodeID(:Terminal, 0, 0), NodeID(:Terminal, 0, 0))) for node_name in keys(u0) if hasfield(Parameters, node_name) @@ -1172,3 +1170,12 @@ function min_low_user_demand_level_factor( one(T) end end + +function mean_input_flows_subnetwork(p::Parameters, subnetwork_id::Int32) + (; mean_input_flows, subnetwork_ids) = p.allocation + subnetwork_idx = searchsortedfirst(subnetwork_ids, subnetwork_id) + return mean_input_flows[subnetwork_idx] +end + +source_edges_subnetwork(p::Parameters, subnetwork_id::Int32) = + keys(mean_input_flows_subnetwork(p, subnetwork_id)) diff --git a/core/src/validation.jl b/core/src/validation.jl index cfed36c6c..81a372421 100644 --- a/core/src/validation.jl +++ b/core/src/validation.jl @@ -182,6 +182,14 @@ function valid_nodes(db::DB)::Bool return !errors end +function database_warning(db::DB)::Nothing + cols = SQLite.columns(db, "Edge") + if "subnetwork_id" in cols.name + @warn "The 'subnetwork_id' column in the 'Edge' table is deprecated since ribasim v2024.12." + end + return nothing +end + """ Test for each node given its node type whether the nodes that # are downstream ('down-edge') of this node are of an allowed type @@ -626,53 +634,6 @@ function valid_discrete_control(p::Parameters, config::Config)::Bool return !errors end -""" -An allocation source edge is valid if either: - - The edge connects the main network to a subnetwork - - The edge comes from a source node -""" -function valid_sources( - p::Parameters, - capacity::JuMP.Containers.SparseAxisArray{Float64, 2, Tuple{NodeID, NodeID}}, - subnetwork_id::Int32, -)::Bool - (; graph) = p - - errors = false - - # Loop over edges that were assigned a capacity - for edge in keys(capacity.data) - - # For an edge (id_a, id_b) in the physical model - # the reverse (id_b, id_a) can exist in the allocation subnetwork - if !haskey(graph, edge...) - edge = reverse(edge) - end - - (id_source, id_dst) = edge - - # Whether the current edge is a source for the current subnetwork - if graph[edge...].subnetwork_id_source == subnetwork_id - from_source_node = id_source.type in allocation_source_nodetypes - - if is_main_network(subnetwork_id) - if !from_source_node - errors = true - @error "The source node of source edge $edge in the main network must be one of $allocation_source_nodetypes." - end - else - from_main_network = is_main_network(graph[id_source].subnetwork_id) - - if !from_source_node && !from_main_network - errors = true - @error "The source node of source edge $edge for subnetwork $subnetwork_id is neither a source node nor is it coming from the main network." - end - end - end - end - return !errors -end - function valid_priorities(priorities::Vector{Int32}, use_allocation::Bool)::Bool if use_allocation && any(iszero, priorities) return false diff --git a/core/test/allocation_test.jl b/core/test/allocation_test.jl index 6da2fec36..3dd6c7ba0 100644 --- a/core/test/allocation_test.jl +++ b/core/test/allocation_test.jl @@ -14,7 +14,8 @@ (; graph, allocation) = p - allocation.mean_input_flows[(NodeID(:FlowBoundary, 1, p), NodeID(:Basin, 2, p))] = 4.5 + allocation.mean_input_flows[1][(NodeID(:FlowBoundary, 1, p), NodeID(:Basin, 2, p))] = + 4.5 allocation_model = p.allocation.allocation_models[1] (; flow) = allocation_model u = ComponentVector(; storage = zeros(length(p.basin.node_id))) @@ -105,11 +106,11 @@ end # In each subnetwork, the connection from the main network to the subnetwork is # interpreted as a source - @test Ribasim.get_allocation_model(p, Int32(3)).problem[:source].axes[1] == + @test Ribasim.get_allocation_model(p, Int32(3)).problem[:source_main_network].axes[1] == [(NodeID(:Basin, 2, p), NodeID(:Pump, 11, p))] - @test Ribasim.get_allocation_model(p, Int32(5)).problem[:source].axes[1] == + @test Ribasim.get_allocation_model(p, Int32(5)).problem[:source_main_network].axes[1] == [(NodeID(:Basin, 6, p), NodeID(:Pump, 24, p))] - @test Ribasim.get_allocation_model(p, Int32(7)).problem[:source].axes[1] == + @test Ribasim.get_allocation_model(p, Int32(7)).problem[:source_main_network].axes[1] == [(NodeID(:Basin, 10, p), NodeID(:Pump, 38, p))] end @@ -174,7 +175,7 @@ end # Running full allocation algorithm (; Δt_allocation) = allocation_models[1] - mean_input_flows[(NodeID(:FlowBoundary, 1, p), NodeID(:Basin, 2, p))] = + mean_input_flows[1][(NodeID(:FlowBoundary, 1, p), NodeID(:Basin, 2, p))] = 4.5 * Δt_allocation Ribasim.update_allocation!(model.integrator) @@ -224,9 +225,9 @@ end allocation t = 0.0 - # Set flows of sources in - mean_input_flows[(NodeID(:FlowBoundary, 58, p), NodeID(:Basin, 16, p))] = 1.0 - mean_input_flows[(NodeID(:FlowBoundary, 59, p), NodeID(:Basin, 44, p))] = 1e-3 + # Set flows of sources in subnetworks + mean_input_flows[2][(NodeID(:FlowBoundary, 58, p), NodeID(:Basin, 16, p))] = 1.0 + mean_input_flows[4][(NodeID(:FlowBoundary, 59, p), NodeID(:Basin, 44, p))] = 1e-3 # Collecting demands u = ComponentVector(; storage = zeros(length(basin.node_id))) @@ -251,24 +252,24 @@ end du = get_du(model.integrator) Ribasim.formulate_storages!(current_storage, du, u, p, t) - current_storage ≈ Float32[ - 1.0346e6, - 1.0346e6, - 1.0346e6, - 1.0346e6, - 1.0346e6, - 13.83, - 40.10, - 6.029e5, - 4641, - 2402, - 6.03995, - 928.8, - 8.017, - 10417, - 5.619, - 10417, - 4.057, + @test current_storage ≈ Float32[ + 1.0346908f6, + 1.03469f6, + 1.0346894f6, + 1.034689f6, + 1.0346888f6, + 13.833241, + 40.109993, + 187761.73, + 4641.365, + 2402.6687, + 6.039952, + 928.84283, + 8.0175905, + 10419.247, + 5.619053, + 10419.156, + 4.057502, ] end @@ -287,7 +288,7 @@ end (; user_demand, graph, allocation, basin, level_demand) = p # Initial "integrated" vertical flux - @test allocation.mean_input_flows[(NodeID(:Basin, 2, p), NodeID(:Basin, 2, p))] ≈ 1e2 + @test allocation.mean_input_flows[1][(NodeID(:Basin, 2, p), NodeID(:Basin, 2, p))] ≈ 1e2 Ribasim.solve!(model) @@ -614,5 +615,5 @@ end # Given a max_level of Inf, the basin capacity is 0.0 because it is not possible for the basin level to be > Inf @test Ribasim.get_basin_capacity(allocation_models[1], u, p, t, basin.node_id[1]) == 0.0 @test Ribasim.get_basin_capacity(allocation_models[1], u, p, t, basin.node_id[2]) == 0.0 - @test Ribasim.get_basin_capacity(allocation_models[1], u, p, t, basin.node_id[3]) == 0.0 + @test Ribasim.get_basin_capacity(allocation_models[2], u, p, t, basin.node_id[3]) == 0.0 end diff --git a/core/test/validation_test.jl b/core/test/validation_test.jl index 547a11ff4..e46c0ef74 100644 --- a/core/test/validation_test.jl +++ b/core/test/validation_test.jl @@ -86,12 +86,7 @@ end graph[NodeID(:Pump, 6, 1)] = NodeMetadata(:pump, 9) function set_edge_metadata!(id_1, id_2, edge_type) - graph[id_1, id_2] = EdgeMetadata(; - id = 0, - type = edge_type, - subnetwork_id_source = 0, - edge = (id_1, id_2), - ) + graph[id_1, id_2] = EdgeMetadata(; id = 0, type = edge_type, edge = (id_1, id_2)) return nothing end @@ -148,12 +143,7 @@ end graph[NodeID(:Basin, 7, 1)] = NodeMetadata(:basin, 0) function set_edge_metadata!(id_1, id_2, edge_type) - graph[id_1, id_2] = EdgeMetadata(; - id = 0, - type = edge_type, - subnetwork_id_source = 0, - edge = (id_1, id_2), - ) + graph[id_1, id_2] = EdgeMetadata(; id = 0, type = edge_type, edge = (id_1, id_2)) return nothing end diff --git a/docs/changelog.qmd b/docs/changelog.qmd index d9c04350e..373ab6cc8 100644 --- a/docs/changelog.qmd +++ b/docs/changelog.qmd @@ -5,6 +5,9 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), ## [Unreleased] +### Changed +- The Edge table no longer supports `subnetwork_id`; this is automatically inferred. [#1956](https://github.com/Deltares/Ribasim/pull/1956) + ## [v2024.11.0] - 2024-10-08 This major new release contains many improvements. diff --git a/docs/concept/allocation.qmd b/docs/concept/allocation.qmd index f90a83cec..91318fd84 100644 --- a/docs/concept/allocation.qmd +++ b/docs/concept/allocation.qmd @@ -42,9 +42,9 @@ The allocation problem is solved per subnetwork, which is given by a subset $S \ Sources are indicated by a set of edges in the subnetwork $$ -E_S^\text{source} \subset E. +E_S^\text{source} \subset E, $$ -That is, if $(i,j) \in E_S^\text{source}$, then the average over the last allocation interval $\Delta t_{\text{alloc}}$ of the of the flow over this edge +which are automatically inferred as all edges that point out of LevelBoundary or FlowBoundary nodes. That is, if $(i,j) \in E_S^\text{source}$, then the average over the last allocation interval $\Delta t_{\text{alloc}}$ of the of the flow over this edge $$ \frac{1}{\Delta t_{\text{alloc}}}\int_{t - \Delta t_{\text{alloc}}}^tQ_{ij}(t') dt' $$ diff --git a/docs/dev/callstacks.qmd b/docs/dev/callstacks.qmd index 5bd4fe372..a36c75cdf 100644 --- a/docs/dev/callstacks.qmd +++ b/docs/dev/callstacks.qmd @@ -74,10 +74,12 @@ toml_path = normpath( @__DIR__, "../../generated_testmodels/main_network_with_subnetworks/ribasim.toml", ) -config = Ribasim.Config(toml_path; allocation_use_allocation = false) +config = Ribasim.Config(toml_path) db_path = Ribasim.database_path(config) db = SQLite.DB(db_path) p = Ribasim.Parameters(db, config) +empty!(p.allocation.subnetwork_ids) +empty!(p.allocation.main_network_connections) graph, verts = tracecall((Ribasim,), Ribasim.initialize_allocation!, (p, config)) plot_graph(graph) ``` diff --git a/docs/guide/examples.ipynb b/docs/guide/examples.ipynb index 31e899dd6..3e09e5986 100644 --- a/docs/guide/examples.ipynb +++ b/docs/guide/examples.ipynb @@ -1203,7 +1203,7 @@ "metadata": {}, "outputs": [], "source": [ - "model.edge.add(model.flow_boundary[1], model.basin[2], subnetwork_id=1)\n", + "model.edge.add(model.flow_boundary[1], model.basin[2])\n", "model.edge.add(model.basin[2], model.user_demand[3])\n", "model.edge.add(model.basin[2], model.linear_resistance[4])\n", "model.edge.add(model.linear_resistance[4], model.basin[5])\n", @@ -1499,7 +1499,7 @@ "metadata": {}, "outputs": [], "source": [ - "model.edge.add(model.flow_boundary[1], model.basin[2], subnetwork_id=2)\n", + "model.edge.add(model.flow_boundary[1], model.basin[2])\n", "model.edge.add(model.basin[2], model.user_demand[3])\n", "model.edge.add(model.level_demand[4], model.basin[2])\n", "model.edge.add(model.user_demand[3], model.basin[5])\n", diff --git a/docs/reference/usage.qmd b/docs/reference/usage.qmd index 63eeedd99..afe2f1529 100644 --- a/docs/reference/usage.qmd +++ b/docs/reference/usage.qmd @@ -200,7 +200,6 @@ to_node_id | Int32 | - edge_type | String | must be "flow" or "control" geom | LineString or MultiLineString | (optional) name | String | (optional, does not have to be unique) -subnetwork_id | Int32 | (optional, denotes source in allocation network) Similarly to the node table, you can use a geometry to visualize the connections between the nodes in QGIS. For instance, you can draw a line connecting the two node coordinates. diff --git a/python/ribasim/ribasim/geometry/edge.py b/python/ribasim/ribasim/geometry/edge.py index dc8f6e202..bc473458f 100644 --- a/python/ribasim/ribasim/geometry/edge.py +++ b/python/ribasim/ribasim/geometry/edge.py @@ -47,7 +47,6 @@ class EdgeSchema(_GeoBaseSchema): from_node_id: Series[Int32] = pa.Field(default=0) to_node_id: Series[Int32] = pa.Field(default=0) edge_type: Series[str] = pa.Field(default="flow") - subnetwork_id: Series[pd.Int32Dtype] = pa.Field(default=pd.NA, nullable=True) geometry: GeoSeries[LineString] = pa.Field(default=None, nullable=True) @classmethod @@ -73,7 +72,6 @@ def add( to_node: NodeData, geometry: LineString | MultiLineString | None = None, name: str = "", - subnetwork_id: int | None = None, edge_id: Optional[NonNegativeInt] = None, **kwargs, ): @@ -92,9 +90,8 @@ def add( The geometry of a line. If not supplied, it creates a straight line between the nodes. name : str An optional name for the edge. - subnetwork_id : int | None - An optional subnetwork id for the edge. This edge indicates a source for - the allocation algorithm, and should thus not be set for every edge in a subnetwork. + id : int + An optional non-negative edge ID. If not supplied, it will be automatically generated. **kwargs : Dict """ if not can_connect(from_node.node_type, to_node.node_type): @@ -125,7 +122,6 @@ def add( "to_node_id": [to_node.node_id], "edge_type": [edge_type], "name": [name], - "subnetwork_id": [subnetwork_id], **kwargs, }, geometry=geometry_to_append, diff --git a/python/ribasim/ribasim/migrations.py b/python/ribasim/ribasim/migrations.py index fc51eff52..31a326e83 100644 --- a/python/ribasim/ribasim/migrations.py +++ b/python/ribasim/ribasim/migrations.py @@ -27,6 +27,9 @@ def edgeschema_migration(gdf: GeoDataFrame, schema_version: int) -> GeoDataFrame warnings.warn("Migrating outdated Edge table.", UserWarning) assert gdf["edge_id"].is_unique, "Edge IDs have to be unique." gdf.set_index("edge_id", inplace=True) + if "subnetwork_id" in gdf.columns: + warnings.warn("Migrating outdated Edge table.", UserWarning) + gdf.drop(columns="subnetwork_id", inplace=True, errors="ignore") return gdf diff --git a/python/ribasim/tests/test_model.py b/python/ribasim/tests/test_model.py index ba703d5fa..f5e341561 100644 --- a/python/ribasim/tests/test_model.py +++ b/python/ribasim/tests/test_model.py @@ -145,7 +145,6 @@ def test_edge_table(basic): df = model.edge.df assert df.geometry.is_unique assert df.from_node_id.dtype == np.int32 - assert df.subnetwork_id.dtype == pd.Int32Dtype() assert df.crs == CRS.from_epsg(28992) diff --git a/python/ribasim_testmodels/ribasim_testmodels/allocation.py b/python/ribasim_testmodels/ribasim_testmodels/allocation.py index dce660924..dce319344 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/allocation.py +++ b/python/ribasim_testmodels/ribasim_testmodels/allocation.py @@ -153,7 +153,7 @@ def subnetwork_model() -> Model: ) model.outlet.add(Node(13, Point(2, 4), subnetwork_id=2), outlet_data) - model.edge.add(model.flow_boundary[1], model.basin[2], subnetwork_id=2) + model.edge.add(model.flow_boundary[1], model.basin[2]) model.edge.add(model.basin[2], model.outlet[3]) model.edge.add(model.outlet[3], model.terminal[4]) model.edge.add(model.basin[2], model.user_demand[10]) @@ -268,7 +268,7 @@ def looped_subnetwork_model() -> Model: ], ) - model.edge.add(model.flow_boundary[5], model.basin[2], subnetwork_id=2) + model.edge.add(model.flow_boundary[5], model.basin[2]) model.edge.add(model.basin[2], model.outlet[3]) model.edge.add(model.outlet[3], model.terminal[4]) model.edge.add(model.basin[2], model.user_demand[1]) @@ -346,7 +346,7 @@ def minimal_subnetwork_model() -> Model: ], ) - model.edge.add(model.flow_boundary[1], model.basin[2], subnetwork_id=2) + model.edge.add(model.flow_boundary[1], model.basin[2]) model.edge.add(model.basin[2], model.pump[3]) model.edge.add(model.pump[3], model.basin[4]) model.edge.add(model.basin[4], model.user_demand[5]) @@ -406,7 +406,7 @@ def allocation_example_model() -> Model: ) model.terminal.add(Node(8, Point(5, 0), subnetwork_id=2)) - model.edge.add(model.flow_boundary[1], model.basin[2], subnetwork_id=2) + model.edge.add(model.flow_boundary[1], model.basin[2]) model.edge.add(model.basin[2], model.user_demand[3]) model.edge.add(model.basin[2], model.linear_resistance[4]) model.edge.add(model.linear_resistance[4], model.basin[5]) @@ -628,7 +628,7 @@ def main_network_with_subnetworks_model() -> Model: [user_demand.Static(return_factor=[0.9], priority=2, min_level=0.0)], ) - model.edge.add(model.flow_boundary[1], model.basin[2], subnetwork_id=1) + model.edge.add(model.flow_boundary[1], model.basin[2]) model.edge.add(model.basin[2], model.linear_resistance[3]) model.edge.add(model.linear_resistance[3], model.basin[4]) model.edge.add(model.basin[4], model.linear_resistance[5]) @@ -686,9 +686,9 @@ def main_network_with_subnetworks_model() -> Model: model.edge.add(model.user_demand[53], model.basin[50]) model.edge.add(model.basin[44], model.user_demand[57]) model.edge.add(model.user_demand[57], model.basin[44]) - model.edge.add(model.basin[2], model.pump[11], subnetwork_id=3) - model.edge.add(model.basin[6], model.pump[24], subnetwork_id=5) - model.edge.add(model.basin[10], model.pump[38], subnetwork_id=7) + model.edge.add(model.basin[2], model.pump[11]) + model.edge.add(model.basin[6], model.pump[24]) + model.edge.add(model.basin[10], model.pump[38]) model.edge.add(model.basin[8], model.user_demand[60]) model.edge.add(model.user_demand[60], model.basin[8]) @@ -708,8 +708,8 @@ def subnetworks_with_sources_model() -> Model: [flow_boundary.Static(flow_rate=[0.003])], ) - model.edge.add(model.flow_boundary[58], model.basin[16], subnetwork_id=3) - model.edge.add(model.flow_boundary[59], model.basin[44], subnetwork_id=7) + model.edge.add(model.flow_boundary[58], model.basin[16]) + model.edge.add(model.flow_boundary[59], model.basin[44]) return model @@ -763,7 +763,7 @@ def level_demand_model() -> Model: [basin.Profile(area=1000.0, level=[0.0, 1.0]), basin.State(level=[2.0])], ) - model.edge.add(model.flow_boundary[1], model.basin[2], subnetwork_id=2) + model.edge.add(model.flow_boundary[1], model.basin[2]) model.edge.add(model.basin[2], model.user_demand[3]) model.edge.add(model.level_demand[4], model.basin[2]) model.edge.add(model.user_demand[3], model.basin[5]) @@ -835,7 +835,6 @@ def flow_demand_model() -> Model: model.edge.add( model.level_boundary[1], model.tabulated_rating_curve[2], - subnetwork_id=2, ) model.edge.add(model.tabulated_rating_curve[2], model.basin[3]) model.edge.add(model.basin[3], model.user_demand[4]) @@ -877,7 +876,7 @@ def linear_resistance_demand_model(): [flow_demand.Static(priority=[1], demand=2.0)], ) - model.edge.add(model.basin[1], model.linear_resistance[2], subnetwork_id=1) + model.edge.add(model.basin[1], model.linear_resistance[2]) model.edge.add(model.linear_resistance[2], model.basin[3]) model.edge.add(model.flow_demand[4], model.linear_resistance[2]) @@ -966,7 +965,7 @@ def fair_distribution_model(): ], ) - model.edge.add(model.level_boundary[1], model.pump[2], subnetwork_id=1) + model.edge.add(model.level_boundary[1], model.pump[2]) model.edge.add(model.pump[2], model.basin[3]) model.edge.add(model.basin[3], model.linear_resistance[4]) model.edge.add(model.linear_resistance[4], model.basin[5]) diff --git a/python/ribasim_testmodels/ribasim_testmodels/invalid.py b/python/ribasim_testmodels/ribasim_testmodels/invalid.py index 049562e6d..1e0657f6e 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/invalid.py +++ b/python/ribasim_testmodels/ribasim_testmodels/invalid.py @@ -214,7 +214,6 @@ def invalid_priorities_model() -> Model: model.edge.add( model.level_boundary[1], model.tabulated_rating_curve[2], - subnetwork_id=2, ) model.edge.add(model.tabulated_rating_curve[2], model.basin[3]) model.edge.add(model.basin[3], model.user_demand[4]) diff --git a/ribasim_qgis/core/nodes.py b/ribasim_qgis/core/nodes.py index 0f6eb5e7e..35d8f3970 100644 --- a/ribasim_qgis/core/nodes.py +++ b/ribasim_qgis/core/nodes.py @@ -254,7 +254,6 @@ def attributes(cls) -> list[QgsField]: QgsField("from_node_id", QVariant.Int), QgsField("to_node_id", QVariant.Int), QgsField("edge_type", QVariant.String), - QgsField("subnetwork_id", QVariant.Int), ] @classmethod