Merge branch 'main' into remove_unconstrained_sections

core/src/parameter.jl

@@ -162,6 +162,8 @@ end
 """
 struct Basin{T, C, V1, V2, V3} <: AbstractParameterNode
     node_id::Indices{NodeID}
+    inflow_ids::Vector{Vector{NodeID}}
+    outflow_ids::Vector{Vector{NodeID}}
     # Vertical fluxes
     vertical_flux_from_input::V1
     vertical_flux::V2
@@ -182,6 +184,8 @@ struct Basin{T, C, V1, V2, V3} <: AbstractParameterNode
 
     function Basin(
         node_id,
+        inflow_ids,
+        outflow_ids,
         vertical_flux_from_input::V1,
         vertical_flux::V2,
         vertical_flux_prev::V3,
@@ -199,6 +203,8 @@ struct Basin{T, C, V1, V2, V3} <: AbstractParameterNode
         is_valid || error("Invalid Basin / profile table.")
         return new{T, C, V1, V2, V3}(
             node_id,
+            inflow_ids,
+            outflow_ids,
             vertical_flux_from_input,
             vertical_flux,
             vertical_flux_prev,

core/src/read.jl

@@ -93,7 +93,7 @@ function parse_static_and_time(
 
     errors = false
     t_end = seconds_since(config.endtime, config.starttime)
-    trivial_timespan = [nextfloat(-Inf), prevfloat(Inf)]
+    trivial_timespan = [0.0, prevfloat(Inf)]
 
     for (node_idx, node_id) in enumerate(node_ids)
         if node_id in static_node_ids
@@ -490,7 +490,7 @@ function Terminal(db::DB, config::Config)::Terminal
     return Terminal(NodeID.(NodeType.Terminal, static.node_id))
 end
 
-function Basin(db::DB, config::Config, chunk_sizes::Vector{Int})::Basin
+function Basin(db::DB, config::Config, graph::MetaGraph, chunk_sizes::Vector{Int})::Basin
     node_id = get_ids(db, "Basin")
     n = length(node_id)
     current_level = zeros(n)
@@ -533,8 +533,12 @@ function Basin(db::DB, config::Config, chunk_sizes::Vector{Int})::Basin
 
     demand = zeros(length(node_id))
 
+    node_id = NodeID.(NodeType.Basin, node_id)
+
     return Basin(
-        Indices(NodeID.(NodeType.Basin, node_id)),
+        Indices(node_id),
+        [collect(inflow_ids(graph, id)) for id in node_id],
+        [collect(outflow_ids(graph, id)) for id in node_id],
         vertical_flux_from_input,
         vertical_flux,
         vertical_flux_prev,
@@ -820,7 +824,7 @@ function UserDemand(db::DB, config::Config)::UserDemand
     realized_bmi = zeros(n_user)
     demand = zeros(n_user, n_priority)
     demand_reduced = zeros(n_user, n_priority)
-    trivial_timespan = [nextfloat(-Inf), prevfloat(Inf)]
+    trivial_timespan = [0.0, prevfloat(Inf)]
     demand_itp = [
         [LinearInterpolation(zeros(2), trivial_timespan) for i in eachindex(priorities)] for j in eachindex(node_ids)
     ]
@@ -885,6 +889,7 @@ function LevelDemand(db::DB, config::Config)::LevelDemand
         static,
         time,
         time_interpolatables = [:min_level, :max_level],
+        defaults = (; min_level = -Inf, max_level = Inf),
     )
 
     if !valid
@@ -1052,7 +1057,7 @@ function Parameters(db::DB, config::Config)::Parameters
     level_demand = LevelDemand(db, config)
     flow_demand = FlowDemand(db, config)
 
-    basin = Basin(db, config, chunk_sizes)
+    basin = Basin(db, config, graph, chunk_sizes)
     subgrid_level = Subgrid(db, config, basin)
 
     p = Parameters(

core/src/schema.jl

@@ -250,16 +250,16 @@ end
 
 @version LevelDemandStaticV1 begin
     node_id::Int32
-    min_level::Float64
-    max_level::Float64
+    min_level::Union{Missing, Float64}
+    max_level::Union{Missing, Float64}
     priority::Int32
 end
 
 @version LevelDemandTimeV1 begin
     node_id::Int32
     time::DateTime
-    min_level::Float64
-    max_level::Float64
+    min_level::Union{Missing, Float64}
+    max_level::Union{Missing, Float64}
     priority::Int32
 end
 

core/src/solve.jl

@@ -607,10 +607,10 @@ function formulate_du!(
     # subtract all outgoing flows
     # add all ingoing flows
     for (i, basin_id) in enumerate(basin.node_id)
-        for inflow_id in inflow_ids(graph, basin_id)
+        for inflow_id in basin.inflow_ids[i]
             du[i] += get_flow(graph, inflow_id, basin_id, storage)
         end
-        for outflow_id in outflow_ids(graph, basin_id)
+        for outflow_id in basin.outflow_ids[i]
             du[i] -= get_flow(graph, basin_id, outflow_id, storage)
         end
     end

core/src/util.jl

@@ -31,17 +31,10 @@ function get_storage_from_level(basin::Basin, state_idx::Int, level::Float64)::F
     storage_discrete = basin.storage[state_idx]
     area_discrete = basin.area[state_idx]
     level_discrete = basin.level[state_idx]
-    bottom = first(level_discrete)
-
-    if level < bottom
-        node_id = basin.node_id.values[state_idx]
-        @error "The level $level of $node_id is lower than the bottom of this basin; $bottom."
-        return NaN
-    end
 
     level_lower_index = searchsortedlast(level_discrete, level)
 
-    # If the level is equal to the bottom then the storage is 0
+    # If the level is at or below the bottom then the storage is 0
     if level_lower_index == 0
         return 0.0
     end
@@ -77,7 +70,10 @@ function get_storages_from_levels(basin::Basin, levels::Vector)::Vector{Float64}
 
     for (i, level) in enumerate(levels)
         storage = get_storage_from_level(basin, i, level)
-        if isnan(storage)
+        bottom = first(basin.level[i])
+        node_id = basin.node_id.values[i]
+        if level < bottom
+            @error "The initial level ($level) of $node_id is below the bottom ($bottom)."
             errors = true
         end
         storages[i] = storage

core/test/allocation_test.jl

@@ -313,6 +313,8 @@ end
 end
 
 @testitem "Allocation level control" begin
+    import JuMP
+
     toml_path = normpath(@__DIR__, "../../generated_testmodels/level_demand/ribasim.toml")
     @test ispath(toml_path)
     model = Ribasim.run(toml_path)
@@ -377,6 +379,14 @@ end
     stage_6_start_idx = findfirst(stage_6)
     u_stage_6(τ) = storage[stage_6_start_idx]
     @test storage[stage_6] ≈ u_stage_6.(t[stage_6]) rtol = 1e-4
+
+    # Isolated LevelDemand + Basin pair to test optional min_level
+    problem = allocation.allocation_models[2].problem
+    @test JuMP.value(only(problem[:F_basin_in])) == 0.0
+    @test JuMP.value(only(problem[:F_basin_out])) == 0.0
+    q = JuMP.normalized_rhs(only(problem[:basin_outflow]))
+    storage_surplus = 1000.0  # Basin #7 is 1000 m2 and 1 m above LevelDemand max_level
+    @test q ≈ storage_surplus / Δt_allocation
 end
 
 @testitem "Flow demand" begin

core/test/utils_test.jl

@@ -38,6 +38,8 @@ end
     demand = zeros(2)
     basin = Ribasim.Basin(
         Indices(NodeID.(:Basin, [5, 7])),
+        [NodeID[]],
+        [NodeID[]],
         [2.0, 3.0],
         [2.0, 3.0],
         [2.0, 3.0],
@@ -92,6 +94,8 @@ end
     demand = zeros(1)
     basin = Ribasim.Basin(
         Indices(NodeID.(:Basin, [1])),
+        [NodeID[]],
+        [NodeID[]],
         zeros(1),
         zeros(1),
         zeros(1),
@@ -114,14 +118,17 @@ end
     @test length(logger.logs) == 1
     @test logger.logs[1].level == Error
     @test logger.logs[1].message ==
-          "The level -1.0 of Basin #1 is lower than the bottom of this basin; 0.0."
+          "The initial level (-1.0) of Basin #1 is below the bottom (0.0)."
 
     # Converting from storages to levels and back should return the same storages
     storages = range(0.0, 2 * storage[end], 50)
     levels = [Ribasim.get_area_and_level(basin, 1, s)[2] for s in storages]
     storages_ = [Ribasim.get_storage_from_level(basin, 1, l) for l in levels]
-
     @test storages ≈ storages_
+
+    # At or below bottom the storage is 0
+    @test Ribasim.get_storage_from_level(basin, 1, 0.0) == 0.0
+    @test Ribasim.get_storage_from_level(basin, 1, -1.0) == 0.0
 end
 
 @testitem "Expand logic_mapping" begin

docs/core/usage.qmd

@@ -479,14 +479,17 @@ min_level     | Float64  | $m$          | -
 
 A `LevelDemand` node associates a minimum and a maximum level with connected basins to be used by the allocation algorithm.
 Below the minimum level the basin has a demand of the supplied priority,
-above the maximum level the basin acts as a source, used by all nodes with demands in order of priority.
-The same `LevelDemand` node can be used for basins in different subnetworks.
+above the maximum level the basin has a surplus and acts as a source, used by all nodes with demands in order of priority.
+The same `LevelDemand` node can be used for Basins in different subnetworks.
+
+Both `min_level` and `max_level` are optional, to be able to handle only the demand or surplus side.
+If both are missing, `LevelDemand` won't have any effects on allocation.
 
 column        | type    | unit         | restriction
 ------------- | ------- | ------------ | -----------
 node_id       | Int32   | -            | sorted
-min_level     | Float64 | $m$          | -
-max_level     | Float64 | $m$          | -
+min_level     | Float64 | $m$          | (optional, default -Inf)
+max_level     | Float64 | $m$          | (optional, default Inf)
 priority      | Int32   | -            | positive
 
 ## LevelDemand / time

python/ribasim/ribasim/schemas.py

@@ -115,16 +115,16 @@ class LevelBoundaryTimeSchema(_BaseSchema):
 
 class LevelDemandStaticSchema(_BaseSchema):
     node_id: Series[Int32] = pa.Field(nullable=False, default=0)
-    min_level: Series[float] = pa.Field(nullable=False)
-    max_level: Series[float] = pa.Field(nullable=False)
+    min_level: Series[float] = pa.Field(nullable=True)
+    max_level: Series[float] = pa.Field(nullable=True)
     priority: Series[Int32] = pa.Field(nullable=False, default=0)
 
 
 class LevelDemandTimeSchema(_BaseSchema):
     node_id: Series[Int32] = pa.Field(nullable=False, default=0)
     time: Series[Timestamp] = pa.Field(nullable=False)
-    min_level: Series[float] = pa.Field(nullable=False)
-    max_level: Series[float] = pa.Field(nullable=False)
+    min_level: Series[float] = pa.Field(nullable=True)
+    max_level: Series[float] = pa.Field(nullable=True)
     priority: Series[Int32] = pa.Field(nullable=False, default=0)
 
 

python/ribasim_testmodels/ribasim_testmodels/allocation.py

@@ -883,7 +883,7 @@ def subnetworks_with_sources_model() -> Model:
 
 
 def level_demand_model() -> Model:
-    """Small model with a LevelDemand."""
+    """Small model with LevelDemand nodes."""
 
     model = Model(
         starttime="2020-01-01",
@@ -922,12 +922,24 @@ def level_demand_model() -> Model:
         [basin.Profile(area=1000.0, level=[0.0, 1.0]), basin.State(level=[0.5])],
     )
 
+    # Isolated LevelDemand + Basin pair to test optional min_level
+    model.level_demand.add(
+        Node(6, Point(3, -1), subnetwork_id=3),
+        [level_demand.Static(max_level=[1.0], priority=1)],
+    )
+    model.basin.add(
+        Node(7, Point(3, 0), subnetwork_id=3),
+        [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.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])
     model.edge.add(model.level_demand[4], model.basin[5])
 
+    model.edge.add(model.level_demand[6], model.basin[7])
+
     return model