Taming cumulative flows

[visr]

Since #1819 our states are cumulative flows, getting further from 0 all the time. This formulation helped the water balance error and mean flows a lot, but came at the cost of errors occuring after months of simulations. Decreasing the reltol in #1874 helped reduce the issues that came to light, though as mentioned in #1874 (comment) it still suffers from reduced effective tolerance over time.

One option to address this is by resetting cumulative flows. A draft function to reset is in #1864, though the difficult part is coming up with a good approach for when to call it, such that it works independently of saveat settings and BMI usage, see #1874 (comment).

We can also look into using custom error norm functions that are more robust to these cumulative flows, somehow keeping the effective errors constant over time.

#1874 was only tested manually, so as a part of this issue we should include tests to show the reduction in error over time. The manual checks consisted of running HWS with #1856, and the gist linked in #1863 (comment).

[Huite]

Another option to consider is to solve for flows AND storage at the same time, and add them both to the unknowns (maybe in ComponentArray form?). My feeling is that the overhead is negligible and it would likely simplify the implementation -- you can assume there's an up-to-date storage available at any time.

It would also bring us more in line with something like a 1D St. Venant simulator, since as far as I can tell, these solve for level (∝ storage) and flow velocity (∝ flow) at the same time too.

This in turn might allow a natural extension for certain basins / connectors (or parts of the network) that do incorporate a momentum balance.

[SouthEndMusic]

I did an experiment with both storage and flow states before, then it was quite detrimental to performance. That was quite a while ago though

[SouthEndMusic]

Regarding the custom norm: we should look into how the error is estimated again:

It's never clear to me from the docs what the order of operations is, I.e. what is per element of the state and when norms are applied. This formula at least indicates a decreasing error with an increasing state value because of the reltol term, maybe for us it makes more sense if that max is a difference.

[evetion]

I believe the implementation resides here https://github.com/SciML/DiffEqBase.jl/blob/01fd8532d8f0aac2c30995c2a74a29ef0a173edf/src/calculate_residuals.jl#L36. Even when you override the norm with internalnorm to init (take care to provide two methods, one vector, one scalar), with one that becomes less big (did some tests with log2), you'll still end up with an increasing error, just slower. If you fix the norm completely, you're essentially setting an (convoluted) absolute tolerance. Doing that (I believe @SouthEndMusic had similar thoughts) and going back to abstol 1e-5 and setting reltol to zero also works fine (below, but requires more testing accross the ecosystem) and might be more predictable?