Remove support for resetting cumulatives from BMI

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