diff --git a/core/src/allocation_optim.jl b/core/src/allocation_optim.jl index fc45ede27..c02b9e1e6 100644 --- a/core/src/allocation_optim.jl +++ b/core/src/allocation_optim.jl @@ -499,6 +499,7 @@ to the capacity of the outflow source. function adjust_capacities_returnflow!( allocation_model::AllocationModel, p::Parameters, + t::Float64, )::Nothing (; graph, user_demand) = p (; problem) = allocation_model @@ -509,7 +510,7 @@ function adjust_capacities_returnflow!( constraint = constraints_outflow[node_id] capacity = JuMP.normalized_rhs(constraint) + - user_demand.return_factor[node_id.idx] * + user_demand.return_factor[node_id.idx](t) * JuMP.value(F[(inflow_id(graph, node_id), node_id)]) JuMP.set_normalized_rhs(constraint, capacity) @@ -973,7 +974,7 @@ function optimize_priority!( adjust_capacities_edge!(allocation_model) adjust_capacities_basin!(allocation_model) adjust_capacities_buffer!(allocation_model) - adjust_capacities_returnflow!(allocation_model, p) + adjust_capacities_returnflow!(allocation_model, p, t) # Adjust demands for next optimization (in case of internal_sources -> collect_demands) for parameter in propertynames(p) diff --git a/core/src/parameter.jl b/core/src/parameter.jl index b3297e1a9..634de4e15 100644 --- a/core/src/parameter.jl +++ b/core/src/parameter.jl @@ -686,7 +686,7 @@ min_level: The level of the source basin below which the UserDemand does not abs demand_itp::Vector{Vector{ScalarInterpolation}} demand_from_timeseries::Vector{Bool} allocated::Matrix{Float64} - return_factor::Vector{Float64} + return_factor::Vector{ScalarInterpolation} min_level::Vector{Float64} end diff --git a/core/src/read.jl b/core/src/read.jl index 96231cfa7..a5d656b8a 100644 --- a/core/src/read.jl +++ b/core/src/read.jl @@ -921,7 +921,7 @@ function user_demand_static!( active::Vector{Bool}, demand::Matrix{Float64}, demand_itp::Vector{Vector{ScalarInterpolation}}, - return_factor::Vector{Float64}, + return_factor::Vector{ScalarInterpolation}, min_level::Vector{Float64}, static::StructVector{UserDemandStaticV1}, ids::Vector{Int32}, @@ -932,7 +932,12 @@ function user_demand_static!( user_demand_idx = searchsortedfirst(ids, first_row.node_id) active[user_demand_idx] = coalesce(first_row.active, true) - return_factor[user_demand_idx] = first_row.return_factor + return_factor_old = return_factor[user_demand_idx] + return_factor[user_demand_idx] = LinearInterpolation( + fill(first_row.return_factor, 2), + return_factor_old.t; + extrapolate = true, + ) min_level[user_demand_idx] = first_row.min_level for row in group @@ -955,7 +960,7 @@ function user_demand_time!( demand::Matrix{Float64}, demand_itp::Vector{Vector{ScalarInterpolation}}, demand_from_timeseries::Vector{Bool}, - return_factor::Vector{Float64}, + return_factor::Vector{ScalarInterpolation}, min_level::Vector{Float64}, time::StructVector{UserDemandTimeV1}, ids::Vector{Int32}, @@ -971,11 +976,24 @@ function user_demand_time!( active[user_demand_idx] = true demand_from_timeseries[user_demand_idx] = true - return_factor[user_demand_idx] = first_row.return_factor + return_factor_itp, is_valid_return = get_scalar_interpolation( + config.starttime, + t_end, + StructVector(group), + NodeID(:UserDemand, first_row.node_id, 0), + :return_factor; + ) + if is_valid_return + return_factor[user_demand_idx] = return_factor_itp + else + @error "The return_factor(t) relationship for UserDemand $(first_row.node_id) from the time table has repeated timestamps, this can not be interpolated." + errors = true + end + min_level[user_demand_idx] = first_row.min_level priority_idx = findsorted(priorities, first_row.priority) - demand_p_itp, is_valid = get_scalar_interpolation( + demand_p_itp, is_valid_demand = get_scalar_interpolation( config.starttime, t_end, StructVector(group), @@ -985,10 +1003,10 @@ function user_demand_time!( ) demand[user_demand_idx, priority_idx] = demand_p_itp(0.0) - if is_valid + if is_valid_demand demand_itp[user_demand_idx][priority_idx] = demand_p_itp else - @error "The demand(t) relationship for UserDemand $node_id of priority $p from the time table has repeated timestamps, this can not be interpolated." + @error "The demand(t) relationship for UserDemand $(first_row.node_id) of priority $(first_row.priority_idx) from the time table has repeated timestamps, this can not be interpolated." errors = true end end @@ -1022,7 +1040,8 @@ function UserDemand(db::DB, config::Config, graph::MetaGraph)::UserDemand ] demand_from_timeseries = fill(false, n_user) allocated = fill(Inf, n_user, n_priority) - return_factor = zeros(n_user) + return_factor = + [LinearInterpolation(zeros(2), trivial_timespan) for i in eachindex(node_ids)] min_level = zeros(n_user) # Process static table diff --git a/core/src/solve.jl b/core/src/solve.jl index adfe45c62..cc941705c 100644 --- a/core/src/solve.jl +++ b/core/src/solve.jl @@ -273,7 +273,7 @@ function formulate_flow!( set_flow!(graph, inflow_edge, q) # Return flow is immediate - set_flow!(graph, outflow_edge, q * return_factor) + set_flow!(graph, outflow_edge, q * return_factor(t)) end return nothing end diff --git a/core/test/run_models_test.jl b/core/test/run_models_test.jl index 32ec8426a..6e911cf02 100644 --- a/core/test/run_models_test.jl +++ b/core/test/run_models_test.jl @@ -385,18 +385,34 @@ end @testitem "UserDemand" begin using SciMLBase: successful_retcode + using Dates + using DataFrames: DataFrame toml_path = normpath(@__DIR__, "../../generated_testmodels/user_demand/ribasim.toml") @test ispath(toml_path) model = Ribasim.run(toml_path) @test successful_retcode(model) - day = 86400.0 - @test only(model.integrator.sol(0day)) == 1000.0 - # constant UserDemand withdraws to 0.9m/900m3 - @test only(model.integrator.sol(150day)) ≈ 900 atol = 5 - # dynamic UserDemand withdraws to 0.5m/509m3 - @test only(model.integrator.sol(180day)) ≈ 509 atol = 1 + seconds_in_day = 86400.0 + @test only(model.integrator.sol(0seconds_in_day)) == 1000.0 + # constant UserDemand withdraws to 0.9m or 900m3 due to min level = 0.9 + @test only(model.integrator.sol(150seconds_in_day)) ≈ 900 atol = 5 + # dynamic UserDemand withdraws to 0.5m or 500m3 due to min level = 0.5 + @test only(model.integrator.sol(220seconds_in_day)) ≈ 500 atol = 1 + + # Trasient return factor + flow = DataFrame(Ribasim.flow_table(model)) + return_factor_itp = model.integrator.p.user_demand.return_factor[3] + flow_in = + filter([:from_node_id, :to_node_id] => (from, to) -> (from, to) == (1, 4), flow) + flow_out = + filter([:from_node_id, :to_node_id] => (from, to) -> (from, to) == (4, 5), flow) + time_seconds = Ribasim.seconds_since.(flow_in.time, model.config.starttime) + @test isapprox( + flow_out.flow_rate, + return_factor_itp.(time_seconds) .* flow_in.flow_rate, + rtol = 1e-1, + ) end @testitem "ManningResistance" begin diff --git a/docs/guide/examples.ipynb b/docs/guide/examples.ipynb index 6bc12fefb..abfb25aeb 100644 --- a/docs/guide/examples.ipynb +++ b/docs/guide/examples.ipynb @@ -189,7 +189,7 @@ "source": [ "q = 10 / 86400 # 10 m³/day\n", "tabulated_rating_curve4 = model.tabulated_rating_curve.add(\n", - " Node(8, Point(4.0, 0.0)),\n", + " Node(8, Point(3.0, -1.0)),\n", " [\n", " tabulated_rating_curve.Static(\n", " level=[0.0, 1.0],\n", @@ -207,7 +207,7 @@ " ],\n", ")\n", "tabulated_rating_curve8 = model.tabulated_rating_curve.add(\n", - " Node(4, Point(3.0, -1.0)),\n", + " Node(4, Point(4.0, 0.0)),\n", " [\n", " tabulated_rating_curve.Static(\n", " level=[0.0, 1.0],\n", @@ -2339,7 +2339,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.12.4" + "version": "3.12.5" } }, "nbformat": 4, diff --git a/python/ribasim_api/tests/test_bmi.py b/python/ribasim_api/tests/test_bmi.py index 23aa3e772..bdd43b236 100644 --- a/python/ribasim_api/tests/test_bmi.py +++ b/python/ribasim_api/tests/test_bmi.py @@ -161,7 +161,7 @@ def test_get_value_ptr_allocation(libribasim, user_demand, tmp_path): # Demand actual_demand = libribasim.get_value_ptr("user_demand.demand") - expected_demand = np.array([1e-4, 0.0]) + expected_demand = np.array([1e-4, 0.0, 0.0]) assert_array_almost_equal(actual_demand, expected_demand) # Run model for a while to build up realized @@ -169,7 +169,7 @@ def test_get_value_ptr_allocation(libribasim, user_demand, tmp_path): # Realized actual_realized = libribasim.get_value_ptr("user_demand.realized") - expected_realized = np.array([1e-4 * 60.0, 0.0]) + expected_realized = np.array([1e-4 * 60.0, 0.0, 0.0]) assert_array_almost_equal(actual_realized, expected_realized) diff --git a/python/ribasim_testmodels/ribasim_testmodels/allocation.py b/python/ribasim_testmodels/ribasim_testmodels/allocation.py index bce1b660f..d49171bc3 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/allocation.py +++ b/python/ribasim_testmodels/ribasim_testmodels/allocation.py @@ -59,12 +59,31 @@ def user_demand_model() -> Model: ) ], ) - model.terminal.add(Node(4, Point(2, 0))) + model.user_demand.add( + Node(4, Point(1, 0)), + [ + user_demand.Time( + time=[ + "2020-08-01 00:00:00", + "2020-09-01 00:00:00", + "2020-10-01 00:00:00", + "2020-11-01 00:00:00", + ], + min_level=0.0, + demand=[0.0, 1e-4, 2e-4, 0.0], + return_factor=[0.0, 0.1, 0.2, 0.3], + priority=1, + ) + ], + ) + model.terminal.add(Node(5, Point(2, 0))) model.edge.add(model.basin[1], model.user_demand[2]) model.edge.add(model.basin[1], model.user_demand[3]) - model.edge.add(model.user_demand[2], model.terminal[4]) - model.edge.add(model.user_demand[3], model.terminal[4]) + model.edge.add(model.basin[1], model.user_demand[4]) + model.edge.add(model.user_demand[2], model.terminal[5]) + model.edge.add(model.user_demand[3], model.terminal[5]) + model.edge.add(model.user_demand[4], model.terminal[5]) return model