Add max_flow_rate to LinearResistance

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/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

@@ -31,6 +31,7 @@
     # 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])
+    # TODO adapt for max_flow_rate of 6e-5
     storage_analytic =
         @. limit_storage + (storage[1] - limit_storage) * exp.(decay_rate * t)
 

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, $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_testmodels/ribasim_testmodels/equations.py

@@ -76,7 +76,9 @@ def linear_resistance_model():
 
     # 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: