diff --git a/core/src/parameter.jl b/core/src/parameter.jl index e3742808e..e066a591f 100644 --- a/core/src/parameter.jl +++ b/core/src/parameter.jl @@ -321,11 +321,15 @@ Requirements: * fraction must be positive. node_id: node ID of the TabulatedRatingCurve node +inflow_id: node ID across the incoming flow edge +outflow_id: node ID across the outgoing flow edge fraction: The fraction in [0,1] of flow the node lets through control_mapping: dictionary from (node_id, control_state) to fraction """ struct FractionalFlow <: AbstractParameterNode node_id::Vector{NodeID} + inflow_id::Vector{NodeID} + outflow_id::Vector{NodeID} fraction::Vector{Float64} control_mapping::Dict{Tuple{NodeID, String}, NamedTuple} end @@ -354,6 +358,8 @@ end """ node_id: node ID of the Pump node +inflow_id: node ID across the incoming flow edge +outflow_ids: node IDs across the outgoing flow edges active: whether this node is active and thus contributes flow flow_rate: target flow rate min_flow_rate: The minimal flow rate of the pump @@ -363,6 +369,8 @@ is_pid_controlled: whether the flow rate of this pump is governed by PID control """ struct Pump{T} <: AbstractParameterNode node_id::Vector{NodeID} + inflow_id::Vector{NodeID} + outflow_ids::Vector{Vector{NodeID}} active::BitVector flow_rate::T min_flow_rate::Vector{Float64} @@ -372,6 +380,8 @@ struct Pump{T} <: AbstractParameterNode function Pump( node_id, + inflow_id, + outflow_ids, active, flow_rate::T, min_flow_rate, @@ -382,6 +392,8 @@ struct Pump{T} <: AbstractParameterNode if valid_flow_rates(node_id, get_tmp(flow_rate, 0), control_mapping) return new{T}( node_id, + inflow_id, + outflow_ids, active, flow_rate, min_flow_rate, @@ -397,6 +409,8 @@ end """ node_id: node ID of the Outlet node +inflow_id: node ID across the incoming flow edge +outflow_ids: node IDs across the outgoing flow edges active: whether this node is active and thus contributes flow flow_rate: target flow rate min_flow_rate: The minimal flow rate of the outlet @@ -406,6 +420,8 @@ is_pid_controlled: whether the flow rate of this outlet is governed by PID contr """ struct Outlet{T} <: AbstractParameterNode node_id::Vector{NodeID} + inflow_id::Vector{NodeID} + outflow_ids::Vector{Vector{NodeID}} active::BitVector flow_rate::T min_flow_rate::Vector{Float64} @@ -416,6 +432,8 @@ struct Outlet{T} <: AbstractParameterNode function Outlet( node_id, + inflow_id, + outflow_ids, active, flow_rate::T, min_flow_rate, @@ -427,6 +445,8 @@ struct Outlet{T} <: AbstractParameterNode if valid_flow_rates(node_id, get_tmp(flow_rate, 0), control_mapping) return new{T}( node_id, + inflow_id, + outflow_ids, active, flow_rate, min_flow_rate, @@ -503,6 +523,9 @@ struct PidControl{T} <: AbstractParameterNode end """ +node_id: node ID of the UserDemand node +inflow_id: node ID across the incoming flow edge +outflow_id: node ID across the outgoing flow edge active: whether this node is active and thus demands water realized_bmi: Cumulative inflow volume, for read or reset by BMI only demand: water flux demand of UserDemand per priority over time @@ -518,6 +541,8 @@ min_level: The level of the source basin below which the UserDemand does not abs """ struct UserDemand <: AbstractParameterNode node_id::Vector{NodeID} + inflow_id::Vector{NodeID} + outflow_id::Vector{NodeID} active::BitVector realized_bmi::Vector{Float64} demand::Matrix{Float64} @@ -530,6 +555,8 @@ struct UserDemand <: AbstractParameterNode function UserDemand( node_id, + inflow_id, + outflow_id, active, realized_bmi, demand, @@ -544,6 +571,8 @@ struct UserDemand <: AbstractParameterNode if valid_demand(node_id, demand_itp, priorities) return new( node_id, + inflow_id, + outflow_id, active, realized_bmi, demand, diff --git a/core/src/read.jl b/core/src/read.jl index 00a17b895..ab30db99b 100644 --- a/core/src/read.jl +++ b/core/src/read.jl @@ -349,7 +349,7 @@ function ManningResistance(db::DB, config::Config, graph::MetaGraph)::ManningRes ) end -function FractionalFlow(db::DB, config::Config)::FractionalFlow +function FractionalFlow(db::DB, config::Config, graph::MetaGraph)::FractionalFlow static = load_structvector(db, config, FractionalFlowStaticV1) parsed_parameters, valid = parse_static_and_time(db, config, "FractionalFlow"; static) @@ -357,8 +357,12 @@ function FractionalFlow(db::DB, config::Config)::FractionalFlow error("Errors occurred when parsing FractionalFlow data.") end + node_id = NodeID.(NodeType.FractionalFlow, parsed_parameters.node_id) + return FractionalFlow( - NodeID.(NodeType.FractionalFlow, parsed_parameters.node_id), + node_id, + inflow_id.(Ref(graph), node_id), + outflow_id.(Ref(graph), node_id), parsed_parameters.fraction, parsed_parameters.control_mapping, ) @@ -427,7 +431,7 @@ function FlowBoundary(db::DB, config::Config)::FlowBoundary return FlowBoundary(node_ids, parsed_parameters.active, parsed_parameters.flow_rate) end -function Pump(db::DB, config::Config, chunk_sizes::Vector{Int})::Pump +function Pump(db::DB, config::Config, graph::MetaGraph, chunk_sizes::Vector{Int})::Pump static = load_structvector(db, config, PumpStaticV1) defaults = (; min_flow_rate = 0.0, max_flow_rate = Inf, active = true) parsed_parameters, valid = parse_static_and_time(db, config, "Pump"; static, defaults) @@ -444,8 +448,12 @@ function Pump(db::DB, config::Config, chunk_sizes::Vector{Int})::Pump parsed_parameters.flow_rate end + node_id = NodeID.(NodeType.Pump, parsed_parameters.node_id) + return Pump( - NodeID.(NodeType.Pump, parsed_parameters.node_id), + node_id, + inflow_id.(Ref(graph), node_id), + [collect(outflow_ids(graph, id)) for id in node_id], BitVector(parsed_parameters.active), flow_rate, parsed_parameters.min_flow_rate, @@ -455,7 +463,7 @@ function Pump(db::DB, config::Config, chunk_sizes::Vector{Int})::Pump ) end -function Outlet(db::DB, config::Config, chunk_sizes::Vector{Int})::Outlet +function Outlet(db::DB, config::Config, graph::MetaGraph, chunk_sizes::Vector{Int})::Outlet static = load_structvector(db, config, OutletStaticV1) defaults = (; min_flow_rate = 0.0, max_flow_rate = Inf, min_crest_level = -Inf, active = true) @@ -473,8 +481,12 @@ function Outlet(db::DB, config::Config, chunk_sizes::Vector{Int})::Outlet parsed_parameters.flow_rate end + node_id = NodeID.(NodeType.Outlet, parsed_parameters.node_id) + return Outlet( - NodeID.(NodeType.Outlet, parsed_parameters.node_id), + node_id, + inflow_id.(Ref(graph), node_id), + [collect(outflow_ids(graph, id)) for id in node_id], BitVector(parsed_parameters.active), flow_rate, parsed_parameters.min_flow_rate, @@ -806,7 +818,7 @@ function user_demand_time!( return errors end -function UserDemand(db::DB, config::Config)::UserDemand +function UserDemand(db::DB, config::Config, graph::MetaGraph)::UserDemand static = load_structvector(db, config, UserDemandStaticV1) time = load_structvector(db, config, UserDemandTimeV1) @@ -865,6 +877,8 @@ function UserDemand(db::DB, config::Config)::UserDemand return UserDemand( node_ids, + inflow_id.(Ref(graph), node_ids), + outflow_id.(Ref(graph), node_ids), active, realized_bmi, demand, @@ -1045,15 +1059,15 @@ function Parameters(db::DB, config::Config)::Parameters linear_resistance = LinearResistance(db, config, graph) manning_resistance = ManningResistance(db, config, graph) tabulated_rating_curve = TabulatedRatingCurve(db, config) - fractional_flow = FractionalFlow(db, config) + fractional_flow = FractionalFlow(db, config, graph) level_boundary = LevelBoundary(db, config) flow_boundary = FlowBoundary(db, config) - pump = Pump(db, config, chunk_sizes) - outlet = Outlet(db, config, chunk_sizes) + pump = Pump(db, config, graph, chunk_sizes) + outlet = Outlet(db, config, graph, chunk_sizes) terminal = Terminal(db, config) discrete_control = DiscreteControl(db, config) pid_control = PidControl(db, config, chunk_sizes) - user_demand = UserDemand(db, config) + user_demand = UserDemand(db, config, graph) level_demand = LevelDemand(db, config) flow_demand = FlowDemand(db, config) diff --git a/core/src/solve.jl b/core/src/solve.jl index 96ec3b6f5..4d20d44f0 100644 --- a/core/src/solve.jl +++ b/core/src/solve.jl @@ -271,22 +271,29 @@ function formulate_flow!( t::Number, )::Nothing (; graph, basin, allocation) = p - (; + i = 0 + for ( node_id, - allocated, + inflow_id, + outflow_id, active, demand_itp, - demand, return_factor, min_level, demand_from_timeseries, - ) = user_demand - - for (i, id) in enumerate(node_id) - src_id = inflow_id(graph, id) - dst_id = outflow_id(graph, id) - - if !active[i] + ) in zip( + user_demand.node_id, + user_demand.inflow_id, + user_demand.outflow_id, + user_demand.active, + user_demand.demand_itp, + user_demand.return_factor, + user_demand.min_level, + user_demand.demand_from_timeseries, + ) + i += 1 + + if !active continue end @@ -297,31 +304,31 @@ function formulate_flow!( # If allocation is not optimized then allocated = Inf, so the result is always # effectively allocated = demand. for priority_idx in eachindex(allocation.priorities) - alloc_prio = allocated[i, priority_idx] - demand_prio = if demand_from_timeseries[i] - demand_itp[i][priority_idx](t) + alloc_prio = user_demand.allocated[i, priority_idx] + demand_prio = if demand_from_timeseries + demand_itp[priority_idx](t) else - demand[i, priority_idx] + user_demand.demand[i, priority_idx] end alloc = min(alloc_prio, demand_prio) q += alloc end # Smoothly let abstraction go to 0 as the source basin dries out - factor_basin = low_storage_factor(storage, basin.node_id, src_id, 10.0) + factor_basin = low_storage_factor(storage, basin.node_id, inflow_id, 10.0) q *= factor_basin # Smoothly let abstraction go to 0 as the source basin # level reaches its minimum level - source_level = get_level(p, src_id, t; storage) - Δsource_level = source_level - min_level[i] + source_level = get_level(p, inflow_id, t; storage) + Δsource_level = source_level - min_level factor_level = reduction_factor(Δsource_level, 0.1) q *= factor_level - set_flow!(graph, src_id, id, q) + set_flow!(graph, inflow_id, node_id, q) # Return flow is immediate - set_flow!(graph, id, dst_id, q * return_factor[i]) + set_flow!(graph, node_id, outflow_id, q * return_factor) end return nothing end @@ -493,15 +500,17 @@ function formulate_flow!( t::Number, )::Nothing (; graph) = p - (; node_id, fraction) = fractional_flow - for (i, id) in enumerate(node_id) - downstream_id = outflow_id(graph, id) - upstream_id = inflow_id(graph, id) + for (node_id, inflow_id, outflow_id, fraction) in zip( + fractional_flow.node_id, + fractional_flow.inflow_id, + fractional_flow.outflow_id, + fractional_flow.fraction, + ) # overwrite the inflow such that flow is conserved over the FractionalFlow - outflow = get_flow(graph, upstream_id, id, storage) * fraction[i] - set_flow!(graph, upstream_id, id, outflow) - set_flow!(graph, id, downstream_id, outflow) + outflow = get_flow(graph, inflow_id, node_id, storage) * fraction + set_flow!(graph, inflow_id, node_id, outflow) + set_flow!(graph, node_id, outflow_id, outflow) end return nothing end @@ -537,22 +546,27 @@ function formulate_flow!( t::Number, )::Nothing (; graph, basin) = p - (; node_id, active, flow_rate, is_pid_controlled) = pump - flow_rate = get_tmp(flow_rate, storage) - for (id, isactive, rate, pid_controlled) in - zip(node_id, active, flow_rate, is_pid_controlled) - src_id = inflow_id(graph, id) - dst_id = outflow_id(graph, id) - - if !isactive || pid_controlled + + for (node_id, inflow_id, outflow_ids, active, flow_rate, is_pid_controlled) in zip( + pump.node_id, + pump.inflow_id, + pump.outflow_ids, + pump.active, + get_tmp(pump.flow_rate, storage), + pump.is_pid_controlled, + ) + if !active || is_pid_controlled continue end - factor = low_storage_factor(storage, basin.node_id, src_id, 10.0) - q = rate * factor + factor = low_storage_factor(storage, basin.node_id, inflow_id, 10.0) + q = flow_rate * factor - set_flow!(graph, src_id, id, q) - set_flow!(graph, id, dst_id, q) + set_flow!(graph, inflow_id, node_id, q) + + for outflow_id in outflow_ids + set_flow!(graph, node_id, outflow_id, q) + end end return nothing end @@ -564,22 +578,36 @@ function formulate_flow!( t::Number, )::Nothing (; graph, basin) = p - (; node_id, active, flow_rate, is_pid_controlled, min_crest_level) = outlet - flow_rate = get_tmp(flow_rate, storage) - for (i, id) in enumerate(node_id) - src_id = inflow_id(graph, id) - dst_id = outflow_id(graph, id) - if !active[i] || is_pid_controlled[i] + for ( + node_id, + inflow_id, + outflow_ids, + active, + flow_rate, + is_pid_controlled, + min_crest_level, + ) in zip( + outlet.node_id, + outlet.inflow_id, + outlet.outflow_ids, + outlet.active, + get_tmp(outlet.flow_rate, storage), + outlet.is_pid_controlled, + outlet.min_crest_level, + ) + if !active || is_pid_controlled continue end - q = flow_rate[i] - q *= low_storage_factor(storage, basin.node_id, src_id, 10.0) + q = flow_rate + q *= low_storage_factor(storage, basin.node_id, inflow_id, 10.0) # No flow of outlet if source level is lower than target level - src_level = get_level(p, src_id, t; storage) - dst_level = get_level(p, dst_id, t; storage) + src_level = get_level(p, inflow_id, t; storage) + # TODO support multiple outflows to FractionalFlow, or refactor FractionalFlow + outflow_id = only(outflow_ids) + dst_level = get_level(p, outflow_id, t; storage) if src_level !== nothing && dst_level !== nothing Δlevel = src_level - dst_level @@ -588,11 +616,14 @@ function formulate_flow!( # No flow out outlet if source level is lower than minimum crest level if src_level !== nothing - q *= reduction_factor(src_level - min_crest_level[i], 0.1) + q *= reduction_factor(src_level - min_crest_level, 0.1) end - set_flow!(graph, src_id, id, q) - set_flow!(graph, id, dst_id, q) + set_flow!(graph, inflow_id, node_id, q) + + for outflow_id in outflow_ids + set_flow!(graph, node_id, outflow_id, q) + end end return nothing end