From 4a468de55141f6600bf409987d0ec25a9db64819 Mon Sep 17 00:00:00 2001
From: Martijn Visser <mgvisser@gmail.com>
Date: Fri, 9 Feb 2024 16:55:48 +0100
Subject: [PATCH] Update test values

---
 core/src/parameter.jl                         |  1 +
 core/test/equations_test.jl                   | 28 +++++++++++--------
 .../ribasim_testmodels/equations.py           |  8 +++---
 3 files changed, 21 insertions(+), 16 deletions(-)

diff --git a/core/src/parameter.jl b/core/src/parameter.jl
index b4f0fcb1b..da492b245 100644
--- a/core/src/parameter.jl
+++ b/core/src/parameter.jl
@@ -205,6 +205,7 @@ 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
 
diff --git a/core/test/equations_test.jl b/core/test/equations_test.jl
index 59a17fe65..29db57890 100644
--- a/core/test/equations_test.jl
+++ b/core/test/equations_test.jl
@@ -18,24 +18,28 @@
 @testitem "LinearResistance" begin
     using SciMLBase: successful_retcode
 
-    toml_path =
-        normpath(@__DIR__, "../../generated_testmodels/linear_resistance/ribasim.toml")
+    toml_path = normpath(@__DIR__, "generated_testmodels/linear_resistance/ribasim.toml")
     @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])
-    # TODO adapt for max_flow_rate of 6e-5
-    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,
diff --git a/python/ribasim_testmodels/ribasim_testmodels/equations.py b/python/ribasim_testmodels/ribasim_testmodels/equations.py
index 734fb1082..a801c7684 100644
--- a/python/ribasim_testmodels/ribasim_testmodels/equations.py
+++ b/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,7 +68,7 @@ def linear_resistance_model():
     state = pd.DataFrame(
         data={
             "node_id": [1],
-            "level": [10.5],
+            "level": [10.0],
         }
     )