From 6cf30c8f5d002b9a715fabc10474e698cbdea50a Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Wed, 4 Sep 2024 07:39:38 +0200 Subject: [PATCH 1/4] Add reduction factors to ManningResistance --- core/src/solve.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/core/src/solve.jl b/core/src/solve.jl index 9738b834f..8a4dc18ef 100644 --- a/core/src/solve.jl +++ b/core/src/solve.jl @@ -454,6 +454,9 @@ function formulate_flow!( q = A / n * ∛(R_h^2) * relaxed_root(Δh / L, 1e-3) + q *= low_storage_factor(u.storage, inflow_id, 10.0) + q *= low_storage_factor(u.storage, outflow_id, 10.0) + set_flow!(graph, inflow_edge, q, du) set_flow!(graph, outflow_edge, q, du) end From ca8b7c8c0fce5fbcefde81cf2f1d24237fbd81c9 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Wed, 4 Sep 2024 08:43:03 +0200 Subject: [PATCH 2/4] Update reduction factors of resistance nodes + temporary fix of FindFirstFunctions --- Manifest.toml | 4 ++-- core/Project.toml | 1 + core/src/Ribasim.jl | 3 +++ core/src/solve.jl | 12 +++--------- core/src/util.jl | 43 +++++++++++++++++++++++++++++++++++++++++++ 5 files changed, 52 insertions(+), 11 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index b5d4db371..a76bfac8e 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 = "ce9b0fd5a3e37477853fac84d64bb99101f27d9c" [[deps.ADTypes]] git-tree-sha1 = "99a6f5d0ce1c7c6afdb759daa30226f71c54f6b0" @@ -1342,7 +1342,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/core/Project.toml b/core/Project.toml index 521d86e91..1268682c6 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" diff --git a/core/src/Ribasim.jl b/core/src/Ribasim.jl index 06c1bf7c3..ee2a202e6 100644 --- a/core/src/Ribasim.jl +++ b/core/src/Ribasim.jl @@ -24,6 +24,9 @@ using OrdinaryDiffEqCore: using OrdinaryDiffEqNonlinearSolve: OrdinaryDiffEqNonlinearSolve, relax!, _compute_rhs! using LineSearches: BackTracking +# TODO: Remove after https://github.com/SciML/FindFirstFunctions.jl/pull/26 +using FindFirstFunctions: FindFirstFunctions + # Interface for defining and solving the ODE problem of the physical layer. using SciMLBase: init, diff --git a/core/src/solve.jl b/core/src/solve.jl index 8a4dc18ef..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,9 +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(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) diff --git a/core/src/util.jl b/core/src/util.jl index 8eee9fcd0..1ce5dc8a8 100644 --- a/core/src/util.jl +++ b/core/src/util.jl @@ -65,6 +65,30 @@ function get_level_from_storage(basin::Basin, state_idx::Int, storage) end end +# TODO: Remove after https://github.com/SciML/FindFirstFunctions.jl/pull/26 +function (g::FindFirstFunctions.Guesser)(x) + (; v, idx_prev, linear_lookup) = g + if linear_lookup + f = (x - first(v)) / (last(v) - first(v)) + if isinf(f) + f > 0 ? lastindex(v) : firstindex(v) + else + i_0, i_f = firstindex(v), lastindex(v) + i_approx = f * (i_f - i_0) + i_0 + target_type = typeof(firstindex(v)) + if i_approx >= typemax(target_type) + lastindex(v) + 1 + elseif i_approx <= typemin(target_type) + firstindex(v) - 1 + else + round(target_type, i_approx) + end + end + else + idx_prev[] + end +end + """ For an element `id` and a vector of elements `ids`, get the range of indices of the last consecutive block of `id`. @@ -384,6 +408,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 +928,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 From bc32447b290da05007f87cae945dd9d03cb7ab50 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Wed, 4 Sep 2024 08:58:34 +0200 Subject: [PATCH 3/4] Make aqua happy --- core/Project.toml | 1 + core/test/aqua_test.jl | 9 ++++++++- 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/core/Project.toml b/core/Project.toml index 1268682c6..b5ac7df39 100644 --- a/core/Project.toml +++ b/core/Project.toml @@ -70,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/test/aqua_test.jl b/core/test/aqua_test.jl index 13a0c0f51..5d2ff1fd2 100644 --- a/core/test/aqua_test.jl +++ b/core/test/aqua_test.jl @@ -1,4 +1,11 @@ @testitem "Aqua" begin import Aqua - Aqua.test_all(Ribasim; ambiguities = false, persistent_tasks = false) + using FindFirstFunctions: Guesser + Aqua.test_all( + Ribasim; + ambiguities = false, + persistent_tasks = false, + # TODO: Remove after https://github.com/SciML/FindFirstFunctions.jl/pull/26 + piracies = (treat_as_own = [Guesser],), + ) end From 21454c33b63f1ed20ef395d9c01af67e84fe4b96 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Wed, 4 Sep 2024 13:22:49 +0200 Subject: [PATCH 4/4] Comments adressed --- Manifest.toml | 6 ++++-- Project.toml | 1 + core/src/Ribasim.jl | 3 --- core/src/util.jl | 24 ------------------------ core/test/aqua_test.jl | 9 +-------- core/test/utils_test.jl | 26 +++++++++++++++++++++++++- 6 files changed, 31 insertions(+), 38 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index a76bfac8e..aedaba779 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.10.5" manifest_format = "2.0" -project_hash = "ce9b0fd5a3e37477853fac84d64bb99101f27d9c" +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" 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/src/Ribasim.jl b/core/src/Ribasim.jl index ee2a202e6..06c1bf7c3 100644 --- a/core/src/Ribasim.jl +++ b/core/src/Ribasim.jl @@ -24,9 +24,6 @@ using OrdinaryDiffEqCore: using OrdinaryDiffEqNonlinearSolve: OrdinaryDiffEqNonlinearSolve, relax!, _compute_rhs! using LineSearches: BackTracking -# TODO: Remove after https://github.com/SciML/FindFirstFunctions.jl/pull/26 -using FindFirstFunctions: FindFirstFunctions - # Interface for defining and solving the ODE problem of the physical layer. using SciMLBase: init, diff --git a/core/src/util.jl b/core/src/util.jl index 1ce5dc8a8..4e4c70666 100644 --- a/core/src/util.jl +++ b/core/src/util.jl @@ -65,30 +65,6 @@ function get_level_from_storage(basin::Basin, state_idx::Int, storage) end end -# TODO: Remove after https://github.com/SciML/FindFirstFunctions.jl/pull/26 -function (g::FindFirstFunctions.Guesser)(x) - (; v, idx_prev, linear_lookup) = g - if linear_lookup - f = (x - first(v)) / (last(v) - first(v)) - if isinf(f) - f > 0 ? lastindex(v) : firstindex(v) - else - i_0, i_f = firstindex(v), lastindex(v) - i_approx = f * (i_f - i_0) + i_0 - target_type = typeof(firstindex(v)) - if i_approx >= typemax(target_type) - lastindex(v) + 1 - elseif i_approx <= typemin(target_type) - firstindex(v) - 1 - else - round(target_type, i_approx) - end - end - else - idx_prev[] - end -end - """ For an element `id` and a vector of elements `ids`, get the range of indices of the last consecutive block of `id`. diff --git a/core/test/aqua_test.jl b/core/test/aqua_test.jl index 5d2ff1fd2..13a0c0f51 100644 --- a/core/test/aqua_test.jl +++ b/core/test/aqua_test.jl @@ -1,11 +1,4 @@ @testitem "Aqua" begin import Aqua - using FindFirstFunctions: Guesser - Aqua.test_all( - Ribasim; - ambiguities = false, - persistent_tasks = false, - # TODO: Remove after https://github.com/SciML/FindFirstFunctions.jl/pull/26 - piracies = (treat_as_own = [Guesser],), - ) + Aqua.test_all(Ribasim; ambiguities = false, persistent_tasks = false) end 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