Update reduction factors of resistance nodes

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