From 51ed6ed35f042464a39baba6d86021d14feff637 Mon Sep 17 00:00:00 2001 From: Bart de Koning <74617371+SouthEndMusic@users.noreply.github.com> Date: Wed, 28 Feb 2024 11:20:00 +0100 Subject: [PATCH] Mean output flows (#1159) Fixes https://github.com/Deltares/Ribasim/issues/935. --------- Co-authored-by: Martijn Visser --- core/src/callback.jl | 82 ++++++++++++++--- core/src/config.jl | 15 +++- core/src/graph.jl | 15 +++- core/src/model.jl | 17 ++-- core/src/parameter.jl | 5 ++ core/src/solve.jl | 14 --- core/src/write.jl | 9 +- core/test/allocation_test.jl | 2 +- core/test/config_test.jl | 2 +- core/test/control_test.jl | 36 ++++---- core/test/equations_test.jl | 18 ++-- core/test/run_models_test.jl | 88 ++++++++++++++++--- docs/core/numerics.qmd | 2 +- docs/core/usage.qmd | 7 +- .../ribasim_testmodels/discrete_control.py | 26 +++--- 15 files changed, 241 insertions(+), 97 deletions(-) diff --git a/core/src/callback.jl b/core/src/callback.jl index 9403416aa..7b312f409 100644 --- a/core/src/callback.jl +++ b/core/src/callback.jl @@ -7,7 +7,7 @@ function set_initial_discrete_controlled_parameters!( storage0::Vector{Float64}, )::Nothing (; p) = integrator - (; basin, discrete_control) = p + (; discrete_control) = p n_conditions = length(discrete_control.condition_value) condition_diffs = zeros(Float64, n_conditions) @@ -30,7 +30,7 @@ Returns the CallbackSet and the SavedValues for flow. """ function create_callbacks( parameters::Parameters, - config::Config; + config::Config, saveat, )::Tuple{CallbackSet, SavedResults} (; starttime, basin, tabulated_rating_curve, discrete_control) = parameters @@ -40,6 +40,9 @@ function create_callbacks( basin_cb = PresetTimeCallback(tstops, update_basin) 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, update_tabulated_rating_curve!) push!(callbacks, tabulated_rating_curve_cb) @@ -55,7 +58,14 @@ function create_callbacks( # 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; saveat, save_start = false) + 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, + ) push!(callbacks, save_flow_cb) # interpolate the levels @@ -87,6 +97,37 @@ function create_callbacks( return callback, saved end +""" +Integrate flows over the last timestep +""" +function integrate_flows!(u, t, integrator)::Nothing + (; p, dt) = integrator + (; graph) = p + (; + flow, + flow_vertical, + flow_prev, + flow_vertical_prev, + flow_integrated, + flow_vertical_integrated, + ) = graph[] + flow = get_tmp(flow, 0) + flow_vertical = get_tmp(flow_vertical, 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 + + copyto!(flow_prev, flow) + copyto!(flow_vertical_prev, flow_vertical) + return nothing +end + """ Listens for changes in condition truths. """ @@ -272,8 +313,7 @@ function discrete_control_affect!( # What the local control state is # TODO: Check time elapsed since control change - control_state_now, control_state_start = - discrete_control.control_state[discrete_control_node_id] + control_state_now, _ = discrete_control.control_state[discrete_control_node_id] control_state_change = false @@ -377,12 +417,34 @@ function set_control_params!(p::Parameters, node_id::NodeID, control_state::Stri end end -"Copy the current flow to the SavedValues" +"Compute the average flows over the last saveat interval and write +them to SavedValues" function save_flow(u, t, integrator) - vcat( - get_tmp(integrator.p.graph[].flow_vertical, 0.0), - get_tmp(integrator.p.graph[].flow, 0.0), - ) + (; 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 diff --git a/core/src/config.jl b/core/src/config.jl index 07dc6afa4..831c82106 100644 --- a/core/src/config.jl +++ b/core/src/config.jl @@ -233,19 +233,26 @@ end "Convert the saveat Float64 from our Config to SciML's saveat" function convert_saveat(saveat::Float64, t_end::Float64)::Union{Float64, Vector{Float64}} + errors = false if iszero(saveat) # every step - Float64[] + saveat = Float64[] elseif saveat == Inf # only the start and end - [0.0, t_end] + saveat = [0.0, t_end] elseif isfinite(saveat) # every saveat seconds - saveat + if saveat !== round(saveat) + errors = true + @error "A finite saveat must be an integer number of seconds." saveat + end else + errors = true @error "Invalid saveat" saveat - error("Invalid saveat") end + + errors && error("Invalid saveat") + return saveat end "Convert the dt from our Config to SciML stepsize control arguments" diff --git a/core/src/graph.jl b/core/src/graph.jl index 8ac7868f0..6ebf2a8a4 100644 --- a/core/src/graph.jl +++ b/core/src/graph.jl @@ -6,7 +6,8 @@ and data of edges (EdgeMetadata): [`EdgeMetadata`](@ref) """ function create_graph(db::DB, config::Config, chunk_sizes::Vector{Int})::MetaGraph - node_rows = execute(db, "SELECT node_id, node_type, subnetwork_id FROM Node ORDER BY fid") + node_rows = + execute(db, "SELECT node_id, node_type, subnetwork_id FROM Node ORDER BY fid") edge_rows = execute( db, "SELECT fid, from_node_type, from_node_id, to_node_type, to_node_id, edge_type, subnetwork_id FROM Edge ORDER BY fid", @@ -44,7 +45,8 @@ function create_graph(db::DB, config::Config, chunk_sizes::Vector{Int})::MetaGra end push!(node_ids[allocation_network_id], node_id) end - graph[node_id] = NodeMetadata(Symbol(snake_case(row.node_type)), allocation_network_id) + 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 @@ -90,7 +92,11 @@ function create_graph(db::DB, config::Config, chunk_sizes::Vector{Int})::MetaGra end 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) @@ -101,8 +107,13 @@ function create_graph(db::DB, config::Config, chunk_sizes::Vector{Int})::MetaGra edges_source, flow_dict, 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 diff --git a/core/src/model.jl b/core/src/model.jl index a7899bfd6..009cd80c2 100644 --- a/core/src/model.jl +++ b/core/src/model.jl @@ -76,10 +76,10 @@ function Model(config::Config)::Model # tell the solver to stop when new data comes in # TODO add all time tables here time_flow_boundary = load_structvector(db, config, FlowBoundaryTimeV1) - tstops_flow_boundary = get_tstops(time_flow_boundary.time, config.starttime) + tstops = Vector{Float64}[] + push!(tstops, get_tstops(time_flow_boundary.time, config.starttime)) time_user_demand = load_structvector(db, config, UserDemandTimeV1) - tstops_user_demand = get_tstops(time_user_demand.time, config.starttime) - tstops = sort(unique(vcat(tstops_flow_boundary, tstops_user_demand))) + push!(tstops, get_tstops(time_user_demand.time, config.starttime)) # use state state = load_structvector(db, config, BasinStateV1) @@ -105,7 +105,10 @@ function Model(config::Config)::Model @assert eps(t_end) < 3600 "Simulation time too long" t0 = zero(t_end) timespan = (t0, t_end) + saveat = convert_saveat(config.solver.saveat, t_end) + saveat isa Float64 && push!(tstops, range(0, t_end; step = saveat)) + tstops = sort(unique(vcat(tstops...))) adaptive, dt = convert_dt(config.solver.dt) jac_prototype = config.solver.sparse ? get_jac_prototype(parameters) : nothing @@ -116,7 +119,7 @@ function Model(config::Config)::Model end @debug "Setup ODEProblem." - callback, saved = create_callbacks(parameters, config; saveat) + callback, saved = create_callbacks(parameters, config, saveat) @debug "Created callbacks." # Initialize the integrator, providing all solver options as described in @@ -154,17 +157,17 @@ function Model(config::Config)::Model end "Get all saved times in seconds since start" -timesteps(model::Model)::Vector{Float64} = model.integrator.sol.t +tsaves(model::Model)::Vector{Float64} = model.integrator.sol.t "Get all saved times as a Vector{DateTime}" function datetimes(model::Model)::Vector{DateTime} - return datetime_since.(timesteps(model), model.config.starttime) + return datetime_since.(tsaves(model), model.config.starttime) end function Base.show(io::IO, model::Model) (; config, integrator) = model t = datetime_since(integrator.t, config.starttime) - nsaved = length(timesteps(model)) + nsaved = length(tsaves(model)) println(io, "Model(ts: $nsaved, t: $t)") end diff --git a/core/src/parameter.jl b/core/src/parameter.jl index fe36b0070..be333a8a9 100644 --- a/core/src/parameter.jl +++ b/core/src/parameter.jl @@ -558,8 +558,13 @@ struct Parameters{T, C1, C2} edges_source::Dict{Int, Set{EdgeMetadata}}, flow_dict::Dict{Tuple{NodeID, NodeID}, Int}, 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, diff --git a/core/src/solve.jl b/core/src/solve.jl index 01fee6c6c..b090dfa26 100644 --- a/core/src/solve.jl +++ b/core/src/solve.jl @@ -622,8 +622,6 @@ function formulate_du!( basin::Basin, storage::AbstractVector, )::Nothing - (; flow_vertical) = graph[] - flow_vertical = get_tmp(flow_vertical, storage) # loop over basins # subtract all outgoing flows # add all ingoing flows @@ -665,15 +663,3 @@ function formulate_flows!(p::Parameters, storage::AbstractVector, t::Number)::No formulate_flow!(level_boundary, p, storage, t) formulate_flow!(terminal, p, storage, t) end - -function track_waterbalance!(u, t, integrator)::Nothing - (; p, tprev, uprev) = integrator - dt = t - tprev - du = u - uprev - p.storage_diff .+= du - p.precipitation.total .+= p.precipitation.value .* dt - p.evaporation.total .+= p.evaporation.value .* dt - p.infiltration.total .+= p.infiltration.value .* dt - p.drainage.total .+= p.drainage.value .* dt - return nothing -end diff --git a/core/src/write.jl b/core/src/write.jl index 17f21a8bf..b94ef0c36 100644 --- a/core/src/write.jl +++ b/core/src/write.jl @@ -64,7 +64,7 @@ function get_storages_and_levels( (; sol, p) = integrator node_id = p.basin.node_id.values::Vector{NodeID} - tsteps = datetime_since.(timesteps(model), config.starttime) + tsteps = datetime_since.(tsaves(model), config.starttime) storage = hcat([collect(u_.storage) for u_ in sol.u]...) level = zero(storage) @@ -148,7 +148,12 @@ function flow_table( nflow = length(unique_edge_ids_flow) ntsteps = length(t) - time = repeat(datetime_since.(t, config.starttime); inner = nflow) + # the timestamp should represent the start of the period, not the end + t_starts = circshift(t, 1) + if !isempty(t) + t_starts[1] = 0.0 + end + time = repeat(datetime_since.(t_starts, config.starttime); inner = nflow) edge_id = repeat(unique_edge_ids_flow; outer = ntsteps) from_node_type = repeat(from_node_type; outer = ntsteps) from_node_id = repeat(from_node_id; outer = ntsteps) diff --git a/core/test/allocation_test.jl b/core/test/allocation_test.jl index 1c44a7eb4..dd361397e 100644 --- a/core/test/allocation_test.jl +++ b/core/test/allocation_test.jl @@ -341,7 +341,7 @@ end model = Ribasim.run(toml_path) storage = Ribasim.get_storages_and_levels(model).storage[1, :] - t = Ribasim.timesteps(model) + t = Ribasim.tsaves(model) p = model.integrator.p (; user_demand, graph, allocation, basin, level_demand) = p diff --git a/core/test/config_test.jl b/core/test/config_test.jl index d07170705..1c4aacf47 100644 --- a/core/test/config_test.jl +++ b/core/test/config_test.jl @@ -68,9 +68,9 @@ end @test convert_saveat(0.0, t_end) == Float64[] @test convert_saveat(60.0, t_end) == 60.0 @test convert_saveat(Inf, t_end) == [0.0, t_end] - @test convert_saveat(Inf, t_end) == [0.0, t_end] @test_throws ErrorException convert_saveat(-Inf, t_end) @test_throws ErrorException convert_saveat(NaN, t_end) + @test_throws ErrorException convert_saveat(3.1415, t_end) @test convert_dt(nothing) == (true, 0.0) @test convert_dt(360.0) == (false, 360.0) diff --git a/core/test/control_test.jl b/core/test/control_test.jl index bd35c147f..ae572e45e 100644 --- a/core/test/control_test.jl +++ b/core/test/control_test.jl @@ -32,15 +32,15 @@ ["off", "active", "on", "off", "inactive"] level = Ribasim.get_storages_and_levels(model).level - timesteps = Ribasim.timesteps(model) + t = Ribasim.tsaves(model) # Control times t_1 = discrete_control.record.time[3] - t_1_index = findfirst(timesteps .≈ t_1) + t_1_index = findfirst(t .≈ t_1) @test level[1, t_1_index] ≈ discrete_control.greater_than[1] t_2 = discrete_control.record.time[4] - t_2_index = findfirst(timesteps .≈ t_2) + t_2_index = findfirst(t .≈ t_2) @test level[2, t_2_index] ≈ discrete_control.greater_than[2] flow = get_tmp(graph[].flow, 0) @@ -56,9 +56,9 @@ end Δt = discrete_control.look_ahead[1] - timesteps = Ribasim.timesteps(model) + t = Ribasim.tsaves(model) t_control = discrete_control.record.time[2] - t_control_index = searchsortedfirst(timesteps, t_control) + t_control_index = searchsortedfirst(t, t_control) greater_than = discrete_control.greater_than[1] flow_t_control = flow_boundary.flow_rate[1](t_control) @@ -80,9 +80,9 @@ end Δt = discrete_control.look_ahead[1] - timesteps = Ribasim.timesteps(model) + t = Ribasim.tsaves(model) t_control = discrete_control.record.time[2] - t_control_index = searchsortedfirst(timesteps, t_control) + t_control_index = searchsortedfirst(t, t_control) greater_than = discrete_control.greater_than[1] level_t_control = level_boundary.level[1](t_control) @@ -100,11 +100,11 @@ end (; basin, pid_control, flow_boundary) = p level = Ribasim.get_storages_and_levels(model).level[1, :] - timesteps = Ribasim.timesteps(model) + t = Ribasim.tsaves(model) target_itp = pid_control.target[1] t_target_change = target_itp.t[2] - idx_target_change = searchsortedlast(timesteps, t_target_change) + idx_target_change = searchsortedlast(t, t_target_change) K_p, K_i, _ = pid_control.pid_params[2](0) level_demand = pid_control.target[2](0) @@ -119,15 +119,13 @@ end phi = atan(du0 / (A * Δlevel) - alpha) / omega a = abs(Δlevel / cos(phi)) # This bound is the exact envelope of the analytical solution - bound = @. a * exp(alpha * timesteps[1:idx_target_change]) + bound = @. a * exp(alpha * t[1:idx_target_change]) eps = 5e-3 # Initial convergence to target level @test all(@. abs(level[1:idx_target_change] - level_demand) < bound + eps) # Later closeness to target level @test all( - @. abs( - level[idx_target_change:end] - target_itp(timesteps[idx_target_change:end]), - ) < 5e-2 + @. abs(level[idx_target_change:end] - target_itp(t[idx_target_change:end])) < 5e-2 ) end @@ -163,7 +161,7 @@ end (; discrete_control) = p (; record, greater_than) = discrete_control level = Ribasim.get_storages_and_levels(model).level[1, :] - timesteps = Ribasim.timesteps(model) + t = Ribasim.tsaves(model) t_none_1 = discrete_control.record.time[2] t_in = discrete_control.record.time[3] @@ -172,9 +170,9 @@ end level_min = greater_than[1] setpoint = greater_than[2] - t_1_none_index = findfirst(timesteps .≈ t_none_1) - t_in_index = findfirst(timesteps .≈ t_in) - t_2_none_index = findfirst(timesteps .≈ t_none_2) + t_1_none_index = findfirst(t .≈ t_none_1) + t_in_index = findfirst(t .≈ t_in) + t_2_none_index = findfirst(t .≈ t_none_2) @test record.control_state == ["out", "none", "in", "none"] @test level[t_1_none_index] ≈ setpoint @@ -194,7 +192,7 @@ end p = model.integrator.p (; discrete_control, pid_control) = p - timesteps = Ribasim.timesteps(model) + t = Ribasim.tsaves(model) level = Ribasim.get_storages_and_levels(model).level[1, :] target_high = @@ -203,7 +201,7 @@ end pid_control.control_mapping[(NodeID(:PidControl, 6), "target_low")].target.u[1] t_target_jump = discrete_control.record.time[2] - t_idx_target_jump = searchsortedlast(timesteps, t_target_jump) + t_idx_target_jump = searchsortedlast(t, t_target_jump) @test isapprox(level[t_idx_target_jump], target_high, atol = 1e-1) @test isapprox(level[end], target_low, atol = 1e-1) diff --git a/core/test/equations_test.jl b/core/test/equations_test.jl index 085dc31ac..95ca06d03 100644 --- a/core/test/equations_test.jl +++ b/core/test/equations_test.jl @@ -25,7 +25,7 @@ @test successful_retcode(model) (; p) = model.integrator - t = Ribasim.timesteps(model) + t = Ribasim.tsaves(model) storage = Ribasim.get_storages_and_levels(model).storage[1, :] A = p.basin.area[1][2] # needs to be constant u0 = A * 10.0 @@ -57,7 +57,7 @@ end @test successful_retcode(model) p = model.integrator.p - t = Ribasim.timesteps(model) + t = Ribasim.tsaves(model) storage = Ribasim.get_storages_and_levels(model).storage[1, :] basin_area = p.basin.area[1][2] storage_min = 50.005 @@ -93,7 +93,7 @@ end p = model.integrator.p (; manning_resistance) = p - t = Ribasim.timesteps(model) + t = Ribasim.tsaves(model) storage_both = Ribasim.get_storages_and_levels(model).storage storage = storage_both[1, :] storage_min = 50.005 @@ -131,7 +131,7 @@ end (; basin, pid_control) = p storage = Ribasim.get_storages_and_levels(model).storage[:] - t = Ribasim.timesteps(model) + t = Ribasim.tsaves(model) SP = pid_control.target[1](0) K_p, K_i, K_d = pid_control.pid_params[1](0) @@ -168,12 +168,14 @@ end @testitem "MiscellaneousNodes" begin using PreallocationTools: get_tmp using SciMLBase: successful_retcode - using Ribasim: timesteps, get_storages_and_levels + using Ribasim: tsaves, get_storages_and_levels toml_path = normpath(@__DIR__, "../../generated_testmodels/misc_nodes/ribasim.toml") @test ispath(toml_path) config = Ribasim.Config(toml_path) - model = Ribasim.run(config) + model = Ribasim.Model(toml_path) + @test config.solver.dt === model.integrator.dt + Ribasim.solve!(model) @test successful_retcode(model) p = model.integrator.p (; flow_boundary, fractional_flow, pump) = p @@ -184,10 +186,8 @@ end frac = fractional_flow.fraction[1] storage_both = get_storages_and_levels(model).storage - t = timesteps(model) + t = tsaves(model) tspan = model.integrator.sol.prob.tspan - - @test config.solver.dt === model.integrator.dt @test t ≈ range(tspan...; step = config.solver.saveat) @test storage_both[1, :] ≈ @. storage_both[1, 1] + t * (frac * q_boundary - q_pump) @test storage_both[2, :] ≈ @. storage_both[2, 1] + t * q_pump diff --git a/core/test/run_models_test.jl b/core/test/run_models_test.jl index a5f615c65..8a1414275 100644 --- a/core/test/run_models_test.jl +++ b/core/test/run_models_test.jl @@ -4,7 +4,7 @@ using Tables.DataAPI: nrow using Dates: DateTime import Arrow - using Ribasim: timesteps + using Ribasim: get_tstops, tsaves toml_path = normpath(@__DIR__, "../../generated_testmodels/trivial/ribasim.toml") @test ispath(toml_path) @@ -84,7 +84,7 @@ end @testset "Results size" begin - nsaved = length(timesteps(model)) + 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 @@ -95,7 +95,7 @@ end @testset "Results values" begin - @test flow.time[1] > DateTime(2020) + @test flow.time[1] == DateTime(2020) @test coalesce.(flow.edge_id[1:4], -1) == [-1, -1, 9, 11] @test flow.from_node_id[1:4] == [6, typemax(Int), 0, 6] @test flow.to_node_id[1:4] == [6, typemax(Int), typemax(Int), 0] @@ -189,18 +189,18 @@ end @test all(isconcretetype, fieldtypes(typeof(p))) @test successful_retcode(model) - @test allunique(Ribasim.timesteps(model)) + @test allunique(Ribasim.tsaves(model)) @test model.integrator.sol.u[end] ≈ Float32[519.8817, 519.8798, 339.3959, 1418.4331] skip = Sys.isapple() atol = 1.5 - @test length(logger.logs) == 8 + @test length(logger.logs) == 11 @test logger.logs[1].level == Debug @test logger.logs[1].message == "Read database into memory." table = Ribasim.flow_table(model) # flows are recorded at the end of each period, and are undefined at the start - @test unique(table.time) == Ribasim.datetimes(model)[2:end] + @test unique(table.time) == Ribasim.datetimes(model)[1:(end - 1)] # inflow = outflow over FractionalFlow t = table.time[1] @@ -355,7 +355,7 @@ end (; level) = level_boundary level = level[1] - timesteps = model.saved.flow.t + t = model.saved.flow.t flow = DataFrame(Ribasim.flow_table(model)) outlet_flow = filter([:from_node_id, :to_node_id] => (from, to) -> from === 2 && to === 3, flow) @@ -364,14 +364,14 @@ end level.t[2] * (outlet.min_crest_level[1] - level.u[1]) / (level.u[2] - level.u[1]) # No outlet flow when upstream level is below minimum crest level - @test all(@. outlet_flow.flow_rate[timesteps <= t_min_crest_level] == 0) + @test all(@. outlet_flow.flow_rate[t <= t_min_crest_level] == 0) - timesteps = Ribasim.timesteps(model) + t = Ribasim.tsaves(model) t_maximum_level = level.t[2] level_basin = Ribasim.get_storages_and_levels(model).level[:] # Basin level converges to stable level boundary level - all(isapprox.(level_basin[timesteps .>= t_maximum_level], level.u[3], atol = 5e-2)) + all(isapprox.(level_basin[t .>= t_maximum_level], level.u[3], atol = 5e-2)) end @testitem "UserDemand" begin @@ -476,3 +476,71 @@ end 0, ) ≈ 5.0 atol = 0.001 skip = Sys.isapple() end + +@testitem "mean_flow" begin + using DataFrames: DataFrame + + toml_path = + normpath(@__DIR__, "../../generated_testmodels/flow_boundary_time/ribasim.toml") + @test ispath(toml_path) + function get_flow(solver_dt::Union{Float64, Nothing}, solver_saveat::Float64) + config = Ribasim.Config(toml_path; solver_dt, solver_saveat) + model = Ribasim.run(config) + df = DataFrame(Ribasim.flow_table(model)) + flow = + filter( + [:from_node_id, :to_node_id] => (from, to) -> from === 3 && to === 2, + df, + ).flow_rate + flow, Ribasim.tsaves(model) + end + + Δt = 24 * 24 * 60.0 + t_end = 3.16224e7 # 366 days + + # t_end % saveat = 0 + saveat = 86400.0 + flow, tstops = get_flow(nothing, saveat) + @test all(flow .≈ 1.0) + @test length(flow) == t_end / saveat + @test length(tstops) == t_end / saveat + 1 + + flow, tstops = get_flow(Δt, saveat) + @test all(flow .≈ 1.0) + @test length(flow) == t_end / saveat + @test length(tstops) == t_end / saveat + 1 + + # t_end % saveat != 0 + saveat = round(10000 * π) + flow, tstops = get_flow(nothing, saveat) + @test all(flow .≈ 1.0) + @test length(flow) == ceil(t_end / saveat) + @test length(tstops) == ceil(t_end / saveat) + 1 + + flow, tstops = get_flow(Δt, saveat) + @test all(flow .≈ 1.0) + @test length(flow) == ceil(t_end / saveat) + @test length(tstops) == ceil(t_end / saveat) + 1 + + # Only save average over all flows in tspan + saveat = Inf + flow, tstops = get_flow(nothing, saveat) + @test all(flow .≈ 1.0) + @test length(flow) == 1 + @test length(tstops) == 2 + + flow, tstops = get_flow(Δt, saveat) + @test all(flow .≈ 1.0) + @test length(flow) == 1 + @test length(tstops) == 2 + + # Save all flows + saveat = 0.0 + flow, tstops = get_flow(nothing, saveat) + @test all(flow .≈ 1.0) + @test length(flow) == length(tstops) - 1 + + flow, tstops = get_flow(Δt, saveat) + @test all(flow .≈ 1.0) + @test length(flow) == length(tstops) - 1 +end diff --git a/docs/core/numerics.qmd b/docs/core/numerics.qmd index e54f69928..ca25ca183 100644 --- a/docs/core/numerics.qmd +++ b/docs/core/numerics.qmd @@ -106,4 +106,4 @@ Rather, the basin gets a very small storage. The reduction factors help with performance, but are also an important tool to avoid getting negative storage in basins. Negative storage needs to be avoided since it is not a real solution, and would introduce water into the model that doesn't exist. Another tool used to avoid negative storage is the [`isoutoutofdomain`](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/) option, which Ribasim makes use of. -This reject timesteps that lead to negative storage, instead retrying with a smaller timestep. +This rejects timesteps that lead to negative storage, instead retrying with a smaller timestep. diff --git a/docs/core/usage.qmd b/docs/core/usage.qmd index 98110cc70..f40220e24 100644 --- a/docs/core/usage.qmd +++ b/docs/core/usage.qmd @@ -267,8 +267,7 @@ to have a reasonable and safe default, a value must be provided in the static ta This table is the transient form of the `Basin` table. The only difference is that a time column is added. The table must by sorted by time, and per time it must be sorted by `node_id`. -A linear interpolation between the given timesteps is currently done if the -solver takes timesteps between the given data points. More options will be available later. +At the given timestamps the values are set in the simulation, such that the timeseries can be seen as forward filled. ## Basin / state @@ -676,8 +675,8 @@ The table is sorted by time, and per time it is sorted by `node_id`. ## Flow - `flow.arrow` -The flow table contains results of the flow on every edge in the model, for each solver -timestep. +The flow table contains calculated mean flows for every flow edge in the model. +In the time column the start of the period is indicated. column | type | unit -------------- | ------------------- | ---- diff --git a/python/ribasim_testmodels/ribasim_testmodels/discrete_control.py b/python/ribasim_testmodels/ribasim_testmodels/discrete_control.py index 9cc53e10e..70bc4f29f 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/discrete_control.py +++ b/python/ribasim_testmodels/ribasim_testmodels/discrete_control.py @@ -69,7 +69,18 @@ def pump_discrete_control_model() -> ribasim.Model: state = pd.DataFrame(data={"node_id": [1, 3], "level": [1.0, 1e-5]}) - basin = ribasim.Basin(profile=profile, state=state) + static = pd.DataFrame( + data={ + "node_id": [3], + "drainage": [0.0], + "potential_evaporation": [0.0], + "infiltration": [0.0], + "precipitation": [1e-9], + "urban_runoff": [0.0], + } + ) + + basin = ribasim.Basin(profile=profile, static=static, state=state) # Setup the discrete control: condition = pd.DataFrame( @@ -325,20 +336,9 @@ def level_boundary_condition_model(): } ) - static = pd.DataFrame( - data={ - "node_id": [3], - "drainage": [0.0], - "potential_evaporation": [0.0], - "infiltration": [0.0], - "precipitation": [0.0], - "urban_runoff": [0.0], - } - ) - state = pd.DataFrame(data={"node_id": [3], "level": [2.5]}) - basin = ribasim.Basin(profile=profile, static=static, state=state) + basin = ribasim.Basin(profile=profile, state=state) # Setup level boundary: level_boundary = ribasim.LevelBoundary(