Farewell, cumulative flow states

[verheem]

Fixes #3037

convergence was removed from flow.nc, only basin.nc now has convergence (breaking change)

Cosmetic storage removal was done in commit 351936b.

[visr]

I just tried a binary from this branch on a large part of the LHM model, and whereas the release finished in 2 hours, I killed this after a night, with an ETA of around 2 days. We'd need to profile to see if anything stands out.

[verheem]

I just tried a binary from this branch on a large part of the LHM model, and whereas the release finished in 2 hours, I killed this after a night, with an ETA of around 2 days. We'd need to profile to see if anything stands out.

Probably due to the level difference threshold of 1e-4. This also made HWS extremely slow. Lets try again with 2cm

[SouthEndMusic]

I realized a few remaining problems:

• The flow integration doesn't work with control by allocation, because the changes in flow due to allocation control are not taken into account in the $f_0$ term of the integration scheme. As a quick fix, at the end of update_allocation!, call water_balance and update flow_rate_prev
• The way flows are cached I think is a bit messy; caching them in p_independent in water_balance! and then copying them to state_and_time_dependent_cache in the callbacks; that is not what the state_and_time_dependent_cache is for. I propose:
  • Introducing a CVector for the different types of flows in state_and_time_dependent_cache and caching flows there when formulating;
  • formulating du at the end of water_balance! based on that vector
• I don't know whether we also did this before, but I think that doing DiscreteControl based on the latest computed instantaneous flows isn't great, because that's highly dependent on with which arguments water_balance! happened to be called last. It would be much better to do DiscreteControl based on the average flow of the last timestep.

I think that with #3113 these things can be handled more cleanly.

[visr]

Thanks for pushing on this. If you think these remaining points are better addressed in or after #3113, that is fine with me.

[SouthEndMusic]

@visr I prefer to do it now, to have a proper state of the main branch. I also don't know when #3113 will be planned, there's some time in the current sprint to work on this.

[SouthEndMusic]

Some more notes:

• I'm refactoring some things around multiple inflows for UserDemand nodes. That doesn't work properly at the moment because the individual inflows aren't stored. At the flow integration equal distribution is assumed
• The flow correction in the flow integration I think is not needed. If I understand correctly, this is about taking into account that we can integrate flow boundary inflow and infiltration exactly. The reason we integrated as many flows exactly as possible in the flow formulation is that it reduced the number of states, but with the storage formulation (and with Combined storage and cumulative flow states: a highly efficient solver algorithm independent approach #3113) that is no longer a concern. The integration scheme in this PR is probably accurate enough. I also discuss whether we should use exact integration or not in Combined storage and cumulative flow states: a highly efficient solver algorithm independent approach #3113 under the water balance error section.

The multiple inflows for UserDemand feature actually adds a bit of complexity for #3113 because with that feature some flows do depend on more than (at most) an upstream and downstream storage.

[SouthEndMusic]

I refactored a lot, and I was able to bring the number of failing tests down to a small amount. Still, I am not happy with this direction. Fundamentally we are back at the same problem we had before the cumulative flow states; the post-solver flow quadrature via Simpson's rule doesn't close the water balance sufficiently (at least to avoid negative tracer masses), and any subsequent step to correct that is either rather expensive (i.e. find a minimal correction of the flows that respects things like non-negativity constraints and return factors) or doesn't take parts of the physics into account (e.g. the correction invents non-existent infiltration to close the water balance).

#3113 really is the way to go; we should not have to guesstimate what flows the solver (implicitly) used internally.

I think that in the new formulation all Basin forcings should be treated as flow states again. I don't think it would significantly drop performance with the custom Jacobian and linear solve (which was the original reason to do exact integration), it cleans up the code, and most importantly it gets rid of all remaining headaches of matching exact integration with the solver results. If we really care about matching the exact integration we can inlcude it in isoutofdomain.