Skip to content

Commit

Permalink
Add max_flow_rate to LinearResistance (#1100)
Browse files Browse the repository at this point in the history
This is needed to be able to model structures that exchange flow based
on a water level difference, but only up to a certain capacity.

Ref Deltares/Ribasim-NL#62

The updated test model now looks like this:


![image](https://github.com/Deltares/Ribasim/assets/4471859/c2c88db2-51ce-48c3-a248-51e038711720)
  • Loading branch information
visr authored Feb 9, 2024
1 parent 88f6fde commit 3565646
Show file tree
Hide file tree
Showing 9 changed files with 42 additions and 27 deletions.
4 changes: 3 additions & 1 deletion core/src/parameter.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
5 changes: 4 additions & 1 deletion core/src/read.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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(
Expand All @@ -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
Expand Down
1 change: 1 addition & 0 deletions core/src/schema.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
11 changes: 5 additions & 6 deletions core/src/solve.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
24 changes: 15 additions & 9 deletions core/test/equations_test.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
10 changes: 5 additions & 5 deletions docs/core/equations.qmd
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
1 change: 1 addition & 0 deletions docs/core/usage.qmd
Original file line number Diff line number Diff line change
Expand Up @@ -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}

Expand Down
1 change: 1 addition & 0 deletions python/ribasim/ribasim/schemas.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)


Expand Down
12 changes: 7 additions & 5 deletions python/ribasim_testmodels/ribasim_testmodels/equations.py
Original file line number Diff line number Diff line change
Expand Up @@ -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],
}
)

Expand All @@ -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:
Expand Down

0 comments on commit 3565646

Please sign in to comment.