A systematic look at interpolation types

[SouthEndMusic]

An overview of considered interpolation types:

Interpolation type	goes through data points	differentiable in the data points	other pros	other cons
Block interpolation	yes	no	area under the curve easy to compute	-
Smoothed block interpolation	no	yes	same area under the curve as block interpolation	the area under the curve does not exactly match in the data points
Linear interpolation	yes	no	has analytically invertible antiderivative	-
Smoothed linear interpolation	no (close)	yes	never overshoots the data	-
PCHIP/Akima interpolation	yes	yes	-	can overshoot the data

Are there other criteria we or modelers care about?

Time interpolations

Interpolations of timeseries do not have to be continuous or differentiable in the time points because of the tstops mechanism. However, smoother interpolations are easier on the solver. For forcings we especially care about the area under the curve because this dictates the volume that the model interchanges with its environment. This is why the smoothed block interpolation was designed in such a way that it has the same area under the curve as block interpolation, be it with a slight delay.

Interpolations that take state dependent input

Interpolations that are part of the RHS of our ODE system and take state dependent input (like storages or levels) affect the Jacobian of our RHS, and thus the convergence of our non-linear solver. That the smoothness of these interpolations matters for performance was shown for example in #2446.

Type stable configurable interpolation types in the core

Currently the type of interpolations is fixed by the type in the corresponding field in one of the AbstractParameterNode object (with a recent change where FlowBoundary has a type parameter to support multiple interpolation types). We could however refactor this data structure to make things more statically typed and flexible. I propose to add a interpolations field to ParametersIndependent of type

@kwdef struct Interpolations
   linear::Vector{LinearInterpolation} = []
   ...
end

and the fields in the AbstractParameterNode objects refer to this with a Tuple{InterpolationType.T, Int}, which is similar to what is already done for TabulatedRatingCurve. Then we evaluate with

function evaluate_interpolation(interpolations::Interpolations, input::T, index::Tuple{InterpolationType.T, Int})::T
    return if index[1] == InterpolationType.linear
      interpolations.linear[index[2]](input)
   elseif ...
   end
end

This also paves the way for node-specific interpolation types were we to go that route (and gives a place to put timeseries we only use as a control input, I couldn't find an existing issue for that).

Related issues

#1249
#1919
#1920