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