From 426e61df0785eed1bacb4177437ca999aa059a1b Mon Sep 17 00:00:00 2001 From: Martijn Visser Date: Mon, 30 Sep 2024 14:52:07 +0200 Subject: [PATCH 1/3] Remove support for resetting cumulatives from BMI --- core/src/bmi.jl | 25 ++++++++++++------------ core/src/callback.jl | 38 +++++++++++++----------------------- core/src/graph.jl | 7 +++---- core/src/parameter.jl | 11 ++++------- core/src/read.jl | 13 ++---------- core/src/solve.jl | 19 +++++++++--------- core/test/run_models_test.jl | 10 +++++----- core/test/time_test.jl | 2 +- 8 files changed, 51 insertions(+), 74 deletions(-) diff --git a/core/src/bmi.jl b/core/src/bmi.jl index 0c12c251e..6a583e8cf 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.infiltration_integrated" # TODO rename to basin.cumulative_infiltration + u.infiltration + elseif name == "basin.drainage_integrated" # TODO rename to 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.realized" # TODO rename to 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..d9b1815d1 100644 --- a/core/src/callback.jl +++ b/core/src/callback.jl @@ -91,30 +91,22 @@ Update with the latest timestep: """ 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 +156,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 +564,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 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..6bb57084e 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)] @@ -101,10 +101,9 @@ function set_current_basin_properties!( 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/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 From 4af27f2dfeb5173b0e13d6cc15524010d46b84ea Mon Sep 17 00:00:00 2001 From: Martijn Visser Date: Mon, 30 Sep 2024 15:29:11 +0200 Subject: [PATCH 2/3] Rename variables in get_value_ptr --- core/src/bmi.jl | 6 +++--- core/src/callback.jl | 13 +++++-------- core/src/solve.jl | 2 +- core/test/bmi_test.jl | 18 +++++++++--------- 4 files changed, 18 insertions(+), 21 deletions(-) diff --git a/core/src/bmi.jl b/core/src/bmi.jl index 6a583e8cf..e81d0ed89 100644 --- a/core/src/bmi.jl +++ b/core/src/bmi.jl @@ -40,15 +40,15 @@ function BMI.get_value_ptr(model::Model, name::AbstractString)::AbstractVector{F p.basin.vertical_flux.infiltration elseif name == "basin.drainage" p.basin.vertical_flux.drainage - elseif name == "basin.infiltration_integrated" # TODO rename to basin.cumulative_infiltration + elseif name == "basin.cumulative_infiltration" u.infiltration - elseif name == "basin.drainage_integrated" # TODO rename to basin.cumulative_drainage + elseif name == "basin.cumulative_drainage" p.basin.cumulative_drainage elseif name == "basin.subgrid_level" p.subgrid.level elseif name == "user_demand.demand" vec(p.user_demand.demand) - elseif name == "user_demand.realized" # TODO rename to user_demand.inflow + elseif name == "user_demand.inflow" u.user_demand_inflow else error("Unknown variable $name") diff --git a/core/src/callback.jl b/core/src/callback.jl index d9b1815d1..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,7 +84,6 @@ 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 @@ -602,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/solve.jl b/core/src/solve.jl index 6bb57084e..2b0ca0b4e 100644 --- a/core/src/solve.jl +++ b/core/src/solve.jl @@ -95,7 +95,7 @@ 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] 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 From 4dc25263a021bec1aec2d863d624dbcb274e7a12 Mon Sep 17 00:00:00 2001 From: Martijn Visser Date: Mon, 30 Sep 2024 18:37:45 +0200 Subject: [PATCH 3/3] Update Python BMI tests --- python/ribasim_api/tests/test_bmi.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) 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):