cumulative flow states

[SouthEndMusic]

Fixes #1808

[SouthEndMusic]

Finally all tests are green again after the refactor. It was kind of a pain to make all the book-keeping work out (quite a bit changed in that regard).

Findings

These findings are based on running the HWS model as taken from the regression tests.

Water balance

Let's start with the best result: for the HWS model the largest absolute water balance error is 8.83354278613524e-10.

This means that per basin per day, less than a μL is not accounted for.

Note that this balance error is independent of the solver tolerances, the water balance is inherently closed due to the problem formulation. The errors left are just (accumulated) floating point errors.

Profile

The refactor didn't introduce any glaring performance problems in the code:

States and Jacobian sparsity

The state vector is still a component vector, with these components:

julia> only(getfield(model.integrator.u, :axes))
Axis(tabulated_rating_curve = 1:11, pump = 12:40, outlet = 41:101, user_demand_inflow = 102:127, user_demand_outflow = 128:153, linear_resistance = 154:153, manning_resistance = 154:239, evaporation = 240:376, infiltration = 377:513, integral = 514:522)

i.e. sorted by node type. This means that we can find convergence bottlenecks in a more fine-grained way than before.

Considerations for the states:

• Each conservative node and not exactly integrable boundary node has a state. So FlowBoundary doesn't need a state, and UserDemand needs states for inflow and outflow separately due to the time-dependent return factor.
• Also precipitation and drainage don't need their own states, as these are exactly integrable (precipitation uses a fixed area).
• evaporation and infiltration now have separate states, a bit more on that later.

The Jacobian sparsity looks like this:

Which admittedly looks quite cool, but there are now many more states (going from 137 to 522) and more non-zeros in the Jacobian. To give an indication on how the Jacobian sparsity comes to be: say we have the following model section.

The state for the manning flow depends on the forcing states of the connected basins, but also the flow states of all nodes connected to these basins. That is because the manning flow depends on the basin levels, which depend on the basin storages, which depend on the cumulative flows over all flow edges connected to the basins.

Performance

As expected, the performance takes a hit. The runtime of the HWS model:

• On main: 104.484030 seconds (63.46 M allocations: 10.354 GiB, 0.65% gc time, 0.01% compilation time) (I thought the HWS model was much faster on main?)
• On this branch: 266.583831 seconds (441.16 M allocations: 104.220 GiB, 1.53% gc time)

I had hoped I could reduce the number of states by clumping evaporation and infliltration together into a single state. The problem with that is that the contribution of these forcings over a timestep can then not exactly be split into an evaporation and an infiltration contribution. The contribution of this state of combined forcings over a timestep $\Delta t = t_1 - t_0$ would be:

$$ \Delta V = \text{infiltration} \cdot \int_{t_0}^{t_1} \phi\left(d(u(t')); 0.1\right) \text{d}t' + \text{evaporation} \cdot \int_{t_0}^{t_1} A\left(u(t')\right)\phi\left(d(u(t');0.1\right)\text{d}t'. $$

where $\phi$ is the reduction factor. Both integrals are unknowns (except when the basin is not in the 'reduction factor regime'), so it cannot be computed exactly how $\Delta V$ is split up into the infiltration and evaporation contributions. If we find the performance more important than this exact split, we could approximate it. Note that this approximation does not affect the water balance, because we still know exactly how much water is leaving the basin. Maybe that is not so bad as the reduction factor is a coarse approximation of what happens when a basin dries out anyway.

[visr]

I made a plot of the "Test Delwaq coupling" artifact "delwaq_map.nc" variable "ribasim_Continuity", which should be as close to 1 as possible, comparing this branch (new) to the current main (old).

So that is very encouraging. (note the +1 on the y axis, so 0 is 1)

import pandas as pd
import xarray as xr

ds_old = xr.open_dataset(r"main\delwaq_map.nc")
da_old = ds_old["ribasim_Continuity"]

ds_new = xr.open_dataset(r"integrated-flow-states\delwaq_map.nc")
da_new = ds_new["ribasim_Continuity"]

df_old = da_old.sel(ribasim_nNodes=3).to_dataframe()["ribasim_Continuity"].rename("old")
df_new = da_new.sel(ribasim_nNodes=3).to_dataframe()["ribasim_Continuity"].rename("new")

df = pd.concat([df_old, df_new], axis=1)
df.plot()

[Huite]

Looks very promising! The robustness + simplicity with regards to the water balances makes me think this is well worth it.

Which admittedly looks quite cool, but there are now many more states (going from 137 to 522) and more non-zeros in the Jacobian.

You can reduce bandwidth with algorithms like: https://en.wikipedia.org/wiki/Cuthill%E2%80%93McKee_algorithm

I'm not sure this is relevant: maybe the solver machinery already does things like this. Second, given the model sizes, a lot of data may fit wholly within CPU caches, so better bandwidth may have no meaningful effects on memory access patterns.

[SouthEndMusic]

You can reduce bandwidth with algorithms like: https://en.wikipedia.org/wiki/Cuthill%E2%80%93McKee_algorithm

@Huite I have dabbled in that before. As far as I'm aware, the goal of reducing the bandwidth is reducing fill-in in the Jacobian factorization. I looked a little bit in what OrdindaryDiffEq.jl does there and the factors I found were quite sparse, so there's probably already some hidden cleverness there

EDIT: This is already part of the matrix factorization: https://dl.acm.org/doi/10.1145/992200.992206

[SouthEndMusic]

I fixed the flaky tests, turns out there was a dictionary involved in setting up the edge data in the flow table which messed up the ordering.