Spline extrapolation on the right side fails
Bug report
Extrapolating a spline from dumux/common/spline.hh
on the right side fails.
What happened / Problem description:
Spline extrapolation failed:
dumux/dumux/common/splinecommon_.hh:636: Dumux::SplineCommon_<ScalarT, ImplementationT>::Scalar Dumux::SplineCommon_<ScalarT, ImplementationT>::h_(int) const [with ScalarT = double; ImplementationT = Dumux::VariableLengthSpline_<double>; Dumux::SplineCommon_<ScalarT, ImplementationT>::Scalar = double]: Assertion `x_(i) > x_(i-1)' failed.
When calling spline.eval(x, /*interpolate=*/true)
, the derivative is evaluated as a forward difference, which doesn't work on the right end.
It looks like the segment index for interpolation on the right is just 1 to high. It is numSamples_() - 1``, but should probably be
numSamples_() - 2``.
What you expected to happen:
I expected spline extrapolation to not fail.
How to reproduce it (as minimally and precisely as possible):
#include <vector>
#include <dumux/common/spline.hh>
using namespace Dumux;
using namespace std;
int main() {
vector<double> x = {0.0, 0.2, 0.4, 0.6, 0.8, 1.0};
vector<double> y = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
Spline<double, -1> spline;
spline.setXYContainers(x, y);
// interpolation
cout << "interpolation: " << spline.eval(0.5) << endl; // returns 1
// lower end extrapolation
cout << "lower end extrapolation: " << spline.eval(-0.1, true) << endl; // returns 1
// upper end extrapolation
cout << "upper end extrapolation: " << spline.eval(1.1, true) << endl; // fails
}
Anything else we need to know?:
There seems to be another small bug in eval
(splinecommon_.hh:141):
Scalar eval(Scalar x, bool extrapolate=false) const
{
assert(extrapolate || applies(x));
// handle extrapolation
if (extrapolate) {
if (x < xMin()) {
Scalar m = evalDerivative(xMin(), /*segmentIdx=*/0);
Scalar y0 = y_(0);
return y0 + m*(x - xMin());
}
else if (x > xMax()) {
Scalar m = evalDerivative(xMax(), /*segmentIdx=*/numSamples_()-1);
Scalar y0 = y_(numSamples_()-1);
return y0 + m*(x - xMax());
}
}
return eval_(x, segmentIdx_(x));
}
I think the calls to evalDerivative
should be to evalDerivative_
(underscore at the end).
Environment:
- DuMux version: master