Add max_flow_rate to LinearResistance

core/src/parameter.jl

@@ -197,13 +197,15 @@ Requirements:
 
 node_id: node ID of the LinearResistance node
 active: whether this node is active and thus contributes flows
-resistance: the resistance to flow; Q = Δh/resistance
+resistance: the resistance to flow; `Q_unlimited = Δh/resistance`
+max_flow_rate: the maximum flow rate allowed through the node; `Q = clamp(Q_unlimited, -max_flow_rate, max_flow_rate)`
 control_mapping: dictionary from (node_id, control_state) to resistance and/or active state
 """
 struct LinearResistance <: AbstractParameterNode
     node_id::Vector{NodeID}
     active::BitVector
     resistance::Vector{Float64}
+    max_flow_rate::Vector{Float64}
     control_mapping::Dict{Tuple{NodeID, String}, NamedTuple}
 end
 

core/src/read.jl

@@ -226,7 +226,9 @@ end
 
 function LinearResistance(db::DB, config::Config)::LinearResistance
     static = load_structvector(db, config, LinearResistanceStaticV1)
-    parsed_parameters, valid = parse_static_and_time(db, config, "LinearResistance"; static)
+    defaults = (; max_flow_rate = Inf, active = true)
+    parsed_parameters, valid =
+        parse_static_and_time(db, config, "LinearResistance"; static, defaults)
 
     if !valid
         error(
@@ -238,6 +240,7 @@ function LinearResistance(db::DB, config::Config)::LinearResistance
         NodeID.(parsed_parameters.node_id),
         BitVector(parsed_parameters.active),
         parsed_parameters.resistance,
+        parsed_parameters.max_flow_rate,
         parsed_parameters.control_mapping,
     )
 end

core/src/schema.jl

@@ -164,6 +164,7 @@ end
     node_id::Int
     active::Union{Missing, Bool}
     resistance::Float64
+    max_flow_rate::Union{Missing, Float64}
     control_state::Union{Missing, String}
 end
 

core/src/solve.jl

@@ -315,17 +315,16 @@ function formulate_flow!(
     t::Number,
 )::Nothing
     (; graph) = p
-    (; node_id, active, resistance) = linear_resistance
+    (; node_id, active, resistance, max_flow_rate) = linear_resistance
     for (i, id) in enumerate(node_id)
         basin_a_id = inflow_id(graph, id)
         basin_b_id = outflow_id(graph, id)
 
         if active[i]
-            q =
-                (
-                    get_level(p, basin_a_id, t; storage) -
-                    get_level(p, basin_b_id, t; storage)
-                ) / resistance[i]
+            h_a = get_level(p, basin_a_id, t; storage)
+            h_b = get_level(p, basin_b_id, t; storage)
+            q_unlimited = (h_a - h_b) / resistance[i]
+            q = clamp(q_unlimited, -max_flow_rate[i], max_flow_rate[i])
 
             # add reduction_factor on highest level
             if q > 0

core/test/equations_test.jl

@@ -23,18 +23,24 @@
     @test ispath(toml_path)
     model = Ribasim.run(toml_path)
     @test successful_retcode(model)
-    p = model.integrator.p
+    (; p) = model.integrator
 
     t = Ribasim.timesteps(model)
     storage = Ribasim.get_storages_and_levels(model).storage[1, :]
-    basin_area = p.basin.area[1][2] # Considered constant
-    # The storage of the basin when it has the same level as the level boundary
-    limit_storage = 450.0
-    decay_rate = -1 / (basin_area * p.linear_resistance.resistance[1])
-    storage_analytic =
-        @. limit_storage + (storage[1] - limit_storage) * exp.(decay_rate * t)
-
-    @test all(isapprox.(storage, storage_analytic; rtol = 0.005)) # Fails with '≈'
+    A = p.basin.area[1][2]  # needs to be constant
+    u0 = A * 10.0
+    L = p.level_boundary.level[1].u[1]
+    R = p.linear_resistance.resistance[1]
+    Q_max = p.linear_resistance.max_flow_rate[1]
+
+    # derivation in https://github.com/Deltares/Ribasim/pull/1100#issuecomment-1934799342
+    t_shift = (u0 - A * (L + R * Q_max)) / Q_max
+    pre_shift = t .< t_shift
+    u_pre(t) = u0 - Q_max * t
+    u_post(t) = A * L + A * R * Q_max * exp(-(t - t_shift) / (A * R))
+
+    @test all(isapprox.(storage[pre_shift], u_pre.(t[pre_shift]); rtol = 1e-4))
+    @test all(isapprox.(storage[.~pre_shift], u_post.(t[.~pre_shift]); rtol = 1e-4))
 end
 
 # Equation: storage' = -Q(level(storage)), storage(t0) = storage0,

docs/core/equations.qmd

@@ -271,15 +271,15 @@ of outlets, or even reversal of flows for high precipitation events.
 ### LinearResistance
 
 A `LinearResistance` connects two basins together. The flow between the two basins
-is determined by a linear relationship:
+is determined by a linear relationship, up to an optional maximum flow rate:
 
 $$
-    Q = \frac{h_a - h_b}{R}
+    Q = \mathrm{clamp}(\frac{h_a - h_b}{R}, -Q_{\max}, Q_{\max})
 $$ {#eq-basinflow}
 
-Here $h_a$ is the water level in the first basin, $h_b$ is the water level in
-the second basin, and $R$ is the resistance of the link. A `LinearResistance` makes
-no assumptions about the direction of the flow: water flows from high to low.
+Here $h_a$ is the water level in the first basin and $h_b$ is the water level in the second basin.
+$R$ is the resistance of the link, and $Q_{\max}$ is the maximum flow rate.
+A `LinearResistance` makes no assumptions about the direction of the flow: water flows from high to low.
 
 ### Terminal
 

docs/core/usage.qmd

@@ -580,6 +580,7 @@ node_id       | Int     | -            | sorted
 control_state | String  | -            | (optional) sorted per node_id
 active        | Bool    | -            | (optional, default true)
 resistance    | Float64 | $sm^{-2}$    | -
+max_flow_rate | Float64 | $m^3 s^{-1}$ | non-negative
 
 # ManningResistance {#sec-manning-resistance}
 

python/ribasim/ribasim/schemas.py

@@ -105,6 +105,7 @@ class LinearResistanceStaticSchema(_BaseSchema):
     node_id: Series[int] = pa.Field(nullable=False)
     active: Series[pa.BOOL] = pa.Field(nullable=True)
     resistance: Series[float] = pa.Field(nullable=False)
+    max_flow_rate: Series[float] = pa.Field(nullable=True)
     control_state: Series[str] = pa.Field(nullable=True)
 
 

python/ribasim_testmodels/ribasim_testmodels/equations.py

@@ -48,9 +48,9 @@ def linear_resistance_model():
     # Setup the basins:
     profile = pd.DataFrame(
         data={
-            "node_id": [1, 1, 1],
-            "area": [0.01, 100.0, 100.0],
-            "level": [0.0, 1.0, 2.0],
+            "node_id": [1, 1],
+            "area": [100.0, 100.0],
+            "level": [0.0, 10.0],
         }
     )
 
@@ -68,15 +68,17 @@ def linear_resistance_model():
     state = pd.DataFrame(
         data={
             "node_id": [1],
-            "level": [10.5],
+            "level": [10.0],
         }
     )
 
     basin = ribasim.Basin(profile=profile, static=static, state=state)
 
     # setup linear resistance:
     linear_resistance = ribasim.LinearResistance(
-        static=pd.DataFrame(data={"node_id": [2], "resistance": [5e4]})
+        static=pd.DataFrame(
+            data={"node_id": [2], "resistance": [5e4], "max_flow_rate": [6e-5]}
+        )
     )
 
     # Setup level boundary: