add realized user demand and expose via BMI (#1190)

Fixes #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 <[email protected]>

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

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

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,

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,

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
 

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,11 +58,36 @@ 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)
         # get_value_ptr does not copy
         @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

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

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],

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,15 +72,15 @@ 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"
 
 
 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 "