Mean output flows

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)
@@ -40,6 +40,9 @@ function create_callbacks(
     basin_cb = PresetTimeCallback(tstops, update_basin)
     push!(callbacks, basin_cb)
 
+    integrating_flows_cb = FunctionCallingCallback(integrate_flows!)
+    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)
@@ -87,6 +90,29 @@ function create_callbacks(
     return callback, saved
 end
 
+function integrate_flows!(u, t, integrator)::Nothing
+    (; p, tprev) = 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)
+    Δt = t - tprev
+
+    @. flow_integrated += 0.5 * (flow + flow_prev) * Δt
+    @. flow_vertical_integrated += 0.5 * (flow_vertical + flow_vertical_prev) * Δt
+
+    copyto!(flow_prev, flow)
+    copyto!(flow_vertical_prev, flow_vertical)
+    return nothing
+end
+
 """
 Listens for changes in condition truths.
 """
@@ -379,10 +405,17 @@ end
 
 "Copy the current flow to the 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),
-    )
+    (; graph) = integrator.p
+    (; flow_integrated, flow_vertical_integrated, tprev_flow_save) = graph[]
+    Δt = t - tprev_flow_save[1]
+    tprev_flow_save[1] = t
+
+    mean_flow_vertical = flow_vertical_integrated / Δt
+    mean_flow = flow_integrated / Δt
+
+    fill!(flow_vertical_integrated, 0.0)
+    fill!(flow_integrated, 0.0)
+    return vcat(mean_flow_vertical, mean_flow)
 end
 
 function update_subgrid_level!(integrator)::Nothing

core/src/graph.jl

@@ -90,19 +90,29 @@ function create_graph(db::DB, config::Config, chunk_sizes::Vector{Int})::MetaGra
     end
 
     flow = zeros(flow_counter)
+    flow_prev = zeros(flow_counter)
+    flow_integrated = zeros(flow_counter)
     flow_vertical = zeros(flow_vertical_counter)
+    flow_vertical_prev = zeros(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
+    tprev_flow_save = zeros(1)
     graph_data = (;
         node_ids,
         edge_ids,
         edges_source,
         flow_dict,
         flow,
+        flow_prev,
+        flow_integrated,
         flow_vertical_dict,
         flow_vertical,
+        flow_vertical_prev,
+        flow_vertical_integrated,
+        tprev_flow_save,
     )
     graph = @set graph.graph_data = graph_data
 

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},
+            tprev_flow_save::Vector{Float64},
         },
         MetaGraphsNext.var"#11#13",
         Float64,

core/test/run_models_test.jl

@@ -476,3 +476,16 @@ end
         0,
     ) ≈ 5.0 atol = 0.001 skip = Sys.isapple()
 end
+
+@testitem "mean_flow" begin
+    toml_path =
+        normpath(@__DIR__, "../../generated_testmodels/linear_resistance/ribasim.toml")
+    @test ispath(toml_path)
+    t_end = 3.16224e7
+
+    config = Ribasim.Config(toml_path; solver_saveat = t_end)
+    model = Ribasim.run(config)
+    mean_flow = model.saved.flow.saveval[1]
+    storage = Ribasim.get_storages_and_levels(model).storage[1, :]
+    @test storage[2] ≈ storage[1] - t_end * mean_flow[3] atol = 1
+end

core/test/time_test.jl

@@ -29,5 +29,5 @@
     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_added_1, flow_expected, rtol = 0.01)
 end