Take discontinuities in flow introduced through BMI into account

core/src/callback.jl

@@ -139,7 +139,6 @@ function integrate_flows!(u, t, integrator)::Nothing
                 0.5 * (get_flow(graph, edge..., 0) + get_flow_prev(graph, edge..., 0)) * dt
         end
     end
-
     copyto!(flow_prev, flow)
     copyto!(vertical_flux_prev, vertical_flux)
     return nothing

core/src/model.jl

@@ -213,6 +213,7 @@ function SciMLBase.step!(model::Model, dt::Float64)::Model
     if round(ntimes) ≈ ntimes
         update_allocation!(integrator)
     end
+    set_previous_flows!(integrator)
     step!(integrator, dt, true)
     return model
 end

core/src/util.jl

@@ -749,3 +749,20 @@ end
 function basin_areas(basin::Basin, state_idx::Int)
     return basin.level_to_area[state_idx].u
 end
+
+"""
+Just before setting a timestep, call water_balance! again
+to get a correct value of the flows for integrating
+"""
+function set_previous_flows!(integrator)
+    (; p, u, t) = integrator
+    (; flow, flow_prev) = p.graph[]
+    (; vertical_flux, vertical_flux_prev) = p.basin
+    du = get_du(integrator)
+    water_balance!(du, u, p, t)
+
+    flow = get_tmp(flow, 0)
+    vertical_flux = get_tmp(vertical_flux, 0)
+    copyto!(flow_prev, flow)
+    copyto!(vertical_flux_prev, vertical_flux)
+end

core/test/bmi_test.jl

@@ -83,15 +83,21 @@ end
     toml_path =
         normpath(@__DIR__, "../../generated_testmodels/minimal_subnetwork/ribasim.toml")
     @test ispath(toml_path)
-    model = BMI.initialize(Ribasim.Model, toml_path)
+    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")
+    # 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, 0.4day)
-    demand = 0.001  # for both users at the start
-    @test all(isapprox.(realized, demand * 0.4day; rtol = 1e-3))
-    BMI.update_until(model, 0.6day)
-    @test all(isapprox.(realized, demand * 0.6day; rtol = 1e-3))
+    BMI.update_until(model, day)
+    @test realized ≈ [demand_start * day, demand_start * day + 0.5 * slope * day^2]
+    demand_later = 2e-3
+    demand[1] = demand_later
+    BMI.update_until(model, 2day)
+    @test realized[1] == demand_start * day + demand_later * day
 end
 
 @testitem "vertical basin flux" begin