From 93366cbda3ca688092ceb0578be3e4b6652e9108 Mon Sep 17 00:00:00 2001 From: Bart de Koning <74617371+SouthEndMusic@users.noreply.github.com> Date: Thu, 29 Feb 2024 19:30:00 +0100 Subject: [PATCH] add realized user demand and expose via BMI (#1190) Fixes https://github.com/Deltares/Ribasim/issues/1098. @visr @HendrikKok would it be OK if `realized` is not an average but a cumulative volume and resetting it to `0.0` is the responsibility of the coupler? Otherwise we need a discussion about where in the Ribasim code these things must happen in such a way that dividing by the time interval happens before `realized` is accessed and resetting afterwards. Also good to know that `demand` is per user per priority, but `realized` is only per user. **For BMI variable name upgrade instructions see:** https://github.com/Deltares/Ribasim/pull/1190/files#diff-88d699c34d26c16edb268bdacbda2afb8dad704b8db2931a83e04b9ca52dc31a. --------- Co-authored-by: Martijn Visser --- core/src/bmi.jl | 16 +++++++------- core/src/callback.jl | 9 +++++++- core/src/parameter.jl | 4 ++++ core/src/read.jl | 2 ++ core/src/solve.jl | 7 +++++-- core/test/bmi_test.jl | 31 +++++++++++++++++++++++++--- core/test/libribasim_test.jl | 2 +- core/test/validation_test.jl | 1 + python/ribasim_api/tests/test_bmi.py | 12 +++++------ 9 files changed, 64 insertions(+), 20 deletions(-) diff --git a/core/src/bmi.jl b/core/src/bmi.jl index af965c996..39a77afb7 100644 --- a/core/src/bmi.jl +++ b/core/src/bmi.jl @@ -32,18 +32,20 @@ function BMI.update_until(model::Model, time)::Model end function BMI.get_value_ptr(model::Model, name::AbstractString) - if name == "volume" + if name == "basin.storage" model.integrator.u.storage - elseif name == "level" + elseif name == "basin.level" get_tmp(model.integrator.p.basin.current_level, 0) - elseif name == "infiltration" + elseif name == "basin.infiltration" model.integrator.p.basin.infiltration - elseif name == "drainage" + elseif name == "basin.drainage" model.integrator.p.basin.drainage - elseif name == "subgrid_level" + elseif name == "basin.subgrid_level" model.integrator.p.subgrid.level - elseif name == "demand" - model.integrator.p.demand + elseif name == "user_demand.demand" + model.integrator.p.user_demand.demand + elseif name == "user_demand.realized" + model.integrator.p.user_demand.realized_bmi else error("Unknown variable $name") end diff --git a/core/src/callback.jl b/core/src/callback.jl index 7b312f409..3fb4107a3 100644 --- a/core/src/callback.jl +++ b/core/src/callback.jl @@ -102,9 +102,10 @@ Integrate flows over the last timestep """ function integrate_flows!(u, t, integrator)::Nothing (; p, dt) = integrator - (; graph) = p + (; graph, user_demand) = p (; flow, + flow_dict, flow_vertical, flow_prev, flow_vertical_prev, @@ -123,6 +124,12 @@ function integrate_flows!(u, t, integrator)::Nothing @. flow_integrated += 0.5 * (flow + flow_prev) * dt @. flow_vertical_integrated += 0.5 * (flow_vertical + flow_vertical_prev) * dt + for (i, id) in enumerate(user_demand.node_id) + src_id = inflow_id(graph, id) + flow_idx = flow_dict[src_id, id] + user_demand.realized_bmi[i] += 0.5 * (flow[flow_idx] + flow_prev[flow_idx]) * dt + end + copyto!(flow_prev, flow) copyto!(flow_vertical_prev, flow_vertical) return nothing diff --git a/core/src/parameter.jl b/core/src/parameter.jl index be333a8a9..6d02cbceb 100644 --- a/core/src/parameter.jl +++ b/core/src/parameter.jl @@ -480,6 +480,7 @@ end demand: water flux demand of UserDemand per priority over time. Each UserDemand has a demand for all priorities, which is 0.0 if it is not provided explicitly. +realized_bmi: Cumulative inflow volume, for read or reset by BMI only. active: whether this node is active and thus demands water allocated: water flux currently allocated to UserDemand per priority return_factor: the factor in [0,1] of how much of the abstracted water is given back to the system @@ -488,6 +489,7 @@ min_level: The level of the source basin below which the UserDemand does not abs struct UserDemand <: AbstractParameterNode node_id::Vector{NodeID} active::BitVector + realized_bmi::Vector{Float64} demand::Vector{Float64} demand_itp::Vector{Vector{ScalarInterpolation}} demand_from_timeseries::BitVector @@ -498,6 +500,7 @@ struct UserDemand <: AbstractParameterNode function UserDemand( node_id, active, + realized_bmi, demand, demand_itp, demand_from_timeseries, @@ -510,6 +513,7 @@ struct UserDemand <: AbstractParameterNode return new( node_id, active, + realized_bmi, demand, demand_itp, demand_from_timeseries, diff --git a/core/src/read.jl b/core/src/read.jl index 71abd2f38..4e4bdcf55 100644 --- a/core/src/read.jl +++ b/core/src/read.jl @@ -742,11 +742,13 @@ function UserDemand(db::DB, config::Config)::UserDemand error("Errors occurred when parsing UserDemand data.") end + realized_bmi = zeros(length(node_ids)) allocated = [fill(Inf, length(priorities)) for id in node_ids] return UserDemand( node_ids, active, + realized_bmi, demand, demand_itp, demand_from_timeseries, diff --git a/core/src/solve.jl b/core/src/solve.jl index b090dfa26..df9f09d4d 100644 --- a/core/src/solve.jl +++ b/core/src/solve.jl @@ -264,7 +264,8 @@ function formulate_flow!( t::Number, )::Nothing (; graph, basin) = p - (; node_id, allocated, active, demand_itp, return_factor, min_level) = user_demand + (; node_id, allocated, active, demand, demand_itp, return_factor, min_level) = + user_demand for (i, id) in enumerate(node_id) src_id = inflow_id(graph, id) @@ -281,7 +282,9 @@ function formulate_flow!( # If allocation is not optimized then allocated = Inf, so the result is always # effectively allocated = demand. for priority_idx in eachindex(allocated[i]) - alloc = min(allocated[i][priority_idx], demand_itp[i][priority_idx](t)) + alloc_prio = allocated[i][priority_idx] + demand_prio = demand_itp[i][priority_idx](t) + alloc = min(alloc_prio, demand_prio) q += alloc end diff --git a/core/test/bmi_test.jl b/core/test/bmi_test.jl index cffe4c8cb..fd33ed312 100644 --- a/core/test/bmi_test.jl +++ b/core/test/bmi_test.jl @@ -43,11 +43,11 @@ end toml_path = normpath(@__DIR__, "../../generated_testmodels/basic/ribasim.toml") model = BMI.initialize(Ribasim.Model, toml_path) - storage0 = BMI.get_value_ptr(model, "volume") + storage0 = BMI.get_value_ptr(model, "basin.storage") @test storage0 ≈ ones(4) @test_throws "Unknown variable foo" BMI.get_value_ptr(model, "foo") BMI.update_until(model, 86400.0) - storage = BMI.get_value_ptr(model, "volume") + storage = BMI.get_value_ptr(model, "basin.storage") # get_value_ptr does not copy @test storage0 === storage != ones(4) end @@ -58,7 +58,15 @@ end toml_path = normpath(@__DIR__, "../../generated_testmodels/basic/ribasim.toml") model = BMI.initialize(Ribasim.Model, toml_path) - for name in ["volume", "level", "infiltration", "drainage"] + for name in [ + "basin.storage", + "basin.level", + "basin.infiltration", + "basin.drainage", + "basin.subgrid_level", + "user_demand.demand", + "user_demand.realized", + ] value_first = BMI.get_value_ptr(model, name) BMI.update_until(model, 86400.0) value_second = BMI.get_value_ptr(model, name) @@ -66,3 +74,20 @@ end @test value_first === value_second end end + +@testitem "realized_user_demand" begin + import BasicModelInterface as BMI + + toml_path = + normpath(@__DIR__, "../../generated_testmodels/minimal_subnetwork/ribasim.toml") + @test ispath(toml_path) + model = BMI.initialize(Ribasim.Model, toml_path) + demand = BMI.get_value_ptr(model, "user_demand.demand") + realized = BMI.get_value_ptr(model, "user_demand.realized") + day = 86400.0 + BMI.update_until(model, 0.4day) + demand = 0.001 # for both users at the start + @test all(isapprox.(realized, demand * 0.4day; rtol = 1e-3)) + BMI.update_until(model, 0.6day) + @test all(isapprox.(realized, demand * 0.6day; rtol = 1e-3)) +end diff --git a/core/test/libribasim_test.jl b/core/test/libribasim_test.jl index 1eb42f431..0d611a42e 100644 --- a/core/test/libribasim_test.jl +++ b/core/test/libribasim_test.jl @@ -3,7 +3,7 @@ # data from which we create pointers for libribasim time = [-1.0] - var_name = "volume" + var_name = "basin.storage" type = ones(UInt8, 8) GC.@preserve time var_name value type toml_path begin diff --git a/core/test/validation_test.jl b/core/test/validation_test.jl index c8f175741..2020c27e9 100644 --- a/core/test/validation_test.jl +++ b/core/test/validation_test.jl @@ -417,6 +417,7 @@ end [NodeID(:UserDemand, 1)], [true], [0.0], + [0.0], [[LinearInterpolation([-5.0, -5.0], [-1.8, 1.8])]], [true], [0.0, -0.0], diff --git a/python/ribasim_api/tests/test_bmi.py b/python/ribasim_api/tests/test_bmi.py index 9d8664d6d..101bf7ed5 100644 --- a/python/ribasim_api/tests/test_bmi.py +++ b/python/ribasim_api/tests/test_bmi.py @@ -62,7 +62,7 @@ def test_update_subgrid_level(libribasim, basic, tmp_path): config_file = str(tmp_path / "ribasim.toml") libribasim.initialize(config_file) libribasim.update_subgrid_level() - level = libribasim.get_value_ptr("subgrid_level") + level = libribasim.get_value_ptr("basin.subgrid_level") # The subgrid levels are initialized with NaN. # After calling update, they should have regular values. assert np.isfinite(level).all() @@ -72,7 +72,7 @@ def test_get_var_type(libribasim, basic, tmp_path): basic.write(tmp_path / "ribasim.toml") config_file = str(tmp_path / "ribasim.toml") libribasim.initialize(config_file) - var_type = libribasim.get_var_type("volume") + var_type = libribasim.get_var_type("basin.storage") assert var_type == "double" @@ -80,7 +80,7 @@ def test_get_var_rank(libribasim, basic, tmp_path): basic.write(tmp_path / "ribasim.toml") config_file = str(tmp_path / "ribasim.toml") libribasim.initialize(config_file) - actual_rank = libribasim.get_var_rank("volume") + actual_rank = libribasim.get_var_rank("basin.storage") expected_rank = 1 assert_array_almost_equal(actual_rank, expected_rank) @@ -89,7 +89,7 @@ def test_get_var_shape(libribasim, basic, tmp_path): basic.write(tmp_path / "ribasim.toml") config_file = str(tmp_path / "ribasim.toml") libribasim.initialize(config_file) - actual_shape = libribasim.get_var_shape("volume") + actual_shape = libribasim.get_var_shape("basin.storage") expected_shape = np.array([4]) assert_array_almost_equal(actual_shape, expected_shape) @@ -98,7 +98,7 @@ def test_get_value_ptr(libribasim, basic, tmp_path): basic.write(tmp_path / "ribasim.toml") config_file = str(tmp_path / "ribasim.toml") libribasim.initialize(config_file) - actual_volume = libribasim.get_value_ptr("volume") + actual_volume = libribasim.get_value_ptr("basin.storage") expected_volume = np.array([1.0, 1.0, 1.0, 1.0]) assert_array_almost_equal(actual_volume, expected_volume) @@ -112,7 +112,7 @@ def test_err_unknown_var(libribasim, basic, tmp_path): config_file = str(tmp_path / "ribasim.toml") libribasim.initialize(config_file) - variable_name = "var-that-does-not-exist" + variable_name = "unknown_node.unknown_variable" error_message = re.escape( f"BMI exception in get_var_type (for variable {variable_name}):" " Message from Ribasim "