From f0a874700d1e1846b4f21f52b749430454e3b2e8 Mon Sep 17 00:00:00 2001 From: Bart de Koning <74617371+SouthEndMusic@users.noreply.github.com> Date: Wed, 4 Sep 2024 14:29:27 +0200 Subject: [PATCH] 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: https://github.com/SciML/FindFirstFunctions.jl/pull/26. --- Manifest.toml | 8 +++++--- Project.toml | 1 + core/Project.toml | 2 ++ core/src/solve.jl | 9 +++------ core/src/util.jl | 19 +++++++++++++++++++ core/test/utils_test.jl | 26 +++++++++++++++++++++++++- 6 files changed, 55 insertions(+), 10 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index b5d4db371..aedaba779 100644 --- a/Manifest.toml +++ b/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" diff --git a/Project.toml b/Project.toml index d75250e4f..f8fd65409 100644 --- a/Project.toml +++ b/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" diff --git a/core/Project.toml b/core/Project.toml index 521d86e91..b5ac7df39 100644 --- a/core/Project.toml +++ b/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" diff --git a/core/src/solve.jl b/core/src/solve.jl index 9738b834f..25dbe56c4 100644 --- a/core/src/solve.jl +++ b/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) diff --git a/core/src/util.jl b/core/src/util.jl index 8eee9fcd0..4e4c70666 100644 --- a/core/src/util.jl +++ b/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 diff --git a/core/test/utils_test.jl b/core/test/utils_test.jl index 60172477f..d05840f52 100644 --- a/core/test/utils_test.jl +++ b/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