Convergence of large outflows from small Basins

[visr]

If you have a small Basin with e.g. a large infiltration, the model is not very stable. We need to build a test model for this case to make sure it is able to handle such a case well. In principle we use a reduction_factor to reduce the infiltration, but for sufficiently large outflows it will step over this smooth landing within one timestep. We need to look into this case to understand the behavior, and how it behaves with both adaptive and fixed timestepping, under different solvers. Perhaps we need to use better tolerances or error norms to handle these cases better.

Related to #1406.

[SouthEndMusic]

This PR proposes a method to handle basins drying out using VectorContinuousCallback, which we should test with a large model containing basins that are drying out: https://github.com/Deltares/Ribasim/issues/1575

[evetion]

Commenting here instead of on the PR, it's good to read https://docs.sciml.ai/DiffEqDocs/stable/basics/faq/#My-ODE-goes-negative-but-should-stay-positive,-what-tools-can-help? and related the PositiveDomain callback, which specifically states its a more performant alternative to isoutofdomain.

However while isoutofdomain is "blind" (requiring many iterations before arriving at a very small timestep), using a callback might not solve it, as I can imagine rootfinding is difficult when we have a incredibly smooth reduction factor applied to q (from some arbitrary threshold onwards). This leads to the question best answered by @visr, why/when did we introduce the reduction factor (for basins) and more generally can we do without (and clamp or something when under a certain threshold/level)?

[evetion]

PS. When someone implements a quick (e.g. GPU) D-Flow FM alternative, there's someone who always asks (not me), but can it handle "droogval" (waterlevel becoming zero)?! It's a harder and more common problem than you'd think.

[visr]

In #639 I commented about PositiveDomain:

One alternative is to use the PositiveDomain callback, but this only works for adaptive=false, and is not easily adapted to enforce only the storage part of the state to be positive.

Though I think also isoutofdomain only really makes sense for adaptive timesteps, since you can reject those. We had to implement an erroring fallback in case we still end up with negative storages in #1406.

All reduction factors that we have are non-physical and only meant for performance. So I'd be happy to remove any of them in favor of another approach as long as it leads to better convergence for large models. They are mainly meant as a way to smoothen discontinuities as suggested in https://discourse.julialang.org/t/handling-instability-when-solving-ode-problems/9019/5.

[Huite]

How is implicit versus explicit behaving? ImplicitEuler might not converge, but it does remain stable until non-convergence, right?

The NLNewton solver allows for relaxation via the relax keyword parameter. Numerical values of relax must lie in the half open unit interval [0,1) where 0 corresponds to no relaxation. Alternatively, relax may be set to a line search algorithm from LineSearches.jl in order to include a line search relaxation step in the Newton iterations. For example, the well known Newton-Armijo iterative scheme can be employed by setting relax=BackTracking() where BackTracking is provided by LineSearches.

https://docs.sciml.ai/DiffEqDocs/stable/features/linear_nonlinear/#Nonlinear-Solvers:-nlsolve-Specification

We haven't enabled this yet, right? Backtracking should really help.

I think the nasty thing about thresholds is that they're essentially arbitrary. Basins may be huge or tiny, so volumes don't cut it; water levels seem reasonable; but there are probably cases where 1 cm is negligible, and cases where 1 cm is very important.

At the end of formulation of rhs, we have du right? In principle du/u tells you something about the rate of change relative to the storage. Do we have acces to dt somewhere? In principle, we could identify "panicking" basins, where the (negative) du * dt exceeds u, or when it nearly exceeds u. You can then dynamically scale down du; better is to let this be a criterion in time step selection.

[SouthEndMusic]

Here's a distillation of a discussion on slack (with @DaniGlez):

Let basins be in either a 'depleted' or 'non-depleted' state. When a basin is non-depleted, it handles inflows and outflows as normal. When it is depleted, it distributes incoming flow proportionally over the outgoing flows (and we can enforce du = 0 to make sure the water balance is closed).

A basin goes into depleted state when u = 0 (or u = eps) which is detected by VectorContinuousCallback. It goes into non-depleted state when du is positive (with another VectorContinuousCallback? Not completely sure about this). An advantage of this approach is that no thresholds have to be chosen.

This can in some situations be less of a shock to the solver than enforcing that there is no outflow at all when the basin is depleted. Plus, the used callback types are especially designed to handle discontinuities in the system.

@Huite we indeed haven't played around with the settings of NLNewton yet. Also if I understand you correctly, you want to enforce that the courant number is less than 1 by scaling down du. I would however argue for the sake of explainability that we modify individual flows (as we also have these as output). I think we can quite easily dig dt out of the integrator, but we have to be careful about at which point of the adaptive time stepping we capture it. Backtracking looks interesting and easy enough to try so I'll do that first.

In response to https://docs.sciml.ai/DiffEqDocs/stable/basics/faq/#My-ODE-goes-negative-but-should-stay-positive,-what-tools-can-help, and as pointed out by @DaniGlez, non-negative storages are not automatically a result of our ODE problem, we have to modify the problem to achieve it (as a 'modelling strategy'). I was trying to get at that if you enforce in the solver that storages are non-negative but you don't modify the ODE problem to reflect that, you have an ill-posed algorithm.

[DaniGlez]

Let me just add that 1) I'm not a hydrologist and 2) other distribution rules other than "all outflows are reduced in the same proportion" might be possible in principle, but if you have cyclic graphs then I imagine there could be issues with unique flow determination under other rules.

[SouthEndMusic]

Completed by https://github.com/Deltares/Ribasim/pull/1697