Merge branch 'main' into quick

.teamcity/Ribasim_Windows/RibasimWindowsProject.kt

@@ -89,7 +89,7 @@ object Windows_TestDelwaqCoupling : BuildType({
     templates(WindowsAgent, GithubCommitStatusIntegration, TestDelwaqCouplingWindows)
     name = "Test Delwaq coupling"
 
-    artifactRules = "ribasim/coupling/delwaq/model"
+    artifactRules = "ribasim/python/ribasim/ribasim/delwaq/model"
 
     triggers {
         vcs {

.teamcity/Templates/RegressionTest.kt

@@ -59,7 +59,7 @@ open class RegressionTest (platformOs: String) : Template() {
                 scriptContent = header +
                         """
                 pixi run python utils/get_benchmark.py %MiniO_credential_token% "benchmark/"
-                pixi run python utils/get_benchmark.py %MiniO_credential_token% "hws_2024_7_0/"
+                pixi run python utils/get_benchmark.py %MiniO_credential_token% "hws_migration_test/"
                 pixi run test-ribasim-regression
                 """.trimIndent()
             }
@@ -73,4 +73,4 @@ open class RegressionTest (platformOs: String) : Template() {
 }
 
 object RegressionTestWindows : RegressionTest("Windows")
-object RegressionTestLinux : RegressionTest("Linux")
+object RegressionTestLinux : RegressionTest("Linux")

core/integration_test/hws_integration_test.jl

@@ -17,7 +17,7 @@
 
     @testset "Results values" begin
         @test basin.node_id == basin_bench.node_id
-        @test all(q -> abs(q) < 0.01, basin.level - basin_bench.level)
+        @test all(q -> abs(q) < 0.06, basin.level - basin_bench.level)
     end
 
     timed = @timed Ribasim.run(toml_path)
@@ -27,6 +27,7 @@
     performance_diff =
         round((timed.time - benchmark_runtime) / benchmark_runtime * 100; digits = 2)
     if performance_diff < 0.0
+        performance_diff = abs(performance_diff)
         @info "Runtime is $(timed.time) and it is $performance_diff % faster than benchmark"
     elseif performance_diff > 0.0 && performance_diff < 0.2
         @info "Runtime is $(timed.time) and it is $performance_diff % slower than benchmark"

core/src/model.jl

@@ -62,17 +62,21 @@ function Model(config::Config)::Model
             error("Invalid discrete control state definition(s).")
         end
 
-        (; pid_control, graph, outlet, basin, tabulated_rating_curve) = parameters
+        (; pid_control, graph, outlet, pump, basin, tabulated_rating_curve) = parameters
         if !valid_pid_connectivity(pid_control.node_id, pid_control.listen_node_id, graph)
             error("Invalid PidControl connectivity.")
         end
 
-        if !valid_outlet_crest_level!(graph, outlet, basin)
-            error("Invalid minimum crest level of outlet")
+        if !valid_min_upstream_level!(graph, outlet, basin)
+            error("Invalid minimum upstream level of Outlet.")
+        end
+
+        if !valid_min_upstream_level!(graph, pump, basin)
+            error("Invalid minimum upstream level of Pump.")
         end
 
         if !valid_tabulated_curve_level(graph, tabulated_rating_curve, basin)
-            error("Invalid level of tabulated rating curve")
+            error("Invalid level of TabulatedRatingCurve.")
         end
 
         # tell the solver to stop when new data comes in

core/src/read.jl

@@ -298,6 +298,7 @@ function TabulatedRatingCurve(
     interpolations = ScalarInterpolation[]
     control_mapping = Dict{Tuple{NodeID, String}, ControlStateUpdate}()
     active = Bool[]
+    max_downstream_level = Float64[]
     errors = false
 
     for node_id in node_ids
@@ -317,12 +318,15 @@ function TabulatedRatingCurve(
                 IterTools.groupby(row -> coalesce(row.control_state, nothing), static_id)
                 control_state = first(group).control_state
                 is_active = coalesce(first(group).active, true)
+                max_level = coalesce(first(group).max_downstream_level, Inf)
                 table = StructVector(group)
-                if !valid_tabulated_rating_curve(node_id, table)
+                rowrange =
+                    findlastgroup(node_id, NodeID.(node_id.type, table.node_id, Ref(0)))
+                if !valid_tabulated_rating_curve(node_id, table, rowrange)
                     errors = true
                 end
                 interpolation = try
-                    qh_interpolation(node_id, table)
+                    qh_interpolation(table, rowrange)
                 catch
                     LinearInterpolation(Float64[], Float64[])
                 end
@@ -339,17 +343,23 @@ function TabulatedRatingCurve(
             end
             push!(interpolations, interpolation)
             push!(active, is_active)
+            push!(max_downstream_level, max_level)
         elseif node_id in time_node_ids
             source = "time"
             # get the timestamp that applies to the model starttime
             idx_starttime = searchsortedlast(time.time, config.starttime)
             pre_table = view(time, 1:idx_starttime)
-            if !valid_tabulated_rating_curve(node_id, pre_table)
+            rowrange =
+                findlastgroup(node_id, NodeID.(node_id.type, pre_table.node_id, Ref(0)))
+
+            if !valid_tabulated_rating_curve(node_id, pre_table, rowrange)
                 errors = true
             end
-            interpolation = qh_interpolation(node_id, pre_table)
+            interpolation = qh_interpolation(pre_table, rowrange)
+            max_level = coalesce(pre_table.max_downstream_level[rowrange][begin], Inf)
             push!(interpolations, interpolation)
             push!(active, true)
+            push!(max_downstream_level, max_level)
         else
             @error "$node_id data not in any table."
             errors = true
@@ -364,6 +374,7 @@ function TabulatedRatingCurve(
         inflow_edge = inflow_edge.(Ref(graph), node_ids),
         outflow_edge = outflow_edge.(Ref(graph), node_ids),
         active,
+        max_downstream_level,
         table = interpolations,
         time,
         control_mapping,

core/src/schema.jl

@@ -210,6 +210,7 @@ end
     active::Union{Missing, Bool}
     level::Float64
     flow_rate::Float64
+    max_downstream_level::Union{Missing, Float64}
     control_state::Union{Missing, String}
 end
 
@@ -218,6 +219,7 @@ end
     time::DateTime
     level::Float64
     flow_rate::Float64
+    max_downstream_level::Union{Missing, Float64}
 end
 
 @version DiscreteControlVariableV1 begin

core/src/solve.jl

@@ -331,22 +331,30 @@ function formulate_flow!(
 )::Nothing
     (; graph) = p
     all_nodes_active = p.all_nodes_active[]
-    (; node_id, active, table, inflow_edge, outflow_edge) = tabulated_rating_curve
+    (; node_id, active, table) = tabulated_rating_curve
 
     for id in node_id
-        upstream_edge = inflow_edge[id.idx]
-        downstream_edge = outflow_edge[id.idx]
-        upstream_basin_id = upstream_edge.edge[1]
+        inflow_edge = tabulated_rating_curve.inflow_edge[id.idx]
+        outflow_edge = tabulated_rating_curve.outflow_edge[id.idx]
+        inflow_id = inflow_edge.edge[1]
+        outflow_id = outflow_edge.edge[2]
+        max_downstream_level = tabulated_rating_curve.max_downstream_level[id.idx]
+
+        h_a = get_level(p, inflow_id, t, du)
+        h_b = get_level(p, outflow_id, t, du)
+        Δh = h_a - h_b
 
         if active[id.idx] || all_nodes_active
-            factor = low_storage_factor(u.storage, upstream_basin_id, 10.0)
-            q = factor * table[id.idx](get_level(p, upstream_basin_id, t, du))
+            factor = low_storage_factor(u.storage, inflow_id, 10.0)
+            q = factor * table[id.idx](h_a)
+            q *= reduction_factor(Δh, 0.02)
+            q *= reduction_factor(max_downstream_level - h_b, 0.02)
         else
             q = 0.0
         end
 
-        set_flow!(graph, upstream_edge, q, du)
-        set_flow!(graph, downstream_edge, q, du)
+        set_flow!(graph, inflow_edge, q, du)
+        set_flow!(graph, outflow_edge, q, du)
     end
     return nothing
 end

core/src/util.jl

@@ -123,8 +123,10 @@ end
 From a table with columns node_id, flow_rate (Q) and level (h),
 create a ScalarInterpolation from level to flow rate for a given node_id.
 """
-function qh_interpolation(node_id::NodeID, table::StructVector)::ScalarInterpolation
-    rowrange = findlastgroup(node_id, NodeID.(node_id.type, table.node_id, Ref(0)))
+function qh_interpolation(
+    table::StructVector,
+    rowrange::UnitRange{Int},
+)::ScalarInterpolation
     level = table.level[rowrange]
     flow_rate = table.flow_rate[rowrange]
 

core/src/validation.jl

@@ -357,18 +357,22 @@ function valid_demand(
 end
 
 """
-Validate Outlet crest level and fill in default values
+Validate Outlet or Pump `min_upstream_level` and fill in default values
 """
-function valid_outlet_crest_level!(graph::MetaGraph, outlet::Outlet, basin::Basin)::Bool
+function valid_min_upstream_level!(
+    graph::MetaGraph,
+    node::Union{Outlet, Pump},
+    basin::Basin,
+)::Bool
     errors = false
-    for (id, crest) in zip(outlet.node_id, outlet.min_upstream_level)
+    for (id, min_upstream_level) in zip(node.node_id, node.min_upstream_level)
         id_in = inflow_id(graph, id)
         if id_in.type == NodeType.Basin
             basin_bottom_level = basin_bottom(basin, id_in)[2]
-            if crest == -Inf
-                outlet.min_upstream_level[id.idx] = basin_bottom_level
-            elseif crest < basin_bottom_level
-                @error "Minimum crest level of $id is lower than bottom of upstream $id_in" crest basin_bottom_level
+            if min_upstream_level == -Inf
+                node.min_upstream_level[id.idx] = basin_bottom_level
+            elseif min_upstream_level < basin_bottom_level
+                @error "Minimum upstream level of $id is lower than bottom of upstream $id_in" min_upstream_level basin_bottom_level
                 errors = true
             end
         end
@@ -397,10 +401,13 @@ function valid_tabulated_curve_level(
     return !errors
 end
 
-function valid_tabulated_rating_curve(node_id::NodeID, table::StructVector)::Bool
+function valid_tabulated_rating_curve(
+    node_id::NodeID,
+    table::StructVector,
+    rowrange::UnitRange{Int},
+)::Bool
     errors = false
 
-    rowrange = findlastgroup(node_id, NodeID.(node_id.type, table.node_id, Ref(0)))
     level = table.level[rowrange]
     flow_rate = table.flow_rate[rowrange]
 

core/test/allocation_test.jl

@@ -465,7 +465,7 @@ end
     model = Ribasim.run(toml_path)
     record_demand = DataFrame(model.integrator.p.allocation.record_demand)
     df_rating_curve_2 = record_demand[record_demand.node_id .== 2, :]
-    @test all(df_rating_curve_2.realized .≈ 2e-3)
+    @test all(df_rating_curve_2.realized .≈ 0.002)
 
     @testset "Results" begin
         allocation_bytes = read(normpath(dirname(toml_path), "results/allocation.arrow"))

core/test/run_models_test.jl

@@ -198,7 +198,7 @@ end
 
     @test successful_retcode(model)
     @test allunique(Ribasim.tsaves(model))
-    @test model.integrator.u ≈ Float32[519.8817, 519.8798, 339.3959, 1418.4331] skip =
+    @test model.integrator.u ≈ Float32[803.7093, 803.68274, 495.241, 1318.3053] skip =
         Sys.isapple() atol = 1.5
 
     @test length(logger.logs) > 10
@@ -236,8 +236,8 @@ end
             model.integrator.p.basin.vertical_flux[Float64[]],
         ).precipitation
     @test length(precipitation) == 4
-    @test model.integrator.u ≈ Float32[472.06555, 472.06366, 367.23883, 1427.9957] atol =
-        2.0 skip = Sys.isapple()
+    @test model.integrator.u ≈ Float32[698.22736, 698.2014, 421.20447, 1334.4354] atol = 2.0 skip =
+        Sys.isapple()
 end
 
 @testitem "Allocation example model" begin
@@ -290,7 +290,7 @@ end
     model = Ribasim.run(toml_path)
     @test model isa Ribasim.Model
     @test successful_retcode(model)
-    @test model.integrator.u ≈ Float32[7.783636, 726.16394] skip = Sys.isapple()
+    @test model.integrator.u ≈ Float32[368.31558, 365.68442] 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.table[end].t[end] == 1.2
 end
@@ -375,7 +375,7 @@ end
     t_min_upstream_level =
         level.t[2] * (outlet.min_upstream_level[1] - level.u[1]) / (level.u[2] - level.u[1])
 
-    # No outlet flow when upstream level is below minimum crest level
+    # No outlet flow when upstream level is below minimum upstream level
     @test all(@. outlet_flow.flow_rate[t <= t_min_upstream_level] == 0)
 
     t = Ribasim.tsaves(model)

core/test/utils_test.jl

@@ -180,7 +180,7 @@ end
 
     jac_prototype_expected = spzeros(Bool, 4, 4)
     jac_prototype_expected[1:2, 1:2] .= true
-    jac_prototype_expected[3:4, 2:4] .= true
+    jac_prototype_expected[2:4, 2:4] .= true
     @test jac_prototype == jac_prototype_expected
 
     toml_path = normpath(@__DIR__, "../../generated_testmodels/pid_control/ribasim.toml")

core/test/validation_test.jl

@@ -23,7 +23,7 @@
     @test logger.logs[3].message == "Basin #1 profile cannot have decreasing areas."
 
     table = StructVector(; flow_rate = [0.0, 0.1], level = [1.0, 2.0], node_id = [5, 5])
-    itp = qh_interpolation(NodeID(:TabulatedRatingCurve, 5, 1), table)
+    itp = qh_interpolation(table, 1:2)
     # constant extrapolation at the bottom end, linear extrapolation at the top end
     itp(0.0) ≈ 0.0
     itp(1.0) ≈ 0.0
@@ -402,7 +402,7 @@ end
 end
 
 @testitem "Outlet upstream level validation" begin
-    using Ribasim: valid_outlet_crest_level!
+    using Ribasim: valid_min_upstream_level!
     using Logging
 
     toml_path = normpath(
@@ -422,13 +422,13 @@ end
 
     logger = TestLogger()
     with_logger(logger) do
-        @test !Ribasim.valid_outlet_crest_level!(graph, outlet, basin)
+        @test !Ribasim.valid_min_upstream_level!(graph, outlet, basin)
     end
 
     @test length(logger.logs) == 1
     @test logger.logs[1].level == Error
     @test logger.logs[1].message ==
-          "Minimum crest level of Outlet #4 is lower than bottom of upstream Basin #3"
+          "Minimum upstream level of Outlet #4 is lower than bottom of upstream Basin #3"
 end
 
 @testitem "Convergence bottleneck" begin