Add water balance output

[visr]

The level and flow output we have right now are samples in a continuous signal.

To get a good sense of the total water balance, and to check if it closes, we need to integrate these over time, to get output in m3 over the past timestep that we can sum.

I think we should implement this using the IntegratingCallback.

Also we should clean up some leftover code from similar old functionality:
https://github.com/Deltares/Ribasim/blob/v2023.09.0/core/src/bmi.jl#L200
https://github.com/Deltares/Ribasim/blob/v2023.09.0/core/src/solve.jl#L1194

EDIT: This issue is now only about water balance calculation and output. Mean flows are #935. This can be done separately, e.g. the water balance below was calculated using sampled flows, ignoring the error for that.

[visr]

The IntegratingCallback doesn't seem to really match our use case, though perhaps some of the code could be used: the Gaussian quadrature that interpolates between the timesteps. The main issue seems to be that it seems to want to integrate over all parameters (integral = allocate_zeros(integrator.p) fails for us). Though in the end these are all DiscreteCallback, so we can build them ourselves.

#644 takes an important first step by making sure the flow output table is complete, by including the gains and losses from non-conservative nodes. Here is a first quick script calculating the water balance error over the basic_transient test model.

using Plots
using DataFrames, Dates
using Ribasim

toml_path = normpath(@__DIR__, "generated_testmodels/basic_transient/basic_transient.toml")
model = Ribasim.run(toml_path)
# use unique since some timesteps are saved multiple times
basin = unique(DataFrame(Ribasim.basin_table(model)), [:time, :node_id])
flow = DataFrame(Ribasim.flow_table(model))
time = unique(Ribasim.datetimes(model))

function water_balance_error(basin, flow, time)
    wbal_errors_m3s = Float64[]
    wbal_errors_m3 = Float64[]
    storages = Float64[]
    for (t1, t2) in zip(time[begin:(end - 1)], time[(begin + 1):end])
        # we have to calculate the storage difference over time ourselves
        dt = 0.001 * Dates.value(t2 - t1)

        s_t1 = sum(filter(:time => t -> t == t1, basin).storage)
        s_t2 = sum(filter(:time => t -> t == t2, basin).storage)
        ds = s_t2 - s_t1
        dsdt = ds / dt

        # for the period (t1,t2) flow is defined on t2
        bc_flows = filter([:time, :edge_id] => (t, e) -> t == t2 && ismissing(e), flow)
        bc_flow = sum(bc_flows.flow)

        wbal_error_m3s = dsdt - bc_flow
        wbal_error_m3 = wbal_error_m3s * dt
        push!(wbal_errors_m3s, wbal_error_m3s)
        push!(wbal_errors_m3, wbal_error_m3)
        push!(storages, s_t2)
    end
    return (; wbal_errors_m3s, wbal_errors_m3, storages)
end

(; wbal_errors_m3s, wbal_errors_m3, storages) = water_balance_error(basin, flow, time)
plot(time[2:end], wbal_errors_m3s; title = "error [m3/s]")
plot(time[2:end], cumsum(wbal_errors_m3); title = "cumulative error [m3]")
plot(time[2:end], storages; title = "total storage [m3]")

The error is not huge, but seems high. Perhaps the default tolerances should be a bit lower. If I reduce them both 1000x, the cumulative error at the end becomes ~1.2 m3 instead of ~120 m3.

We should probably integrate some form of this into Ribasim, to have the error as part of our output. This just does dS/dt = I - O over the whole model. It would be nice to split the terms per node type for more insight, e.g. in the style of MODFLOW / Ribasim 7.

     CUMULATIVE VOLUME      L**3       RATES FOR THIS TIME STEP      L**3/T          PACKAGE NAME
     ------------------                 ------------------------                     ----------------

           IN:                                      IN:
           ---                                      ---
              STO-SS =      136438.0189                STO-SS =          10.7687     STORAGE
              STO-SY =      182397.4547                STO-SY =        6424.1892     STORAGE
                 WEL =           0.0000                   WEL =           0.0000     WEL
                 RIV =      929004.1292                   RIV =       18699.1633     RIV
                 GHB =           0.0000                   GHB =           0.0000     GHB-TIDAL
                 RCH =     1055750.0000                   RCH =       51250.0000     RCH-ZONE_1
                 RCH =      218400.0000                   RCH =        6400.0000     RCH-ZONE_2
                 RCH =      737100.0000                   RCH =       23400.0000     RCH-ZONE_3
                 EVT =           0.0000                   EVT =           0.0000     EVT

            TOTAL IN =     3259089.6027              TOTAL IN =      106184.1213

          OUT:                                     OUT:
          ----                                     ----
              STO-SS =      137364.2498                STO-SS =       11404.8487     STORAGE
              STO-SY =      641573.7206                STO-SY =       23973.7320     STORAGE
                 WEL =       32200.0000                   WEL =         670.0000     WEL
                 RIV =      120735.7289                   RIV =        8577.0472     RIV
                 GHB =     2319223.3807                   GHB =       61295.9057     GHB-TIDAL
                 RCH =           0.0000                   RCH =           0.0000     RCH-ZONE_1
                 RCH =           0.0000                   RCH =           0.0000     RCH-ZONE_2
                 RCH =           0.0000                   RCH =           0.0000     RCH-ZONE_3
                 EVT =        7992.5228                   EVT =         262.5877     EVT

           TOTAL OUT =     3259089.6027             TOTAL OUT =      106184.1213

            IN - OUT =      -6.7656E-06              IN - OUT =      -8.0561E-07

 PERCENT DISCREPANCY =          -0.00     PERCENT DISCREPANCY =          -0.00

[visr]

First step: add the percentage error to stdout.

[visr]

Fixed by #1367