From f6401e2ca88298c2a49b1321227e2e516a247dca Mon Sep 17 00:00:00 2001
From: Bart de Koning <bart.dekoning@deltares.nl>
Date: Fri, 25 Aug 2023 14:15:11 +0200
Subject: [PATCH 01/22] Implement hard stop of outlet flow

---
 core/src/jac.jl   | 13 +++++++++++++
 core/src/solve.jl | 40 ++++++++++++++++++++++++++++++++--------
 2 files changed, 45 insertions(+), 8 deletions(-)

diff --git a/core/src/jac.jl b/core/src/jac.jl
index c474bd2cd..9f0497db0 100644
--- a/core/src/jac.jl
+++ b/core/src/jac.jl
@@ -388,6 +388,19 @@ function formulate_jac!(
         controlled_node_id = only(outneighbors(graph_control, id))
         controls_pump = insorted(controlled_node_id, pump.node_id)
 
+        # No flow of outlet if source level is lower than target level
+        if !controls_pump
+            src_id = only(inneighbors(graph_flow, controlled_node_id))
+            dst_id = only(outneighbors(graph_flow, controlled_node_id))
+
+            src_level = get_level(p, src_id, t)
+            dst_level = get_level(p, dst_id, t)
+
+            if (src_level < dst_level)
+                continue
+            end
+        end
+
         if !controls_pump
             if !insorted(controlled_node_id, outlet.node_id)
                 error(
diff --git a/core/src/solve.jl b/core/src/solve.jl
index e61e68698..146896075 100644
--- a/core/src/solve.jl
+++ b/core/src/solve.jl
@@ -557,7 +557,7 @@ function continuous_control!(
         if !active[i]
             du.integral[i] = 0.0
             u.integral[i] = 0.0
-            return
+            continue
         end
 
         du.integral[i] = error[i]
@@ -568,6 +568,19 @@ function continuous_control!(
         controlled_node_id = only(outneighbors(graph_control, id))
         controls_pump = (controlled_node_id in pump.node_id)
 
+        # No flow of outlet if source level is lower than target level
+        if !controls_pump
+            src_id = only(inneighbors(graph_flow, controlled_node_id))
+            dst_id = only(outneighbors(graph_flow, controlled_node_id))
+
+            src_level = get_level(p, src_id, t)
+            dst_level = get_level(p, dst_id, t)
+
+            if (src_level < dst_level)
+                continue
+            end
+        end
+
         if controls_pump
             controlled_node_idx = findsorted(pump.node_id, controlled_node_id)
 
@@ -858,22 +871,33 @@ function formulate!(
     node::Union{Pump, Outlet},
     p::Parameters,
     storage::AbstractVector{Float64},
+    t::Float64,
 )::Nothing
     (; connectivity, basin) = p
     (; graph_flow, flow) = connectivity
     (; node_id, active, flow_rate, is_pid_controlled) = node
+    is_outlet = nameof(typeof(node)) == :Outlet
+
     for (id, isactive, rate, pid_controlled) in
         zip(node_id, active, flow_rate, is_pid_controlled)
         src_id = only(inneighbors(graph_flow, id))
         dst_id = only(outneighbors(graph_flow, id))
 
-        if !isactive
-            flow[src_id, id] = 0.0
-            flow[id, dst_id] = 0.0
-            continue
+        # No flow of outlet if source level is lower than target level
+        non_flowing_outlet = false
+        if is_outlet
+            src_level = get_level(p, src_id, t)
+            dst_level = get_level(p, dst_id, t)
+
+            if (src_level < dst_level)
+                non_flowing_outlet = true
+            end
         end
 
-        if pid_controlled
+        # No flow of outlet if source level is lower than target level
+        if !isactive || pid_controlled || non_flowing_outlet
+            flow[src_id, id] = 0.0
+            flow[id, dst_id] = 0.0
             continue
         end
 
@@ -935,8 +959,8 @@ function formulate_flows!(
     formulate!(tabulated_rating_curve, p, t)
     formulate!(flow_boundary, p, t)
     formulate!(fractional_flow, p)
-    formulate!(pump, p, storage)
-    formulate!(outlet, p, storage)
+    formulate!(pump, p, storage, t)
+    formulate!(outlet, p, storage, t)
 
     return nothing
 end

From 6ad246088c676b0681d7d93da315d06c03e24ff2 Mon Sep 17 00:00:00 2001
From: Bart de Koning <bart.dekoning@deltares.nl>
Date: Fri, 1 Sep 2023 14:55:41 +0200
Subject: [PATCH 02/22] Add smooth reduction factor

---
 core/src/solve.jl | 94 +++++++++++++++++++++++++++++++++--------------
 core/src/utils.jl | 21 ++++++++++-
 2 files changed, 86 insertions(+), 29 deletions(-)

diff --git a/core/src/solve.jl b/core/src/solve.jl
index 146896075..712f40a5a 100644
--- a/core/src/solve.jl
+++ b/core/src/solve.jl
@@ -509,12 +509,12 @@ function formulate!(
         bottom = basin.level[i][1]
         fixed_area = basin.area[i][end]
         depth = max(level - bottom, 0.0)
-        reduction_factor = min(depth, 0.1) / 0.1
+        reduction_factor_basin = reduction_factor(depth, 0.1)
 
         precipitation = fixed_area * basin.precipitation[i]
-        evaporation = area * reduction_factor * basin.potential_evaporation[i]
+        evaporation = area * reduction_factor_basin * basin.potential_evaporation[i]
         drainage = basin.drainage[i]
-        infiltration = reduction_factor * basin.infiltration[i]
+        infiltration = reduction_factor_basin * basin.infiltration[i]
 
         du.storage[i] += precipitation - evaporation + drainage - infiltration
     end
@@ -576,16 +576,21 @@ function continuous_control!(
             src_level = get_level(p, src_id, t)
             dst_level = get_level(p, dst_id, t)
 
-            if (src_level < dst_level)
-                continue
+            if !isnothing(src_level) && !isnothing(dst_level)
+                Δlevel = src_level - dst_level
+                reduction_factor_outlet = reduction_factor(Δlevel, 0.1)
+            else
+                reduction_factor_outlet = 1.0
             end
+        else
+            reduction_factor_outlet = 1.0
         end
 
         if controls_pump
             controlled_node_idx = findsorted(pump.node_id, controlled_node_id)
 
             listened_basin_storage = u.storage[listened_node_idx]
-            reduction_factor = min(listened_basin_storage, 10.0) / 10.0
+            reduction_factor_basin = reduction_factor(listened_basin_storage, 10.0)
         else
             controlled_node_idx = findsorted(outlet.node_id, controlled_node_id)
 
@@ -594,12 +599,13 @@ function continuous_control!(
             has_index, upstream_basin_idx = id_index(basin.node_id, upstream_node_id)
             if has_index
                 upstream_basin_storage = u.storage[upstream_basin_idx]
-                reduction_factor = min(upstream_basin_storage, 10.0) / 10.0
+                reduction_factor_basin = reduction_factor(upstream_basin_storage, 10.0)
             else
-                reduction_factor = 1.0
+                reduction_factor_basin = 1.0
             end
         end
 
+        reduction_factor = reduction_factor_basin * reduction_factor_outlet
         flow_rate = 0.0
 
         K_p, K_i, K_d = pid_params[i](t)
@@ -868,34 +874,59 @@ function formulate!(flow_boundary::FlowBoundary, p::Parameters, t::Float64)::Not
 end
 
 function formulate!(
-    node::Union{Pump, Outlet},
+    pump::Pump,
     p::Parameters,
     storage::AbstractVector{Float64},
     t::Float64,
 )::Nothing
     (; connectivity, basin) = p
     (; graph_flow, flow) = connectivity
-    (; node_id, active, flow_rate, is_pid_controlled) = node
-    is_outlet = nameof(typeof(node)) == :Outlet
+    (; node_id, active, flow_rate, is_pid_controlled) = pump
 
     for (id, isactive, rate, pid_controlled) in
         zip(node_id, active, flow_rate, is_pid_controlled)
         src_id = only(inneighbors(graph_flow, id))
         dst_id = only(outneighbors(graph_flow, id))
 
-        # No flow of outlet if source level is lower than target level
-        non_flowing_outlet = false
-        if is_outlet
-            src_level = get_level(p, src_id, t)
-            dst_level = get_level(p, dst_id, t)
+        if !isactive || pid_controlled
+            flow[src_id, id] = 0.0
+            flow[id, dst_id] = 0.0
+            continue
+        end
 
-            if (src_level < dst_level)
-                non_flowing_outlet = true
-            end
+        hasindex, basin_idx = id_index(basin.node_id, src_id)
+
+        q = rate
+
+        if hasindex
+            # Pumping from basin
+            s = storage[basin_idx]
+            reduction_factor_basin = reduction_factor(s, 10.0)
+            q = reduction_factor_basin * rate
         end
 
-        # No flow of outlet if source level is lower than target level
-        if !isactive || pid_controlled || non_flowing_outlet
+        flow[src_id, id] = q
+        flow[id, dst_id] = q
+    end
+    return nothing
+end
+
+function formulate!(
+    outlet::Outlet,
+    p::Parameters,
+    storage::AbstractVector{Float64},
+    t::Float64,
+)::Nothing
+    (; connectivity, basin) = p
+    (; graph_flow, flow) = connectivity
+    (; node_id, active, flow_rate, is_pid_controlled) = outlet
+
+    for (id, isactive, rate, pid_controlled) in
+        zip(node_id, active, flow_rate, is_pid_controlled)
+        src_id = only(inneighbors(graph_flow, id))
+        dst_id = only(outneighbors(graph_flow, id))
+
+        if !isactive || pid_controlled
             flow[src_id, id] = 0.0
             flow[id, dst_id] = 0.0
             continue
@@ -903,14 +934,23 @@ function formulate!(
 
         hasindex, basin_idx = id_index(basin.node_id, src_id)
 
+        q = rate
+
         if hasindex
-            # Pumping from basin
+            # Flowing from basin
             s = storage[basin_idx]
-            reduction_factor = min(s, 10.0) / 10.0
-            q = reduction_factor * rate
-        else
-            # Pumping from level boundary
-            q = rate
+            reduction_factor_basin = reduction_factor(s, 10.0)
+            q *= reduction_factor_basin
+        end
+
+        # No flow of outlet if source level is lower than target level
+        src_level = get_level(p, src_id, t)
+        dst_level = get_level(p, dst_id, t)
+
+        if !isnothing(src_level) && !isnothing(dst_level)
+            Δlevel = src_level - dst_level
+            reduction_factor_outlet = reduction_factor(Δlevel, 0.1)
+            q *= reduction_factor_outlet
         end
 
         flow[src_id, id] = q
diff --git a/core/src/utils.jl b/core/src/utils.jl
index cdb7b4360..3b8793210 100644
--- a/core/src/utils.jl
+++ b/core/src/utils.jl
@@ -410,7 +410,7 @@ end
 Get the current water level of a node ID.
 The ID can belong to either a Basin or a LevelBoundary.
 """
-function get_level(p::Parameters, node_id::Int, t::Float64)::Float64
+function get_level(p::Parameters, node_id::Int, t::Float64)::Union{Float64, Nothing}
     (; basin, level_boundary) = p
     # since the node_id fields are already Indices, Dictionary creation is instant
     basin = Dictionary(basin.node_id, basin.current_level)
@@ -419,7 +419,8 @@ function get_level(p::Parameters, node_id::Int, t::Float64)::Float64
         gettokenvalue(basin, token)
     else
         boundary = Dictionary(level_boundary.node_id, level_boundary.level)
-        boundary[node_id](t)
+        hasindex, token = gettoken(boundary, node_id, nothing)
+        return hasindex ? gettokenvalue(boundary, token)(t) : nothing
     end
 end
 
@@ -899,3 +900,19 @@ function Base.getindex(fv::FlatVector, i::Int)
     v = fv.v[d + 1]
     return v[r + 1]
 end
+
+"""
+Function that goes smoothly from 0 to 1 in the interval [0,threshold],
+and is constant outside this interval.
+"""
+function reduction_factor(x::Float64, threshold::Float64)::Float64
+    return if x < 0
+        0.0
+    elseif x < threshold
+        x_scaled = x / threshold
+        x_scaled_sq = x_scaled^2
+        -2 * x_scaled_sq * x_scaled + 3 * x_scaled_sq
+    else
+        1.0
+    end
+end

From f4a53b292e334f7a963524d43c206a005505e766 Mon Sep 17 00:00:00 2001
From: Bart de Koning <bart.dekoning@deltares.nl>
Date: Fri, 1 Sep 2023 16:59:27 +0200
Subject: [PATCH 03/22] Fix merge issues

---
 core/src/solve.jl | 28 +++++++++++++++-------------
 core/src/utils.jl | 14 +++++++++-----
 2 files changed, 24 insertions(+), 18 deletions(-)

diff --git a/core/src/solve.jl b/core/src/solve.jl
index a161e5450..1bf04f63a 100644
--- a/core/src/solve.jl
+++ b/core/src/solve.jl
@@ -595,8 +595,8 @@ function continuous_control!(
             src_id = only(inneighbors(graph_flow, controlled_node_id))
             dst_id = only(outneighbors(graph_flow, controlled_node_id))
 
-            src_level = get_level(p, src_id, t)
-            dst_level = get_level(p, dst_id, t)
+            src_level = get_level(p, src_id, current_level, t)
+            dst_level = get_level(p, dst_id, current_level, t)
 
             if !isnothing(src_level) && !isnothing(dst_level)
                 Δlevel = src_level - dst_level
@@ -627,7 +627,7 @@ function continuous_control!(
             end
         end
 
-        reduction_factor = reduction_factor_basin * reduction_factor_outlet
+        reduction_factor_all = reduction_factor_basin * reduction_factor_outlet
         flow_rate = 0.0
 
         K_p, K_i, K_d = pid_params[i](t)
@@ -635,17 +635,17 @@ function continuous_control!(
         if !iszero(K_d)
             # dlevel/dstorage = 1/area
             area = current_area[listened_node_idx]
-            D = 1.0 - K_d * reduction_factor / area
+            D = 1.0 - K_d * reduction_factor_all / area
         else
             D = 1.0
         end
 
         if !iszero(K_p)
-            flow_rate += reduction_factor * K_p * pid_error[i] / D
+            flow_rate += reduction_factor_all * K_p * pid_error[i] / D
         end
 
         if !iszero(K_i)
-            flow_rate += reduction_factor * K_i * integral_value[i] / D
+            flow_rate += reduction_factor_all * K_i * integral_value[i] / D
         end
 
         if !iszero(K_d)
@@ -914,7 +914,7 @@ end
 function formulate!(pump::Pump, p::Parameters, flow, storage::AbstractVector)::Nothing
     (; connectivity, basin) = p
     (; graph_flow) = connectivity
-    (; node_id, active, flow_rate, is_pid_controlled) = node
+    (; 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)
@@ -947,13 +947,15 @@ end
 function formulate!(
     outlet::Outlet,
     p::Parameters,
-    storage::AbstractVector{Float64},
+    flow,
+    current_level,
+    storage::AbstractVector,
     t::Float64,
 )::Nothing
     (; connectivity, basin) = p
-    (; graph_flow, flow) = connectivity
+    (; graph_flow) = connectivity
     (; node_id, active, flow_rate, is_pid_controlled) = outlet
-
+    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 = only(inneighbors(graph_flow, id))
@@ -977,8 +979,8 @@ function formulate!(
         end
 
         # No flow of outlet if source level is lower than target level
-        src_level = get_level(p, src_id, t)
-        dst_level = get_level(p, dst_id, t)
+        src_level = get_level(p, src_id, current_level, t)
+        dst_level = get_level(p, dst_id, current_level, t)
 
         if !isnothing(src_level) && !isnothing(dst_level)
             Δlevel = src_level - dst_level
@@ -1036,7 +1038,7 @@ function formulate_flows!(
     formulate!(flow_boundary, p, flow, t)
     formulate!(fractional_flow, flow, p)
     formulate!(pump, p, flow, storage)
-    formulate!(outlet, p, flow, storage)
+    formulate!(outlet, p, flow, current_level, storage, t)
 
     return nothing
 end
diff --git a/core/src/utils.jl b/core/src/utils.jl
index 2f68af8a6..99d178635 100644
--- a/core/src/utils.jl
+++ b/core/src/utils.jl
@@ -406,7 +406,12 @@ end
 Get the current water level of a node ID.
 The ID can belong to either a Basin or a LevelBoundary.
 """
-function get_level(p::Parameters, node_id::Int, current_level, t::Float64)::Real
+function get_level(
+    p::Parameters,
+    node_id::Int,
+    current_level,
+    t::Float64,
+)::Union{Real, Nothing}
     (; basin, level_boundary) = p
     # since the node_id fields are already Indices, Dictionary creation is instant
     basin = Dictionary(basin.node_id, current_level)
@@ -415,7 +420,7 @@ function get_level(p::Parameters, node_id::Int, current_level, t::Float64)::Real
         gettokenvalue(basin, token)
     else
         boundary = Dictionary(level_boundary.node_id, level_boundary.level)
-        hasindex, token = gettoken(boundary, node_id, nothing)
+        hasindex, token = gettoken(boundary, node_id)
         return hasindex ? gettokenvalue(boundary, token)(t) : nothing
     end
 end
@@ -901,13 +906,12 @@ end
 Function that goes smoothly from 0 to 1 in the interval [0,threshold],
 and is constant outside this interval.
 """
-function reduction_factor(x::Float64, threshold::Float64)::Float64
+function reduction_factor(x::Real, threshold::Real)::Real
     return if x < 0
         0.0
     elseif x < threshold
         x_scaled = x / threshold
-        x_scaled_sq = x_scaled^2
-        -2 * x_scaled_sq * x_scaled + 3 * x_scaled_sq
+        (-2 * x_scaled + 3) * x_scaled^2
     else
         1.0
     end

From b67112cd0f148bbdb4215e58f64a1c557bc89086 Mon Sep 17 00:00:00 2001
From: Bart de Koning <bart.dekoning@deltares.nl>
Date: Fri, 1 Sep 2023 17:12:55 +0200
Subject: [PATCH 04/22] Adjust some failing tests

---
 core/test/control.jl                                        | 4 ++--
 python/ribasim_testmodels/ribasim_testmodels/pid_control.py | 2 +-
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/core/test/control.jl b/core/test/control.jl
index 4480099fe..51a30414a 100644
--- a/core/test/control.jl
+++ b/core/test/control.jl
@@ -188,6 +188,6 @@ end
     t_target_jump = discrete_control.record.time[2]
     t_idx_target_jump = searchsortedlast(timesteps, t_target_jump)
 
-    @test isapprox(level[t_idx_target_jump], target_high, atol = 1e-4)
-    @test isapprox(level[end], target_low, atol = 1e-2)
+    @test isapprox(level[t_idx_target_jump], target_high, atol = 1e-1)
+    @test isapprox(level[end], target_low, atol = 1e-1)
 end
diff --git a/python/ribasim_testmodels/ribasim_testmodels/pid_control.py b/python/ribasim_testmodels/ribasim_testmodels/pid_control.py
index 7e7082435..333f4536c 100644
--- a/python/ribasim_testmodels/ribasim_testmodels/pid_control.py
+++ b/python/ribasim_testmodels/ribasim_testmodels/pid_control.py
@@ -290,7 +290,7 @@ def discrete_control_of_pid_control_model():
                 "node_id": [6, 6],
                 "control_state": ["target_high", "target_low"],
                 "listen_node_id": [3, 3],
-                "target": [6.0, 4.0],
+                "target": [5.0, 3.0],
                 "proportional": 2 * [1e-2],
                 "integral": 2 * [1e-8],
                 "derivative": 2 * [-1e-1],

From d8b0951e2a89288769edb9fbee6e13cb6082b4dd Mon Sep 17 00:00:00 2001
From: Bart de Koning <bart.dekoning@deltares.nl>
Date: Fri, 1 Sep 2023 17:17:29 +0200
Subject: [PATCH 05/22] Adjust other failing test

---
 python/ribasim_testmodels/ribasim_testmodels/pid_control.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/python/ribasim_testmodels/ribasim_testmodels/pid_control.py b/python/ribasim_testmodels/ribasim_testmodels/pid_control.py
index 333f4536c..b8fc718d9 100644
--- a/python/ribasim_testmodels/ribasim_testmodels/pid_control.py
+++ b/python/ribasim_testmodels/ribasim_testmodels/pid_control.py
@@ -113,7 +113,7 @@ def pid_control_model():
         static=pd.DataFrame(
             data={
                 "node_id": [4],
-                "level": [1.0],  # Not relevant
+                "level": [5.0],  # Not relevant
             }
         )
     )

From 0e5006e5f90c8f25b2a925bae38f945327f2be50 Mon Sep 17 00:00:00 2001
From: Bart de Koning <bart.dekoning@deltares.nl>
Date: Fri, 1 Sep 2023 17:19:12 +0200
Subject: [PATCH 06/22] Remove PID control Jacobian contributions section

---
 docs/core/equations.qmd | 75 -----------------------------------------
 1 file changed, 75 deletions(-)

diff --git a/docs/core/equations.qmd b/docs/core/equations.qmd
index 1c29c6ffc..f67d3a494 100644
--- a/docs/core/equations.qmd
+++ b/docs/core/equations.qmd
@@ -465,81 +465,6 @@ Q_\text{pump/outlet} = \text{clip}\left(\phi(u_\text{us})\frac{K_pe + K_iI + K_d
 $$
 where the clipping is again done last. Note that to compute this, $\hat{f}_\text{PID}$ has to be known first, meaning that the PID controlled pump/outlet flow rate has to be computed after all other contributions to the PID controlled basin's storage are known.
 
-## Jacobian contributions
-
-For convenience here we denote a time derivative of a variable as a dot over the symbol.
-
-For the integral state we simply get the contribution
-$$
-\frac{\partial \dot{I}}{\partial u_\text{PID}} = - \frac{1}{A(u_\text{PID})}.
-$$
-
-Note that when the calculated pump flow lies outside $[Q_\min, Q_\max]$, the pump flow is locally constant and thus all partial derivatives of $Q_\text{pump/outlet}$ with respect to states are $0$. Assuming that it lies in $[Q_\min, Q_\max]$, we derive the Jacobian contributions. Define the enumerator and denominator of the fraction in $Q_\text{pump/outlet}$:
-
-$$
-\begin{align}
-    E &= K_pe + K_iI + K_d \left(\dot{\text{SP}}-\frac{\hat{f}_\text{PID}}{A(u_\text{PID})}\right) \\
-    D &= 1 \pm \phi(u_\text{us})\frac{K_d}{A(u_\text{PID})}
-\end{align}
-$$
-
-Then
-
-$$
-\frac{\partial Q_\text{pump/outlet}}{\partial I} = K_i\frac{\phi(u_\text{us})}{D}.
-$$
-
-For the derivative of $Q_\text{pump/outlet}$ with respect to storages we need the derivatives of the enumerator and denominator.
-
-For the denominator we distinguish 2 possibilities:
-
-- If the controlled node is a pump, then $u_\text{us} \equiv u_\text{PID}$ and so $D$ is only a function of $u_\text{PID}$:
-$$
-\frac{\text{d}D}{\text{d}u_\text{PID}} = - K_d\left[\frac{\phi'(u_\text{PID})}{A(u_\text{PID})}- \phi(u_\text{PID})\frac{A'(u_\text{PID})}{A(u_\text{PID})^2}\right]
-$$
-- If the controlled node is an outlet, then $D$ is a function of two distinct storages (assuming the upstream node is indeed a basin):
-$$
-\begin{align}
-    \frac{\partial D}{\partial u_\text{PID}} &= \phi(u_\text{us})K_d\frac{A'(u_\text{PID})}{A(u_\text{PID})^2}, \\
-    \frac{\partial D}{\partial u_\text{us}} &= -\phi'(u_\text{us})K_d\frac{1}{A(u_\text{PID})}.
-\end{align}
-$$
-
-For the enumerator there distinguish 2 cases:
-
-- The partial derivative with respect to $u_\text{PID}$:
-$$
-\frac{\partial E}{\partial u_\text{PID}} = -\frac{K_p}{A(u_\text{PID})} - K_d\frac{A(u_\text{PID})\frac{\partial\hat{f}_\text{PID}}{\partial u_\text{PID}} - \hat{f}_\text{PID}A'(u_\text{PID})}{A(u_\text{PID})^2}.
-$$
-- The partial derivative with respect to a different storage $u_i$:
-$$
-\frac{\partial E}{\partial u_i} = - \frac{K_d}{A(u_\text{PID})}\frac{\partial \hat{f}_\text{PID}}{\partial u_i}.
-$$
-
-For computational efficiency we exploit that $\phi'$ is mostly $0$. Note that
-
-$$
-\frac{\partial\hat{f}_\text{PID}}{\partial u_i} = \hat{J}[i,\text{PID}],
-$$
-\emph{i.e.} the Jacobian term of the dependence of $u_\text{PID}$ on $u_i$ without the contribution of the PID controlled pump/outlet. Note thus that also the Jacobian contribution of the PID controlled pump can only be computed after the other contributions to the Jacobian are known.
-
-The partial derivatives of $Q_\text{pump}$ are then given by
-$$
-\begin{align}
-    \frac{\partial Q_\text{pump}}{\partial u_\text{PID}} &= \phi'(u_\text{PID})\frac{E}{D} + \frac{\phi(u_\text{PID})}{D^2}\left[D\frac{\partial E}{\partial u_\text{PID}} - E\frac{\text{d}D}{\text{d}u_\text{PID}}\right], \\
-    \frac{\partial Q_\text{pump}}{\partial u_i} &= \frac{\phi(u_\text{PID})}{D}\frac{\partial E}{\partial u_i} \qquad (i \ne \text{PID}).
-\end{align}
-$$
-
-The partial derivatives of $Q_\text{outlet}$ are given by
-$$
-\begin{align}
-    \frac{\partial Q_\text{outlet}}{\partial u_\text{PID}} &= \frac{\phi(u_\text{us})}{D^2}\left[D\frac{\partial E}{\partial u_\text{PID}} - E\frac{\partial D}{\partial u_\text{PID}}\right], \\
-    \frac{\partial Q_\text{outlet}}{\partial u_\text{us}} &= \phi'(u_\text{us})\frac{E}{D} + \frac{\phi(u_\text{us})}{D^2}\left[D\frac{\partial E}{\partial u_\text{us}} - E\frac{\partial D}{\partial u_\text{us}}\right], \\
-    \frac{\partial Q_\text{outlet}}{\partial u_i} &= \frac{\phi(u_\text{us})}{D}\frac{\partial E}{\partial u_i} \qquad (i \ne \text{PID, us}).
-\end{align}
-$$
-
 ## The sign of the parameters
 
 Note by @eq-error that the error is positive if the setpoint is larger than the basin level and negative if the setpoint is smaller than the basin level.

From 66c7cebf258028915103855599d3f2ee164788b4 Mon Sep 17 00:00:00 2001
From: Bart de Koning <bart.dekoning@deltares.nl>
Date: Mon, 4 Sep 2023 13:46:38 +0200
Subject: [PATCH 07/22] Add min crest level

---
 core/src/create.jl     | 4 +++-
 core/src/validation.jl | 1 +
 2 files changed, 4 insertions(+), 1 deletion(-)

diff --git a/core/src/create.jl b/core/src/create.jl
index 8221f4504..05e4dccce 100644
--- a/core/src/create.jl
+++ b/core/src/create.jl
@@ -432,7 +432,8 @@ end
 
 function Outlet(db::DB, config::Config, chunk_size::Int)::Outlet
     static = load_structvector(db, config, OutletStaticV1)
-    defaults = (; min_flow_rate = 0.0, max_flow_rate = NaN, active = true)
+    defaults =
+        (; min_flow_rate = 0.0, max_flow_rate = NaN, min_crest_level = NaN, active = true)
     parsed_parameters, valid = parse_static_and_time(db, config, "Outlet"; static, defaults)
     is_pid_controlled = falses(length(parsed_parameters.node_id))
 
@@ -453,6 +454,7 @@ function Outlet(db::DB, config::Config, chunk_size::Int)::Outlet
         flow_rate,
         parsed_parameters.min_flow_rate,
         parsed_parameters.max_flow_rate,
+        parsed_parameters.min_crest_level,
         parsed_parameters.control_mapping,
         is_pid_controlled,
     )
diff --git a/core/src/validation.jl b/core/src/validation.jl
index 2eb1919d4..cafabd4db 100644
--- a/core/src/validation.jl
+++ b/core/src/validation.jl
@@ -141,6 +141,7 @@ end
     flow_rate::Float64
     min_flow_rate::Union{Missing, Float64}
     max_flow_rate::Union{Missing, Float64}
+    min_crest_level::Union{Missing, Float64}
     control_state::Union{Missing, String}
 end
 

From 491bbd408d25f6ee1866ea08f010727d443f4540 Mon Sep 17 00:00:00 2001
From: Bart de Koning <bart.dekoning@deltares.nl>
Date: Mon, 4 Sep 2023 13:47:56 +0200
Subject: [PATCH 08/22] Process min crest level

---
 core/src/solve.jl                | 19 +++++++++++++------
 docs/core/usage.qmd              | 20 +++++++++++++++++++-
 python/ribasim/ribasim/models.py |  3 ++-
 3 files changed, 34 insertions(+), 8 deletions(-)

diff --git a/core/src/solve.jl b/core/src/solve.jl
index 1bf04f63a..78e00c50e 100644
--- a/core/src/solve.jl
+++ b/core/src/solve.jl
@@ -318,6 +318,7 @@ struct Outlet{T} <: AbstractParameterNode
     flow_rate::T
     min_flow_rate::Vector{Float64}
     max_flow_rate::Vector{Float64}
+    min_crest_level::Vector{Float64}
     control_mapping::Dict{Tuple{Int, String}, NamedTuple}
     is_pid_controlled::BitVector
 
@@ -501,7 +502,7 @@ function set_current_basin_properties!(
 end
 
 """
-Linearize the evaporation flux when at small water depths
+Smoothly let the evaporation flux go to 0 when at small water depths
 Currently at less than 0.1 m.
 """
 function formulate!(
@@ -954,14 +955,13 @@ function formulate!(
 )::Nothing
     (; connectivity, basin) = p
     (; graph_flow) = connectivity
-    (; node_id, active, flow_rate, is_pid_controlled) = outlet
+    (; node_id, active, flow_rate, is_pid_controlled, min_crest_level) = outlet
     flow_rate = get_tmp(flow_rate, storage)
-    for (id, isactive, rate, pid_controlled) in
-        zip(node_id, active, flow_rate, is_pid_controlled)
+    for (i, id) in enumerate(node_id)
         src_id = only(inneighbors(graph_flow, id))
         dst_id = only(outneighbors(graph_flow, id))
 
-        if !isactive || pid_controlled
+        if !active[i] || is_pid_controlled[i]
             flow[src_id, id] = 0.0
             flow[id, dst_id] = 0.0
             continue
@@ -988,6 +988,13 @@ function formulate!(
             q *= reduction_factor_outlet
         end
 
+        # No flow out outlet if source level is lower than minimum crest level
+        if !isnothing(src_level) && !isnan(min_crest_level[i])
+            reduction_factor_min_crest_level =
+                reduction_factor(src_level - min_crest_level[i], 0.1)
+            q *= reduction_factor_min_crest_level
+        end
+
         flow[src_id, id] = q
         flow[id, dst_id] = q
     end
@@ -1091,7 +1098,7 @@ function water_balance!(
         t,
     )
 
-    # Negative storage musn't decrease, based on Shampine's et. al. advice
+    # Negative storage must not decrease, based on Shampine's et. al. advice
     # https://docs.sciml.ai/DiffEqCallbacks/stable/step_control/#DiffEqCallbacks.PositiveDomain
     for i in eachindex(u.storage)
         if u.storage[i] < 0
diff --git a/docs/core/usage.qmd b/docs/core/usage.qmd
index 83afbdc67..2f6bae071 100644
--- a/docs/core/usage.qmd
+++ b/docs/core/usage.qmd
@@ -151,6 +151,8 @@ name it must have in the GeoPackage if it is stored there.
   - `TabulatedRatingCurve / time`: dynamic rating curve
 - Pump: pump water from a source node to a destination node
   - `Pump / static`: flow rate
+- Outlet: let water flow with a prescribed flux under the conditions of positive head difference and the upstream level being higher than the minimum crest level
+  - `Outlet / static`: flow rate, minimum crest level
 - Terminal: Water sink without state or properties
   - `Terminal / static`: - (only node IDs)
 - DiscreteControl: Set parameters of other nodes based on model state conditions (e.g. basin level)
@@ -311,10 +313,26 @@ discharge | Float64  | $m^3 s^{-1}$ | non-negative
 
 ### Pump and outlet
 
-Pump/conduct water from a source node to a destination node.
+Pump water from a source node to a destination node.
 The set flow rate will be pumped unless the intake storage is less than $10~m^3$,
 in which case the flow rate will be linearly reduced to $0~m^3/s$.
 The intake must be either a Basin or LevelBoundary.
+When PID controlled, the pump must point away from the controlled basin in terms of edges.
+
+column                | type    | unit         | restriction
+---------             | ------- | ------------ | -----------
+node_id               | Int     | -            | sorted
+active                | Bool    | -            | (optional, default true)
+flow_rate             | Float64 | $m^3 s^{-1}$ | non-negative
+min_flow_rate         | Float64 | $m^3 s^{-1}$ | (optional, default 0.0)
+max_flow_rate         | Float64 | $m^3 s^{-1}$ | (optional)
+min_crest_level       | Float64 | $m$          | (optional)
+control_state         | String  | -            | (optional)
+
+### Outlet
+
+The outlet is very similar to the pump. The outlet has two additional physical constraints: water only flows trough the outlet when the head difference is positive (i.e. water flows down by gravity) and the upstream level must be above the minimum crest level, if the upstream level is defined.
+When PID controlled, the outlet must point towards the controlled basin in terms of edges.
 
 column                | type    | unit         | restriction
 ---------             | ------- | ------------ | -----------
diff --git a/python/ribasim/ribasim/models.py b/python/ribasim/ribasim/models.py
index da6f31259..9fb2ab244 100644
--- a/python/ribasim/ribasim/models.py
+++ b/python/ribasim/ribasim/models.py
@@ -1,6 +1,6 @@
 # generated by datamodel-codegen:
 #   filename:  root.schema.json
-#   timestamp: 2023-08-21T14:05:19+00:00
+#   timestamp: 2023-09-04T11:45:02+00:00
 
 from __future__ import annotations
 
@@ -130,6 +130,7 @@ class OutletStatic(BaseModel):
     max_flow_rate: Optional[float] = Field(None, description="max_flow_rate")
     remarks: Optional[str] = Field("", description="a hack for pandera")
     active: Optional[bool] = Field(None, description="active")
+    min_crest_level: Optional[float] = Field(None, description="min_crest_level")
     flow_rate: float = Field(..., description="flow_rate")
     node_id: int = Field(..., description="node_id")
     control_state: Optional[str] = Field(None, description="control_state")

From 8b19210fdfd95c2d0ba29e047cc1f2425011dd70 Mon Sep 17 00:00:00 2001
From: Bart de Koning <bart.dekoning@deltares.nl>
Date: Mon, 4 Sep 2023 13:48:42 +0200
Subject: [PATCH 09/22] Update equation docs

---
 docs/core/equations.qmd              | 63 +++++++++++++++++-----------
 docs/schema/OutletStatic.schema.json |  7 ++++
 2 files changed, 45 insertions(+), 25 deletions(-)

diff --git a/docs/core/equations.qmd b/docs/core/equations.qmd
index f67d3a494..894c7a469 100644
--- a/docs/core/equations.qmd
+++ b/docs/core/equations.qmd
@@ -90,6 +90,25 @@ The presence of division by the basin area means that areas of size zero are not
 :::
 
 # Natural water balance terms
+
+## The reduction factor {#sec-reduction_factor}
+At several points in the equations below a *reduction factor* is used. This is a term that makes certain transitions more smooth, for instance when a pump stops providing water when its source basin dries up. The reduction factor is given by
+
+$$
+    \phi(x; p) =
+    \begin{align}
+        \begin{cases}
+        0 &\text{if}\quad x < 0 \\
+            -2 \left(\frac{x}{p}\right)^3 + 3\left(\frac{x}{p}\right)^2 &\text{if}\quad 0 \le x \le p \\
+            1 &\text{if}\quad x > p
+        \end{cases}
+    \end{align},
+$$
+
+where $p > 0$ is the threshold value which determines the interval $[0,p]$ of the smooth transition between $0$ and $1$, see the plot below.
+
+<!-- Add reduction factor plot here-->
+
 ## Precipitation
 
 The precipitation term is given by
@@ -107,14 +126,10 @@ surface may change as well, depending on the slopes of the surface waters.
 The evaporation term is given by
 
 $$
-    Q_E = E_\text{pot} \cdot A(S) \cdot r.
+    Q_E = E_\text{pot} \cdot A(S) \cdot \phi(d;0.1).
 $$ {#eq-evap}
 
-Here $E_\text{pot} = E_\text{pot}(t)$ is the potential evaporation rate and $A$ is the wetted area. $r$ is the reduction factor which depends on the depth $d$ and is given by
-$$
-    r = \frac{\min\{d,0.1\}}{0.1} \le 1.
-$$ {#eq-reduc}
-It provides a smooth gradient as $S \rightarrow 0$ (but not at $d=0.1$).
+Here $E_\text{pot} = E_\text{pot}(t)$ is the potential evaporation rate and $A$ is the wetted area. $\phi$ is the [reduction factor](equations.qmd#sec-reduction_factor) which depends on the depth $d$. It provides a smooth gradient as $S \rightarrow 0$.
 
 A straightforward formulation $Q_E = \mathrm{max}(E_\text{pot} A(S),
 0)$ is unsuitable, as $\frac{\mathrm{d}Q_E}{\mathrm{d}S}(S=0)$ is then not well-defined.
@@ -183,26 +198,24 @@ basins (external nodes) in a network.
 2. `LevelBoundary`
 3. `FlowBoundary`
 
-### Pump and outlet
+### Pump {#sec-pump}
 
-The behaviour of pumps and outlets is very straight forward if these nodes are not PID controlled. Their flow is given by a fixed flow rate $q$, multiplied by a reduction factor:
+The behaviour of pumps is very straight forward if these nodes are not PID controlled. Their flow is given by a fixed flow rate $q$, multiplied by a reduction factor:
 $$
-Q_\text{pump/outlet} = \phi(u)q
+Q_\text{pump} = \phi(u; 10.0)q
 $$
 
-Here $u$ is the storage of the upstream basin. The reduction factor
+Here $u$ is the storage of the upstream basin. The [reduction factor](equations.qmd#sec-reduction_factor) $\phi$ makes sure that the flow of the pump goes smootly to $0$ as the upstream basin dries out.
+
+### Outlet {#sec-outlet}
 
+The outlet is very similar to the pump, but it has a few extra [reduction factors](equations.qmd#sec-reduction_factor) for physical constraints:
 $$
-\phi(u) =
-\begin{align}
-    \begin{cases}
-        \frac{u}{10} &\text{if}& 0 \leq u \le 10 \\
-        1 &\text{if}& u \geq 10
-    \end{cases}
-\end{align}
-$$ {#eq-pumpreductionfactor}
+Q_\text{outlet} = \phi(u_a; 10.0)\phi(\Delta h; 0.1) \phi(h_a-h_\text{min};0.1)q.
+$$
+The subscript $a$ denotes the upstream node and $b$ the downstream node. The first reduction factor is equivalent to the one for the pump. The second one makes sure that the outlet flow goes to zero as the head difference $\Delta h = h_a - h_b$ goes to zero. The last one makes sure that the outlet only produces flow when the upstream level is above the minimum chrest level $h_\text{min}$.
 
-makes sure that the flow of the pump goes to $0$ as the upstream basin dries out.
+Not all node types upstream or downstream of the outlet have a defined level. If this is the case, and therefore the reduction factor cannot be computed, it is defined to be $1.0$.
 
 ### TabulatedRatingCurve
 
@@ -415,19 +428,19 @@ $$ {#eq-PIDflow}
 
 for given constant parameters $K_p,K_i,K_d$. The pump or outlet can have associated minimum and maximum flow rates $Q_\min, Q_\max$, and so
 $$
-Q_\text{pump/outlet} = \text{clip}(\phi(u_\text{us})Q_\text{PID}; Q_\min, Q_\max).
+Q_\text{pump/outlet} = \text{clip}(\Phi Q_\text{PID}; Q_\min, Q_\max).
 $$
 
-Here $u_\text{us}$ is the storage of the basin upstream of the pump or outlet, $\phi$ is the reduction factor (see @eq-pumpreductionfactor) and
+Here $u_\text{us}$ is the storage of the basin upstream of the pump or outlet, $\Phi$ is the product of [reduction factors](equations.qmd#sec-reduction_factor) associated with the [pump](equations.qmd#sec-pump) or [outlet](equations.qmd#sec-outlet) and
 
 $$
 \text{clip}(Q; Q_\min, Q_\max) =
 \begin{align}
     \begin{cases}
-        Q_\min & \text{if} & Q < Q_\min, \\
-        Q & \text{if} & Q_\min \leq Q \leq Q_\max, \\
-        Q_\max & \text{if} & Q > Q_\max.
-    \end{cases}
+        Q_\min & \text{if} & Q < Q_\min \\
+        Q & \text{if} & Q_\min \leq Q \leq Q_\max \\
+        Q_\max & \text{if} & Q > Q_\max
+    \end{cases}.
 \end{align}
 $$
 
diff --git a/docs/schema/OutletStatic.schema.json b/docs/schema/OutletStatic.schema.json
index d307169fc..af03e5539 100644
--- a/docs/schema/OutletStatic.schema.json
+++ b/docs/schema/OutletStatic.schema.json
@@ -21,6 +21,13 @@
                 "boolean"
             ]
         },
+        "min_crest_level": {
+            "format": "default",
+            "description": "min_crest_level",
+            "type": [
+                "number"
+            ]
+        },
         "flow_rate": {
             "format": "double",
             "description": "flow_rate",

From efeb4c44c2d96be5e282435f1155ac6cfbe204c6 Mon Sep 17 00:00:00 2001
From: Bart de Koning <bart.dekoning@deltares.nl>
Date: Mon, 4 Sep 2023 14:00:31 +0200
Subject: [PATCH 10/22] Pass all tests

---
 core/src/solve.jl       | 4 +++-
 core/test/validation.jl | 1 +
 2 files changed, 4 insertions(+), 1 deletion(-)

diff --git a/core/src/solve.jl b/core/src/solve.jl
index 78e00c50e..1adcf4c68 100644
--- a/core/src/solve.jl
+++ b/core/src/solve.jl
@@ -328,6 +328,7 @@ struct Outlet{T} <: AbstractParameterNode
         flow_rate::T,
         min_flow_rate,
         max_flow_rate,
+        min_crest_level,
         control_mapping,
         is_pid_controlled,
     ) where {T}
@@ -338,6 +339,7 @@ struct Outlet{T} <: AbstractParameterNode
                 flow_rate,
                 min_flow_rate,
                 max_flow_rate,
+                min_crest_level,
                 control_mapping,
                 is_pid_controlled,
             )
@@ -969,7 +971,7 @@ function formulate!(
 
         hasindex, basin_idx = id_index(basin.node_id, src_id)
 
-        q = rate
+        q = flow_rate[i]
 
         if hasindex
             # Flowing from basin
diff --git a/core/test/validation.jl b/core/test/validation.jl
index 64720ea41..f34c66184 100644
--- a/core/test/validation.jl
+++ b/core/test/validation.jl
@@ -230,6 +230,7 @@ end
             [-1.0],
             [NaN],
             [NaN],
+            [NaN],
             Dict{Tuple{Int, String}, NamedTuple}(),
             [false],
         )

From 5e36dc488d2f2ff0484736c4e96cc962c17452d1 Mon Sep 17 00:00:00 2001
From: Bart de Koning <bart.dekoning@deltares.nl>
Date: Mon, 4 Sep 2023 19:59:38 +0200
Subject: [PATCH 11/22] Add tests

---
 core/test/run_models.jl                       |  27 +++++
 docs/core/usage.qmd                           |   6 +-
 .../ribasim_testmodels/__init__.py            |   2 +
 .../ribasim_testmodels/basic.py               | 107 +++++++++++++++++-
 4 files changed, 138 insertions(+), 4 deletions(-)

diff --git a/core/test/run_models.jl b/core/test/run_models.jl
index e95197086..411e5c035 100644
--- a/core/test/run_models.jl
+++ b/core/test/run_models.jl
@@ -133,6 +133,33 @@ end
     @test A ≈ 10 * h
 end
 
+@testset "Outlet constraints" begin
+    toml_path = normpath(@__DIR__, "../../data/outlet/outlet.toml")
+    @test ispath(toml_path)
+
+    model = Ribasim.run(toml_path)
+    p = model.integrator.p
+    (; level_boundary, outlet) = p
+    (; level) = level_boundary
+    level = level[1]
+
+    timesteps = model.saved_flow.t
+    outlet_flow = [saveval[1] for saveval in model.saved_flow.saveval]
+
+    t_min_crest_level =
+        level.t[2] * (outlet.min_crest_level[1] - level.u[1]) / (level.u[2] - level.u[1])
+
+    # No outlet flow when upstream level is below minimum crest level
+    @test all(@. outlet_flow[timesteps <= t_min_crest_level] == 0)
+
+    timesteps = Ribasim.timesteps(model)
+    t_maximum_level = level.t[2]
+    level_basin = Ribasim.get_storages_and_levels(model).level[:]
+
+    # Basin level converges to stable level boundary level
+    all(isapprox.(level_basin[timesteps .>= t_maximum_level], level.u[3], atol = 5e-2))
+end
+
 @testset "ManningResistance" begin
     """
     Apply the "standard step method" finite difference method to find a
diff --git a/docs/core/usage.qmd b/docs/core/usage.qmd
index 2f6bae071..526103d24 100644
--- a/docs/core/usage.qmd
+++ b/docs/core/usage.qmd
@@ -311,7 +311,7 @@ level     | Float64  | $m$          | -
 discharge | Float64  | $m^3 s^{-1}$ | non-negative
 
 
-### Pump and outlet
+### Pump
 
 Pump water from a source node to a destination node.
 The set flow rate will be pumped unless the intake storage is less than $10~m^3$,
@@ -326,12 +326,11 @@ active                | Bool    | -            | (optional, default true)
 flow_rate             | Float64 | $m^3 s^{-1}$ | non-negative
 min_flow_rate         | Float64 | $m^3 s^{-1}$ | (optional, default 0.0)
 max_flow_rate         | Float64 | $m^3 s^{-1}$ | (optional)
-min_crest_level       | Float64 | $m$          | (optional)
 control_state         | String  | -            | (optional)
 
 ### Outlet
 
-The outlet is very similar to the pump. The outlet has two additional physical constraints: water only flows trough the outlet when the head difference is positive (i.e. water flows down by gravity) and the upstream level must be above the minimum crest level, if the upstream level is defined.
+The outlet is very similar to the pump. The outlet has two additional physical constraints: water only flows trough the outlet when the head difference is positive (i.e. water flows down by gravity), and the upstream level must be above the minimum crest level if the upstream level is defined.
 When PID controlled, the outlet must point towards the controlled basin in terms of edges.
 
 column                | type    | unit         | restriction
@@ -341,6 +340,7 @@ active                | Bool    | -            | (optional, default true)
 flow_rate             | Float64 | $m^3 s^{-1}$ | non-negative
 min_flow_rate         | Float64 | $m^3 s^{-1}$ | (optional, default 0.0)
 max_flow_rate         | Float64 | $m^3 s^{-1}$ | (optional)
+min_crest_level       | Float64 | $m$          | (optional)
 control_state         | String  | -            | (optional)
 
 ### LevelBoundary
diff --git a/python/ribasim_testmodels/ribasim_testmodels/__init__.py b/python/ribasim_testmodels/ribasim_testmodels/__init__.py
index 49cad3bdf..2f15f4f77 100644
--- a/python/ribasim_testmodels/ribasim_testmodels/__init__.py
+++ b/python/ribasim_testmodels/ribasim_testmodels/__init__.py
@@ -4,6 +4,7 @@
 from ribasim_testmodels.basic import (
     basic_model,
     basic_transient_model,
+    outlet_model,
     tabulated_rating_curve_model,
 )
 from ribasim_testmodels.bucket import bucket_model
@@ -62,4 +63,5 @@
     "invalid_edge_types_model",
     "discrete_control_of_pid_control_model",
     "level_boundary_condition_model",
+    "outlet_model",
 ]
diff --git a/python/ribasim_testmodels/ribasim_testmodels/basic.py b/python/ribasim_testmodels/ribasim_testmodels/basic.py
index 8bbb72ccc..9da4892c2 100644
--- a/python/ribasim_testmodels/ribasim_testmodels/basic.py
+++ b/python/ribasim_testmodels/ribasim_testmodels/basic.py
@@ -5,7 +5,7 @@
 
 
 def basic_model() -> ribasim.Model:
-    """Set up a basic model with all node types and static forcing"""
+    """Set up a basic model with most node types and static forcing"""
 
     # Setup the basins:
     profile = pd.DataFrame(
@@ -379,3 +379,108 @@ def tabulated_rating_curve_model() -> ribasim.Model:
     )
 
     return model
+
+
+def outlet_model():
+    """Set up a basic model with an outlet that encounters various physical constraints."""
+
+    # Set up the nodes:
+    xy = np.array(
+        [
+            (0.0, 0.0),  # 1: LevelBoundary
+            (1.0, 0.0),  # 2: Outlet
+            (2.0, 0.0),  # 3: Basin
+        ]
+    )
+    node_xy = gpd.points_from_xy(x=xy[:, 0], y=xy[:, 1])
+
+    node_type = ["LevelBoundary", "Outlet", "Basin"]
+
+    # Make sure the feature id starts at 1: explicitly give an index.
+    node = ribasim.Node(
+        static=gpd.GeoDataFrame(
+            data={"type": node_type},
+            index=pd.Index(np.arange(len(xy)) + 1, name="fid"),
+            geometry=node_xy,
+            crs="EPSG:28992",
+        )
+    )
+
+    # Setup the edges:
+    from_id = np.array([1, 2], dtype=np.int64)
+    to_id = np.array([2, 3], dtype=np.int64)
+    lines = ribasim.utils.geometry_from_connectivity(node, from_id, to_id)
+    edge = ribasim.Edge(
+        static=gpd.GeoDataFrame(
+            data={
+                "from_node_id": from_id,
+                "to_node_id": to_id,
+                "edge_type": len(from_id) * ["flow"],
+            },
+            geometry=lines,
+            crs="EPSG:28992",
+        )
+    )
+
+    # Setup the basins:
+    profile = pd.DataFrame(
+        data={
+            "node_id": 3,
+            "area": 1000.0,
+            "level": [0.0, 1.0],
+        }
+    )
+
+    static = pd.DataFrame(
+        data={
+            "node_id": [3],
+            "drainage": 0.0,
+            "potential_evaporation": 0.0,
+            "infiltration": 0.0,
+            "precipitation": 0.0,
+            "urban_runoff": 0.0,
+        }
+    )
+
+    state = pd.DataFrame(data={"node_id": [3], "level": 1e-3})
+
+    basin = ribasim.Basin(profile=profile, static=static, state=state)
+
+    # Setup the level boundary:
+    level_boundary = ribasim.LevelBoundary(
+        time=pd.DataFrame(
+            data={
+                "node_id": 1,
+                "time": [
+                    "2020-01-01 00:00:00",
+                    "2020-06-01 00:00:00",
+                    "2021-01-01 00:00:00",
+                ],
+                "level": [1.0, 3.0, 3.0],
+            }
+        )
+    )
+
+    # Setup the outlet
+    outlet = ribasim.Outlet(
+        static=pd.DataFrame(
+            data={
+                "node_id": [2],
+                "flow_rate": 1e-3,
+                "min_crest_level": 2.0,
+            }
+        )
+    )
+
+    model = ribasim.Model(
+        modelname="outlet",
+        node=node,
+        edge=edge,
+        basin=basin,
+        outlet=outlet,
+        level_boundary=level_boundary,
+        starttime="2020-01-01 00:00:00",
+        endtime="2021-01-01 00:00:00",
+    )
+
+    return model

From 65ea1680697522e19e99d549f659a45e8d3c4ce3 Mon Sep 17 00:00:00 2001
From: Bart de Koning <bart.dekoning@deltares.nl>
Date: Tue, 5 Sep 2023 13:54:49 +0200
Subject: [PATCH 12/22] Add plot to docs with source code

---
 docs/core/equations.qmd | 39 ++++++++++++++++++++++++++++++++++++++-
 1 file changed, 38 insertions(+), 1 deletion(-)

diff --git a/docs/core/equations.qmd b/docs/core/equations.qmd
index 894c7a469..f2c884959 100644
--- a/docs/core/equations.qmd
+++ b/docs/core/equations.qmd
@@ -107,7 +107,44 @@ $$
 
 where $p > 0$ is the threshold value which determines the interval $[0,p]$ of the smooth transition between $0$ and $1$, see the plot below.
 
-<!-- Add reduction factor plot here-->
+```{python}
+# | code-fold: true
+import numpy as np
+import matplotlib.pyplot as plt
+
+def f(x, p = 3):
+    x_scaled = x / p
+    phi = (-2 * x_scaled + 3) * x_scaled**2
+    phi = np.where(x < 0, 0, phi)
+    phi = np.where(x > p, 1, phi)
+
+    return phi
+
+fontsize = 15
+p = 3
+N = 100
+x_min = -1
+x_max = 4
+x = np.linspace(x_min,x_max,N)
+phi = f(x,p)
+
+fig,ax = plt.subplots(dpi=80)
+ax.plot(x,phi)
+
+y_lim = ax.get_ylim()
+
+ax.set_xticks([0,p], [0,"$p$"], fontsize=fontsize)
+ax.set_yticks([0,1], [0,1], fontsize=fontsize)
+ax.hlines([0,1],x_min,x_max, color = "k", ls = ":", zorder=-1)
+ax.vlines([0,p], *y_lim, color = "k", ls = ":")
+ax.set_xlim(x_min,x_max)
+ax.set_xlabel("$x$", fontsize=fontsize)
+ax.set_ylabel("$\phi(x;p)$", fontsize=fontsize)
+ax.set_ylim(y_lim)
+
+fig.tight_layout()
+plt.show()
+```
 
 ## Precipitation
 

From ebd9dbb71e1806aaf654ec69d0adb688ed22574b Mon Sep 17 00:00:00 2001
From: Martijn Visser <mgvisser@gmail.com>
Date: Tue, 12 Sep 2023 17:06:37 +0200
Subject: [PATCH 13/22] use formulate_flow! function name

---
 core/src/solve.jl | 9 +++++++--
 1 file changed, 7 insertions(+), 2 deletions(-)

diff --git a/core/src/solve.jl b/core/src/solve.jl
index 8d9e53842..77d3e42a5 100644
--- a/core/src/solve.jl
+++ b/core/src/solve.jl
@@ -930,7 +930,12 @@ function formulate_flow!(
     end
 end
 
-function formulate!(pump::Pump, p::Parameters, flow::AbstractMatrix, storage::AbstractVector)::Nothing
+function formulate_flow!(
+    pump::Pump,
+    p::Parameters,
+    flow::AbstractMatrix,
+    storage::AbstractVector,
+)::Nothing
     (; connectivity, basin) = p
     (; graph_flow) = connectivity
     (; node_id, active, flow_rate, is_pid_controlled) = pump
@@ -963,7 +968,7 @@ function formulate!(pump::Pump, p::Parameters, flow::AbstractMatrix, storage::Ab
     return nothing
 end
 
-function formulate!(
+function formulate_flow!(
     outlet::Outlet,
     p::Parameters,
     flow::AbstractMatrix,

From e4c67ff741ed747f9f0358772a251d924bf77b13 Mon Sep 17 00:00:00 2001
From: Martijn Visser <mgvisser@gmail.com>
Date: Wed, 13 Sep 2023 12:04:53 +0200
Subject: [PATCH 14/22] apply reduction factor with *= like other functions

---
 core/src/solve.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/core/src/solve.jl b/core/src/solve.jl
index 77d3e42a5..decc6d841 100644
--- a/core/src/solve.jl
+++ b/core/src/solve.jl
@@ -959,7 +959,7 @@ function formulate_flow!(
             # Pumping from basin
             s = storage[basin_idx]
             reduction_factor_basin = reduction_factor(s, 10.0)
-            q = reduction_factor_basin * rate
+            q *= reduction_factor_basin
         end
 
         flow[src_id, id] = q

From 464f33901f2d8b5e81ebf1f1763632817b13326e Mon Sep 17 00:00:00 2001
From: Martijn Visser <mgvisser@gmail.com>
Date: Wed, 13 Sep 2023 12:10:33 +0200
Subject: [PATCH 15/22] convert isnothing(x) to x === nothing

This is easier on the compiler, generally recommended for performance.
---
 build/libribasim/src/libribasim.jl |  6 +++---
 core/src/bmi.jl                    |  6 +++---
 core/src/config.jl                 |  8 ++++----
 core/src/create.jl                 |  4 ++--
 core/src/io.jl                     |  6 +++---
 core/src/solve.jl                  | 10 +++++-----
 core/src/utils.jl                  | 10 +++++-----
 core/test/libribasim.jl            |  6 +++---
 8 files changed, 28 insertions(+), 28 deletions(-)

diff --git a/build/libribasim/src/libribasim.jl b/build/libribasim/src/libribasim.jl
index d4551e159..ccaf3fb54 100644
--- a/build/libribasim/src/libribasim.jl
+++ b/build/libribasim/src/libribasim.jl
@@ -27,7 +27,7 @@ This expands to
 ```
 try
     global model
-    isnothing(model) && error("Model not initialized")
+    model === nothing && error("Model not initialized")
     BMI.update(model)
 catch
     Base.invokelatest(Base.display_error, current_exceptions())
@@ -40,7 +40,7 @@ macro try_c(ex)
     return quote
         try
             global model
-            isnothing(model) && error("Model not initialized")
+            model === nothing && error("Model not initialized")
             $(esc(ex))
         catch e
             global last_error_message
@@ -81,7 +81,7 @@ end
 
 Base.@ccallable function finalize()::Cint
     @try_c_uninitialized begin
-        if !isnothing(model)
+        if model !== nothing
             BMI.finalize(model)
         end
         model = nothing
diff --git a/core/src/bmi.jl b/core/src/bmi.jl
index 2890ead2c..7aa8b5497 100644
--- a/core/src/bmi.jl
+++ b/core/src/bmi.jl
@@ -53,7 +53,7 @@ function BMI.initialize(T::Type{Model}, config::Config)::Model
         for id in pid_control.node_id
             id_controlled = only(outneighbors(connectivity.graph_control, id))
             pump_idx = findsorted(pump.node_id, id_controlled)
-            if isnothing(pump_idx)
+            if pump_idx === nothing
                 outlet_idx = findsorted(outlet.node_id, id_controlled)
                 outlet.is_pid_controlled[outlet_idx] = true
             else
@@ -247,7 +247,7 @@ function get_value(
 
         if hasindex_basin
             _, level = get_area_and_level(basin, basin_idx, u[basin_idx])
-        elseif !isnothing(level_boundary_idx)
+        elseif level_boundary_idx !== nothing
             level = level_boundary.level[level_boundary_idx](t + Δt)
         else
             error(
@@ -260,7 +260,7 @@ function get_value(
     elseif variable == "flow_rate"
         flow_boundary_idx = findsorted(flow_boundary.node_id, feature_id)
 
-        if isnothing(flow_boundary_idx)
+        if flow_boundary_idx === nothing
             error("Flow condition node #$feature_id is not a flow boundary.")
         end
 
diff --git a/core/src/config.jl b/core/src/config.jl
index fee738f45..6073a2187 100644
--- a/core/src/config.jl
+++ b/core/src/config.jl
@@ -84,7 +84,7 @@ end
 
 function Base.convert(::Type{Compression}, str::AbstractString)
     i = findfirst(==(Symbol(str)) ∘ Symbol, instances(Compression))
-    if isnothing(i)
+    if i === nothing
         throw(
             ArgumentError(
                 "Compression algorithm $str not supported, choose one of: $(join(instances(Compression), " ")).",
@@ -149,7 +149,7 @@ function Base.show(io::IO, c::Config)
     println(io, "Ribasim Config")
     for field in fieldnames(typeof(c))
         f = getfield(c, field)
-        isnothing(f) || println(io, "\t$field\t= $f")
+        f === nothing || println(io, "\t$field\t= $f")
     end
 end
 
@@ -157,7 +157,7 @@ function Base.show(io::IO, c::TableOption)
     first = true
     for field in fieldnames(typeof(c))
         f = getfield(c, field)
-        if !isnothing(f)
+        if f !== nothing
             first && (first = false; println(io))
             println(io, "\t\t$field\t= $f")
         end
@@ -180,7 +180,7 @@ const algorithms = Dict{String, Type}(
 "Create an OrdinaryDiffEqAlgorithm from solver config"
 function algorithm(solver::Solver)::OrdinaryDiffEqAlgorithm
     algotype = get(algorithms, solver.algorithm, nothing)
-    if isnothing(algotype)
+    if algotype === nothing
         options = join(keys(algorithms), ", ")
         error("Given solver algorithm $(solver.algorithm) not supported.\n\
             Available options are: ($(options)).")
diff --git a/core/src/create.jl b/core/src/create.jl
index 94e1907d6..7988fb6f3 100644
--- a/core/src/create.jl
+++ b/core/src/create.jl
@@ -74,14 +74,14 @@ function parse_static_and_time(
     end
 
     # Get node IDs of static nodes if the static table exists
-    static_node_ids = if isnothing(static)
+    static_node_ids = if static === nothing
         Set{Int}()
     else
         Set(static.node_id)
     end
 
     # Get node IDs of transient nodes if the time table exists
-    time_node_ids = if isnothing(time)
+    time_node_ids = if time === nothing
         Set{Int}()
     else
         Set(time.node_id)
diff --git a/core/src/io.jl b/core/src/io.jl
index 4be67f2db..afe50db34 100644
--- a/core/src/io.jl
+++ b/core/src/io.jl
@@ -50,7 +50,7 @@ function load_data(
     path = getfield(getfield(config, snake_case(node)), kind)
     sqltable = tablename(schema)
 
-    table = if !isnothing(path)
+    table = if path !== nothing
         table_path = input_path(config, path)
         Table(read(table_path))
     elseif exists(db, sqltable)
@@ -76,7 +76,7 @@ function load_structvector(
 )::StructVector{T} where {T <: AbstractRow}
     table = load_data(db, config, T)
 
-    if isnothing(table)
+    if table === nothing
         return StructVector{T}(undef, 0)
     end
 
@@ -98,7 +98,7 @@ function load_structvector(
     table = StructVector{T}(nt)
     sv = Legolas._schema_version_from_record_type(T)
     tableschema = Tables.schema(table)
-    if declared(sv) && !isnothing(tableschema)
+    if declared(sv) && tableschema !== nothing
         validate(tableschema, sv)
         # R = Legolas.record_type(sv)
         # foreach(R, Tables.rows(table))  # construct each row
diff --git a/core/src/solve.jl b/core/src/solve.jl
index decc6d841..baa514d31 100644
--- a/core/src/solve.jl
+++ b/core/src/solve.jl
@@ -601,11 +601,11 @@ function continuous_control!(
             src_level = get_level(p, src_id, current_level, t)
             dst_level = get_level(p, dst_id, current_level, t)
 
-            if !isnothing(src_level) && !isnothing(dst_level)
+            if src_level === nothing || dst_level === nothing
+                reduction_factor_outlet = 1.0
+            else
                 Δlevel = src_level - dst_level
                 reduction_factor_outlet = reduction_factor(Δlevel, 0.1)
-            else
-                reduction_factor_outlet = 1.0
             end
         else
             reduction_factor_outlet = 1.0
@@ -1005,14 +1005,14 @@ function formulate_flow!(
         src_level = get_level(p, src_id, current_level, t)
         dst_level = get_level(p, dst_id, current_level, t)
 
-        if !isnothing(src_level) && !isnothing(dst_level)
+        if src_level !== nothing && dst_level !== nothing
             Δlevel = src_level - dst_level
             reduction_factor_outlet = reduction_factor(Δlevel, 0.1)
             q *= reduction_factor_outlet
         end
 
         # No flow out outlet if source level is lower than minimum crest level
-        if !isnothing(src_level) && !isnan(min_crest_level[i])
+        if src_level !== nothing && !isnan(min_crest_level[i])
             reduction_factor_min_crest_level =
                 reduction_factor(src_level - min_crest_level[i], 0.1)
             q *= reduction_factor_min_crest_level
diff --git a/core/src/utils.jl b/core/src/utils.jl
index e84219e34..155de2b70 100644
--- a/core/src/utils.jl
+++ b/core/src/utils.jl
@@ -222,11 +222,11 @@ Ribasim.findlastgroup(2, [5,4,2,2,5,2,2,2,1])
 """
 function findlastgroup(id::Int, ids::AbstractVector{Int})::UnitRange{Int}
     idx_block_end = findlast(==(id), ids)
-    if isnothing(idx_block_end)
+    if idx_block_end === nothing
         return 1:0
     end
     idx_block_begin = findprev(!=(id), ids, idx_block_end)
-    idx_block_begin = if isnothing(idx_block_begin)
+    idx_block_begin = if idx_block_begin === nothing
         1
     else
         # can happen if that id is the only ID in ids
@@ -353,7 +353,7 @@ function set_static_value!(
 )::NamedTuple
     for (i, id) in enumerate(node_id)
         idx = findsorted(static.node_id, id)
-        isnothing(idx) && continue
+        idx === nothing && continue
         row = static[idx]
         set_table_row!(table, row, i)
     end
@@ -380,7 +380,7 @@ function set_current_value!(
                 row -> row.node_id == id && !isnan(getproperty(row, symbol)),
                 pre_table,
             )
-            if !isnothing(idx)
+            if idx !== nothing
                 vector[i] = getproperty(pre_table, symbol)[idx]
             end
         end
@@ -456,7 +456,7 @@ function basin_bottoms(
 )::Tuple{Float64, Float64}
     bottom_a = basin_bottom(basin, basin_a_id)
     bottom_b = basin_bottom(basin, basin_b_id)
-    if isnothing(bottom_a) && isnothing(bottom_b)
+    if bottom_a === bottom_b === nothing
         error(lazy"No bottom defined on either side of $id")
     end
     bottom_a = something(bottom_a, bottom_b)
diff --git a/core/test/libribasim.jl b/core/test/libribasim.jl
index 81ad77772..bb604b92d 100644
--- a/core/test/libribasim.jl
+++ b/core/test/libribasim.jl
@@ -19,9 +19,9 @@ toml_path = normpath(@__DIR__, "../../data/basic/basic.toml")
         toml_path_ptr = Base.unsafe_convert(Cstring, toml_path)
 
         # safe to finalize uninitialized model
-        @test isnothing(libribasim.model)
+        @test libribasim.model === nothing
         @test libribasim.finalize() == 0
-        @test isnothing(libribasim.model)
+        @test libribasim.model === nothing
 
         # cannot get time of uninitialized model
         @test libribasim.last_error_message == ""
@@ -41,6 +41,6 @@ toml_path = normpath(@__DIR__, "../../data/basic/basic.toml")
         @test unsafe_string(type_ptr) == "double"
 
         @test libribasim.finalize() == 0
-        @test isnothing(libribasim.model)
+        @test libribasim.model === nothing
     end
 end

From 328d8d9b77e834661698ddc79cf2bb03f86d6503 Mon Sep 17 00:00:00 2001
From: Martijn Visser <mgvisser@gmail.com>
Date: Wed, 13 Sep 2023 13:22:51 +0200
Subject: [PATCH 16/22] keeping it brief

---
 core/src/solve.jl | 48 +++++++++++++++++++----------------------------
 1 file changed, 19 insertions(+), 29 deletions(-)

diff --git a/core/src/solve.jl b/core/src/solve.jl
index baa514d31..634acdab0 100644
--- a/core/src/solve.jl
+++ b/core/src/solve.jl
@@ -523,12 +523,12 @@ function formulate!(
         bottom = basin.level[i][1]
         fixed_area = basin.area[i][end]
         depth = max(level - bottom, 0.0)
-        reduction_factor_basin = reduction_factor(depth, 0.1)
+        factor = reduction_factor(depth, 0.1)
 
         precipitation = fixed_area * basin.precipitation[i]
-        evaporation = area * reduction_factor_basin * basin.potential_evaporation[i]
+        evaporation = area * factor * basin.potential_evaporation[i]
         drainage = basin.drainage[i]
-        infiltration = reduction_factor_basin * basin.infiltration[i]
+        infiltration = factor * basin.infiltration[i]
 
         du.storage[i] += precipitation - evaporation + drainage - infiltration
     end
@@ -602,20 +602,20 @@ function continuous_control!(
             dst_level = get_level(p, dst_id, current_level, t)
 
             if src_level === nothing || dst_level === nothing
-                reduction_factor_outlet = 1.0
+                factor_outlet = 1.0
             else
                 Δlevel = src_level - dst_level
-                reduction_factor_outlet = reduction_factor(Δlevel, 0.1)
+                factor_outlet = reduction_factor(Δlevel, 0.1)
             end
         else
-            reduction_factor_outlet = 1.0
+            factor_outlet = 1.0
         end
 
         if controls_pump
             controlled_node_idx = findsorted(pump.node_id, controlled_node_id)
 
             listened_basin_storage = u.storage[listened_node_idx]
-            reduction_factor_basin = reduction_factor(listened_basin_storage, 10.0)
+            factor_basin = reduction_factor(listened_basin_storage, 10.0)
         else
             controlled_node_idx = findsorted(outlet.node_id, controlled_node_id)
 
@@ -624,13 +624,13 @@ function continuous_control!(
             has_index, upstream_basin_idx = id_index(basin.node_id, upstream_node_id)
             if has_index
                 upstream_basin_storage = u.storage[upstream_basin_idx]
-                reduction_factor_basin = reduction_factor(upstream_basin_storage, 10.0)
+                factor_basin = reduction_factor(upstream_basin_storage, 10.0)
             else
-                reduction_factor_basin = 1.0
+                factor_basin = 1.0
             end
         end
 
-        reduction_factor_all = reduction_factor_basin * reduction_factor_outlet
+        factor = factor_basin * factor_outlet
         flow_rate = 0.0
 
         K_p, K_i, K_d = pid_params[i](t)
@@ -638,17 +638,17 @@ function continuous_control!(
         if !iszero(K_d)
             # dlevel/dstorage = 1/area
             area = current_area[listened_node_idx]
-            D = 1.0 - K_d * reduction_factor_all / area
+            D = 1.0 - K_d * factor / area
         else
             D = 1.0
         end
 
         if !iszero(K_p)
-            flow_rate += reduction_factor_all * K_p * pid_error[i] / D
+            flow_rate += factor * K_p * pid_error[i] / D
         end
 
         if !iszero(K_i)
-            flow_rate += reduction_factor_all * K_i * integral_value[i] / D
+            flow_rate += factor * K_i * integral_value[i] / D
         end
 
         if !iszero(K_d)
@@ -773,11 +773,8 @@ function formulate_flow!(
         if active[i]
             hasindex, basin_idx = id_index(basin.node_id, upstream_basin_id)
             @assert hasindex "TabulatedRatingCurve must be downstream of a Basin"
-            s = storage[basin_idx]
-            reduction_factor = clamp(s, 0.0, 10.0) / 10.0
-            q =
-                reduction_factor *
-                tables[i](get_level(p, upstream_basin_id, current_level, t))
+            factor = reduction_factor(storage[basin_idx], 10.0)
+            q = factor * tables[i](get_level(p, upstream_basin_id, current_level, t))
         else
             q = 0.0
         end
@@ -957,9 +954,7 @@ function formulate_flow!(
 
         if hasindex
             # Pumping from basin
-            s = storage[basin_idx]
-            reduction_factor_basin = reduction_factor(s, 10.0)
-            q *= reduction_factor_basin
+            q *= reduction_factor(storage[basin_idx], 10.0)
         end
 
         flow[src_id, id] = q
@@ -996,9 +991,7 @@ function formulate_flow!(
 
         if hasindex
             # Flowing from basin
-            s = storage[basin_idx]
-            reduction_factor_basin = reduction_factor(s, 10.0)
-            q *= reduction_factor_basin
+            q *= reduction_factor(storage[basin_idx], 10.0)
         end
 
         # No flow of outlet if source level is lower than target level
@@ -1007,15 +1000,12 @@ function formulate_flow!(
 
         if src_level !== nothing && dst_level !== nothing
             Δlevel = src_level - dst_level
-            reduction_factor_outlet = reduction_factor(Δlevel, 0.1)
-            q *= reduction_factor_outlet
+            q *= reduction_factor(Δlevel, 0.1)
         end
 
         # No flow out outlet if source level is lower than minimum crest level
         if src_level !== nothing && !isnan(min_crest_level[i])
-            reduction_factor_min_crest_level =
-                reduction_factor(src_level - min_crest_level[i], 0.1)
-            q *= reduction_factor_min_crest_level
+            q *= reduction_factor(src_level - min_crest_level[i], 0.1)
         end
 
         flow[src_id, id] = q

From 54a57823bd77bea999f1e02c4d9a475c9cf9b010 Mon Sep 17 00:00:00 2001
From: Martijn Visser <mgvisser@gmail.com>
Date: Wed, 13 Sep 2023 13:32:06 +0200
Subject: [PATCH 17/22] make reduction_factor type stable, add tests

---
 core/src/utils.jl  |  6 +++---
 core/test/utils.jl | 10 ++++++++++
 2 files changed, 13 insertions(+), 3 deletions(-)

diff --git a/core/src/utils.jl b/core/src/utils.jl
index 155de2b70..65cb8372c 100644
--- a/core/src/utils.jl
+++ b/core/src/utils.jl
@@ -906,13 +906,13 @@ end
 Function that goes smoothly from 0 to 1 in the interval [0,threshold],
 and is constant outside this interval.
 """
-function reduction_factor(x::Real, threshold::Real)::Real
+function reduction_factor(x::T, threshold::Real)::T where {T <: Real}
     return if x < 0
-        0.0
+        zero(T)
     elseif x < threshold
         x_scaled = x / threshold
         (-2 * x_scaled + 3) * x_scaled^2
     else
-        1.0
+        one(T)
     end
 end
diff --git a/core/test/utils.jl b/core/test/utils.jl
index a0ae96bd7..0544cbe99 100644
--- a/core/test/utils.jl
+++ b/core/test/utils.jl
@@ -206,3 +206,13 @@ end
     @test isempty(fv)
     @test length(fv) == 0
 end
+
+@testset "reduction_factor" begin
+    @test Ribasim.reduction_factor(-2.0, 2.0) === 0.0
+    @test Ribasim.reduction_factor(0.0f0, 2.0) === 0.0f0
+    @test Ribasim.reduction_factor(0.0, 2.0) === 0.0
+    @test Ribasim.reduction_factor(1.0f0, 2.0) === 0.5f0
+    @test Ribasim.reduction_factor(1.0, 2.0) === 0.5
+    @test Ribasim.reduction_factor(3.0f0, 2.0) === 1.0f0
+    @test Ribasim.reduction_factor(3.0, 2.0) === 1.0
+end

From 079b31a6a2a15ddf0c7d89d4ebcbc6b285dd5907 Mon Sep 17 00:00:00 2001
From: Martijn Visser <mgvisser@gmail.com>
Date: Wed, 13 Sep 2023 13:47:17 +0200
Subject: [PATCH 18/22] update state for new reduction_factor

---
 core/test/run_models.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/core/test/run_models.jl b/core/test/run_models.jl
index 6847eaab1..ad800deef 100644
--- a/core/test/run_models.jl
+++ b/core/test/run_models.jl
@@ -86,7 +86,7 @@ end
     model = Ribasim.run(toml_path)
     @test model isa Ribasim.Model
     @test successful_retcode(model)
-    @test model.integrator.sol.u[end] ≈ Float32[8.41143, 725.5068] skip = Sys.isapple()
+    @test model.integrator.sol.u[end] ≈ Float32[7.783636, 726.16394] skip = Sys.isapple()
     # the highest level in the dynamic table is updated to 1.2 from the callback
     @test model.integrator.p.tabulated_rating_curve.tables[end].t[end] == 1.2
 end

From db28e8f26abc1812002f6691ea4cff17ba902d11 Mon Sep 17 00:00:00 2001
From: Martijn Visser <mgvisser@gmail.com>
Date: Wed, 13 Sep 2023 14:00:18 +0200
Subject: [PATCH 19/22] fix QGIS class LevelBoundaryTime(Input)

---
 qgis/core/nodes.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/qgis/core/nodes.py b/qgis/core/nodes.py
index 5f666e405..94d5624ea 100644
--- a/qgis/core/nodes.py
+++ b/qgis/core/nodes.py
@@ -397,9 +397,10 @@ class LevelBoundaryStatic(Input):
 
 
 class LevelBoundaryTime(Input):
-    input_type = "LevelBoundary / static"
+    input_type = "LevelBoundary / time"
     geometry_type = "No Geometry"
     attributes = [
+        QgsField("time", QVariant.DateTime),
         QgsField("node_id", QVariant.Int),
         QgsField("time", QVariant.DateTime),
         QgsField("level", QVariant.Double),

From 5b576355e8f3ce7348673540af5102d724a43456 Mon Sep 17 00:00:00 2001
From: Martijn Visser <mgvisser@gmail.com>
Date: Wed, 13 Sep 2023 14:19:21 +0200
Subject: [PATCH 20/22] make Outlet green and add min_crest_level to QGIS

I prefer to avoid yellow since it is the color of selected nodes in QGIS.
---
 python/ribasim/ribasim/geometry/node.py | 2 +-
 qgis/core/nodes.py                      | 3 ++-
 2 files changed, 3 insertions(+), 2 deletions(-)

diff --git a/python/ribasim/ribasim/geometry/node.py b/python/ribasim/ribasim/geometry/node.py
index 998fdaeda..a61f9aa4d 100644
--- a/python/ribasim/ribasim/geometry/node.py
+++ b/python/ribasim/ribasim/geometry/node.py
@@ -147,7 +147,7 @@ def plot(self, ax=None, zorder=None) -> Any:
             "ManningResistance": "r",
             "TabulatedRatingCurve": "g",
             "Pump": "0.5",  # grayscale level
-            "Outlet": "y",
+            "Outlet": "g",
             "Terminal": "m",
             "FlowBoundary": "m",
             "DiscreteControl": "k",
diff --git a/qgis/core/nodes.py b/qgis/core/nodes.py
index 94d5624ea..a12cb068a 100644
--- a/qgis/core/nodes.py
+++ b/qgis/core/nodes.py
@@ -177,7 +177,7 @@ def renderer(self) -> QgsCategorizedSymbolRenderer:
             "LevelBoundary": (QColor("green"), "LevelBoundary", shape.Circle),
             "FlowBoundary": (QColor("purple"), "FlowBoundary", shape.Hexagon),
             "Pump": (QColor("gray"), "Pump", shape.Hexagon),
-            "Outlet": (QColor("yellow"), "Outlet", shape.Hexagon),
+            "Outlet": (QColor("green"), "Outlet", shape.Hexagon),
             "ManningResistance": (QColor("red"), "ManningResistance", shape.Diamond),
             "Terminal": (QColor("purple"), "Terminal", shape.Square),
             "DiscreteControl": (QColor("black"), "DiscreteControl", shape.Star),
@@ -429,6 +429,7 @@ class OutletStatic(Input):
         QgsField("flow_rate", QVariant.Double),
         QgsField("min_flow_rate", QVariant.Double),
         QgsField("max_flow_rate", QVariant.Double),
+        QgsField("min_crest_level", QVariant.Double),
         QgsField("control_state", QVariant.String),
     ]
 

From 953e915bfa34b22d63851b9c4eeaa00fb54f73f9 Mon Sep 17 00:00:00 2001
From: Martijn Visser <mgvisser@gmail.com>
Date: Wed, 13 Sep 2023 15:32:16 +0200
Subject: [PATCH 21/22] forgot to commit, so the last layer was not in
 gpkg_contents

Making it not show up in QGIS.
---
 python/ribasim/ribasim/model.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/python/ribasim/ribasim/model.py b/python/ribasim/ribasim/model.py
index 66c2e3bdc..4b271865f 100644
--- a/python/ribasim/ribasim/model.py
+++ b/python/ribasim/ribasim/model.py
@@ -200,6 +200,7 @@ def _write_tables(self, directory: FilePath) -> None:
                 )
                 if isinstance(input_entry, TableModel) and not is_geometry:
                     input_entry.write_table(connection)
+            connection.commit()
 
         shutil.move(temp_path, gpkg_path)
         return

From a642f9f12dedaf932a24cf19870b3e4685b9bd5f Mon Sep 17 00:00:00 2001
From: Martijn Visser <mgvisser@gmail.com>
Date: Wed, 13 Sep 2023 16:36:49 +0200
Subject: [PATCH 22/22] add note on Outlet to example model

---
 docs/python/examples.ipynb | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/docs/python/examples.ipynb b/docs/python/examples.ipynb
index addae4e7a..acb00ca7b 100644
--- a/docs/python/examples.ipynb
+++ b/docs/python/examples.ipynb
@@ -635,7 +635,7 @@
    "source": [
     "# Model with discrete control\n",
     "\n",
-    "The model constructed below consists of a single basin which slowly drains trough a `TabulatedRatingCurve`, but is held within a range around a target level (setpoint) by two connected pumps."
+    "The model constructed below consists of a single basin which slowly drains trough a `TabulatedRatingCurve`, but is held within a range around a target level (setpoint) by two connected pumps. These two pumps behave like a reversible pump. When pumping can be done in only one direction, and the other direction is only possible under gravity, use an Outlet for that direction."
    ]
   },
   {
@@ -660,7 +660,7 @@
     "        (2.0, 0.0),  # 4: LevelBoundary\n",
     "        (-1.0, 0.0),  # 5: TabulatedRatingCurve\n",
     "        (-2.0, 0.0),  # 6: Terminal\n",
-    "        (0.0, 1.5),  # 7: DiscreteControl\n",
+    "        (1.0, 0.0),  # 7: DiscreteControl\n",
     "    ]\n",
     ")\n",
     "\n",
@@ -1449,7 +1449,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.4"
+   "version": "3.10.12"
   }
  },
  "nbformat": 4,