From f7eb8952a0e6fbed7f2c6913a73c759f9abdbd63 Mon Sep 17 00:00:00 2001 From: Martijn Visser Date: Tue, 1 Oct 2024 11:08:10 +0200 Subject: [PATCH] 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 | 25 +++++++------- core/src/callback.jl | 51 +++++++++++----------------- core/src/graph.jl | 7 ++-- core/src/parameter.jl | 11 +++--- core/src/read.jl | 13 ++----- core/src/solve.jl | 21 ++++++------ core/test/bmi_test.jl | 18 +++++----- core/test/run_models_test.jl | 10 +++--- core/test/time_test.jl | 2 +- python/ribasim_api/tests/test_bmi.py | 16 ++++----- 10 files changed, 74 insertions(+), 100 deletions(-) diff --git a/core/src/bmi.jl b/core/src/bmi.jl index 0c12c251e..e81d0ed89 100644 --- a/core/src/bmi.jl +++ b/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 diff --git a/core/src/callback.jl b/core/src/callback.jl index 8990d86f8..0cdd0d23b 100644 --- a/core/src/callback.jl +++ b/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 diff --git a/core/src/graph.jl b/core/src/graph.jl index d948511ec..d6c0669ba 100644 --- a/core/src/graph.jl +++ b/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 diff --git a/core/src/parameter.jl b/core/src/parameter.jl index 6525363f3..576c12769 100644 --- a/core/src/parameter.jl +++ b/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} diff --git a/core/src/read.jl b/core/src/read.jl index f57689956..7983a3e5e 100644 --- a/core/src/read.jl +++ b/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, diff --git a/core/src/solve.jl b/core/src/solve.jl index 398adf5c6..2b0ca0b4e 100644 --- a/core/src/solve.jl +++ b/core/src/solve.jl @@ -87,7 +87,7 @@ 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)] @@ -95,16 +95,15 @@ function set_current_basin_properties!( 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] diff --git a/core/test/bmi_test.jl b/core/test/bmi_test.jl index 6e0a51800..a9d714874 100644 --- a/core/test/bmi_test.jl +++ b/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 diff --git a/core/test/run_models_test.jl b/core/test/run_models_test.jl index 8e60ef289..065df785b 100644 --- a/core/test/run_models_test.jl +++ b/core/test/run_models_test.jl @@ -131,8 +131,8 @@ end @test model isa Ribasim.Model (; basin) = model.integrator.p @test basin.current_storage[Float64[]] ≈ [1000] - @test basin.vertical_flux_from_input.precipitation == [0.0] - @test basin.vertical_flux_from_input.drainage == [0.0] + @test basin.vertical_flux.precipitation == [0.0] + @test basin.vertical_flux.drainage == [0.0] du = get_du(model.integrator) @test du.evaporation == [0.0] @test du.infiltration == [0.0] @@ -154,9 +154,9 @@ end (; u, p, t) = integrator Ribasim.water_balance!(du, u, p, t) stor = integrator.p.basin.current_storage[parent(du)] - prec = p.basin.vertical_flux_from_input.precipitation + prec = p.basin.vertical_flux.precipitation evap = du.evaporation - drng = p.basin.vertical_flux_from_input.drainage + drng = p.basin.vertical_flux.drainage infl = du.infiltration # The dynamic data has missings, but these are not set. @test prec == [0.0] @@ -239,7 +239,7 @@ end @test successful_retcode(model) @test allunique(Ribasim.tsaves(model)) du = get_du(model.integrator) - precipitation = model.integrator.p.basin.vertical_flux_from_input.precipitation + precipitation = model.integrator.p.basin.vertical_flux.precipitation @test length(precipitation) == 4 @test model.integrator.p.basin.current_storage[parent(du)] ≈ Float32[697.30591, 697.2799, 419.19034, 1334.3859] atol = 2.0 skip = Sys.isapple() diff --git a/core/test/time_test.jl b/core/test/time_test.jl index 14b535117..68886c208 100644 --- a/core/test/time_test.jl +++ b/core/test/time_test.jl @@ -98,7 +98,7 @@ end model = Ribasim.Model(config) (; basin) = model.integrator.p starting_precipitation = - basin.vertical_flux_from_input.precipitation[1] * Ribasim.basin_areas(basin, 1)[end] + basin.vertical_flux.precipitation[1] * Ribasim.basin_areas(basin, 1)[end] BMI.update_until(model, saveat) mean_precipitation = only(model.saved.flow.saveval).precipitation[1] # Given that precipitation stops after 15 of the 20 days diff --git a/python/ribasim_api/tests/test_bmi.py b/python/ribasim_api/tests/test_bmi.py index bdd43b236..bfbd8e3e8 100644 --- a/python/ribasim_api/tests/test_bmi.py +++ b/python/ribasim_api/tests/test_bmi.py @@ -143,13 +143,13 @@ def test_get_value_ptr_basin_forcing(libribasim, leaky_bucket, tmp_path): # Run model for a while to build up integrated forcing libribasim.update_until(60.0) - # Infiltration (integrated) - actual_infiltration = libribasim.get_value_ptr("basin.infiltration_integrated") + # Cumulative infiltration + actual_infiltration = libribasim.get_value_ptr("basin.cumulative_infiltration") expected_infiltration = np.array([0.0]) assert_array_almost_equal(actual_infiltration, expected_infiltration) - # Drainage (integrated) - actual_drainage = libribasim.get_value_ptr("basin.drainage_integrated") + # Cumulative drainage + actual_drainage = libribasim.get_value_ptr("basin.cumulative_drainage") expected_drainage = np.array([0.003 * 60.0]) assert_array_almost_equal(actual_drainage, expected_drainage) @@ -167,10 +167,10 @@ def test_get_value_ptr_allocation(libribasim, user_demand, tmp_path): # Run model for a while to build up realized libribasim.update_until(60.0) - # Realized - actual_realized = libribasim.get_value_ptr("user_demand.realized") - expected_realized = np.array([1e-4 * 60.0, 0.0, 0.0]) - assert_array_almost_equal(actual_realized, expected_realized) + # Realized inflow + actual_inflow = libribasim.get_value_ptr("user_demand.inflow") + expected_inflow = np.array([1e-4 * 60.0, 0.0, 0.0]) + assert_array_almost_equal(actual_inflow, expected_inflow) def test_err_unknown_var(libribasim, basic, tmp_path):