Improve continuous control performance

[SouthEndMusic]

Fixes #2445.

The model repro-models/HollandseDelta_parameterized makes heavy use of ContinuousControl. I had the hypothesis that the fact that the functions that define the continuous relations for ContinuousControl are piecewise linear is hard on the solver due to the lack of smoothness. To test this I changed the interpolation type to PCHIPInterpolation (see here), not as a change to merge, but just to see what a smoother interpolation type does. And indeed, with this smoother interpolation I can locally run the first 4 months of the simulation in 9 min.

There are several things we can do to improve the ContinuousControl performance structurally:

• Indeed use a smoother interpolation type than LinearInterpolation. We of course have to use an interpolation type that can be understood by modelers, so this could be a good use case for smoothed linear interpolation
• Treat ContinuousControl more similar to DiscreteControl, in the sense that the continuous control logic is not part of the ODE system, but the control logic is computed at the start of the timestep and then kept fixed over the timestep. This approach probably has the best performance gain, but also results in a loss of accuracy. Alternatively we could make the continuous-but-not-part-of-the-ODE-system-control a new node type (or add a switch in the TOML between callback computed or ODE computed), but then we need to make the subtle difference very clear in the docs.

@visr @gijsber @Huite let me know whether you think this warrants a design session. Maybe in the particular case of the HollandseDelta model the schematization can also be made a bit more solver friendly

As a side node, I think we can expect similar issues with models that have a lot of TabulatedRatingCurve nodes with very sharp transitions in the q(h) relations

[SouthEndMusic]

This model also contains a lot of Manning resistance nodes, so I tried some things there as well without much succes. One thing I tried was relaxing the Manning resistance flow, i.e. computing the flow for a new timestep with a fraction of the flow of the previous timestep.

In general I think that the Manning resistance nodes still cause stability problems because around a head difference of 0 the flow is very sensitive to level fluctuations. I discussed this with Claude and one of the things it suggested that I thought was interesting is introducing an artificial viscosity term to dampen oscillations induced by Manning resistance nodes. I don't yet fully understand how to implement that, @Huite do you have thoughts on this?

[gijsber]

Hi @SouthEndMusic thanks for looking into this.
While I do not have a real objection against fixing the continuous control in a times tep, I would start with the smoothed linear interpolation.
If we then come to the conclusion that performance we still want a change (trading accuracy against performance) we can keep this smooth interpolation while going for your second suggestion.

[visr]

Yeah if we can find a smooth interpolation for ContinuousControl / TabulatedRatingCurve that helps performance and is still quite close to the linear interpolation we could just enable that always.

We could try different interpolation functions and smoothness settings and pick one.

[Huite]

Indeed use a smoother interpolation type than LinearInterpolation. We of course have to use an interpolation type that can be understood by modelers, so this could be a good use case for smoothed linear interpolation

Yeah if we can find a smooth interpolation for ContinuousControl / TabulatedRatingCurve that helps performance and is still quite close to the linear interpolation we could just enable that always.

You could consider relaxing the "close to linear interpolation" requirement since reality is generally smooth rather than piecewise linear. The reason we picked linear interpolation -- I think -- is because it's conceptually simple to implement and understand for a modeler, it obviously has no basis in reality. If e.g. a spline (or whatever) is numerically advantageous, I think most would be happy with it as long as it matches the control points exactly and doesn't look too bendy.

Splines are a somewhat common approach in unsaturated zone modeling for the constitutive functions. Something like Mualem-Van Genuchten needs to be computed at every formulation to compute conductivity and moisture content from pressure; it's cheaper to approximate with a spline (and it'll support all kind of tabulated/empirical constitutive functions as well).

In general I think that the Manning resistance nodes still cause stability problems because around a head difference of 0 the flow is very sensitive to level fluctuations. I discussed this with Claude and one of the things it suggested that I thought was interesting is introducing an artificial viscosity term to dampen oscillations induced by Manning resistance nodes. I don't yet fully understand how to implement that, @Huite do you have thoughts on this?

I've heard of artificial viscosity, but our formulation doesn't consider viscosity explicitly. So you have to alter some values/parameters and my guess is that it's effectively a similar approach as the relaxed root. I guess you'd increase viscosity around 0, thereby decreasing flows in the vincinity? You still have to go through 0 though and meet the square root further away.

What's causing the instability exactly? My guess is that it's reversing the direction of the flow, not so much the steep derative?
I'd guess that a one-directional Manning resistance is more stable, but that might be worthwhile to check with an ad hoc experiment.

[SouthEndMusic]

I think most would be happy with it as long as it matches the control points exactly and doesn't look too bendy.

That brings us back to discussion about how we let the modeler see exactly how we interpolate their data. For block and linear that is indeed easily understood, but for more smooth interpolation types it is not and cannot be returned as a finite amount of points (unless you sample a lot). Then the question is whether documenting the interpolation type properly is enough, or whether we should add some plotting functionality to Ribasim Python.

What's causing the instability exactly? My guess is that it's reversing the direction of the flow, not so much the steep derative?
I'd guess that a one-directional Manning resistance is more stable, but that might be worthwhile to check with an ad hoc experiment.

My hunch is that it is the steep derivative around zero. That makes it hard for the solver to find the equilibrium of the model (as is the case from what I understood from recent Ribasim NL results). It might be somewhat similar to the reduction factor: when we made that less steep by applying it over a much larger storage interval it gave a large performance boost.

I'd guess that a one-directional Manning resistance is more stable, but that might be worthwhile to check with an ad hoc experiment.

It is not immediately clear to me how this would be formulated in a way that is both smooth in the flow direction reversal and not a departure from the current formulation in a neighborhood around that point.

[Huite]

It is not immediately clear to me how this would be formulated in a way that is both smooth in the flow direction reversal and not a departure from the current formulation in a neighborhood around that point.

What I mean with ad hoc is: hopefully there's a simple model which triggers the instability. Then just "break" the functionality it and allow flow in only one direction. The steep derivative remains, so then it's the obvious culprit. But on reflection, it does seem more obvious that it's the steepness rather than the flow direction switching. My thinking was that that maybe the flow reversal causes "bouncy" behavior and instability through the network thereby causing the Newton iteration to struggle; but that doesn't really sound very parsimonious when there's also a singularity to blame.

On that topic: singularities are unphysical. So maybe we should consider a conceptual model that has less singular behavior.

[Huite]

Continuing my last point: Manning and Chezy assume turbulent flow. This is the wrong type of physics for near zero flows; at some point we get laminar flow. If you were to approximate with Poisseuille flow, it's numerically excellent: fully linear!

What we could do is compute Reynolds number, and choose the appropriate approximation. ChatGPT-o3 suggests:

• Laminar: Re < 500 or depth < 1 cm in small channels: Poisseuille / Darcy-Weisbach
• Transitional: 500 <= Re <= 4000: Churchill or Swamee-Jain
• Fully turbulent: Re > 4000: existing Manning (or Colebrook-White)

You can use e.g. a continuous weighting function to keep the function smooth.

In retrospect, I obviously should've done this from the start! This is superior in terms of physical representation and numerical characteristics!

[evetion]

Claude Sonnet Thinking 3.7 when asked a similar question, with links to literature

The Manning equation is widely used for modeling flow resistance in hydrological models, but it has important limitations when applied to small, non-turbulent flows.

The standard Manning equation was developed for fully turbulent flows in rough channels. For small flows where Reynolds numbers are low (typically Re < 2000), the flow regime shifts from turbulent to transitional or laminar, and Manning's resistance formulation becomes physically inappropriate.

Several researchers have investigated this problem:

1. Ferguson, R. (2010). "Time to abandon the Manning equation?" Earth Surface Processes and Landforms, 35(15), 1873-1876. https://doi.org/10.1002/esp.2091

   • This paper critically examines Manning's equation limitations and suggests alternative approaches.

2. Chanson, H. (2004). "The Hydraulics of Open Channel Flow: An Introduction" discusses transitional flows and the breakdown of standard resistance equations. https://doi.org/10.1016/B978-0-7506-5978-9.X5000-4

3. Lawrence, D.S.L. (1997). "Macroscale surface roughness and frictional resistance in overland flow." Earth Surface Processes and Landforms, 22(4), 365-382. https://doi.org/10.1002/(SICI)1096-9837(199704)22:4<365::AID-ESP695>3.0.CO;2-6

   • Examines resistance in shallow overland flows.

For proper handling of small, non-turbulent flows in a ManningResistance node, consider:

1. Reynolds Number-Based Transition: Implement a smooth transition between resistance laws based on the Reynolds number. As suggested by Cheng (2008) in "Formulas for friction factor in transitional regimes" https://doi.org/10.1061/(ASCE)0733-9429(2008)134:9(1357)

2. Alternative Resistance Equations: For very low flows, use Hagen-Poiseuille for laminar flow resistance instead of Manning's equation.

3. Darcy-Weisbach Formulation: This approach can handle both turbulent and laminar regimes more consistently than Manning's.

4. Minimum Flow Threshold: Implement a minimum flow threshold below which alternative formulations are used.

The most theoretically sound approach is to use a composite resistance law that handles the full spectrum of flow conditions, as described by the VPE (Variable Power Equation) in Ferguson's 2007 paper "Flow resistance equations for gravel- and boulder-bed streams" https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2006wr005422.

[SouthEndMusic]

I now use SmoothedLinearInterpolation, see SciML/DataInterpolations.jl#441. Here's an example of the first ContinuousControl function of the Hollandse Delta model:

[evetion]

I now use SmoothedLinearInterpolation, see SciML/DataInterpolations.jl#441. Here's an example of the first ContinuousControl function of the Hollandse Delta model:

I understood this was even faster, now running a significant portion of HD in minutes? I'd say how CC works internally is non-breaking, so I don't mind too much that it doesn't hit it's data points exactly. That said, it would be good to document any large changes in behavior.

[visr]

I think this is very reasonable, and great that it provides such a speedup. We can always revisit interpolation choices in #2452.

Should we directly do the same for TabulatedRatingCurve?

[Huite]

Should we directly do the same for TabulatedRatingCurve?

Yes, I think that makes a lot of sense!