From b327a006641c07be72e34f5331b8e6a64c983750 Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Thu, 16 May 2024 17:15:12 +0200 Subject: [PATCH] Bugfixes --- core/src/bmi.jl | 2 +- core/src/callback.jl | 17 +++++++++++------ core/src/util.jl | 11 +++++------ core/src/write.jl | 11 +++++------ core/test/bmi_test.jl | 2 ++ core/test/utils_test.jl | 9 +++------ 6 files changed, 27 insertions(+), 25 deletions(-) diff --git a/core/src/bmi.jl b/core/src/bmi.jl index 06a5e7cd4..f05722297 100644 --- a/core/src/bmi.jl +++ b/core/src/bmi.jl @@ -48,7 +48,7 @@ function BMI.get_value_ptr(model::Model, name::AbstractString)::AbstractVector{F elseif name == "user_demand.demand" vec(model.integrator.p.user_demand.demand) elseif name == "user_demand.realized" - model.saved.stored_for_bmi.integrand.user_realized + model.saved.integrated_bmi_data.user_realized else error("Unknown variable $name") end diff --git a/core/src/callback.jl b/core/src/callback.jl index 16f8ab422..6068dc030 100644 --- a/core/src/callback.jl +++ b/core/src/callback.jl @@ -17,7 +17,11 @@ function TrapezoidIntegrationCallback(integrand_func!, integral_value)::Discrete cache, integral_value, ) - return DiscreteCallback((u, t, integrator) -> t != 0, affect!) + return DiscreteCallback( + (u, t, integrator) -> t != 0, + affect!; + save_positions = (false, false), + ) end function (affect!::TrapezoidIntegrationAffect)(integrator)::Nothing @@ -27,10 +31,11 @@ function (affect!::TrapezoidIntegrationAffect)(integrator)::Nothing copyto!(integrand_value_prev, integrand_value) integrand_func!(integrand_value, integrator.p) - cache += integrand_value_prev - cache += integrand_value - cache *= 0.5 * dt - integral += cache + cache .= integrand_value_prev + cache .+= integrand_value + cache .*= 0.5 * dt + + integral .+= cache return nothing end @@ -164,7 +169,7 @@ function save_flow(u, t, integrator) flow_mean = copy(integrated_flow) flow_mean ./= Δt fill!(integrated_flow, 0.0) - inflow_mean, outflow_mean = compute_mean_inoutflows(flow_mean.flow, graph, basin) + inflow_mean, outflow_mean = compute_mean_inoutflows(flow_mean.flow, basin) return SavedFlow(; flow = flow_mean, inflow = inflow_mean, outflow = outflow_mean) end diff --git a/core/src/util.jl b/core/src/util.jl index f48d3c7ae..44cd77fb9 100644 --- a/core/src/util.jl +++ b/core/src/util.jl @@ -763,26 +763,25 @@ end function compute_mean_inoutflows( flow_mean::AbstractVector, - graph::MetaGraph, basin::Basin, )::Tuple{Vector{Float64}, Vector{Float64}} - (; node_id) = basin + (; node_id, inflow_edges, outflow_edges) = basin # Divide the flows over edges to Basin inflow and outflow, regardless of edge direction. inflow_mean = zeros(length(node_id)) outflow_mean = zeros(length(node_id)) for (i, basin_id) in enumerate(node_id) - for inflow_edge in basin.inflow_edges[i] - q = get_flow(graph, inflow_edge, 0) + for inflow_edge in inflow_edges[i] + q = flow_mean[inflow_edge.flow_idx] if q > 0 inflow_mean[i] += q else outflow_mean[i] -= q end end - for outflow_edge in basin.outflow_edges[i] - q = get_flow(graph, outflow_edge, 0) + for outflow_edge in outflow_edges[i] + q = flow_mean[outflow_edge.flow_idx] if q > 0 inflow_mean[i] += q else diff --git a/core/src/write.jl b/core/src/write.jl index 2002c7d6e..c821c23ce 100644 --- a/core/src/write.jl +++ b/core/src/write.jl @@ -114,8 +114,6 @@ function basin_table( balance_error = zeros(nrows) relative_error = zeros(nrows) - @show length(model.integrator.sol.t) - idx_row = 0 for cvec in saved.flow.saveval for (precipitation_, evaporation_, drainage_, infiltration_) in zip( @@ -132,9 +130,6 @@ function basin_table( end end - @show ntsteps - @show length(saved.flow.saveval) - time = repeat(data.time[begin:(end - 1)]; inner = nbasin) Δtime_seconds = seconds.(diff(data.time)) Δtime = repeat(Δtime_seconds; inner = nbasin) @@ -146,6 +141,9 @@ function basin_table( storage_increase = max(storage_flow, 0.0) storage_decrease = max(-storage_flow, 0.0) + inflow_rate[i] + precipitation[i] + total_in = inflow_rate[i] + precipitation[i] + drainage[i] - storage_increase total_out = outflow_rate[i] + evaporation[i] + infiltration[i] - storage_decrease balance_error[i] = total_in - total_out @@ -224,7 +222,8 @@ function flow_table( from_node_id = repeat(from_node_id; outer = ntsteps) to_node_type = repeat(to_node_type; outer = ntsteps) to_node_id = repeat(to_node_id; outer = ntsteps) - flow_rate = average_over_saveats(saved.flow, :flow, saveats) + flow_rate = + FlatVector([collect(saved_flows.flow.flow) for saved_flows in saved.flow.saveval]) return (; time, diff --git a/core/test/bmi_test.jl b/core/test/bmi_test.jl index 3e4847c27..055562580 100644 --- a/core/test/bmi_test.jl +++ b/core/test/bmi_test.jl @@ -86,6 +86,8 @@ end model = BMI.initialize(Ribasim.Model, toml_path) demand = BMI.get_value_ptr(model, "user_demand.demand") realized = BMI.get_value_ptr(model, "user_demand.realized") + @test realized === + model.integrator.opts.callback.discrete_callbacks[4].affect!.integral.user_realized day = 86400.0 BMI.update_until(model, 0.4day) demand = 0.001 # for both users at the start diff --git a/core/test/utils_test.jl b/core/test/utils_test.jl index 50bfe9076..5eeb8b81f 100644 --- a/core/test/utils_test.jl +++ b/core/test/utils_test.jl @@ -64,7 +64,7 @@ end using Dictionaries: Indices using StructArrays: StructVector using Logging - using Ribasim: NodeID + using Ribasim: EdgeMetadata, NodeID level = [ 0.0, @@ -94,11 +94,8 @@ end demand = zeros(1) basin = Ribasim.Basin( Indices(NodeID.(:Basin, [1])), - [NodeID[]], - [NodeID[]], - zeros(1), - zeros(1), - zeros(1), + [EdgeMetadata[]], + [EdgeMetadata[]], zeros(1), zeros(1), zeros(1),