From db46d1e45aca519fa772ba8438579a968ab60744 Mon Sep 17 00:00:00 2001 From: Bart de Koning <74617371+SouthEndMusic@users.noreply.github.com> Date: Thu, 28 Mar 2024 16:18:26 +0100 Subject: [PATCH] Split vertical fluxes on basin (#1300) Fixes https://github.com/Deltares/Ribasim/issues/661. --------- Co-authored-by: Martijn Visser --- core/src/allocation_optim.jl | 6 +- core/src/bmi.jl | 8 ++- core/src/callback.jl | 126 +++++++++++++++++------------------ core/src/graph.jl | 42 ------------ core/src/model.jl | 3 +- core/src/parameter.jl | 42 ++++++------ core/src/read.jl | 41 ++++++++---- core/src/solve.jl | 80 +++++++++------------- core/src/util.jl | 38 +++++++++++ core/src/write.jl | 56 +++++++++++----- core/test/bmi_test.jl | 2 + core/test/run_models_test.jl | 53 ++++++++++----- core/test/time_test.jl | 75 ++++++++++++++++++--- core/test/utils_test.jl | 2 + docs/core/usage.qmd | 32 +++++---- 15 files changed, 360 insertions(+), 246 deletions(-) diff --git a/core/src/allocation_optim.jl b/core/src/allocation_optim.jl index 17e9eab37..6f3606d4b 100644 --- a/core/src/allocation_optim.jl +++ b/core/src/allocation_optim.jl @@ -321,9 +321,13 @@ function get_basin_data( node_id::NodeID, ) (; graph, basin, level_demand) = p + (; vertical_flux) = basin (; Δt_allocation) = allocation_model @assert node_id.type == NodeType.Basin - influx = get_flow(graph, node_id, 0.0) + vertical_flux = get_tmp(vertical_flux, 0) + _, basin_idx = id_index(basin.node_id, node_id) + # NOTE: Instantaneous + influx = get_influx(basin, node_id) _, basin_idx = id_index(basin.node_id, node_id) storage_basin = u.storage[basin_idx] control_inneighbors = inneighbor_labels_type(graph, node_id, EdgeType.control) diff --git a/core/src/bmi.jl b/core/src/bmi.jl index 39a77afb7..7414d3dfe 100644 --- a/core/src/bmi.jl +++ b/core/src/bmi.jl @@ -37,9 +37,13 @@ function BMI.get_value_ptr(model::Model, name::AbstractString) elseif name == "basin.level" get_tmp(model.integrator.p.basin.current_level, 0) elseif name == "basin.infiltration" - model.integrator.p.basin.infiltration + get_tmp(model.integrator.p.basin.vertical_flux, 0).infiltration elseif name == "basin.drainage" - model.integrator.p.basin.drainage + get_tmp(model.integrator.p.basin.vertical_flux, 0).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 elseif name == "basin.subgrid_level" model.integrator.p.subgrid.level elseif name == "user_demand.demand" diff --git a/core/src/callback.jl b/core/src/callback.jl index 9f3db5c1a..8aa34f519 100644 --- a/core/src/callback.jl +++ b/core/src/callback.jl @@ -36,13 +36,13 @@ function create_callbacks( (; starttime, basin, tabulated_rating_curve, discrete_control) = parameters callbacks = SciMLBase.DECallback[] + integrating_flows_cb = FunctionCallingCallback(integrate_flows!; func_start = false) + push!(callbacks, integrating_flows_cb) + tstops = get_tstops(basin.time.time, starttime) basin_cb = PresetTimeCallback(tstops, update_basin; save_positions = (false, false)) push!(callbacks, basin_cb) - integrating_flows_cb = FunctionCallingCallback(integrate_flows!; func_start = false) - push!(callbacks, integrating_flows_cb) - tstops = get_tstops(tabulated_rating_curve.time.time, starttime) tabulated_rating_curve_cb = PresetTimeCallback( tstops, @@ -61,16 +61,17 @@ function create_callbacks( push!(callbacks, allocation_cb) end + # If saveat is a vector which contains 0.0 this callback will still be called + # at t = 0.0 despite save_start = false + saveat = saveat isa Vector ? filter(x -> x != 0.0, saveat) : saveat + saved_vertical_flux = SavedValues(Float64, typeof(basin.vertical_flux_integrated)) + save_vertical_flux_cb = + SavingCallback(save_vertical_flux, saved_vertical_flux; saveat, save_start = false) + push!(callbacks, save_vertical_flux_cb) + # save the flows over time, as a Vector of the nonzeros(flow) saved_flow = SavedValues(Float64, Vector{Float64}) - save_flow_cb = SavingCallback( - save_flow, - saved_flow; - # If saveat is a vector which contains 0.0 this callback will still be called - # at t = 0.0 despite save_start = false - saveat = saveat isa Vector ? filter(x -> x != 0.0, saveat) : saveat, - save_start = false, - ) + save_flow_cb = SavingCallback(save_flow, saved_flow; saveat, save_start = false) push!(callbacks, save_flow_cb) # interpolate the levels @@ -85,7 +86,7 @@ function create_callbacks( push!(callbacks, export_cb) end - saved = SavedResults(saved_flow, saved_subgrid_level) + saved = SavedResults(saved_flow, saved_vertical_flux, saved_subgrid_level) n_conditions = length(discrete_control.node_id) if n_conditions > 0 @@ -108,27 +109,20 @@ Integrate flows over the last timestep """ function integrate_flows!(u, t, integrator)::Nothing (; p, dt) = integrator - (; graph, user_demand) = p - (; - flow, - flow_dict, - flow_vertical, - flow_prev, - flow_vertical_prev, - flow_integrated, - flow_vertical_integrated, - ) = graph[] + (; graph, user_demand, basin) = p + (; flow, flow_dict, flow_prev, flow_integrated) = graph[] + (; vertical_flux, vertical_flux_prev, vertical_flux_integrated, vertical_flux_bmi) = + basin flow = get_tmp(flow, 0) - flow_vertical = get_tmp(flow_vertical, 0) - + vertical_flux = get_tmp(vertical_flux, 0) if !isempty(flow_prev) && isnan(flow_prev[1]) # If flow_prev is not populated yet copyto!(flow_prev, flow) - copyto!(flow_vertical_prev, flow_vertical) end @. flow_integrated += 0.5 * (flow + flow_prev) * dt - @. flow_vertical_integrated += 0.5 * (flow_vertical + flow_vertical_prev) * dt + @. vertical_flux_integrated += 0.5 * (vertical_flux + vertical_flux_prev) * dt + @. vertical_flux_bmi += 0.5 * (vertical_flux + vertical_flux_prev) * dt for (i, id) in enumerate(user_demand.node_id) src_id = inflow_id(graph, id) @@ -137,10 +131,37 @@ function integrate_flows!(u, t, integrator)::Nothing end copyto!(flow_prev, flow) - copyto!(flow_vertical_prev, flow_vertical) + copyto!(vertical_flux_prev, vertical_flux) return nothing end +"Compute the average flows over the last saveat interval and write +them to SavedValues" +function save_flow(u, t, integrator) + (; flow_integrated) = integrator.p.graph[] + + Δt = get_Δt(integrator) + flow_mean = copy(flow_integrated) + flow_mean ./= Δt + fill!(flow_integrated, 0.0) + + return flow_mean +end + +"Compute the average vertical fluxes over the last saveat interval and write +them to SavedValues" +function save_vertical_flux(u, t, integrator) + (; basin) = integrator.p + (; vertical_flux_integrated) = basin + + Δt = get_Δt(integrator) + vertical_flux_mean = copy(vertical_flux_integrated) + vertical_flux_mean ./= Δt + fill!(vertical_flux_integrated, 0.0) + + return vertical_flux_mean +end + """ Listens for changes in condition truths. """ @@ -430,36 +451,6 @@ function set_control_params!(p::Parameters, node_id::NodeID, control_state::Stri end end -"Compute the average flows over the last saveat interval and write -them to SavedValues" -function save_flow(u, t, integrator) - (; dt, p) = integrator - (; graph) = p - (; flow_integrated, flow_vertical_integrated, saveat) = graph[] - - Δt = if iszero(saveat) - dt - elseif isinf(saveat) - t - else - t_end = integrator.sol.prob.tspan[2] - if t_end - t > saveat - saveat - else - # The last interval might be shorter than saveat - rem = t % saveat - iszero(rem) ? saveat : rem - end - end - - mean_flow_all = vcat(flow_vertical_integrated, flow_integrated) - mean_flow_all ./= Δt - fill!(flow_vertical_integrated, 0.0) - fill!(flow_integrated, 0.0) - - return mean_flow_all -end - function update_subgrid_level!(integrator)::Nothing basin_level = get_tmp(integrator.p.basin.current_level, 0) subgrid = integrator.p.subgrid @@ -476,18 +467,21 @@ end "Load updates from 'Basin / time' into the parameters" function update_basin(integrator)::Nothing - (; basin) = integrator.p - (; node_id, time) = basin + (; p, u) = integrator + (; basin) = p + (; storage) = u + (; node_id, time, vertical_flux_from_input, vertical_flux, vertical_flux_prev) = basin t = datetime_since(integrator.t, integrator.p.starttime) + vertical_flux = get_tmp(vertical_flux, integrator.u) rows = searchsorted(time.time, t) timeblock = view(time, rows) table = (; - basin.precipitation, - basin.potential_evaporation, - basin.drainage, - basin.infiltration, + vertical_flux_from_input.precipitation, + vertical_flux_from_input.potential_evaporation, + vertical_flux_from_input.drainage, + vertical_flux_from_input.infiltration, ) for row in timeblock @@ -496,6 +490,12 @@ function update_basin(integrator)::Nothing set_table_row!(table, row, i) end + for (i, id) in enumerate(basin.node_id) + update_vertical_flux!(basin, storage, i) + end + + # Forget about vertical fluxes to handle discontinuous forcing from basin_update + copyto!(vertical_flux_prev, vertical_flux) return nothing end diff --git a/core/src/graph.jl b/core/src/graph.jl index 2d19c717d..cb453e88d 100644 --- a/core/src/graph.jl +++ b/core/src/graph.jl @@ -24,10 +24,6 @@ function create_graph(db::DB, config::Config, chunk_sizes::Vector{Int})::MetaGra flow_counter = 0 # Dictionary from flow edge to index in flow vector flow_dict = Dict{Tuple{NodeID, NodeID}, Int}() - # The number of nodes with vertical flow (interaction with outside of model) - flow_vertical_counter = 0 - # Dictionary from node ID to index in vertical flow vector - flow_vertical_dict = Dict{NodeID, Int}() graph = MetaGraph( DiGraph(); label_type = NodeID, @@ -49,10 +45,6 @@ function create_graph(db::DB, config::Config, chunk_sizes::Vector{Int})::MetaGra end graph[node_id] = NodeMetadata(Symbol(snake_case(row.node_type)), allocation_network_id) - if row.node_type in nonconservative_nodetypes - flow_vertical_counter += 1 - flow_vertical_dict[node_id] = flow_vertical_counter - end end errors = false @@ -105,12 +97,8 @@ function create_graph(db::DB, config::Config, chunk_sizes::Vector{Int})::MetaGra flow = zeros(flow_counter) flow_prev = fill(NaN, flow_counter) flow_integrated = zeros(flow_counter) - flow_vertical = zeros(flow_vertical_counter) - flow_vertical_prev = fill(NaN, flow_vertical_counter) - flow_vertical_integrated = zeros(flow_vertical_counter) if config.solver.autodiff flow = DiffCache(flow, chunk_sizes) - flow_vertical = DiffCache(flow_vertical, chunk_sizes) end graph_data = (; node_ids, @@ -120,10 +108,6 @@ function create_graph(db::DB, config::Config, chunk_sizes::Vector{Int})::MetaGra flow, flow_prev, flow_integrated, - flow_vertical_dict, - flow_vertical, - flow_vertical_prev, - flow_vertical_integrated, config.solver.saveat, ) graph = @set graph.graph_data = graph_data @@ -195,15 +179,6 @@ function set_flow!(graph::MetaGraph, id_src::NodeID, id_dst::NodeID, q::Number): return nothing end -""" -Set the given flow q on the horizontal (self-loop) edge from id to id. -""" -function set_flow!(graph::MetaGraph, id::NodeID, q::Number)::Nothing - (; flow_vertical_dict, flow_vertical) = graph[] - get_tmp(flow_vertical, q)[flow_vertical_dict[id]] = q - return nothing -end - """ Add the given flow q to the existing flow over the edge between the given nodes. """ @@ -213,15 +188,6 @@ function add_flow!(graph::MetaGraph, id_src::NodeID, id_dst::NodeID, q::Number): return nothing end -""" -Add the given flow q to the flow over the edge on the horizontal (self-loop) edge from id to id. -""" -function add_flow!(graph::MetaGraph, id::NodeID, q::Number)::Nothing - (; flow_vertical_dict, flow_vertical) = graph[] - get_tmp(flow_vertical, q)[flow_vertical_dict[id]] += q - return nothing -end - """ Get the flow over the given edge (val is needed for get_tmp from ForwardDiff.jl). """ @@ -230,14 +196,6 @@ function get_flow(graph::MetaGraph, id_src::NodeID, id_dst::NodeID, val)::Number return get_tmp(flow, val)[flow_dict[id_src, id_dst]] end -""" -Get the flow over the given horizontal (selfloop) edge (val is needed for get_tmp from ForwardDiff.jl). -""" -function get_flow(graph::MetaGraph, id::NodeID, val)::Number - (; flow_vertical_dict, flow_vertical) = graph[] - return get_tmp(flow_vertical, val)[flow_vertical_dict[id]] -end - """ Get the inneighbor node IDs of the given node ID (label) over the given edge type in the graph. diff --git a/core/src/model.jl b/core/src/model.jl index afa24347e..111407a34 100644 --- a/core/src/model.jl +++ b/core/src/model.jl @@ -1,5 +1,6 @@ -struct SavedResults +struct SavedResults{V1 <: ComponentVector{Float64}} flow::SavedValues{Float64, Vector{Float64}} + vertical_flux::SavedValues{Float64, V1} subgrid_level::SavedValues{Float64, Vector{Float64}} end diff --git a/core/src/parameter.jl b/core/src/parameter.jl index e1d150ce1..214d056a4 100644 --- a/core/src/parameter.jl +++ b/core/src/parameter.jl @@ -151,12 +151,14 @@ else T = Vector{Float64} end """ -struct Basin{T, C} <: AbstractParameterNode +struct Basin{T, C, V1, V2, V3} <: AbstractParameterNode node_id::Indices{NodeID} - precipitation::Vector{Float64} - potential_evaporation::Vector{Float64} - drainage::Vector{Float64} - infiltration::Vector{Float64} + # Vertical fluxes + vertical_flux_from_input::V1 + vertical_flux::V2 + vertical_flux_prev::V3 + vertical_flux_integrated::V3 + vertical_flux_bmi::V3 # Cache this to avoid recomputation current_level::T current_area::T @@ -171,10 +173,11 @@ struct Basin{T, C} <: AbstractParameterNode function Basin( node_id, - precipitation, - potential_evaporation, - drainage, - infiltration, + vertical_flux_from_input::V1, + vertical_flux::V2, + vertical_flux_prev::V3, + vertical_flux_integrated::V3, + vertical_flux_bmi::V3, current_level::T, current_area::T, area, @@ -182,15 +185,16 @@ struct Basin{T, C} <: AbstractParameterNode storage, demand, time::StructVector{BasinTimeV1, C, Int}, - ) where {T, C} + ) where {T, C, V1, V2, V3} is_valid = valid_profiles(node_id, level, area) is_valid || error("Invalid Basin / profile table.") - return new{T, C}( + return new{T, C, V1, V2, V3}( node_id, - precipitation, - potential_evaporation, - drainage, - infiltration, + vertical_flux_from_input, + vertical_flux, + vertical_flux_prev, + vertical_flux_integrated, + vertical_flux_bmi, current_level, current_area, area, @@ -553,7 +557,7 @@ struct Subgrid end # TODO Automatically add all nodetypes here -struct Parameters{T, C1, C2} +struct Parameters{T, C1, C2, V1, V2, V3} starttime::DateTime graph::MetaGraph{ Int64, @@ -569,17 +573,13 @@ struct Parameters{T, C1, C2} flow::T, flow_prev::Vector{Float64}, flow_integrated::Vector{Float64}, - flow_vertical_dict::Dict{NodeID, Int}, - flow_vertical::T, - flow_vertical_prev::Vector{Float64}, - flow_vertical_integrated::Vector{Float64}, saveat::Float64, }, MetaGraphsNext.var"#11#13", Float64, } allocation::Allocation - basin::Basin{T, C1} + basin::Basin{T, C1, V1, V2, V3} linear_resistance::LinearResistance manning_resistance::ManningResistance tabulated_rating_curve::TabulatedRatingCurve{C2} diff --git a/core/src/read.jl b/core/src/read.jl index ccf7e9704..0ae86147d 100644 --- a/core/src/read.jl +++ b/core/src/read.jl @@ -488,15 +488,11 @@ function Basin(db::DB, config::Config, chunk_sizes::Vector{Int})::Basin current_level = zeros(n) current_area = zeros(n) - if config.solver.autodiff - current_level = DiffCache(current_level, chunk_sizes) - current_area = DiffCache(current_area, chunk_sizes) - end - - precipitation = zeros(length(node_id)) - potential_evaporation = zeros(length(node_id)) - drainage = zeros(length(node_id)) - infiltration = zeros(length(node_id)) + precipitation = zeros(n) + potential_evaporation = zeros(n) + evaporation = zeros(n) + drainage = zeros(n) + infiltration = zeros(n) table = (; precipitation, potential_evaporation, drainage, infiltration) area, level, storage = create_storage_tables(db, config) @@ -509,14 +505,33 @@ function Basin(db::DB, config::Config, chunk_sizes::Vector{Int})::Basin set_current_value!(table, node_id, time, config.starttime) check_no_nans(table, "Basin") + vertical_flux_from_input = + ComponentVector(; precipitation, potential_evaporation, drainage, infiltration) + vertical_flux = ComponentVector(; + precipitation = copy(precipitation), + evaporation, + drainage = copy(drainage), + infiltration = copy(infiltration), + ) + vertical_flux_prev = zero(vertical_flux) + vertical_flux_integrated = zero(vertical_flux) + vertical_flux_bmi = zero(vertical_flux) + + if config.solver.autodiff + current_level = DiffCache(current_level, chunk_sizes) + current_area = DiffCache(current_area, chunk_sizes) + vertical_flux = DiffCache(vertical_flux, chunk_sizes) + end + demand = zeros(length(node_id)) return Basin( Indices(NodeID.(NodeType.Basin, node_id)), - precipitation, - potential_evaporation, - drainage, - infiltration, + vertical_flux_from_input, + vertical_flux, + vertical_flux_prev, + vertical_flux_integrated, + vertical_flux_bmi, current_level, current_area, area, diff --git a/core/src/solve.jl b/core/src/solve.jl index df9f09d4d..f581fc591 100644 --- a/core/src/solve.jl +++ b/core/src/solve.jl @@ -14,13 +14,12 @@ function water_balance!( du .= 0.0 get_tmp(graph[].flow, storage) .= 0.0 - get_tmp(graph[].flow_vertical, storage) .= 0.0 # Ensures current_* vectors are current set_current_basin_properties!(basin, storage) # Basin forcings - formulate_basins!(du, basin, graph, storage) + formulate_basins!(du, basin, storage) # First formulate intermediate flows formulate_flows!(p, storage, t) @@ -52,34 +51,42 @@ end 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, storage::AbstractVector, i::Int)::Nothing + (; current_level, current_area, vertical_flux_from_input, vertical_flux) = basin + current_level = get_tmp(current_level, storage) + current_area = get_tmp(current_area, storage) + vertical_flux = get_tmp(vertical_flux, storage) + + level = current_level[i] + area = current_area[i] + + bottom = basin.level[i][1] + fixed_area = basin.area[i][end] + depth = max(level - bottom, 0.0) + factor = reduction_factor(depth, 0.1) + + precipitation = fixed_area * vertical_flux_from_input.precipitation[i] + evaporation = area * factor * vertical_flux_from_input.potential_evaporation[i] + drainage = vertical_flux_from_input.drainage[i] + infiltration = factor * vertical_flux_from_input.infiltration[i] + + vertical_flux.precipitation[i] = precipitation + vertical_flux.evaporation[i] = evaporation + vertical_flux.drainage[i] = drainage + vertical_flux.infiltration[i] = infiltration + + return nothing +end + function formulate_basins!( du::AbstractVector, basin::Basin, - graph::MetaGraph, storage::AbstractVector, )::Nothing - (; node_id, current_level, current_area) = basin - current_level = get_tmp(current_level, storage) - current_area = get_tmp(current_area, storage) - - for (i, id) in enumerate(node_id) + for (i, id) in enumerate(basin.node_id) # add all precipitation that falls within the profile - level = current_level[i] - area = current_area[i] - - bottom = basin.level[i][1] - fixed_area = basin.area[i][end] - depth = max(level - bottom, 0.0) - factor = reduction_factor(depth, 0.1) - - precipitation = fixed_area * basin.precipitation[i] - evaporation = area * factor * basin.potential_evaporation[i] - drainage = basin.drainage[i] - infiltration = factor * basin.infiltration[i] - - influx = precipitation - evaporation + drainage - infiltration - du.storage[i] += influx - set_flow!(graph, id, influx) + update_vertical_flux!(basin, storage, i) + du.storage[i] += get_influx(basin, i) end return nothing end @@ -303,7 +310,6 @@ function formulate_flow!( # Return flow is immediate set_flow!(graph, id, dst_id, q * return_factor[i]) - set_flow!(graph, id, -q * (1 - return_factor[i])) end return nothing end @@ -505,28 +511,6 @@ function formulate_flow!( return nothing end -function formulate_flow!( - level_boundary::LevelBoundary, - p::Parameters, - storage::AbstractVector, - t::Number, -)::Nothing - (; graph) = p - (; node_id) = level_boundary - - for id in node_id - for in_id in inflow_ids(graph, id) - q = get_flow(graph, in_id, id, storage) - add_flow!(graph, id, -q) - end - for out_id in outflow_ids(graph, id) - q = get_flow(graph, id, out_id, storage) - add_flow!(graph, id, q) - end - end - return nothing -end - function formulate_flow!( flow_boundary::FlowBoundary, p::Parameters, @@ -547,7 +531,6 @@ function formulate_flow!( # Adding water is always possible set_flow!(graph, id, dst_id, rate) - set_flow!(graph, id, rate) end end end @@ -663,6 +646,5 @@ function formulate_flows!(p::Parameters, storage::AbstractVector, t::Number)::No # do these last since they rely on formulated input flows formulate_flow!(fractional_flow, p, storage, t) - formulate_flow!(level_boundary, p, storage, t) formulate_flow!(terminal, p, storage, t) end diff --git a/core/src/util.jl b/core/src/util.jl index 3b2131d1e..879f89726 100644 --- a/core/src/util.jl +++ b/core/src/util.jl @@ -707,3 +707,41 @@ function Base.get( nothing end end + +""" +Get the time interval between (flow) saves +""" +function get_Δt(integrator)::Float64 + (; p, t, dt) = integrator + (; saveat) = p.graph[] + if iszero(saveat) + dt + elseif isinf(saveat) + t + else + t_end = integrator.sol.prob.tspan[2] + if t_end - t > saveat + saveat + else + # The last interval might be shorter than saveat + rem = t % saveat + iszero(rem) ? saveat : rem + end + end +end + +function get_influx(basin::Basin, node_id::NodeID)::Float64 + has_index, basin_idx = id_index(basin.node_id, node_id) + if !has_index + error("Sum of vertical fluxes requested for non-basin $id.") + end + return get_influx(basin, basin_idx) +end + +function get_influx(basin::Basin, basin_idx::Int)::Float64 + (; vertical_flux) = basin + vertical_flux = get_tmp(vertical_flux, 0) + (; precipitation, evaporation, drainage, infiltration) = vertical_flux + return precipitation[basin_idx] - evaporation[basin_idx] + drainage[basin_idx] - + infiltration[basin_idx] +end diff --git a/core/src/write.jl b/core/src/write.jl index 9f855e324..4964f1754 100644 --- a/core/src/write.jl +++ b/core/src/write.jl @@ -84,15 +84,53 @@ function basin_table( node_id::Vector{Int32}, storage::Vector{Float64}, level::Vector{Float64}, + precipitation::Vector{Union{Missing, Float64}}, + evaporation::Vector{Union{Missing, Float64}}, + drainage::Vector{Union{Missing, Float64}}, + infiltration::Vector{Union{Missing, Float64}}, } + (; saved) = model + (; vertical_flux) = saved + data = get_storages_and_levels(model) + storage = vec(data.storage) + level = vec(data.level) + nbasin = length(data.node_id) ntsteps = length(data.time) + nrows = nbasin * ntsteps + + precipitation = Vector{Union{Missing, Float64}}(missing, nrows) + evaporation = Vector{Union{Missing, Float64}}(missing, nrows) + drainage = Vector{Union{Missing, Float64}}(missing, nrows) + infiltration = Vector{Union{Missing, Float64}}(missing, nrows) + + idx_row = 0 + + for vec in vertical_flux.saveval + for (precipitation_, evaporation_, drainage_, infiltration_) in + zip(vec.precipitation, vec.evaporation, vec.drainage, vec.infiltration) + idx_row += 1 + precipitation[idx_row] = precipitation_ + evaporation[idx_row] = evaporation_ + drainage[idx_row] = drainage_ + infiltration[idx_row] = infiltration_ + end + end time = repeat(data.time; inner = nbasin) node_id = repeat(Int32.(data.node_id); outer = ntsteps) - return (; time, node_id, storage = vec(data.storage), level = vec(data.level)) + return (; + time, + node_id, + storage, + level, + precipitation, + evaporation, + drainage, + infiltration, + ) end "Create a flow result table from the saved data" @@ -110,28 +148,14 @@ function flow_table( (; config, saved, integrator) = model (; t, saveval) = saved.flow (; graph) = integrator.p - (; flow_dict, flow_vertical_dict) = graph[] + (; flow_dict) = graph[] - # self-loops have no edge ID from_node_type = String[] from_node_id = Int32[] to_node_type = String[] to_node_id = Int32[] unique_edge_ids_flow = Union{Int32, Missing}[] - vertical_flow_node_ids = Vector{NodeID}(undef, length(flow_vertical_dict)) - for (node_id, index) in flow_vertical_dict - vertical_flow_node_ids[index] = node_id - end - - for id in vertical_flow_node_ids - push!(from_node_type, string(id.type)) - push!(from_node_id, id.value) - push!(to_node_type, string(id.type)) - push!(to_node_id, id.value) - push!(unique_edge_ids_flow, missing) - end - flow_edge_ids = Vector{Tuple{NodeID, NodeID}}(undef, length(flow_dict)) for (edge_id, index) in flow_dict flow_edge_ids[index] = edge_id diff --git a/core/test/bmi_test.jl b/core/test/bmi_test.jl index fd33ed312..7b37b2017 100644 --- a/core/test/bmi_test.jl +++ b/core/test/bmi_test.jl @@ -63,6 +63,8 @@ end "basin.level", "basin.infiltration", "basin.drainage", + "basin.infiltration_integrated", + "basin.drainage_integrated", "basin.subgrid_level", "user_demand.demand", "user_demand.realized", diff --git a/core/test/run_models_test.jl b/core/test/run_models_test.jl index 8f285c810..1a41b33b6 100644 --- a/core/test/run_models_test.jl +++ b/core/test/run_models_test.jl @@ -44,8 +44,26 @@ (DateTime, Union{Int32, Missing}, String, Int32, String, Int32, Float64), ) @test Tables.schema(basin) == Tables.Schema( - (:time, :node_id, :storage, :level), - (DateTime, Int32, Float64, Float64), + ( + :time, + :node_id, + :storage, + :level, + :precipitation, + :evaporation, + :drainage, + :infiltration, + ), + ( + DateTime, + Int32, + Float64, + Float64, + Union{Missing, Float64}, + Union{Missing, Float64}, + Union{Missing, Float64}, + Union{Missing, Float64}, + ), ) @test Tables.schema(control) == Tables.Schema( (:time, :control_node_id, :truth_state, :control_state), @@ -88,8 +106,8 @@ @testset "Results size" begin nsaved = length(tsaves(model)) @test nsaved > 10 - # t0 has no flow, 2 flow edges and 2 boundary condition flows - @test nrow(flow) == (nsaved - 1) * 4 + # t0 has no flow, 2 flow edges + @test nrow(flow) == (nsaved - 1) * 2 @test nrow(basin) == nsaved @test nrow(control) == 0 @test nrow(allocation) == 0 @@ -98,9 +116,9 @@ @testset "Results values" begin @test flow.time[1] == DateTime(2020) - @test coalesce.(flow.edge_id[1:4], -1) == [-1, -1, 1, 2] - @test flow.from_node_id[1:4] == [6, 922, 6, 0] - @test flow.to_node_id[1:4] == [6, 922, 0, 922] + @test coalesce.(flow.edge_id[1:2], -1) == [1, 2] + @test flow.from_node_id[1:2] == [6, 0] + @test flow.to_node_id[1:2] == [0, 922] @test basin.storage[1] ≈ 1.0 @test basin.level[1] ≈ 0.044711584 @@ -125,10 +143,11 @@ end model = Ribasim.run(toml_path) @test model isa Ribasim.Model @test model.integrator.u.storage ≈ [1000] - @test model.integrator.p.basin.precipitation == [0.0] - @test model.integrator.p.basin.potential_evaporation == [0.0] - @test model.integrator.p.basin.drainage == [0.0] - @test model.integrator.p.basin.infiltration == [0.0] + vertical_flux = Ribasim.get_tmp(model.integrator.p.basin.vertical_flux, 0) + @test vertical_flux.precipitation == [0.0] + @test vertical_flux.evaporation == [0.0] + @test vertical_flux.drainage == [0.0] + @test vertical_flux.infiltration == [0.0] @test successful_retcode(model) end @@ -142,10 +161,11 @@ end @test model isa Ribasim.Model stor = model.integrator.u.storage - prec = model.integrator.p.basin.precipitation - evap = model.integrator.p.basin.potential_evaporation - drng = model.integrator.p.basin.drainage - infl = model.integrator.p.basin.infiltration + vertical_flux = Ribasim.get_tmp(model.integrator.p.basin.vertical_flux, 0) + prec = vertical_flux.precipitation + evap = vertical_flux.evaporation + drng = vertical_flux.drainage + infl = vertical_flux.infiltration # The dynamic data has missings, but these are not set. @test prec == [0.0] @test evap == [0.0] @@ -235,7 +255,8 @@ end @test model isa Ribasim.Model @test successful_retcode(model) @test allunique(Ribasim.tsaves(model)) - @test length(model.integrator.p.basin.precipitation) == 4 + precipitation = Ribasim.get_tmp(model.integrator.p.basin.vertical_flux, 0).precipitation + @test length(precipitation) == 4 @test model.integrator.sol.u[end] ≈ Float32[472.02444, 472.02252, 367.6387, 1427.981] skip = Sys.isapple() end diff --git a/core/test/time_test.jl b/core/test/time_test.jl index 2410f1933..7977ad705 100644 --- a/core/test/time_test.jl +++ b/core/test/time_test.jl @@ -12,22 +12,77 @@ flow = DataFrame(Ribasim.flow_table(model)) # only from March to September the FlowBoundary varies is_summer(t::DateTime) = 3 <= month(t) < 10 - - flow_added_1 = - filter( - [:time, :from_node_id, :to_node_id] => - (t, from, to) -> is_summer(t) && from == 1 && to == 1, - flow, - ).flow_rate flow_1_to_2 = filter( [:time, :from_node_id, :to_node_id] => (t, from, to) -> is_summer(t) && from == 1 && to == 2, flow, ) - @test flow_added_1 == flow_1_to_2.flow_rate - t = Ribasim.seconds_since.(flow_1_to_2.time, model.config.starttime) flow_expected = @. 1 + sin(0.5 * π * (t - t[1]) / (t[end] - t[1]))^2 # some difference is expected since the modeled flow is for the period up to t - @test isapprox(flow_added_1, flow_expected, rtol = 0.005) + @test isapprox(flow_1_to_2.flow_rate, flow_expected, rtol = 0.005) +end + +@testitem "vertical_flux_means" begin + using DataFrames: DataFrame, transform!, ByRow + + toml_path = + normpath(@__DIR__, "../../generated_testmodels/basic_transient/ribasim.toml") + @test ispath(toml_path) + model = Ribasim.run(toml_path) + basin = model.integrator.p.basin + n_basin = length(basin.node_id) + basin_table = DataFrame(Ribasim.basin_table(model)) + + # No vertical flux data for last saveat + t_end = last(basin_table).time + data_end = filter(:time => t -> t == t_end, basin_table) + @test all(ismissing.(data_end.precipitation)) + @test all(ismissing.(data_end.evaporation)) + @test all(ismissing.(data_end.drainage)) + @test all(ismissing.(data_end.infiltration)) + + time_table = DataFrame(basin.time) + time_table[!, "basin_idx"] = [ + Ribasim.id_index(basin.node_id, node_id)[2] for + node_id in Ribasim.NodeID.(:Basin, time_table.node_id) + ] + time_table[!, "area"] = [ + Ribasim.get_area_and_level(basin, idx, storage)[1] for + (idx, storage) in zip(time_table.basin_idx, basin_table.storage) + ] + # Mean areas are sufficient to compute the mean flows + # (assuming the saveats coincide with the solver timepoints), + # as the potential evaporation is constant over the saveat intervals + time_table[!, "mean_area"] .= 0.0 + n_basins = length(basin.node_id) + n_times = length(unique(time_table.time)) - 1 + for basin_idx in 1:n_basins + for time_idx in 1:n_times + idx_1 = n_basins * (time_idx - 1) + basin_idx + idx_2 = n_basins * time_idx + basin_idx + mean_area = (time_table.area[idx_1] + time_table.area[idx_2]) / 2 + time_table.mean_area[idx_1] = mean_area + end + end + + filter!(:time => t -> t !== t_end, basin_table) + filter!(:time => t -> t !== t_end, time_table) + + @test all( + isapprox( + basin_table.evaporation, + time_table.mean_area .* time_table.potential_evaporation; + rtol = 1e-4, + ), + ) + + fixed_area = Dict( + id.value => basin.area[Ribasim.id_index(basin.node_id, id)[2]][end] for + id in basin.node_id + ) + transform!(time_table, :node_id => ByRow(id -> fixed_area[id]) => :fixed_area) + @test all( + basin_table.precipitation .≈ time_table.fixed_area .* time_table.precipitation, + ) end diff --git a/core/test/utils_test.jl b/core/test/utils_test.jl index b5a9ee0e3..c5e4d7af3 100644 --- a/core/test/utils_test.jl +++ b/core/test/utils_test.jl @@ -43,6 +43,7 @@ end [2.0, 3.0], [2.0, 3.0], [2.0, 3.0], + [2.0, 3.0], darea, area, level, @@ -121,6 +122,7 @@ end zeros(1), zeros(1), zeros(1), + zeros(1), [area], [level], [storage], diff --git a/docs/core/usage.qmd b/docs/core/usage.qmd index e3da5d330..f8fec8b9c 100644 --- a/docs/core/usage.qmd +++ b/docs/core/usage.qmd @@ -705,27 +705,36 @@ derivative | Float64 | - | - ## Basin - `basin.arrow` -The basin table contains results of the storage and level of each basin at every solver -timestep. The initial condition is also written to the file. +The basin table contains: + +- results of the storage and level of each basin, which are instantaneous values; +- results of the vertical fluxes on each basin, which are mean values over the `saveat` intervals. In the time column the start of the period is indicated. This means that for the final time in the table the mean vertical fluxes are not computed, and thus are `missing`. + +The initial condition is also written to the file. + +column | type | unit +------------- | ------------------------ | ---- +time | DateTime | - +node_id | Int32 | - +storage | Float64 | $m^3$ +level | Float64 | $m$ +precipitation | Union{Float64, Missing} | $m^3 s^{-1}$ +evaporation | Union{Float64, Missing} | $m^3 s^{-1}$ +drainage | Union{Float64, Missing} | $m^3 s^{-1}$ +infiltration | Union{Float64, Missing} | $m^3 s^{-1}$ -column | type | unit --------- | -------- | ---- -time | DateTime | - -node_id | Int32 | - -storage | Float64 | $m^3$ -level | Float64 | $m$ The table is sorted by time, and per time it is sorted by `node_id`. ## Flow - `flow.arrow` -The flow table contains calculated mean flows for every flow edge in the model. +The flow table contains calculated mean flows over the `saveat` intervals for every flow edge in the model. In the time column the start of the period is indicated. column | type | unit -------------- | --------------------- | ---- time | DateTime | - -edge_id | Union{Int32, Missing} | - +edge_id | Int32 | - from_node_type | String | - from_node_id | Int32 | - to_node_type | String | - @@ -734,8 +743,7 @@ flow_rate | Float64 | $m^3 s^{-1}$ The table is sorted by time, and per time the same `edge_id` order is used, though not sorted. The `edge_id` value is the same as the `fid` written to the Edge table, and can be used to directly look up the Edge geometry. -Flows that are added to the model at a node, have a missing `edge_id`, and identical `from_node_id` and `to_node_id`. -Flows out of the model always have a negative sign, and additions a positive sign. +Flows from the "from" to the "to" node have a positive sign, and if the flow is reversed it will be negative. ## DiscreteControl - `control.arrow`