From 9eec9a6037c96d74cb27a2f8c9bc1cc95017472f Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Fri, 30 Oct 2020 16:12:22 +0100
Subject: [PATCH] [bugfix][monotonecubicspline] Enable increasing and
decreasing values for both function and argument
---
dumux/common/monotonecubicspline.hh | 115 ++++++------
.../common/spline/test_monotonecubicspline.cc | 167 ++++++++++--------
2 files changed, 155 insertions(+), 127 deletions(-)
diff --git a/dumux/common/monotonecubicspline.hh b/dumux/common/monotonecubicspline.hh
index 9fa462e634..75abcfd61f 100644
--- a/dumux/common/monotonecubicspline.hh
+++ b/dumux/common/monotonecubicspline.hh
@@ -77,7 +77,7 @@ public:
// check some requirements
assert (x.size() == y.size());
assert (x.size() >=2);
- assert (std::is_sorted(x.begin(), x.end()));
+ assert (std::is_sorted(x.begin(), x.end()) || std::is_sorted(x.rbegin(), x.rend()));
// save a copy of the control points
x_ = x;
@@ -87,7 +87,8 @@ public:
numPoints_ = x.size();
// whether we are increasing
- increasing_ = y_.back() > y_.front();
+ increasingX_ = x_.back() > x_.front();
+ increasingY_ = y_.back() > y_.front();
// the slope at every control point
m_.resize(numPoints_);
@@ -114,21 +115,12 @@ public:
*/
Scalar eval(const Scalar x) const
{
- if (x <= x_.front())
+ if ((x <= x_.front() && increasingX_) || (x >= x_.front() && !increasingX_))
return y_.front() + m_.front()*(x - x_.front());
- else if (x > x_.back())
+ else if ((x > x_.back() && increasingX_) || (x < x_.back() && !increasingX_))
return y_.back() + m_.back()*(x - x_.back());
- else
- {
- const auto lookUpIndex = std::distance(x_.begin(), std::lower_bound(x_.begin(), x_.end(), x));
- assert(lookUpIndex != 0);
- // interpolate parametrization parameter t in [0,1]
- const auto h = (x_[lookUpIndex] - x_[lookUpIndex-1]);
- const auto t = (x - x_[lookUpIndex-1])/h;
- return y_[lookUpIndex-1]*Basis::h00(t) + h*m_[lookUpIndex-1]*Basis::h10(t)
- + y_[lookUpIndex]*Basis::h01(t) + h*m_[lookUpIndex]*Basis::h11(t);
- }
+ return eval_(x);
}
/*!
@@ -138,22 +130,12 @@ public:
*/
Scalar evalDerivative(const Scalar x) const
{
- if (x <= x_.front())
+ if ((x <= x_.front() && increasingX_) || (x >= x_.front() && !increasingX_))
return m_.front();
- else if (x > x_.back())
+ else if ((x > x_.back() && increasingX_) || (x < x_.back() && !increasingX_))
return m_.back();
- else
- {
- const auto lookUpIndex = std::distance(x_.begin(), std::lower_bound(x_.begin(), x_.end(), x));
- assert(lookUpIndex != 0);
- // interpolate parametrization parameter t in [0,1]
- const auto h = (x_[lookUpIndex] - x_[lookUpIndex-1]);
- const auto t = (x - x_[lookUpIndex-1])/h;
- const auto dtdx = 1.0/h;
- return y_[lookUpIndex-1]*Basis::dh00(t)*dtdx + m_[lookUpIndex-1]*Basis::dh10(t)
- + y_[lookUpIndex]*Basis::dh01(t)*dtdx + m_[lookUpIndex]*Basis::dh11(t);
- }
+ return evalDerivative_(x);
}
/*!
@@ -164,40 +146,39 @@ public:
*/
Scalar evalInverse(const Scalar y) const
{
- if (increasing_)
- {
- if (y <= y_.front())
- return x_.front() + (y - y_.front())/m_.front();
- else if (y > y_.back())
- return x_.back() + (y - y_.back())/m_.back();
- else
- {
- const auto lookUpIndex = std::distance(y_.begin(), std::lower_bound(y_.begin(), y_.end(), y));
- assert(lookUpIndex != 0);
-
- return evalInverse_(y, lookUpIndex);
- }
- }
+ if ((y <= y_.front() && increasingY_) || (y >= y_.front() && !increasingY_))
+ return x_.front() + (y - y_.front())/m_.front();
+ else if ((y > y_.back() && increasingY_) || (y < y_.back() && !increasingY_))
+ return x_.back() + (y - y_.back())/m_.back();
- else
- {
- if (y >= y_.front())
- return x_.front() + (y - y_.front())/m_.front();
- else if (y < y_.back())
- return x_.back() + (y - y_.back())/m_.back();
- else
- {
- const auto lookUpIndex = y_.size() - std::distance(y_.rbegin(), std::lower_bound(y_.rbegin(), y_.rend(), y));
- assert(lookUpIndex != 0);
-
- return evalInverse_(y, lookUpIndex);
- }
- }
+ return evalInverse_(y);
}
private:
- Scalar evalInverse_(const Scalar y, const std::size_t lookUpIndex) const
+ Scalar eval_(const Scalar x) const
{
+ // interpolate parametrization parameter t in [0,1]
+ const auto lookUpIndex = lookUpIndex_(x_, x, increasingX_);
+ const auto h = (x_[lookUpIndex] - x_[lookUpIndex-1]);
+ const auto t = (x - x_[lookUpIndex-1])/h;
+ return y_[lookUpIndex-1]*Basis::h00(t) + h*m_[lookUpIndex-1]*Basis::h10(t)
+ + y_[lookUpIndex]*Basis::h01(t) + h*m_[lookUpIndex]*Basis::h11(t);
+ }
+
+ Scalar evalDerivative_(const Scalar x) const
+ {
+ // interpolate parametrization parameter t in [0,1]
+ const auto lookUpIndex = lookUpIndex_(x_, x, increasingX_);
+ const auto h = (x_[lookUpIndex] - x_[lookUpIndex-1]);
+ const auto t = (x - x_[lookUpIndex-1])/h;
+ const auto dtdx = 1.0/h;
+ return y_[lookUpIndex-1]*Basis::dh00(t)*dtdx + m_[lookUpIndex-1]*Basis::dh10(t)
+ + y_[lookUpIndex]*Basis::dh01(t)*dtdx + m_[lookUpIndex]*Basis::dh11(t);
+ }
+
+ Scalar evalInverse_(const Scalar y) const
+ {
+ const auto lookUpIndex = lookUpIndex_(y_, y, increasingY_);
auto localPolynomial = [&](const auto x) {
// interpolate parametrization parameter t in [0,1]
const auto h = (x_[lookUpIndex] - x_[lookUpIndex-1]);
@@ -211,11 +192,31 @@ private:
return findScalarRootBrent(x_[lookUpIndex-1]-eps, x_[lookUpIndex]+eps, localPolynomial);
}
+ auto lookUpIndex_(const std::vector<Scalar>& vec, const Scalar v, bool increasing) const
+ {
+ return increasing ? lookUpIndexIncreasing_(vec, v) : lookUpIndexDecreasing_(vec, v);
+ }
+
+ auto lookUpIndexIncreasing_(const std::vector<Scalar>& vec, const Scalar v) const
+ {
+ const auto lookUpIndex = std::distance(vec.begin(), std::lower_bound(vec.begin(), vec.end(), v));
+ assert(lookUpIndex != 0);
+ return lookUpIndex;
+ }
+
+ auto lookUpIndexDecreasing_(const std::vector<Scalar>& vec, const Scalar v) const
+ {
+ const auto lookUpIndex = vec.size() - std::distance(vec.rbegin(), std::lower_bound(vec.rbegin(), vec.rend(), v));
+ assert(lookUpIndex != 0);
+ return lookUpIndex;
+ }
+
std::vector<Scalar> x_; //!< the x-coordinates
std::vector<Scalar> y_; //!< the y-coordinates
std::vector<Scalar> m_; //!< the slope for each control point
std::size_t numPoints_; //!< the number of control points
- bool increasing_; //!< if we are increasing monotone or not
+ bool increasingX_; //!< if we are increasing monotone or not
+ bool increasingY_; //!< if we are increasing monotone or not
};
} // end namespace Dumux
diff --git a/test/common/spline/test_monotonecubicspline.cc b/test/common/spline/test_monotonecubicspline.cc
index dc59b47701..04596dcb17 100644
--- a/test/common/spline/test_monotonecubicspline.cc
+++ b/test/common/spline/test_monotonecubicspline.cc
@@ -45,82 +45,109 @@ int main(int argc, char** argv)
{
Dune::MPIHelper::instance(argc, argv);
- // we test the spline interpolation against a sample function
- const auto f = [](double x){ return x*x*x; };
- const auto df = [](double x){ return 3*x*x; };
-
- // create some test samples
- const auto testPoints = Dumux::linspace(0.0, 4.0, 1000);
- const auto ref = eval(f, testPoints);
- const auto refDeriv = eval(df, testPoints);
-
- // create the spline sample points
- const auto samplePoints = Dumux::linspace(0.0, 5.0, 10);
- const auto y = eval(f, samplePoints);
-
- // create the spline
- Dumux::MonotoneCubicSpline<double> spline(samplePoints, y);
-
- // evaluate spline and derivative
- const auto result = eval([&](const double x) { return spline.eval(x); }, testPoints);
- const auto resultDeriv = eval([&](const double x) { return spline.evalDerivative(x); }, testPoints);
-
- // compute largest difference
- auto diff = result; auto diffDeriv = result;
- std::transform(result.begin(), result.end(), ref.begin(), diff.begin(), [](auto a, auto b){ return std::abs(a-b); });
- std::transform(resultDeriv.begin(), resultDeriv.end(), refDeriv.begin(), diffDeriv.begin(), [](auto a, auto b){ return std::abs(a-b); });
- const auto maxNorm = std::accumulate(diff.begin(), diff.end(), diff[0], [](auto a, auto b){ return std::max(a, b); })
- /(*std::max_element(ref.begin(), ref.end()));
- const auto maxNormDeriv = std::accumulate(diffDeriv.begin(), diffDeriv.end(), diffDeriv[0], [](auto a, auto b){ return std::max(a, b); })
- /(*std::max_element(refDeriv.begin(), refDeriv.end()));
- std::cout << "Maximum error: " << maxNorm << "\n";
- std::cout << "Maximum error in derivative: " << maxNormDeriv << "\n";
-
- if (maxNorm > 0.0008 || maxNormDeriv > 0.013)
- DUNE_THROW(Dune::Exception, "Maximum error in spline interpolation too large!");
-
- // test inverse by evaluating (x = f^-1(f(x))) for monotonically increasing function
+ const auto test = [](auto f, auto df, const auto& testPoints, const auto& samplePoints, const std::string& prefix)
{
- const auto resultX = eval([&](double x){ return spline.evalInverse(spline.eval(x)); }, testPoints);
- auto diffInverse = resultX;
- std::transform(resultX.begin(), resultX.end(), testPoints.begin(), diffInverse.begin(), [](auto a, auto b){ return std::abs(a-b); });
- const auto maxNormInverse = std::accumulate(diffInverse.begin(), diffInverse.end(), diffInverse[0], [](auto a, auto b){ return std::max(a, b); })
- /(*std::max_element(testPoints.begin(), testPoints.end()));
+ // create some test samples
+ const auto ref = eval(f, testPoints);
+ const auto refDeriv = eval(df, testPoints);
+
+ // create the spline sample points
+ const auto y = eval(f, samplePoints);
+
+ // create the spline
+ Dumux::MonotoneCubicSpline<double> spline(samplePoints, y);
+
+ // evaluate spline and derivative
+ const auto result = eval([&](const double x) { return spline.eval(x); }, testPoints);
+ const auto resultDeriv = eval([&](const double x) { return spline.evalDerivative(x); }, testPoints);
+
+ // compute largest difference
+ auto diff = result; auto diffDeriv = result;
+ std::transform(result.begin(), result.end(), ref.begin(), diff.begin(), [](auto a, auto b){ return std::abs(a-b); });
+ std::transform(resultDeriv.begin(), resultDeriv.end(), refDeriv.begin(), diffDeriv.begin(), [](auto a, auto b){ return std::abs(a-b); });
+ const auto maxNorm = std::accumulate(diff.begin(), diff.end(), diff[0], [](auto a, auto b){ return std::max(a, b); })
+ /std::abs(*std::max_element(ref.begin(), ref.end(), [](auto a, auto b){ return abs(a) < abs(b); }));
+ const auto maxNormDeriv = std::accumulate(diffDeriv.begin(), diffDeriv.end(), diffDeriv[0], [](auto a, auto b){ return std::max(a, b); })
+ /std::abs(*std::max_element(refDeriv.begin(), refDeriv.end(), [](auto a, auto b){ return abs(a) < abs(b); }));
+ std::cout << "Maximum error: " << std::scientific << maxNorm << "\n";
+ std::cout << "Maximum error in derivative: " << std::scientific << maxNormDeriv << "\n";
+
+ if (maxNorm > 0.0008 || maxNormDeriv > 0.013)
+ DUNE_THROW(Dune::Exception, "Maximum error in spline interpolation too large!");
+ // test inverse by evaluating (x = f^-1(f(x))) for monotonically increasing function
+ {
+ const auto resultX = eval([&](double x){ return spline.evalInverse(spline.eval(x)); }, testPoints);
+ auto diffInverse = resultX;
+ std::transform(resultX.begin(), resultX.end(), testPoints.begin(), diffInverse.begin(), [](auto a, auto b){ return std::abs(a-b); });
+ const auto maxNormInverse = std::accumulate(diffInverse.begin(), diffInverse.end(), diffInverse[0], [](auto a, auto b){ return std::max(a, b); })
+ /std::abs(*std::max_element(testPoints.begin(), testPoints.end(), [](auto a, auto b){ return abs(a) < abs(b); }));
+
+
+ std::cout << "Maximum error in identity using the inverse (mon. incr.): " << std::scientific << maxNormInverse << "\n";
+ if (maxNormInverse > 1e-13)
+ DUNE_THROW(Dune::Exception, "Maximum error in spline interpolation too large!");
+ }
+ // test inverse by evaluating (x = f^-1(f(x))) for monotonically decreasing function
+ {
+ auto reverseTest = testPoints;
+ std::reverse(reverseTest.begin(), reverseTest.end());
+ const auto resultX = eval([&](double x){ return spline.evalInverse(spline.eval(x)); }, reverseTest);
+ auto diffInverse = resultX;
+ std::transform(resultX.begin(), resultX.end(), reverseTest.begin(), diffInverse.begin(), [](auto a, auto b){ return std::abs(a-b); });
+ const auto maxNormInverse = std::accumulate(diffInverse.begin(), diffInverse.end(), diffInverse[0], [](auto a, auto b){ return std::max(a, b); })
+ /std::abs(*std::max_element(reverseTest.begin(), reverseTest.end(), [](auto a, auto b){ return abs(a) < abs(b); }));
+
+ std::cout << "Maximum error in identity using the inverse (mon. decr.): " << std::scientific << maxNormInverse << "\n";
+ if (maxNormInverse > 1e-13)
+ DUNE_THROW(Dune::Exception, "Maximum error in spline interpolation too large!");
+ }
+
+ // plot with Gnuplot (plot a bit more so we can see the linear extension)
+ const auto plotPoints = Dumux::linspace(-1.0, 5.0, 1000);
+ const auto refPlot = eval(f, plotPoints);
+ const auto refDerivPlot = eval(df, plotPoints);
+ const auto resultPlot = eval([&](const double x) { return spline.eval(x); }, plotPoints);
+ const auto resultDerivPlot = eval([&](const double x) { return spline.evalDerivative(x); }, plotPoints);
+
+ Dumux::GnuplotInterface<double> gnuplot(/*persist=*/false);
+ gnuplot.setOpenPlotWindow(false);
+ gnuplot.addDataSetToPlot(plotPoints, refPlot, prefix + "exp_reference");
+ gnuplot.addDataSetToPlot(plotPoints, refDerivPlot, prefix + "exp_reference_derivative");
+ gnuplot.addDataSetToPlot(plotPoints, resultPlot, prefix + "monotspline");
+ gnuplot.addDataSetToPlot(plotPoints, resultDerivPlot, prefix + "monotspline_derivative");
+ gnuplot.plot(prefix + "monotspline");
+ };
- std::cout << "Maximum error in identity using the inverse (mon. incr.): " << std::scientific << maxNormInverse << "\n";
- if (maxNormInverse > 1e-13)
- DUNE_THROW(Dune::Exception, "Maximum error in spline interpolation too large!");
- }
- // test inverse by evaluating (x = f^-1(f(x))) for monotonically decreasing function
{
- auto reverseTest = testPoints;
- std::reverse(reverseTest.begin(), reverseTest.end());
- const auto resultX = eval([&](double x){ return spline.evalInverse(spline.eval(x)); }, reverseTest);
- auto diffInverse = resultX;
- std::transform(resultX.begin(), resultX.end(), reverseTest.begin(), diffInverse.begin(), [](auto a, auto b){ return std::abs(a-b); });
- const auto maxNormInverse = std::accumulate(diffInverse.begin(), diffInverse.end(), diffInverse[0], [](auto a, auto b){ return std::max(a, b); })
- /(*std::max_element(reverseTest.begin(), reverseTest.end()));
-
- std::cout << "Maximum error in identity using the inverse (mon. decr.): " << std::scientific << maxNormInverse << "\n";
- if (maxNormInverse > 1e-13)
- DUNE_THROW(Dune::Exception, "Maximum error in spline interpolation too large!");
+ // we test the spline interpolation against a sample function
+ // monotonically increasing function
+ const auto f = [](double x){ return x*x*x; };
+ const auto df = [](double x){ return 3*x*x; };
+
+ const auto testPoints = Dumux::linspace(0.0, 4.0, 1000);
+ const auto samplePoints = Dumux::linspace(0.0, 5.0, 10);
+ test(f, df, testPoints, samplePoints, "x3in");
+
+ const auto testPoints2 = Dumux::linspace(4.0, 0.0, 1000);
+ const auto samplePoints2 = Dumux::linspace(5.0, 0.0, 10);
+ test(f, df, testPoints2, samplePoints2, "x3de");
}
- // plot with Gnuplot (plot a bit more so we can see the linear extension)
- const auto plotPoints = Dumux::linspace(-1.0, 5.0, 1000);
- const auto refPlot = eval(f, plotPoints);
- const auto refDerivPlot = eval(df, plotPoints);
- const auto resultPlot = eval([&](const double x) { return spline.eval(x); }, plotPoints);
- const auto resultDerivPlot = eval([&](const double x) { return spline.evalDerivative(x); }, plotPoints);
-
- Dumux::GnuplotInterface<double> gnuplot(/*persist=*/false);
- gnuplot.setOpenPlotWindow(false);
- gnuplot.addDataSetToPlot(plotPoints, refPlot, "exp_reference");
- gnuplot.addDataSetToPlot(plotPoints, refDerivPlot, "exp_reference_derivative");
- gnuplot.addDataSetToPlot(plotPoints, resultPlot, "monotspline");
- gnuplot.addDataSetToPlot(plotPoints, resultDerivPlot, "monotspline_derivative");
- gnuplot.plot("monotspline");
+ {
+ // we test the spline interpolation against a sample function
+ // monotonically decreasing function
+ const auto f = [](double x){ return -x*x*x; };
+ const auto df = [](double x){ return -3*x*x; };
+
+ const auto testPoints = Dumux::linspace(0.0, 4.0, 1000);
+ const auto samplePoints = Dumux::linspace(0.0, 5.0, 10);
+ test(f, df, testPoints, samplePoints, "mx3in");
+
+ const auto testPoints2 = Dumux::linspace(4.0, 0.0, 1000);
+ const auto samplePoints2 = Dumux::linspace(5.0, 0.0, 10);
+ test(f, df, testPoints2, samplePoints2, "mx3de");
+ }
return 0;
}
--
GitLab