Smooth clamp function

[Huite]

I've noticed clamp is used four times in solve.jl. The problem with clamp is that it is non-smooth so it won't play nice with anything that utilizes derivatives, like Newton Raphson.

A straightforward clamp_smooth is:

function clamp_smooth(x, lo, hi, m)
    d = (hi - lo)
    @assert m <= (0.5 * d)
    x_norm = (x - lo) / d
    a = 1 / (1 - m)
    
    if x_norm < 0
        y = 0
    elseif x_norm < m
        y = (a / (2 * m)) * x_norm^2
    elseif x_norm < (1 - m)
        y = a * x_norm + 0.5 * (1 - a)
    elseif x_norm < 1
        y = 1 - (a / (2 * m)) * (1 - x_norm)^2
    else
        y = 1
    end
    
    return lo + y * d
end

Where m is the smoothing parameter. This is only C1 continuous which is fine for us now but maybe needs some attention if we need Hessians in the future.

In terms of its behavior, this is with a large m (clamp in orange):

Applied to a pump, it means that the pump is a bit "lazy", then "panicky": it pumps slightly less when at half-capacity, and it pumps a little more when over half-capacity. In general, the value for m is chosen very small, so deviations are negligible.

The reduction_factor is similar in idea:

Ribasim/core/src/util.jl

Line 367 in 180f283

function reduction_factor(x::T, threshold::Real)::T where {T <: Real}

It uses a cubic approximation to go from 0 to 1 from 0 to threshold (but it doesn't approximate a piecewise linear function like clamp well).

Although it's more costly to evaluate than a basic clamp, it should be helpful for the solver and give some performance benefits.

[evetion]

it should be helpful for the solver and give some performance benefits.

Lots of maybe here, is there a way to actually measure this (apart from implementing it?). There's always a performance need (Ribasim-NL will only get bigger), but ideally we could order ideas based on their expected benefits.

[Huite]

Lots of maybe here, is there a way to actually measure this (apart from implementing it?)

Not really. It makes the problem better behaved, so the solver will require less iterations. Our runtimes are basically linear with the number of iterations.

We don't have an a priori way of calculating how many iterations are required.

we could order ideas based on their expected benefits.

Sure, but if you make it a cost-benefit, then the cost in this case is pretty small; if the resulting benefits are small but not tiny, then it's not too bad. In this case the cost is small because it's a single function replacement. Some deviation from the 1:1 line is unavoidable if you want preserve the clamping boundaries. But the smoothing parameter directly controls the deviation. A (normalized) value of 0.001 is probably already quite effective, and really barely influences flows.

The biggest time sink is the benchmarking: you might not see the benefit unless the model has sufficient Outlets hitting their maximum flow rate or pumps interacting with ContinuousControl. So the smart thing is likely to prioritize other things until our model zoo has grown a bit; then benchmarking is very simple and benefits should be obvious (or not).

Secondly, from a philosophical point of view: it's almost certainly a good idea to "make the whole RHS of the Ribasim ODE problem C1 smooth" (see also #1919).