Remove support for resetting cumulatives from BMI (#1852)

The coupler is switching to an approach where they don't need to reset
the cumulatives to 0. Since this complicated the Ribasim code, let's get
rid of it.

The reason for the coupler to move away from resetting cumulatives is
that this possibly broke with #1819.

The dedicated BMI arrays like `vertical_flux_bmi` and
`user_demand.realized_bmi` are removed. The `BMI.get_value_ptr` now just
points to the cumulative state, or if it is not a state (drainage,
precipitation) to `basin.cumulative_drainage` or
`basin.cumulative_precipitation`.

This is breaking because I also renamed these BMI variables to be more
consistent with the current code, where we use the term cumulative
internally.

```
user_demand.realized -> user_demand.inflow
basin.drainage_integrated -> basin.cumulative_drainage
basin.infiltration_integrated ->  basin.cumulative_infiltration
```

core/src/bmi.jl

@@ -31,24 +31,25 @@ function BMI.update_until(model::Model, time::Float64)::Model
 end
 
 function BMI.get_value_ptr(model::Model, name::AbstractString)::AbstractVector{Float64}
+    (; u, p) = model.integrator
     if name == "basin.storage"
-        model.integrator.p.basin.current_storage[parent(model.integrator.u)]
+        p.basin.current_storage[parent(u)]
     elseif name == "basin.level"
-        model.integrator.p.basin.current_level[parent(model.integrator.u)]
+        p.basin.current_level[parent(u)]
     elseif name == "basin.infiltration"
-        model.integrator.p.basin.vertical_flux_from_input.infiltration
+        p.basin.vertical_flux.infiltration
     elseif name == "basin.drainage"
-        model.integrator.p.basin.vertical_flux_from_input.drainage
-    elseif name == "basin.infiltration_integrated"
-        model.integrator.p.basin.vertical_flux_bmi.infiltration
-    elseif name == "basin.drainage_integrated"
-        model.integrator.p.basin.vertical_flux_bmi.drainage
+        p.basin.vertical_flux.drainage
+    elseif name == "basin.cumulative_infiltration"
+        u.infiltration
+    elseif name == "basin.cumulative_drainage"
+        p.basin.cumulative_drainage
     elseif name == "basin.subgrid_level"
-        model.integrator.p.subgrid.level
+        p.subgrid.level
     elseif name == "user_demand.demand"
-        vec(model.integrator.p.user_demand.demand)
-    elseif name == "user_demand.realized"
-        model.integrator.p.user_demand.realized_bmi
+        vec(p.user_demand.demand)
+    elseif name == "user_demand.inflow"
+        u.user_demand_inflow
     else
         error("Unknown variable $name")
     end

core/src/callback.jl

@@ -22,10 +22,10 @@ function create_callbacks(
     save_basin_state_cb = SavingCallback(save_basin_state, saved_basin_states; saveat)
     push!(callbacks, save_basin_state_cb)
 
-    # Update cumulative flows (exact integration and for BMI)
-    integrated_flows_cb =
+    # Update cumulative flows (exact integration and for allocation)
+    cumulative_flows_cb =
         FunctionCallingCallback(update_cumulative_flows!; func_start = false)
-    push!(callbacks, integrated_flows_cb)
+    push!(callbacks, cumulative_flows_cb)
 
     # Update Basin forcings
     tstops = get_tstops(basin.time.time, starttime)
@@ -84,37 +84,28 @@ end
 
 """
 Update with the latest timestep:
-- Cumulative flows/forcings which are exposed via the BMI
 - Cumulative flows/forcings which are integrated exactly
 - Cumulative flows/forcings which are input for the allocation algorithm
 - Cumulative flows/forcings which are realized demands in the allocation context
 
 """
 function update_cumulative_flows!(u, t, integrator)::Nothing
-    (; p, uprev, tprev, dt) = integrator
-    (; basin, user_demand, flow_boundary, allocation) = p
-    (; vertical_flux_bmi, vertical_flux_from_input) = basin
+    (; p, tprev, dt) = integrator
+    (; basin, flow_boundary, allocation) = p
+    (; vertical_flux) = basin
 
     # Update tprev
     p.tprev[] = t
 
-    # Update cumulative flows exposed via the BMI
-    @. user_demand.realized_bmi += u.user_demand_inflow - uprev.user_demand_inflow
-    @. vertical_flux_bmi.drainage += vertical_flux_from_input.drainage * dt
-    @. vertical_flux_bmi.evaporation += u.evaporation - uprev.evaporation
-    @. vertical_flux_bmi.infiltration += u.infiltration - uprev.infiltration
-
     # Update cumulative forcings which are integrated exactly
-    @. basin.cumulative_drainage += vertical_flux_from_input.drainage * dt
-    @. basin.cumulative_drainage_saveat += vertical_flux_from_input.drainage * dt
+    @. basin.cumulative_drainage += vertical_flux.drainage * dt
+    @. basin.cumulative_drainage_saveat += vertical_flux.drainage * dt
 
     # Precipitation depends on fixed area
     for node_id in basin.node_id
         fixed_area = basin_areas(basin, node_id.idx)[end]
-        added_precipitation =
-            fixed_area * vertical_flux_from_input.precipitation[node_id.idx] * dt
+        added_precipitation = fixed_area * vertical_flux.precipitation[node_id.idx] * dt
 
-        vertical_flux_bmi.precipitation[node_id.idx] += added_precipitation
         basin.cumulative_precipitation[node_id.idx] += added_precipitation
         basin.cumulative_precipitation_saveat[node_id.idx] += added_precipitation
     end
@@ -164,16 +155,14 @@ function flow_update_on_edge(
 )::Float64
     (; u, uprev, p, t, tprev, dt) = integrator
     (; basin, flow_boundary) = p
-    (; vertical_flux_from_input) = basin
+    (; vertical_flux) = basin
     from_id, to_id = edge_src
     if from_id == to_id
         @assert from_id.type == to_id.type == NodeType.Basin
         idx = from_id.idx
         fixed_area = basin_areas(basin, idx)[end]
-        (
-            fixed_area * vertical_flux_from_input.precipitation[idx] +
-            vertical_flux_from_input.drainage[idx]
-        ) * dt - (u.evaporation[idx] - uprev.evaporation[idx]) -
+        (fixed_area * vertical_flux.precipitation[idx] + vertical_flux.drainage[idx]) * dt -
+        (u.evaporation[idx] - uprev.evaporation[idx]) -
         (u.infiltration[idx] - uprev.infiltration[idx])
     elseif from_id.type == NodeType.FlowBoundary
         if flow_boundary.active[from_id.idx]
@@ -574,17 +563,17 @@ end
 function update_basin!(integrator)::Nothing
     (; p) = integrator
     (; basin) = p
-    (; node_id, time, vertical_flux_from_input) = basin
+    (; node_id, time, vertical_flux) = basin
     t = datetime_since(integrator.t, integrator.p.starttime)
 
     rows = searchsorted(time.time, t)
     timeblock = view(time, rows)
 
     table = (;
-        vertical_flux_from_input.precipitation,
-        vertical_flux_from_input.potential_evaporation,
-        vertical_flux_from_input.drainage,
-        vertical_flux_from_input.infiltration,
+        vertical_flux.precipitation,
+        vertical_flux.potential_evaporation,
+        vertical_flux.drainage,
+        vertical_flux.infiltration,
     )
 
     for row in timeblock
@@ -612,15 +601,13 @@ function update_allocation!(integrator)::Nothing
         return nothing
     end
 
-    # Divide by the allocation Δt to obtain the mean input flows
-    # from the integrated flows
+    # Divide by the allocation Δt to get the mean input flows from the cumulative flows
     (; Δt_allocation) = allocation_models[1]
     for edge in keys(mean_input_flows)
         mean_input_flows[edge] /= Δt_allocation
     end
 
-    # Divide by the allocation Δt to obtain the mean realized flows
-    # from the integrated flows
+    # Divide by the allocation Δt to get the mean realized flows from the cumulative flows
     for edge in keys(mean_realized_flows)
         mean_realized_flows[edge] /= Δt_allocation
     end

core/src/graph.jl

@@ -271,9 +271,8 @@ end
 function get_influx(du::ComponentVector, id::NodeID, p::Parameters)
     @assert id.type == NodeType.Basin
     (; basin) = p
-    (; vertical_flux_from_input) = basin
+    (; vertical_flux) = basin
     fixed_area = basin_areas(basin, id.idx)[end]
-    return fixed_area * vertical_flux_from_input.precipitation[id.idx] +
-           vertical_flux_from_input.drainage[id.idx] - du.evaporation[id.idx] -
-           du.infiltration[id.idx]
+    return fixed_area * vertical_flux.precipitation[id.idx] +
+           vertical_flux.drainage[id.idx] - du.evaporation[id.idx] - du.infiltration[id.idx]
 end

core/src/parameter.jl

@@ -305,13 +305,12 @@ else
     T = Vector{Float64}
 end
 """
-@kwdef struct Basin{C, V1, V2} <: AbstractParameterNode
+@kwdef struct Basin{C, V} <: AbstractParameterNode
     node_id::Vector{NodeID}
     inflow_ids::Vector{Vector{NodeID}} = [NodeID[]]
     outflow_ids::Vector{Vector{NodeID}} = [NodeID[]]
     # Vertical fluxes
-    vertical_flux_from_input::V1 = zeros(length(node_id))
-    vertical_flux_bmi::V2 = zeros(length(node_id))
+    vertical_flux::V = zeros(length(node_id))
     # Initial_storage
     storage0::Vector{Float64} = zeros(length(node_id))
     storage_prev_saveat::Vector{Float64} = zeros(length(node_id))
@@ -731,7 +730,6 @@ inflow_edge: incoming flow edge
 outflow_edge: outgoing flow edge metadata
     The ID of the source node is always the ID of the UserDemand node
 active: whether this node is active and thus demands water
-realized_bmi: Cumulative inflow volume, for read or reset by BMI only
 demand: water flux demand of UserDemand per priority (node_idx, priority_idx)
     Each UserDemand has a demand for all priorities,
     which is 0.0 if it is not provided explicitly.
@@ -748,7 +746,6 @@ min_level: The level of the source basin below which the UserDemand does not abs
     inflow_edge::Vector{EdgeMetadata} = []
     outflow_edge::Vector{EdgeMetadata} = []
     active::Vector{Bool} = fill(true, length(node_id))
-    realized_bmi::Vector{Float64} = zeros(length(node_id))
     demand::Matrix{Float64}
     demand_reduced::Matrix{Float64}
     demand_itp::Vector{Vector{ScalarInterpolation}}
@@ -818,11 +815,11 @@ const ModelGraph = MetaGraph{
     Float64,
 }
 
-@kwdef struct Parameters{C1, C2, C3, C4, C5, V1, V2}
+@kwdef struct Parameters{C1, C2, C3, C4, C5, V}
     starttime::DateTime
     graph::ModelGraph
     allocation::Allocation
-    basin::Basin{C1, V1, V2}
+    basin::Basin{C1, V}
     linear_resistance::LinearResistance
     manning_resistance::ManningResistance
     tabulated_rating_curve::TabulatedRatingCurve{C2}

core/src/read.jl

@@ -577,14 +577,8 @@ function Basin(db::DB, config::Config, graph::MetaGraph)::Basin
     set_current_value!(table, node_id, time, config.starttime)
     check_no_nans(table, "Basin")
 
-    vertical_flux_from_input =
+    vertical_flux =
         ComponentVector(; precipitation, potential_evaporation, drainage, infiltration)
-    vertical_flux_bmi = ComponentVector(;
-        precipitation = zero(precipitation),
-        evaporation = zero(evaporation),
-        drainage = zero(drainage),
-        infiltration = zero(infiltration),
-    )
 
     demand = zeros(length(node_id))
 
@@ -642,8 +636,7 @@ function Basin(db::DB, config::Config, graph::MetaGraph)::Basin
         node_id,
         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_from_input,
-        vertical_flux_bmi,
+        vertical_flux,
         current_level,
         current_area,
         storage_to_level,
@@ -1035,7 +1028,6 @@ function UserDemand(db::DB, config::Config, graph::MetaGraph)::UserDemand
     n_user = length(node_ids)
     n_priority = length(priorities)
     active = fill(true, n_user)
-    realized_bmi = zeros(n_user)
     demand = zeros(n_user, n_priority)
     demand_reduced = zeros(n_user, n_priority)
     trivial_timespan = [0.0, prevfloat(Inf)]
@@ -1088,7 +1080,6 @@ function UserDemand(db::DB, config::Config, graph::MetaGraph)::UserDemand
         inflow_edge = inflow_edge.(Ref(graph), node_ids),
         outflow_edge = outflow_edge.(Ref(graph), node_ids),
         active,
-        realized_bmi,
         demand,
         demand_reduced,
         demand_itp,

core/src/solve.jl

@@ -87,24 +87,23 @@ function set_current_basin_properties!(
         current_cumulative_drainage,
         cumulative_precipitation,
         cumulative_drainage,
-        vertical_flux_from_input,
+        vertical_flux,
     ) = basin
     current_storage = current_storage[parent(du)]
     current_level = current_level[parent(du)]
     current_area = current_area[parent(du)]
     current_cumulative_precipitation = current_cumulative_precipitation[parent(du)]
     current_cumulative_drainage = current_cumulative_drainage[parent(du)]
 
-    # The exactly integrated precipitation and drainage up to the t of this water_balance call
+    # The exact cumulative precipitation and drainage up to the t of this water_balance call
     dt = t - p.tprev[]
     for node_id in basin.node_id
         fixed_area = basin_areas(basin, node_id.idx)[end]
         current_cumulative_precipitation[node_id.idx] =
             cumulative_precipitation[node_id.idx] +
-            fixed_area * vertical_flux_from_input.precipitation[node_id.idx] * dt
+            fixed_area * vertical_flux.precipitation[node_id.idx] * dt
     end
-    @. current_cumulative_drainage =
-        cumulative_drainage + dt * vertical_flux_from_input.drainage
+    @. current_cumulative_drainage = cumulative_drainage + dt * vertical_flux.drainage
 
     formulate_storages!(current_storage, du, u, p, t)
 
@@ -232,7 +231,7 @@ 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_from_input, current_level, current_area) = basin
+    (; vertical_flux, current_level, current_area) = basin
     current_level = current_level[parent(du)]
     current_area = current_area[parent(du)]
 
@@ -244,8 +243,8 @@ function update_vertical_flux!(basin::Basin, du::AbstractVector)::Nothing
         depth = max(level - bottom, 0.0)
         factor = reduction_factor(depth, 0.1)
 
-        evaporation = area * factor * vertical_flux_from_input.potential_evaporation[id.idx]
-        infiltration = factor * vertical_flux_from_input.infiltration[id.idx]
+        evaporation = area * factor * vertical_flux.potential_evaporation[id.idx]
+        infiltration = factor * vertical_flux.infiltration[id.idx]
 
         du.evaporation[id.idx] = evaporation
         du.infiltration[id.idx] = infiltration
@@ -338,7 +337,7 @@ Formulate the time derivative of the storage in a single Basin.
 """
 function formulate_dstorage(du::ComponentVector, p::Parameters, t::Number, node_id::NodeID)
     (; basin) = p
-    (; inflow_ids, outflow_ids, vertical_flux_from_input) = basin
+    (; inflow_ids, outflow_ids, vertical_flux) = basin
     @assert node_id.type == NodeType.Basin
     dstorage = 0.0
     for inflow_id in inflow_ids[node_id.idx]
@@ -349,8 +348,8 @@ function formulate_dstorage(du::ComponentVector, p::Parameters, t::Number, node_
     end
 
     fixed_area = basin_areas(basin, node_id.idx)[end]
-    dstorage += fixed_area * vertical_flux_from_input.precipitation[node_id.idx]
-    dstorage += vertical_flux_from_input.drainage[node_id.idx]
+    dstorage += fixed_area * vertical_flux.precipitation[node_id.idx]
+    dstorage += vertical_flux.drainage[node_id.idx]
     dstorage -= du.evaporation[node_id.idx]
     dstorage -= du.infiltration[node_id.idx]
 

core/test/bmi_test.jl

@@ -63,11 +63,11 @@ end
         "basin.level",
         "basin.infiltration",
         "basin.drainage",
-        "basin.infiltration_integrated",
-        "basin.drainage_integrated",
+        "basin.cumulative_infiltration",
+        "basin.cumulative_drainage",
         "basin.subgrid_level",
         "user_demand.demand",
-        "user_demand.realized",
+        "user_demand.inflow",
     ]
         value_first = BMI.get_value_ptr(model, name)
         BMI.update_until(model, 86400.0)
@@ -77,7 +77,7 @@ end
     end
 end
 
-@testitem "realized_user_demand" begin
+@testitem "UserDemand inflow" begin
     import BasicModelInterface as BMI
 
     toml_path =
@@ -86,19 +86,19 @@ end
     config = Ribasim.Config(toml_path; allocation_use_allocation = false)
     model = Ribasim.Model(config)
     demand = BMI.get_value_ptr(model, "user_demand.demand")
-    realized = BMI.get_value_ptr(model, "user_demand.realized")
+    inflow = BMI.get_value_ptr(model, "user_demand.inflow")
     # One year in seconds
     year = model.integrator.p.user_demand.demand_itp[2][1].t[2]
     demand_start = 1e-3
     slope = 1e-3 / year
     day = 86400.0
     BMI.update_until(model, day)
-    @test realized ≈ [demand_start * day, demand_start * day + 0.5 * slope * day^2] atol =
+    @test inflow ≈ [demand_start * day, demand_start * day + 0.5 * slope * day^2] atol =
         1e-3
     demand_later = 2e-3
     demand[1] = demand_later
     BMI.update_until(model, 2day)
-    @test realized[1] ≈ demand_start * day + demand_later * day atol = 1e-3
+    @test inflow[1] ≈ demand_start * day + demand_later * day atol = 1e-3
 end
 
 @testitem "vertical basin flux" begin
@@ -114,6 +114,6 @@ end
     Δt = 5 * 86400.0
     BMI.update_until(model, Δt)
 
-    drainage_integrated = BMI.get_value_ptr(model, "basin.drainage_integrated")
-    @test drainage_integrated ≈ Δt * drainage_flux
+    cumulative_drainage = BMI.get_value_ptr(model, "basin.cumulative_drainage")
+    @test cumulative_drainage ≈ Δt * drainage_flux
 end