From d6ca23b86731bf437145a5afd80d7326b983cdff Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 15 Oct 2024 11:19:24 +0200 Subject: [PATCH 1/6] Add flow limiter (step_limiter!) --- core/src/config.jl | 29 ++++++-- core/src/solve.jl | 134 +++++++++++++++++++++++++++++------ core/src/util.jl | 9 +++ core/test/run_models_test.jl | 11 ++- 4 files changed, 153 insertions(+), 30 deletions(-) diff --git a/core/src/config.jl b/core/src/config.jl index bfbea84ef..ba6d3e197 100644 --- a/core/src/config.jl +++ b/core/src/config.jl @@ -249,6 +249,18 @@ const algorithms = Dict{String, Type}( "Euler" => Euler, ) +""" +Check whether the given function has a method that accepts the given kwarg. +Note that it is possible that methods exist that accept :a and :b individually, +but not both. +""" +function function_accepts_kwarg(f, kwarg)::Bool + for method in methods(f) + kwarg in Base.kwarg_decl(method) && return true + end + return false +end + "Create an OrdinaryDiffEqAlgorithm from solver config" function algorithm(solver::Solver; u0 = [])::OrdinaryDiffEqAlgorithm algotype = get(algorithms, solver.algorithm, nothing) @@ -258,19 +270,22 @@ function algorithm(solver::Solver; u0 = [])::OrdinaryDiffEqAlgorithm Available options are: ($(options)).") end kwargs = Dict{Symbol, Any}() + if algotype <: OrdinaryDiffEqNewtonAdaptiveAlgorithm kwargs[:nlsolve] = NLNewton(; relax = Ribasim.MonitoredBackTracking(; z_tmp = copy(u0), dz_tmp = copy(u0)), ) end - # not all algorithms support this keyword - kwargs[:autodiff] = solver.autodiff - try - algotype(; kwargs...) - catch - pop!(kwargs, :autodiff) - algotype(; kwargs...) + + if function_accepts_kwarg(algotype, :step_limiter!) + kwargs[:step_limiter!] = Ribasim.limit_flow! end + + if function_accepts_kwarg(algotype, :autodiff) + kwargs[:autodiff] = solver.autodiff + end + + algotype(; kwargs...) end "Convert the saveat Float64 from our Config to SciML's saveat" diff --git a/core/src/solve.jl b/core/src/solve.jl index 9c2555908..da2520c52 100644 --- a/core/src/solve.jl +++ b/core/src/solve.jl @@ -320,29 +320,14 @@ function formulate_flow!( )::Nothing (; allocation) = p all_nodes_active = p.all_nodes_active[] - for ( - id, - inflow_edge, - outflow_edge, - active, - demand_itp, - demand, - allocated, - return_factor, - min_level, - demand_from_timeseries, - ) in zip( + for (id, inflow_edge, outflow_edge, active, allocated, return_factor, min_level) in zip( user_demand.node_id, user_demand.inflow_edge, user_demand.outflow_edge, user_demand.active, - user_demand.demand_itp, - # TODO permute these so the nodes are the last dimension, for performance - eachrow(user_demand.demand), eachrow(user_demand.allocated), user_demand.return_factor, user_demand.min_level, - user_demand.demand_from_timeseries, ) if !(active || all_nodes_active) continue @@ -356,11 +341,7 @@ function formulate_flow!( # effectively allocated = demand. for priority_idx in eachindex(allocation.priorities) alloc_prio = allocated[priority_idx] - demand_prio = if demand_from_timeseries - demand_itp[priority_idx](t) - else - demand[priority_idx] - end + demand_prio = get_demand(user_demand, id, priority_idx, t) alloc = min(alloc_prio, demand_prio) q += alloc end @@ -718,3 +699,114 @@ function formulate_flows!( formulate_flow!(du, user_demand, p, t, current_storage, current_level) end end + +""" +Clamp the cumulative flow states between bounds given my minimum and maximum +flow rates for the last time step if these flow rate bounds are known. +""" +function limit_flow!( + u::ComponentVector, + integrator::DEIntegrator, + p::Parameters, + t::Number, +)::Nothing + (; uprev, dt) = integrator + (; + pump, + outlet, + linear_resistance, + user_demand, + tabulated_rating_curve, + basin, + allocation, + ) = p + + # TabulatedRatingCurve + for (id, active) in zip(tabulated_rating_curve.node_id, tabulated_rating_curve.active) + limit_flow!( + u.tabulated_rating_curve, + uprev.tabulated_rating_curve, + id, + 0.0, + Inf, + active, + dt, + ) + end + + # Pump + for (id, min_flow_rate, max_flow_rate, active) in + zip(pump.node_id, pump.min_flow_rate, pump.max_flow_rate, pump.active) + limit_flow!(u.pump, uprev.pump, id, min_flow_rate, max_flow_rate, active, dt) + end + + # Outlet + for (id, min_flow_rate, max_flow_rate, active) in + zip(outlet.node_id, outlet.min_flow_rate, outlet.max_flow_rate, outlet.active) + limit_flow!(u.outlet, uprev.outlet, id, min_flow_rate, max_flow_rate, active, dt) + end + + # LinearResistance + for (id, max_flow_rate, active) in zip( + linear_resistance.node_id, + linear_resistance.max_flow_rate, + linear_resistance.active, + ) + limit_flow!( + u.linear_resistance, + uprev.linear_resistance, + id, + -max_flow_rate, + max_flow_rate, + active, + dt, + ) + end + + # UserDemand inflow + for (id, active) in zip(user_demand.node_id, user_demand.active) + total_demand = sum( + get_demand(user_demand, id, priority_idx, t) for + priority_idx in eachindex(allocation.priorities) + ) + limit_flow!( + u.user_demand_inflow, + uprev.user_demand_inflow, + id, + 0.0, + total_demand, + active, + dt, + ) + end + + # Evaporation and Infiltration + for (id, infiltration) in zip(basin.node_id, basin.vertical_flux.infiltration) + limit_flow!(u.evaporation, uprev.evaporation, id, 0.0, Inf, true, dt) + limit_flow!(u.infiltration, uprev.infiltration, id, 0.0, infiltration, true, dt) + end + + return nothing +end + +function limit_flow!( + u_component, + uprev_component, + id::NodeID, + min_flow_rate::Number, + max_flow_rate::Number, + active::Bool, + dt::Number, +)::Nothing + u_prev = uprev_component[id.idx] + if active + u_component[id.idx] = clamp( + u_component[id.idx], + u_prev + min_flow_rate * dt, + u_prev + max_flow_rate * dt, + ) + else + u_component[id.idx] = uprev_component[id.idx] + end + return nothing +end diff --git a/core/src/util.jl b/core/src/util.jl index f8f17c8c0..77127e99c 100644 --- a/core/src/util.jl +++ b/core/src/util.jl @@ -1134,3 +1134,12 @@ function isoutofdomain(u, p, t) formulate_storages!(current_storage, u, u, p, t) any(<(0), current_storage) end + +function get_demand(user_demand, id, priority_idx, t)::Float64 + (; demand_from_timeseries, demand_itp, demand) = user_demand + if demand_from_timeseries[id.idx] + demand_itp[id.idx][priority_idx](t) + else + demand[id.idx, priority_idx] + end +end diff --git a/core/test/run_models_test.jl b/core/test/run_models_test.jl index f033b9f29..bfc634ab1 100644 --- a/core/test/run_models_test.jl +++ b/core/test/run_models_test.jl @@ -184,6 +184,7 @@ end using Logging: Debug, with_logger using LoggingExtras using SciMLBase: successful_retcode + using OrdinaryDiffEqBDF: QNDF import Tables using Dates @@ -197,11 +198,17 @@ end end @test model isa Ribasim.Model - p = model.integrator.p + + (; integrator) = model + (; p, alg) = integrator + @test p isa Ribasim.Parameters @test isconcretetype(typeof(p)) @test all(isconcretetype, fieldtypes(typeof(p))) + @test alg isa QNDF + @test alg.step_limiter! == Ribasim.limit_flow! + @test successful_retcode(model) @test length(model.integrator.sol) == 2 # start and end @test model.integrator.p.basin.current_storage[Float64[]] ≈ @@ -242,7 +249,7 @@ end precipitation = model.integrator.p.basin.vertical_flux.precipitation @test length(precipitation) == 4 @test model.integrator.p.basin.current_storage[parent(du)] ≈ - Float32[720.23611, 694.8785, 415.22371, 1334.3859] atol = 2.0 skip = Sys.isapple() + Float32[721.17656, 695.8066, 416.66188, 1334.4879] atol = 2.0 skip = Sys.isapple() end @testitem "Allocation example model" begin From f9d70ef8c6ed76d715926bfc8193aca339e16461 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Mon, 21 Oct 2024 11:43:34 +0200 Subject: [PATCH 2/6] Add a bit more documentation --- core/src/solve.jl | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/core/src/solve.jl b/core/src/solve.jl index da2520c52..3863dc633 100644 --- a/core/src/solve.jl +++ b/core/src/solve.jl @@ -701,7 +701,7 @@ function formulate_flows!( end """ -Clamp the cumulative flow states between bounds given my minimum and maximum +Clamp the cumulative flow states within the minimum and maximum flow rates for the last time step if these flow rate bounds are known. """ function limit_flow!( @@ -721,7 +721,7 @@ function limit_flow!( allocation, ) = p - # TabulatedRatingCurve + # TabulatedRatingCurve flow is in [0, ∞) and can be inactive for (id, active) in zip(tabulated_rating_curve.node_id, tabulated_rating_curve.active) limit_flow!( u.tabulated_rating_curve, @@ -734,19 +734,19 @@ function limit_flow!( ) end - # Pump + # Pump flow is in [min_flow_rate, max_flow_rate] and can be inactive for (id, min_flow_rate, max_flow_rate, active) in zip(pump.node_id, pump.min_flow_rate, pump.max_flow_rate, pump.active) limit_flow!(u.pump, uprev.pump, id, min_flow_rate, max_flow_rate, active, dt) end - # Outlet + # Outlet flow is in [min_flow_rate, max_flow_rate] and can be inactive for (id, min_flow_rate, max_flow_rate, active) in zip(outlet.node_id, outlet.min_flow_rate, outlet.max_flow_rate, outlet.active) limit_flow!(u.outlet, uprev.outlet, id, min_flow_rate, max_flow_rate, active, dt) end - # LinearResistance + # LinearResistance flow is in [-max_flow_rate, max_flow_rate] and can be inactive for (id, max_flow_rate, active) in zip( linear_resistance.node_id, linear_resistance.max_flow_rate, @@ -763,7 +763,7 @@ function limit_flow!( ) end - # UserDemand inflow + # UserDemand inflow is in [0, total_demand] and can be inactive for (id, active) in zip(user_demand.node_id, user_demand.active) total_demand = sum( get_demand(user_demand, id, priority_idx, t) for @@ -780,7 +780,8 @@ function limit_flow!( ) end - # Evaporation and Infiltration + # Evaporation is in [0, ∞) (a stricter upper bound would require also estimating the area) + # Infiltration is in [0, infiltration] for (id, infiltration) in zip(basin.node_id, basin.vertical_flux.infiltration) limit_flow!(u.evaporation, uprev.evaporation, id, 0.0, Inf, true, dt) limit_flow!(u.infiltration, uprev.infiltration, id, 0.0, infiltration, true, dt) From f4b23daac657ce2c6095b518600c11bb8671e291 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Mon, 21 Oct 2024 13:33:09 +0200 Subject: [PATCH 3/6] Refactor reduction factor --- core/src/allocation_optim.jl | 2 +- core/src/bmi.jl | 4 +- core/src/callback.jl | 18 +++--- core/src/parameter.jl | 21 +++++-- core/src/read.jl | 5 -- core/src/solve.jl | 104 ++++++++++++++++++++++++----------- core/src/util.jl | 33 ++++------- core/src/write.jl | 5 +- core/test/io_test.jl | 7 ++- core/test/run_models_test.jl | 14 ++--- core/test/utils_test.jl | 33 ++--------- 11 files changed, 130 insertions(+), 116 deletions(-) diff --git a/core/src/allocation_optim.jl b/core/src/allocation_optim.jl index dc9613688..23b83d8ca 100644 --- a/core/src/allocation_optim.jl +++ b/core/src/allocation_optim.jl @@ -311,7 +311,7 @@ function get_basin_data( (; mean_input_flows) = allocation @assert node_id.type == NodeType.Basin influx = mean_input_flows[(node_id, node_id)][] - storage_basin = basin.current_storage[parent(u)][node_id.idx] + storage_basin = basin.current_properties.current_storage[parent(u)][node_id.idx] control_inneighbors = inneighbor_labels_type(graph, node_id, EdgeType.control) if isempty(control_inneighbors) level_demand_idx = 0 diff --git a/core/src/bmi.jl b/core/src/bmi.jl index 8e5c9960c..c02b27278 100644 --- a/core/src/bmi.jl +++ b/core/src/bmi.jl @@ -33,9 +33,9 @@ end function BMI.get_value_ptr(model::Model, name::AbstractString)::AbstractVector{Float64} (; u, p) = model.integrator if name == "basin.storage" - p.basin.current_storage[parent(u)] + p.basin.current_properties.current_storage[parent(u)] elseif name == "basin.level" - p.basin.current_level[parent(u)] + p.basin.current_properties.current_level[parent(u)] elseif name == "basin.infiltration" p.basin.vertical_flux.infiltration elseif name == "basin.drainage" diff --git a/core/src/callback.jl b/core/src/callback.jl index 3ec332faf..a857fc87d 100644 --- a/core/src/callback.jl +++ b/core/src/callback.jl @@ -312,8 +312,9 @@ function update_cumulative_flows!(u, t, integrator)::Nothing end # Update the Basin concentrations again based on the removed mass - basin.concentration_state .= basin.mass ./ basin.current_storage[parent(u)] - basin.storage_prev .= basin.current_storage[parent(u)] + basin.concentration_state .= + basin.mass ./ basin.current_properties.current_storage[parent(u)] + basin.storage_prev .= basin.current_properties.current_storage[parent(u)] return nothing end @@ -357,8 +358,8 @@ function save_basin_state(u, t, integrator) (; p) = integrator (; basin) = p du = get_du(integrator) - current_storage = basin.current_storage[parent(du)] - current_level = basin.current_level[parent(du)] + current_storage = basin.current_properties.current_storage[parent(du)] + current_level = basin.current_properties.current_level[parent(du)] water_balance!(du, u, p, t) SavedBasinState(; storage = copy(current_storage), level = copy(current_level), t) end @@ -444,7 +445,7 @@ function check_water_balance_error!( (; u, p, t) = integrator (; basin, water_balance_abstol, water_balance_reltol) = p errors = false - current_storage = basin.current_storage[parent(u)] + current_storage = basin.current_properties.current_storage[parent(u)] # The initial storage is irrelevant for the storage rate and can only cause # floating point truncation errors @@ -510,7 +511,8 @@ end function check_negative_storage(u, t, integrator)::Nothing (; basin) = integrator.p - (; node_id, current_storage) = basin + (; node_id, current_properties) = basin + (; current_storage) = current_properties du = get_du(integrator) set_current_basin_properties!(du, u, integrator.p, t) current_storage = current_storage[parent(du)] @@ -733,7 +735,7 @@ end function update_subgrid_level!(integrator)::Nothing (; p) = integrator du = get_du(integrator) - basin_level = p.basin.current_level[parent(du)] + basin_level = p.basin.current_properties.current_level[parent(du)] subgrid = integrator.p.subgrid for (i, (index, interp)) in enumerate(zip(subgrid.basin_index, subgrid.interpolations)) subgrid.level[i] = interp(basin_level[index]) @@ -819,7 +821,7 @@ update_userd_conc!(integrator)::Nothing = function update_allocation!(integrator)::Nothing (; p, t, u) = integrator (; allocation, basin) = p - (; current_storage) = basin + (; current_storage) = basin.current_properties (; allocation_models, mean_input_flows, mean_realized_flows) = allocation # Make sure current storages are up to date diff --git a/core/src/parameter.jl b/core/src/parameter.jl index 8dcf03705..4737d0bcb 100644 --- a/core/src/parameter.jl +++ b/core/src/parameter.jl @@ -299,6 +299,21 @@ In-memory storage of saved instantaneous storages and levels for writing to resu t::Float64 end +""" +Caches of current basin properties +""" +struct CurrentBasinProperties + current_storage::Cache + current_low_storage_factor::Cache + current_level::Cache + current_area::Cache + current_cumulative_precipitation::Cache + current_cumulative_drainage::Cache + function CurrentBasinProperties(n) + new((cache(n) for _ in 1:6)...) + end +end + """ Requirements: @@ -331,11 +346,7 @@ end cumulative_precipitation_saveat::Vector{Float64} = zeros(length(node_id)) cumulative_drainage_saveat::Vector{Float64} = zeros(length(node_id)) # Cache this to avoid recomputation - current_storage::Cache = cache(length(node_id)) - current_level::Cache = cache(length(node_id)) - current_area::Cache = cache(length(node_id)) - current_cumulative_precipitation::Cache = cache(length(node_id)) - current_cumulative_drainage::Cache = cache(length(node_id)) + current_properties::CurrentBasinProperties = CurrentBasinProperties(length(node_id)) # Discrete values for interpolation storage_to_level::Vector{ LinearInterpolationIntInv{ diff --git a/core/src/read.jl b/core/src/read.jl index 2a9ae9bf5..749554c3b 100644 --- a/core/src/read.jl +++ b/core/src/read.jl @@ -574,13 +574,10 @@ end function Basin(db::DB, config::Config, graph::MetaGraph)::Basin node_id = get_ids(db, "Basin") n = length(node_id) - current_level = cache(n) - current_area = cache(n) evaporate_mass = config.solver.evaporate_mass precipitation = zeros(n) potential_evaporation = zeros(n) - evaporation = zeros(n) drainage = zeros(n) infiltration = zeros(n) table = (; precipitation, potential_evaporation, drainage, infiltration) @@ -686,8 +683,6 @@ function Basin(db::DB, config::Config, graph::MetaGraph)::Basin inflow_ids = [collect(inflow_ids(graph, id)) for id in node_id], outflow_ids = [collect(outflow_ids(graph, id)) for id in node_id], vertical_flux, - current_level, - current_area, storage_to_level, level_to_area, demand, diff --git a/core/src/solve.jl b/core/src/solve.jl index 3863dc633..acc6b4806 100644 --- a/core/src/solve.jl +++ b/core/src/solve.jl @@ -8,14 +8,17 @@ function water_balance!( t::Number, )::Nothing (; basin, pid_control) = p + (; current_storage, current_low_storage_factor, current_level) = + basin.current_properties du .= 0.0 # Ensures current_* vectors are current set_current_basin_properties!(du, u, p, t) - current_storage = basin.current_storage[parent(du)] - current_level = basin.current_level[parent(du)] + current_storage = current_storage[parent(du)] + current_low_storage_factor = current_low_storage_factor[parent(du)] + current_level = current_level[parent(du)] # Notes on the ordering of these formulations: # - Continuous control can depend on flows (which are not continuously controlled themselves), @@ -27,7 +30,7 @@ function water_balance!( update_vertical_flux!(basin, du) # Formulate intermediate flows (non continuously controlled) - formulate_flows!(du, p, t, current_storage, current_level) + formulate_flows!(du, p, t, current_storage, current_low_storage_factor, current_level) # Compute continuous control formulate_continuous_control!(du, p, t) @@ -38,6 +41,7 @@ function water_balance!( p, t, current_storage, + current_low_storage_factor, current_level; continuous_control_type = ContinuousControlType.Continuous, ) @@ -51,6 +55,7 @@ function water_balance!( p, t, current_storage, + current_low_storage_factor, current_level; continuous_control_type = ContinuousControlType.PID, ) @@ -79,17 +84,19 @@ function set_current_basin_properties!( t::Number, )::Nothing (; basin) = p + (; current_properties, cumulative_precipitation, cumulative_drainage, vertical_flux) = + basin (; current_storage, + current_low_storage_factor, current_level, current_area, current_cumulative_precipitation, current_cumulative_drainage, - cumulative_precipitation, - cumulative_drainage, - vertical_flux, - ) = basin + ) = current_properties + current_storage = current_storage[parent(du)] + current_low_storage_factor = current_low_storage_factor[parent(du)] current_level = current_level[parent(du)] current_area = current_area[parent(du)] current_cumulative_precipitation = current_cumulative_precipitation[parent(du)] @@ -107,7 +114,9 @@ function set_current_basin_properties!( formulate_storages!(current_storage, du, u, p, t) - for (i, s) in enumerate(current_storage) + for (id, s) in zip(basin.node_id, current_storage) + i = id.idx + current_low_storage_factor[i] = reduction_factor(s, 10.0) current_level[i] = get_level_from_storage(basin, i, s) current_area[i] = basin.level_to_area[i](current_level[i]) end @@ -144,7 +153,8 @@ function formulate_storage!( basin::Basin, du::ComponentVector, ) - (; current_cumulative_precipitation, current_cumulative_drainage) = basin + (; current_cumulative_precipitation, current_cumulative_drainage) = + basin.current_properties current_cumulative_precipitation = current_cumulative_precipitation[parent(du)] current_cumulative_drainage = current_cumulative_drainage[parent(du)] @@ -185,7 +195,8 @@ Smoothly let the evaporation flux go to 0 when at small water depths Currently at less than 0.1 m. """ function update_vertical_flux!(basin::Basin, du::AbstractVector)::Nothing - (; vertical_flux, current_level, current_area) = basin + (; vertical_flux, current_properties) = basin + (; current_level, current_area) = current_properties current_level = current_level[parent(du)] current_area = current_area[parent(du)] @@ -211,7 +222,7 @@ function set_error!(pid_control::PidControl, p::Parameters, du::ComponentVector, (; basin) = p (; listen_node_id, target, error) = pid_control error = error[parent(du)] - current_level = basin.current_level[parent(du)] + current_level = basin.current_properties.current_level[parent(du)] for i in eachindex(listen_node_id) listened_node_id = listen_node_id[i] @@ -229,7 +240,7 @@ function formulate_pid_control!( )::Nothing (; basin) = p (; node_id, active, target, listen_node_id, error) = pid_control - (; current_area) = basin + (; current_area) = basin.current_properties current_area = current_area[parent(du)] error = error[parent(du)] @@ -315,7 +326,7 @@ function formulate_flow!( user_demand::UserDemand, p::Parameters, t::Number, - current_storage::Vector, + current_low_storage_factor::Vector, current_level::Vector, )::Nothing (; allocation) = p @@ -348,7 +359,7 @@ function formulate_flow!( # Smoothly let abstraction go to 0 as the source basin dries out inflow_id = inflow_edge.edge[1] - factor_basin = low_storage_factor(current_storage, inflow_id, 10.0) + factor_basin = get_low_storage_factor(current_low_storage_factor, inflow_id) q *= factor_basin # Smoothly let abstraction go to 0 as the source basin @@ -371,7 +382,7 @@ function formulate_flow!( linear_resistance::LinearResistance, p::Parameters, t::Number, - current_storage::Vector, + current_low_storage_factor::Vector, current_level::Vector, )::Nothing all_nodes_active = p.all_nodes_active[] @@ -389,11 +400,10 @@ function formulate_flow!( 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_resistance_node( - current_storage, + current_low_storage_factor, q, inflow_id, outflow_id, - 10.0, ) du.linear_resistance[id.idx] = q end @@ -409,10 +419,9 @@ function formulate_flow!( tabulated_rating_curve::TabulatedRatingCurve, p::Parameters, t::Number, - current_storage::Vector, + current_low_storage_factor::Vector, current_level::Vector, )::Nothing - (; basin) = p all_nodes_active = p.all_nodes_active[] (; node_id, active, table) = tabulated_rating_curve @@ -428,7 +437,7 @@ function formulate_flow!( Δh = h_a - h_b if active[id.idx] || all_nodes_active - factor = low_storage_factor(current_storage, inflow_id, 10.0) + factor = get_low_storage_factor(current_low_storage_factor, inflow_id) q = factor * table[id.idx](h_a) q *= reduction_factor(Δh, 0.02) q *= reduction_factor(max_downstream_level - h_b, 0.02) @@ -485,7 +494,7 @@ function formulate_flow!( manning_resistance::ManningResistance, p::Parameters, t::Number, - current_storage::Vector, + current_low_storage_factor::Vector, current_level::Vector, )::Nothing (; @@ -540,11 +549,10 @@ function formulate_flow!( q = A / n * ∛(R_h^2) * relaxed_root(Δh / L, 1e-3) q *= low_storage_factor_resistance_node( - current_storage, + current_low_storage_factor, q, inflow_id, outflow_id, - 10.0, ) du.manning_resistance[id.idx] = q end @@ -556,7 +564,7 @@ function formulate_flow!( pump::Pump, p::Parameters, t::Number, - current_storage::Vector, + current_low_storage_factor::Vector, current_level::Vector, continuous_control_type_::ContinuousControlType.T, )::Nothing @@ -594,7 +602,7 @@ function formulate_flow!( src_level = get_level(p, inflow_id, t, current_level) dst_level = get_level(p, outflow_id, t, current_level) - factor = low_storage_factor(current_storage, inflow_id, 10.0) + factor = get_low_storage_factor(current_low_storage_factor, inflow_id) q = flow_rate * factor q *= reduction_factor(src_level - min_upstream_level, 0.02) @@ -611,7 +619,7 @@ function formulate_flow!( outlet::Outlet, p::Parameters, t::Number, - current_storage::Vector, + current_low_storage_factor::Vector, current_level::Vector, continuous_control_type_::ContinuousControlType.T, )::Nothing @@ -650,7 +658,7 @@ function formulate_flow!( dst_level = get_level(p, outflow_id, t, current_level) q = flow_rate - q *= low_storage_factor(current_storage, inflow_id, 10.0) + q *= get_low_storage_factor(current_low_storage_factor, inflow_id) # No flow of outlet if source level is lower than target level Δlevel = src_level - dst_level @@ -669,6 +677,7 @@ function formulate_flows!( p::Parameters, t::Number, current_storage::Vector, + current_low_storage_factor::Vector, current_level::Vector; continuous_control_type::ContinuousControlType.T = ContinuousControlType.None, )::Nothing @@ -681,22 +690,51 @@ function formulate_flows!( user_demand, ) = p - formulate_flow!(du, pump, p, t, current_storage, current_level, continuous_control_type) + formulate_flow!( + du, + pump, + p, + t, + current_low_storage_factor, + current_level, + continuous_control_type, + ) formulate_flow!( du, outlet, p, t, - current_storage, + current_low_storage_factor, current_level, continuous_control_type, ) if continuous_control_type == ContinuousControlType.None - formulate_flow!(du, linear_resistance, p, t, current_storage, current_level) - formulate_flow!(du, manning_resistance, p, t, current_storage, current_level) - formulate_flow!(du, tabulated_rating_curve, p, t, current_storage, current_level) - formulate_flow!(du, user_demand, p, t, current_storage, current_level) + formulate_flow!( + du, + linear_resistance, + p, + t, + current_low_storage_factor, + current_level, + ) + formulate_flow!( + du, + manning_resistance, + p, + t, + current_low_storage_factor, + current_level, + ) + formulate_flow!( + du, + tabulated_rating_curve, + p, + t, + current_low_storage_factor, + current_level, + ) + formulate_flow!(du, user_demand, p, t, current_low_storage_factor, current_level) end end diff --git a/core/src/util.jl b/core/src/util.jl index 77127e99c..536410404 100644 --- a/core/src/util.jl +++ b/core/src/util.jl @@ -375,17 +375,11 @@ function reduction_factor(x::T, threshold::Real)::T where {T <: Real} end end -"If id is a Basin with storage below the threshold, return a reduction factor != 1" -function low_storage_factor( - storage::AbstractVector{T}, +function get_low_storage_factor( + current_low_storage_factor::Vector{T}, id::NodeID, - threshold::Real, -)::T where {T <: Real} - if id.type == NodeType.Basin - reduction_factor(storage[id.idx], threshold) - else - one(T) - end +)::T where {T} + return id.type == NodeType.Basin ? current_low_storage_factor[id.idx] : one(T) end """ @@ -393,16 +387,15 @@ For resistance nodes, give a reduction factor based on the upstream node as defined by the flow direction. """ function low_storage_factor_resistance_node( - current_storage, + current_low_storage_factor, q, inflow_id, outflow_id, - threshold, ) if q > 0 - low_storage_factor(current_storage, inflow_id, threshold) + get_low_storage_factor(current_low_storage_factor, inflow_id) else - low_storage_factor(current_storage, outflow_id, threshold) + get_low_storage_factor(current_low_storage_factor, outflow_id) end end @@ -648,7 +641,7 @@ function get_variable_ref( u = build_state_vector(p) ref = if node_id.type == NodeType.Basin && variable == "level" - PreallocationRef(basin.current_level, node_id.idx) + PreallocationRef(basin.current_properties.current_level, node_id.idx) elseif variable == "flow_rate" && node_id.type != NodeType.FlowBoundary if listen if node_id.type ∉ conservative_nodetypes @@ -861,13 +854,7 @@ end # Custom overloads reduction_factor(x::GradientTracer, threshold::Real) = x -low_storage_factor_resistance_node( - storage, - q::GradientTracer, - inflow_id, - outflow_id, - threshold, -) = q +low_storage_factor_resistance_node(storage, q::GradientTracer, inflow_id, outflow_id) = 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 @@ -1129,7 +1116,7 @@ end Check whether any storages are negative given the state u. """ function isoutofdomain(u, p, t) - (; current_storage) = p.basin + (; current_storage) = p.basin.current_properties current_storage = current_storage[parent(u)] formulate_storages!(current_storage, u, u, p, t) any(<(0), current_storage) diff --git a/core/src/write.jl b/core/src/write.jl index 87afce433..2204771eb 100644 --- a/core/src/write.jl +++ b/core/src/write.jl @@ -105,7 +105,10 @@ function basin_state_table( # ensure the levels are up-to-date set_current_basin_properties!(du, u, p, t) - return (; node_id = Int32.(p.basin.node_id), level = p.basin.current_level[Float64[]]) + return (; + node_id = Int32.(p.basin.node_id), + level = p.basin.current_properties.current_level[Float64[]], + ) end "Create the basin result table from the saved data" diff --git a/core/test/io_test.jl b/core/test/io_test.jl index da482360f..36fe71778 100644 --- a/core/test/io_test.jl +++ b/core/test/io_test.jl @@ -151,9 +151,10 @@ end config = Ribasim.Config(toml_path) model = Ribasim.Model(config) - storage1_begin = copy(model.integrator.p.basin.current_storage[Float64[]]) + storage1_begin = + copy(model.integrator.p.basin.current_properties.current_storage[Float64[]]) solve!(model) - storage1_end = model.integrator.p.basin.current_storage[Float64[]] + storage1_end = model.integrator.p.basin.current_properties.current_storage[Float64[]] @test storage1_begin != storage1_end # copy state results to input @@ -169,6 +170,6 @@ end end model = Ribasim.Model(toml_path) - storage2_begin = model.integrator.p.basin.current_storage[Float64[]] + storage2_begin = model.integrator.p.basin.current_properties.current_storage[Float64[]] @test storage1_end ≈ storage2_begin end diff --git a/core/test/run_models_test.jl b/core/test/run_models_test.jl index bfc634ab1..fbb8486d8 100644 --- a/core/test/run_models_test.jl +++ b/core/test/run_models_test.jl @@ -130,7 +130,7 @@ end model = Ribasim.run(toml_path) @test model isa Ribasim.Model (; basin) = model.integrator.p - @test basin.current_storage[Float64[]] ≈ [1000] + @test basin.current_properties.current_storage[Float64[]] ≈ [1000] @test basin.vertical_flux.precipitation == [0.0] @test basin.vertical_flux.drainage == [0.0] du = get_du(model.integrator) @@ -153,7 +153,7 @@ end du = get_du(integrator) (; u, p, t) = integrator Ribasim.water_balance!(du, u, p, t) - stor = integrator.p.basin.current_storage[parent(du)] + stor = integrator.p.basin.current_properties.current_storage[parent(du)] prec = p.basin.vertical_flux.precipitation evap = du.evaporation drng = p.basin.vertical_flux.drainage @@ -211,7 +211,7 @@ end @test successful_retcode(model) @test length(model.integrator.sol) == 2 # start and end - @test model.integrator.p.basin.current_storage[Float64[]] ≈ + @test model.integrator.p.basin.current_properties.current_storage[Float64[]] ≈ Float32[828.5386, 801.88289, 492.290, 1318.3053] skip = Sys.isapple() atol = 1.5 @test length(logger.logs) > 10 @@ -248,7 +248,7 @@ end du = get_du(model.integrator) precipitation = model.integrator.p.basin.vertical_flux.precipitation @test length(precipitation) == 4 - @test model.integrator.p.basin.current_storage[parent(du)] ≈ + @test model.integrator.p.basin.current_properties.current_storage[parent(du)] ≈ Float32[721.17656, 695.8066, 416.66188, 1334.4879] atol = 2.0 skip = Sys.isapple() end @@ -302,7 +302,7 @@ end model = Ribasim.run(toml_path) @test model isa Ribasim.Model @test successful_retcode(model) - @test model.integrator.p.basin.current_storage[Float64[]] ≈ + @test model.integrator.p.basin.current_properties.current_storage[Float64[]] ≈ Float32[368.31558, 365.68442] skip = Sys.isapple() # the highest level in the dynamic table is updated to 1.2 from the callback @test model.integrator.p.tabulated_rating_curve.table[end].t[end] == 1.2 @@ -411,7 +411,7 @@ end model = Ribasim.Model(toml_path) (; u, p, t, sol) = model.integrator - current_storage = p.basin.current_storage[Float64[]] + current_storage = p.basin.current_properties.current_storage[Float64[]] day = 86400.0 @@ -505,7 +505,7 @@ end du = get_du(model.integrator) (; p, t) = model.integrator - h_actual = p.basin.current_level[parent(du)][1:50] + h_actual = p.basin.current_properties.current_level[parent(du)][1:50] x = collect(10.0:20.0:990.0) h_expected = standard_step_method(x, 5.0, 1.0, 0.04, h_actual[end], 1.0e-6) diff --git a/core/test/utils_test.jl b/core/test/utils_test.jl index f8a3fd12e..74e082f08 100644 --- a/core/test/utils_test.jl +++ b/core/test/utils_test.jl @@ -21,10 +21,6 @@ end level_to_area = LinearInterpolation.(area, level) storage_to_level = invert_integral.(level_to_area) demand = zeros(2) - current_level = cache(2) - current_area = cache(2) - current_level[Float64[]] .= [2.0, 3.0] - current_area[Float64[]] .= [2.0, 3.0] substances = OrderedSet([:test]) concentration_state = zeros(2, 1) @@ -33,8 +29,6 @@ end basin = Ribasim.Basin(; node_id = NodeID.(:Basin, [5, 7], [1, 2]), - current_level, - current_area, storage_to_level, level_to_area, demand, @@ -46,6 +40,11 @@ end concentration_time = StructVector{Ribasim.BasinConcentrationV1}(undef, 0), ) + (; current_level, current_area) = basin.current_properties + + current_level[Float64[]] .= [2.0, 3.0] + current_area[Float64[]] .= [2.0, 3.0] + @test Ribasim.basin_levels(basin, 2)[1] === 4.0 @test Ribasim.basin_bottom(basin, NodeID(:Basin, 5, 1))[2] === 0.0 @test Ribasim.basin_bottom(basin, NodeID(:Basin, 7, 2))[2] === 4.0 @@ -261,28 +260,6 @@ end @test reduction_factor(-Inf, 2.0) === 0.0 end -@testitem "low_storage_factor" begin - 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 - @test low_storage_factor([0.0f0], node_id, 2.0) === 0.0f0 - @test low_storage_factor([0.0], node_id, 2.0) === 0.0 - @test low_storage_factor([1.0f0], node_id, 2.0) === 0.5f0 - @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([3.0, 3.0], 1.0, node_id_1, node_id_2, 2.0) == - 1.0 - @test low_storage_factor_resistance_node([1.0, 3.0], 1.0, node_id_1, node_id_2, 2.0) == - 0.5 - @test low_storage_factor_resistance_node([1.0, 3.0], -1.0, node_id_1, node_id_2, 2.0) == - 1.0 -end - @testitem "constraints_from_nodes" begin using Ribasim: Model, From a27c4619fa18307fe7e5740574ad2c538cfb36db Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Mon, 21 Oct 2024 15:43:12 +0200 Subject: [PATCH 4/6] Estimate lower bounds for user demand inflow and infiltration --- core/src/callback.jl | 1 + core/src/model.jl | 2 ++ core/src/parameter.jl | 5 +++ core/src/solve.jl | 61 ++++++++++++++++++++++++++++-------- core/src/util.jl | 29 +++++++++++++++++ core/test/run_models_test.jl | 3 +- 6 files changed, 87 insertions(+), 14 deletions(-) diff --git a/core/src/callback.jl b/core/src/callback.jl index a857fc87d..f4020f6fc 100644 --- a/core/src/callback.jl +++ b/core/src/callback.jl @@ -315,6 +315,7 @@ function update_cumulative_flows!(u, t, integrator)::Nothing basin.concentration_state .= basin.mass ./ basin.current_properties.current_storage[parent(u)] basin.storage_prev .= basin.current_properties.current_storage[parent(u)] + basin.level_prev .= basin.current_properties.current_level[parent(u)] return nothing end diff --git a/core/src/model.jl b/core/src/model.jl index d4aa9af1e..f8e7ed0cf 100644 --- a/core/src/model.jl +++ b/core/src/model.jl @@ -118,6 +118,8 @@ function Model(config::Config)::Model # Synchronize level with storage set_current_basin_properties!(du0, u0, parameters, t0) + parameters.basin.level_prev .= + parameters.basin.current_properties.current_level[parent(u0)] saveat = convert_saveat(config.solver.saveat, t_end) saveat isa Float64 && push!(tstops, range(0, t_end; step = saveat)) diff --git a/core/src/parameter.jl b/core/src/parameter.jl index 4737d0bcb..e6c6e4bd5 100644 --- a/core/src/parameter.jl +++ b/core/src/parameter.jl @@ -1,3 +1,6 @@ +const LOW_STORAGE_THRESHOLD = 10.0 +const USER_DEMAND_MIN_LEVEL_THRESHOLD = 0.1 + const SolverStats = @NamedTuple{ time::Float64, rhs_calls::Int, @@ -364,6 +367,8 @@ end # Data source for concentration updates concentration_time::StructVector{BasinConcentrationV1, D, Int} + # Level for each Basin at the previous time step + level_prev::Vector{Float64} = zeros(length(node_id)) # Concentrations # Config setting to enable/disable evaporation of mass evaporate_mass::Bool = true diff --git a/core/src/solve.jl b/core/src/solve.jl index acc6b4806..68517a7ee 100644 --- a/core/src/solve.jl +++ b/core/src/solve.jl @@ -116,7 +116,7 @@ function set_current_basin_properties!( for (id, s) in zip(basin.node_id, current_storage) i = id.idx - current_low_storage_factor[i] = reduction_factor(s, 10.0) + current_low_storage_factor[i] = reduction_factor(s, LOW_STORAGE_THRESHOLD) current_level[i] = get_level_from_storage(basin, i, s) current_area[i] = basin.level_to_area[i](current_level[i]) end @@ -366,7 +366,7 @@ function formulate_flow!( # level reaches its minimum level source_level = get_level(p, inflow_id, t, current_level) Δsource_level = source_level - min_level - factor_level = reduction_factor(Δsource_level, 0.1) + factor_level = reduction_factor(Δsource_level, USER_DEMAND_MIN_LEVEL_THRESHOLD) q *= factor_level du.user_demand_inflow[id.idx] = q du.user_demand_outflow[id.idx] = q * return_factor(t) @@ -759,6 +759,11 @@ function limit_flow!( allocation, ) = p + du = get_du(integrator) + set_current_basin_properties!(du, u, p, t) + current_storage = basin.current_properties.current_storage[parent(u)] + current_level = basin.current_properties.current_level[parent(u)] + # TabulatedRatingCurve flow is in [0, ∞) and can be inactive for (id, active) in zip(tabulated_rating_curve.node_id, tabulated_rating_curve.active) limit_flow!( @@ -801,28 +806,58 @@ function limit_flow!( ) end - # UserDemand inflow is in [0, total_demand] and can be inactive - for (id, active) in zip(user_demand.node_id, user_demand.active) - total_demand = sum( - get_demand(user_demand, id, priority_idx, t) for - priority_idx in eachindex(allocation.priorities) - ) + # UserDemand inflow bounds depend on multiple aspects of the simulation + for (id, active, inflow_edge, demand_from_timeseries) in zip( + user_demand.node_id, + user_demand.active, + user_demand.inflow_edge, + user_demand.demand_from_timeseries, + ) + min_flow_rate, max_flow_rate = if demand_from_timeseries + # Bounding the flow rate if the demand comes from a time series is hard + 0, Inf + else + inflow_id = inflow_edge.edge[1] + factor_basin_min = + min_low_storage_factor(current_storage, basin.storage_prev, inflow_id) + factor_level_min = min_low_user_demand_level_factor( + current_level, + basin.level_prev, + user_demand.min_level, + id, + inflow_id, + ) + allocated_total = + is_active(allocation) ? sum(user_demand.allocated[id.idx, :]) : + sum(user_demand.demand[id.idx, :]) + factor_basin_min * factor_level_min * allocated_total, allocated_total + end limit_flow!( u.user_demand_inflow, uprev.user_demand_inflow, id, - 0.0, - total_demand, + min_flow_rate, + max_flow_rate, active, dt, ) end - # Evaporation is in [0, ∞) (a stricter upper bound would require also estimating the area) - # Infiltration is in [0, infiltration] + # Evaporation is in [0, ∞) (stricter bounds would require also estimating the area) + # Infiltration is in [f * infiltration, infiltration] where f is a rough estimate of the smallest low storage factor + # reduction factor value that was attained over the last timestep for (id, infiltration) in zip(basin.node_id, basin.vertical_flux.infiltration) + factor_min = min_low_storage_factor(current_storage, basin.storage_prev, id) limit_flow!(u.evaporation, uprev.evaporation, id, 0.0, Inf, true, dt) - limit_flow!(u.infiltration, uprev.infiltration, id, 0.0, infiltration, true, dt) + limit_flow!( + u.infiltration, + uprev.infiltration, + id, + factor_min * infiltration, + infiltration, + true, + dt, + ) end return nothing diff --git a/core/src/util.jl b/core/src/util.jl index 536410404..1aadbad63 100644 --- a/core/src/util.jl +++ b/core/src/util.jl @@ -1130,3 +1130,32 @@ function get_demand(user_demand, id, priority_idx, t)::Float64 demand[id.idx, priority_idx] end end + +function min_low_storage_factor(storage_now::Vector{T}, storage_prev, id) where {T} + if id.type == NodeType.Basin + reduction_factor( + min(storage_now[id.idx], storage_prev[id.idx]) - 2LOW_STORAGE_THRESHOLD, + LOW_STORAGE_THRESHOLD, + ) + else + one(T) + end +end + +function min_low_user_demand_level_factor( + level_now::Vector{T}, + level_prev, + min_level, + id_user_demand, + id_inflow, +) where {T} + if id_inflow.type == NodeType.Basin + reduction_factor( + min(level_now[id_inflow.idx], level_prev[id_inflow.idx]) - + min_level[id_user_demand.idx] - 2USER_DEMAND_MIN_LEVEL_THRESHOLD, + USER_DEMAND_MIN_LEVEL_THRESHOLD, + ) + else + one(T) + end +end diff --git a/core/test/run_models_test.jl b/core/test/run_models_test.jl index fbb8486d8..a8a1a92a9 100644 --- a/core/test/run_models_test.jl +++ b/core/test/run_models_test.jl @@ -410,7 +410,8 @@ end @test ispath(toml_path) model = Ribasim.Model(toml_path) - (; u, p, t, sol) = model.integrator + (; integrator) = model + (; u, p, t, sol) = integrator current_storage = p.basin.current_properties.current_storage[Float64[]] day = 86400.0 From 500691a3ab3425bb8b83a72b2d90a389fd931c0a Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Mon, 21 Oct 2024 15:54:17 +0200 Subject: [PATCH 5/6] Add some docstrings --- core/src/util.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/core/src/util.jl b/core/src/util.jl index 1aadbad63..a6aad4d96 100644 --- a/core/src/util.jl +++ b/core/src/util.jl @@ -1131,6 +1131,11 @@ function get_demand(user_demand, id, priority_idx, t)::Float64 end end +""" +Estimate the minimum reduction factor achieved over the last time step by +estimating the lowest storage achieved over the last time step. To make sure +it is an underestimate of the minimum, 2LOW_STORAGE_THRESHOLD is subtracted from this lowest storage. +""" function min_low_storage_factor(storage_now::Vector{T}, storage_prev, id) where {T} if id.type == NodeType.Basin reduction_factor( @@ -1142,6 +1147,11 @@ function min_low_storage_factor(storage_now::Vector{T}, storage_prev, id) where end end +""" +Estimate the minimum minimum level reduction factor achieved over the last time step by +estimating the lowest level achieved over the last time step. To make sure +it is an underestimate of the minimum, 2USER_DEMAND_MIN_LEVEL_THRESHOLD is subtracted from this lowest level. +""" function min_low_user_demand_level_factor( level_now::Vector{T}, level_prev, From 47a3c6c2ee965500c23e2390f00cb4e2cc2d18ed Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Thu, 24 Oct 2024 16:24:33 +0200 Subject: [PATCH 6/6] Comments adressed --- core/src/model.jl | 3 +++ core/src/parameter.jl | 6 +++++ core/src/solve.jl | 5 +++++ core/src/util.jl | 4 +++- core/test/run_models_test.jl | 43 ++++++++++++++++++++++++++++++++++++ 5 files changed, 60 insertions(+), 1 deletion(-) diff --git a/core/src/model.jl b/core/src/model.jl index f8e7ed0cf..92d7ac3a2 100644 --- a/core/src/model.jl +++ b/core/src/model.jl @@ -118,6 +118,9 @@ function Model(config::Config)::Model # Synchronize level with storage set_current_basin_properties!(du0, u0, parameters, t0) + + # Previous level is used to estimate the minimum level that was attained during a time step + # in limit_flow! parameters.basin.level_prev .= parameters.basin.current_properties.current_level[parent(u0)] diff --git a/core/src/parameter.jl b/core/src/parameter.jl index d35621605..97bc1cab3 100644 --- a/core/src/parameter.jl +++ b/core/src/parameter.jl @@ -1,4 +1,8 @@ + +# Universal reduction factor threshold for the low storage factor const LOW_STORAGE_THRESHOLD = 10.0 + +# Universal reduction factor threshold for the minimum upstream level of UserDemand nodes const USER_DEMAND_MIN_LEVEL_THRESHOLD = 0.1 const SolverStats = @NamedTuple{ @@ -308,6 +312,8 @@ Caches of current basin properties """ struct CurrentBasinProperties current_storage::Cache + # Low storage factor for reducing flows out of drying basins + # given the current storages current_low_storage_factor::Cache current_level::Cache current_area::Cache diff --git a/core/src/solve.jl b/core/src/solve.jl index 68517a7ee..6d306fdf3 100644 --- a/core/src/solve.jl +++ b/core/src/solve.jl @@ -759,6 +759,9 @@ function limit_flow!( allocation, ) = p + # The current storage and level based on the proposed u are used to estimate the lowest + # storage and level attained in the last time step to estimate whether there was an effect + # of reduction factors du = get_du(integrator) set_current_basin_properties!(du, u, p, t) current_storage = basin.current_properties.current_storage[parent(u)] @@ -817,6 +820,8 @@ function limit_flow!( # Bounding the flow rate if the demand comes from a time series is hard 0, Inf else + # The lower bound is estimated as the lowest inflow given the minimum values + # of the reduction factors involved (with a margin) inflow_id = inflow_edge.edge[1] factor_basin_min = min_low_storage_factor(current_storage, basin.storage_prev, inflow_id) diff --git a/core/src/util.jl b/core/src/util.jl index a6aad4d96..d12498c9b 100644 --- a/core/src/util.jl +++ b/core/src/util.jl @@ -1135,6 +1135,7 @@ end Estimate the minimum reduction factor achieved over the last time step by estimating the lowest storage achieved over the last time step. To make sure it is an underestimate of the minimum, 2LOW_STORAGE_THRESHOLD is subtracted from this lowest storage. +This is done to not be too strict in clamping the flow in the limiter """ function min_low_storage_factor(storage_now::Vector{T}, storage_prev, id) where {T} if id.type == NodeType.Basin @@ -1148,9 +1149,10 @@ function min_low_storage_factor(storage_now::Vector{T}, storage_prev, id) where end """ -Estimate the minimum minimum level reduction factor achieved over the last time step by +Estimate the minimum level reduction factor achieved over the last time step by estimating the lowest level achieved over the last time step. To make sure it is an underestimate of the minimum, 2USER_DEMAND_MIN_LEVEL_THRESHOLD is subtracted from this lowest level. +This is done to not be too strict in clamping the flow in the limiter """ function min_low_user_demand_level_factor( level_now::Vector{T}, diff --git a/core/test/run_models_test.jl b/core/test/run_models_test.jl index a8a1a92a9..fc6c5e66d 100644 --- a/core/test/run_models_test.jl +++ b/core/test/run_models_test.jl @@ -593,3 +593,46 @@ end @test all(flow .≈ 1.0) @test length(flow) == length(tstops) - 1 end + +@testitem "stroboscopic_forcing" begin + import BasicModelInterface as BMI + using SciMLBase: successful_retcode + + toml_path = normpath(@__DIR__, "../../generated_testmodels/bucket/ribasim.toml") + model = BMI.initialize(Ribasim.Model, toml_path) + + drn = BMI.get_value_ptr(model, "basin.drainage") + drn_sum = BMI.get_value_ptr(model, "basin.cumulative_drainage") + inf = BMI.get_value_ptr(model, "basin.infiltration") + inf_sum = BMI.get_value_ptr(model, "basin.cumulative_infiltration") + + nday = 300 + inf_out = fill(NaN, nday) + drn_out = fill(NaN, nday) + + Δt::Float64 = 86400.0 + + for day in 0:(nday - 1) + if iseven(day) + drn .= 25.0 / Δt + inf .= 0.0 + else + drn .= 0.0 + inf .= 25.0 / Δt + end + BMI.update_until(model, day * Δt) + + inf_out[day + 1] = only(inf_sum) + drn_out[day + 1] = only(drn_sum) + end + + @test successful_retcode(model) + + Δdrn = diff(drn_out) + Δinf = diff(inf_out) + + @test all(Δdrn[1:2:end] .== 0.0) + @test all(isapprox.(Δdrn[2:2:end], 25.0; atol = 1e-10)) + @test all(isapprox.(Δinf[1:2:end], 25.0; atol = 1e-10)) + @test all(Δinf[2:2:end] .== 0.0) +end