Update reduction factors of resistance nodes (#1796)

These changes are needed for running the `De Dommel` model for Ribasim
NL and were requested by Neeltje. The main change is that
`ManningResistance` has low storage factors now, which is important to
have if the connected basins have different bottom levels.

I made an unified approach for low storage factors of `LinearResistance`
and `ManningResistance`, where the former already had low storage
factors for both connected basins. That however leads to weird behavior
where water cannot flow into an empty basin via the resistance node, so
now I made it so that the reduction factor is based on the storage of
the upstream basin as defined by the flow direction. This feels a bit
discontinuous but I think it's fine.

I also fixed an upstream issue in `FindFirstFunctions`, for now the
manifest points at the branch with the fix:
SciML/FindFirstFunctions.jl#26.

Manifest.toml

@@ -2,7 +2,7 @@
 
 julia_version = "1.10.5"
 manifest_format = "2.0"
-project_hash = "77c0bd9bc2d2b8efcc64c81eafef58824c0d85ae"
+project_hash = "d694251e41cfe037d14e0d129b2e2e49035c0990"
 
 [[deps.ADTypes]]
 git-tree-sha1 = "99a6f5d0ce1c7c6afdb759daa30226f71c54f6b0"
@@ -557,7 +557,9 @@ version = "1.12.0"
     Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [[deps.FindFirstFunctions]]
-git-tree-sha1 = "fa0ba2042021409deb144f868abafde0b06be8f0"
+git-tree-sha1 = "9d43caa23d1b97b57e3d70defac42cbe1f08b7d2"
+repo-rev = "catch_inexact_rounding_error"
+repo-url = "https://github.com/SouthEndMusic/FindFirstFunctions.jl"
 uuid = "64ca27bc-2ba2-4a57-88aa-44e436879224"
 version = "1.3.0"
 
@@ -1342,7 +1344,7 @@ uuid = "295af30f-e4ad-537b-8983-00126c2a3abe"
 version = "3.5.18"
 
 [[deps.Ribasim]]
-deps = ["Accessors", "Arrow", "BasicModelInterface", "CodecZstd", "ComponentArrays", "Configurations", "DBInterface", "DataInterpolations", "DataStructures", "Dates", "DiffEqCallbacks", "EnumX", "FiniteDiff", "ForwardDiff", "Graphs", "HiGHS", "IterTools", "JuMP", "Legolas", "LineSearches", "LinearSolve", "Logging", "LoggingExtras", "MetaGraphsNext", "OrdinaryDiffEqBDF", "OrdinaryDiffEqCore", "OrdinaryDiffEqLowOrderRK", "OrdinaryDiffEqNonlinearSolve", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqSDIRK", "OrdinaryDiffEqTsit5", "PreallocationTools", "SQLite", "SciMLBase", "SparseArrays", "SparseConnectivityTracer", "StructArrays", "Tables", "TerminalLoggers", "TranscodingStreams"]
+deps = ["Accessors", "Arrow", "BasicModelInterface", "CodecZstd", "ComponentArrays", "Configurations", "DBInterface", "DataInterpolations", "DataStructures", "Dates", "DiffEqCallbacks", "EnumX", "FindFirstFunctions", "FiniteDiff", "ForwardDiff", "Graphs", "HiGHS", "IterTools", "JuMP", "Legolas", "LineSearches", "LinearSolve", "Logging", "LoggingExtras", "MetaGraphsNext", "OrdinaryDiffEqBDF", "OrdinaryDiffEqCore", "OrdinaryDiffEqLowOrderRK", "OrdinaryDiffEqNonlinearSolve", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqSDIRK", "OrdinaryDiffEqTsit5", "PreallocationTools", "SQLite", "SciMLBase", "SparseArrays", "SparseConnectivityTracer", "StructArrays", "Tables", "TerminalLoggers", "TranscodingStreams"]
 path = "core"
 uuid = "aac5e3d9-0b8f-4d4f-8241-b1a7a9632635"
 version = "2024.10.0"

Project.toml

@@ -17,6 +17,7 @@ DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 EnumX = "4e289a0a-7415-4d19-859d-a7e5c4648b56"
+FindFirstFunctions = "64ca27bc-2ba2-4a57-88aa-44e436879224"
 FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"

core/Project.toml

@@ -17,6 +17,7 @@ DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 EnumX = "4e289a0a-7415-4d19-859d-a7e5c4648b56"
+FindFirstFunctions = "64ca27bc-2ba2-4a57-88aa-44e436879224"
 FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
@@ -69,6 +70,7 @@ DataStructures = "0.18"
 Dates = "<0.0.1, 1"
 DiffEqCallbacks = "3.6"
 EnumX = "1.0"
+FindFirstFunctions = "1.3"
 FiniteDiff = "2.21"
 ForwardDiff = "0.10"
 Graphs = "1.9"

core/src/solve.jl

@@ -84,9 +84,7 @@ function set_current_basin_properties!(
     current_level = current_level[parent(du)]
     current_area = current_area[parent(du)]
 
-    storage = u.storage
-
-    for (i, s) in enumerate(storage)
+    for (i, s) in enumerate(u.storage)
         current_level[i] = get_level_from_storage(basin, i, s)
         current_area[i] = basin.level_to_area[i](current_level[i])
     end
@@ -312,9 +310,7 @@ function formulate_flow!(
             h_b = get_level(p, outflow_id, t, du)
             q_unlimited = (h_a - h_b) / resistance[id.idx]
             q = clamp(q_unlimited, -max_flow_rate[id.idx], max_flow_rate[id.idx])
-
-            q *= low_storage_factor(u.storage, inflow_id, 10.0)
-            q *= low_storage_factor(u.storage, outflow_id, 10.0)
+            q *= low_storage_factor_resistance_node(u, q, inflow_id, outflow_id, 10.0)
 
             set_flow!(graph, inflow_edge, q, du)
             set_flow!(graph, outflow_edge, q, du)
@@ -453,6 +449,7 @@ function formulate_flow!(
         Δh = h_a - h_b
 
         q = A / n * ∛(R_h^2) * relaxed_root(Δh / L, 1e-3)
+        q *= low_storage_factor_resistance_node(u, q, inflow_id, outflow_id, 10.0)
 
         set_flow!(graph, inflow_edge, q, du)
         set_flow!(graph, outflow_edge, q, du)

core/src/util.jl

@@ -384,6 +384,18 @@ function low_storage_factor(
     end
 end
 
+"""
+For resistance nodes, give a reduction factor based on the upstream node
+as defined by the flow direction.
+"""
+function low_storage_factor_resistance_node(u, q, inflow_id, outflow_id, threshold)
+    if q > 0
+        low_storage_factor(u.storage, inflow_id, threshold)
+    else
+        low_storage_factor(u.storage, outflow_id, threshold)
+    end
+end
+
 """Whether the given node node is flow constraining by having a maximum flow rate."""
 function is_flow_constraining(type::NodeType.T)::Bool
     type in (NodeType.LinearResistance, NodeType.Pump, NodeType.Outlet)
@@ -892,6 +904,13 @@ end
 
 # Custom overloads
 reduction_factor(x::GradientTracer, threshold::Real) = x
+low_storage_factor_resistance_node(
+    storage::ComponentVector{<:GradientTracer},
+    q,
+    inflow_id,
+    outflow_id,
+    threshold,
+) = q
 relaxed_root(x::GradientTracer, threshold::Real) = x
 get_level_from_storage(basin::Basin, state_idx::Int, storage::GradientTracer) = storage
 stop_declining_negative_storage!(du, u::ComponentVector{<:GradientTracer}) = nothing

core/test/utils_test.jl

@@ -233,7 +233,7 @@ end
 end
 
 @testitem "low_storage_factor" begin
-    using Ribasim: NodeID, low_storage_factor, EdgeMetadata, EdgeType
+    using Ribasim: NodeID, low_storage_factor, low_storage_factor_resistance_node
 
     node_id = NodeID(:Basin, 5, 1)
     @test low_storage_factor([-2.0], node_id, 2.0) === 0.0
@@ -243,6 +243,30 @@ end
     @test low_storage_factor([1.0], node_id, 2.0) === 0.5
     @test low_storage_factor([3.0f0], node_id, 2.0) === 1.0f0
     @test low_storage_factor([3.0], node_id, 2.0) === 1.0
+
+    node_id_1 = NodeID(:Basin, 5, 1)
+    node_id_2 = NodeID(:Basin, 6, 2)
+    @test low_storage_factor_resistance_node(
+        (; storage = [3.0, 3.0]),
+        1.0,
+        node_id_1,
+        node_id_2,
+        2.0,
+    ) == 1.0
+    @test low_storage_factor_resistance_node(
+        (; storage = [1.0, 3.0]),
+        1.0,
+        node_id_1,
+        node_id_2,
+        2.0,
+    ) == 0.5
+    @test low_storage_factor_resistance_node(
+        (; storage = [1.0, 3.0]),
+        -1.0,
+        node_id_1,
+        node_id_2,
+        2.0,
+    ) == 1.0
 end
 
 @testitem "constraints_from_nodes" begin