[bug] Terminals sometimes "leak" water

[evetion]

I hit this during the implementation of #1849, where edges to Terminals have negative flows when checked in update_cumulative_flows! using the flow from flow_update_on_edge. For example in the HWS model I get the following:

┌ Warning: Terminal #5415 has flow of -338.54177272319794 at 6.862248444895375e6 for a timestep of 91.35735275302858
└ @ Ribasim /Users/evetion/code/Ribasim/core/src/callback.jl:205
┌ Warning: Terminal #5415 has flow of -85.64921969175339 at 6.862339802248128e6 for a timestep of 91.35735275302858
└ @ Ribasim /Users/evetion/code/Ribasim/core/src/callback.jl:205
┌ Warning: Terminal #5415 has flow of -8.970432549715042 at 6.862431159600881e6 for a timestep of 91.35735275302858
└ @ Ribasim /Users/evetion/code/Ribasim/core/src/callback.jl:205
┌ Warning: Terminal #5415 has flow of -0.8940373957157135 at 6.862489481604243e6 for a timestep of 58.32200336178765
└ @ Ribasim /Users/evetion/code/Ribasim/core/src/callback.jl:205
┌ Warning: Terminal #5415 has flow of -0.9066413938999176 at 6.862868324295173e6 for a timestep of 378.8426909304897
└ @ Ribasim /Users/evetion/code/Ribasim/core/src/callback.jl:205
┌ Warning: Terminal #5415 has flow of -0.1415431797504425 at 6.863247166986104e6 for a timestep of 378.8426909304897
└ @ Ribasim /Users/evetion/code/Ribasim/core/src/callback.jl:205
┌ Warning: Terminal #5415 has flow of -0.06280365586280823 at 6.8643238840039205e6 for a timestep of 1076.7170178166448
└ @ Ribasim /Users/evetion/code/Ribasim/core/src/callback.jl:205
┌ Warning: Terminal #5415 has flow of -0.009804785251617432 at 6.865400601021737e6 for a timestep of 1076.7170178166448
└ @ Ribasim /Users/evetion/code/Ribasim/core/src/callback.jl:205
┌ Warning: Terminal #5415 has flow of -3.9677118062973022 at 6.882244989985872e6 for a timestep of 0.059130550789141804
└ @ Ribasim /Users/evetion/code/Ribasim/core/src/callback.jl:205
┌ Warning: Terminal #5415 has flow of -0.628738522529602 at 6.882245101221386e6 for a timestep of 0.11123551420152002
└ @ Ribasim /Users/evetion/code/Ribasim/core/src/callback.jl:205
┌ Warning: Terminal #5415 has flow of -1.5169505487529635e6 at 8.982259552166957e6 for a timestep of 3340.4478330416605
└ @ Ribasim /Users/evetion/code/Ribasim/core/src/callback.jl:205
┌ Warning: Terminal #5415 has flow of -2858.7084048986435 at 9.849623705491675e6 for a timestep of 23.70549167560891
└ @ Ribasim /Users/evetion/code/Ribasim/core/src/callback.jl:205
...

This is the Terminal at Herentals in the HWS model with an Outlet(!) before it (the only Terminal one to show such behaviour). I can imagine this happens during smaller timesteps to balance the problem, and you'll never see this on the larger saveat timesteps. Given the rather large flows/timesteps in the above, this should be investigated (could also be a bug in the code to check such flows).

[evetion]

To reproduce, insert the following snippet into update_cumulative_flows or see #1849.

    for (inflow_edge, outflow_edge) in zip(state_inflow_edge, state_outflow_edge)
        from_flow = inflow_edge.edge[1]
        to_flow = outflow_edge.edge[2]
        if from_flow.type == NodeType.Basin
            flow = flow_update_on_edge(integrator, inflow_edge.edge)

[visr]

I cannot get your snippet to work. Can you maybe make such a snippet a PR on main that would fail the HWS integration test by checking for negative flows?

I tried reproducing by running the HWS model with saveat=0:

starttime = 2023-01-20 00:00:00
endtime = 2023-05-30 00:00:00
crs = "EPSG:28992"
input_dir = "."
results_dir = "results"
ribasim_version = "2024.10.0"

[solver]
autodiff = false
saveat=0
water_balance_abstol = 1e30 # optional, default 1e-3
water_balance_reltol = 1e20 # optional, default 1e-2

But I don't find any leaky Terminal nodes:

import ribasim

model = ribasim.Model.read(r"d:\repo\ribasim\Ribasim\models\hws\minio\hws.toml")
terminals = model.node_table().df[model.node_table().df.node_type == "Terminal"].index.to_numpy()

flow = pd.read_feather(
    r"d:\repo\ribasim\Ribasim\models\hws\minio\results\flow.arrow",
    dtype_backend="pyarrow",
)
flow = flow[flow.to_node_id.isin(terminals)]
flow.flow_rate.min()  # 0.0
flow.flow_rate.max()  # 31.31

[evetion]

What about this on main? I'm on HWS version 2024_7_0.

diff --git core/src/callback.jl core/src/callback.jl
index 8990d86f..c3bf1d2c 100644
--- core/src/callback.jl
+++ core/src/callback.jl
@@ -92,7 +92,14 @@ Update with the latest timestep:
 """
 function update_cumulative_flows!(u, t, integrator)::Nothing
     (; p, uprev, tprev, dt) = integrator
-    (; basin, user_demand, flow_boundary, allocation) = p
+    (;
+        basin,
+        user_demand,
+        flow_boundary,
+        allocation,
+        state_inflow_edge,
+        state_outflow_edge,
+    ) = p
     (; vertical_flux_bmi, vertical_flux_from_input) = basin

     # Update tprev
@@ -150,6 +157,25 @@ function update_cumulative_flows!(u, t, integrator)::Nothing
             end
         end
     end
+
+    for (inflow_edge, outflow_edge) in zip(state_inflow_edge, state_outflow_edge)
+        from_flow = inflow_edge.edge[1]
+        to_flow = outflow_edge.edge[2]
+        if from_flow.type == NodeType.Basin
+            flow = flow_update_on_edge(integrator, inflow_edge.edge)
+            if to_flow.type == NodeType.Terminal && flow < 0 && to_flow.value != 0
+                @warn "Unsupported inflow type $(to_flow.type) #$(to_flow.value) with flow $flow"
+            end
+        end
+
+        if to_flow.type == NodeType.Basin
+            flow = flow_update_on_edge(integrator, outflow_edge.edge)
+            if from_flow.type == NodeType.Terminal && flow > 0 && from_flow.value != 0
+                @warn "Unsupported outflow type $(from_flow.type) #$(from_flow.value) with flow $flow"
+            end
+        end
+    end
+
     return nothing
 end

[evetion]

Just a thought, it could be my mac here, producing slightly different stuff than linux/windows. Then again seeing ┌ Warning: Unsupported inflow type Terminal #5415 with flow -2858.7084048986435 is bad.

[visr]

#1856 shows this reproduces on CI and @evetion managed to reproduce locally based on the results as well:

They are very short negative spikes so likely not much impact on the water balance.

[visr]

Earlier I couldn't reproduce just because I didn't run quite long enough to hit the error.
I saw that ImplicitEuler doesn't have this issue, but is taking much much longer, e.g. 30 minutes instead of 10 seconds or so.
It started getting slow in March as the Meuse started surging. Note in the picture above that this Outlet that goes to the Terminal is under DiscreteControl based on Monsin. The HWS still has huge DiscreteControl tables for Rhine and Meuse. This is known to be problematic, and we need to replace these with ContinuousControl instead. I quickly tried removing DiscreteControl altogether, but unfortunately that didn't eliminate the leaky Terminal.