Merge branch 'feature/monotone-cubic-spline' into 'master'

Feature/monotone cubic spline

Closes #712

See merge request !1620

dumux/common/CMakeLists.txt

@@ -9,6 +9,7 @@ boundaryflag.hh
 boundarytypes.hh
 boundingboxtree.hh
 cubicspline.hh
+cubicsplinehermitebasis.hh
 defaultmappertraits.hh
 defaultusagemessage.hh
 deprecated.hh
@@ -23,6 +24,7 @@ intersectionmapper.hh
 intrange.hh
 loggingparametertree.hh
 math.hh
+monotonecubicspline.hh
 numericdifferentiation.hh
 optional.hh
 parameters.hh

dumux/common/cubicsplinehermitebasis.hh

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 3 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup Common
+ * \brief The cubic hermite spline basis
+ */
+#ifndef DUMUX_COMMON_CUBIC_SPLINE_HERMITE_BASIS_HH
+#define DUMUX_COMMON_CUBIC_SPLINE_HERMITE_BASIS_HH
+
+namespace Dumux {
+
+/*!
+ * \ingroup Common
+ * \brief The cubic spline hermite basis
+ * \note see https://en.wikipedia.org/wiki/Cubic_Hermite_spline
+ */
+template<class Scalar = double>
+struct CubicSplineHermiteBasis
+{
+    static constexpr Scalar h00(const Scalar t)
+    { return t*(2.0*t*t - 3.0*t) + 1.0; }
+
+    static constexpr Scalar h10(const Scalar t)
+    { return t*(t*t - 2.0*t + 1.0); }
+
+    static constexpr Scalar h01(const Scalar t)
+    { return t*t*(3.0 - 2.0*t); }
+
+    static constexpr Scalar h11(const Scalar t)
+    { return t*t*(t - 1.0); }
+
+    static constexpr Scalar dh00(const Scalar t)
+    { return 6.0*t*(t - 1.0); }
+
+    static constexpr Scalar dh10(const Scalar t)
+    { return t*(3.0*t - 4.0) + 1.0; }
+
+    static constexpr Scalar dh01(const Scalar t)
+    { return 6.0*t*(1.0 - t); }
+
+    static constexpr Scalar dh11(const Scalar t)
+    { return t*(3.0*t - 2.0); }
+};
+
+} // end namespace Dumux
+
+#endif

dumux/common/monotonecubicspline.hh

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 3 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup Common
+ * \brief A monotone cubic spline
+ */
+#ifndef DUMUX_COMMON_MONOTONE_CUBIC_SPLINE_HH
+#define DUMUX_COMMON_MONOTONE_CUBIC_SPLINE_HH
+
+#include <cmath>
+#include <cassert>
+#include <iterator>
+#include <vector>
+#include <algorithm>
+
+#include <dune/common/float_cmp.hh>
+
+// Hermite basis functions
+#include <dumux/common/cubicsplinehermitebasis.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup Common
+ * \brief A monotone cubic spline
+ * \note Contruction after Fritsch & Butland (1984) (see https://doi.org/10.1137/0905021)
+ * \note The resulting interpolation is globally monotone but only C^1
+ */
+template<class Scalar = double>
+class MonotoneCubicSpline
+{
+    using Basis = CubicSplineHermiteBasis<Scalar>;
+public:
+    /*!
+     * \brief Default constructor
+     */
+    MonotoneCubicSpline() = default;
+
+    /*!
+     * \brief Contruct a monotone cubic spline from the control points (x[i], y[i])
+     * \note if the data set is monotone, monotonicity is preserved
+     * \param x a vector of x-coordinates
+     * \param y a vector of y-coordinates
+     */
+    MonotoneCubicSpline(const std::vector<Scalar>& x, const std::vector<Scalar> y)
+    {
+        // check some requirements
+        assert (x.size() == y.size());
+        assert (x.size() >=2);
+        assert (std::is_sorted(x.begin(), x.end()));
+
+        updatePoints(x, y);
+    }
+
+    /*!
+     * \brief Create a monotone cubic spline from the control points (x[i], y[i])
+     * \param x a vector of x-coordinates
+     * \param y a vector of y-coordinates
+     */
+    void updatePoints(const std::vector<Scalar>& x, const std::vector<Scalar>& y)
+    {
+        // save a copy of the control points
+        x_ = x;
+        y_ = y;
+
+        // the number of control points
+        numPoints_ = x.size();
+
+        // the slope at every control point
+        m_.resize(numPoints_);
+
+        // compute slopes (see Fritsch & Butland (1984), Eq. (5))
+        Scalar deltaX = (x[1]-x[0]);
+        Scalar secant = m_.front() = (y[1]-y[0])/deltaX;
+        Scalar prevDeltaX = deltaX;
+        Scalar prevSecant = secant;
+        for (int i = 1; i < numPoints_-1; ++i, prevSecant = secant, prevDeltaX = deltaX)
+        {
+            deltaX = (x[i+1]-x[i]);
+            secant = (y[i+1]-y[i])/deltaX;
+            const auto alpha = (prevDeltaX + 2*deltaX)/(3*(prevDeltaX + deltaX));
+            m_[i] = prevSecant*secant > 0.0 ? prevSecant*secant/(alpha*secant + (1.0-alpha)*prevSecant) : 0.0;
+        }
+        m_.back() = secant;
+    }
+
+    /*!
+     * \brief Evaluate the y value at a given x value
+     * \param x the x-coordinate
+     * \note We extrapolate linearly if out of bounds
+     */
+    Scalar eval(const Scalar x) const
+    {
+        if (x <= x_.front())
+            return y_.front() + m_.front()*(x - x_.front());
+        else if (x > x_.back())
+            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);
+        }
+    }
+
+    /*!
+     * \brief Evaluate the first derivative dy/dx at a given x value
+     * \param x the x-coordinate
+     * \note We extrapolate linearly if out of bounds
+     */
+    Scalar evalDerivative(const Scalar x) const
+    {
+        if (x <= x_.front())
+            return m_.front();
+        else if (x > x_.back())
+            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);
+        }
+    }
+
+private:
+    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
+};
+
+} // end namespace Dumux
+
+#endif

test/common/spline/CMakeLists.txt

@@ -2,3 +2,5 @@ dune_add_test(SOURCES test_spline.cc
               LABELS unit)
 dune_add_test(SOURCES test_cubicspline.cc
               LABELS unit)
+dune_add_test(SOURCES test_monotonecubicspline.cc
+              LABELS unit)

test/common/spline/test_cubicspline.cc

@@ -23,6 +23,7 @@
 #include <config.h>
 #include <cmath>
 #include <vector>
+#include <numeric>
 #include <algorithm>
 #include <functional>
 

test/common/spline/test_monotonecubicspline.cc

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 3 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Test the simple cubic spline implementation
+ */
+#include <config.h>
+#include <cmath>
+#include <vector>
+#include <numeric>
+#include <algorithm>
+#include <functional>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/parallel/mpihelper.hh>
+#include <dumux/common/monotonecubicspline.hh>
+#include <dumux/io/gnuplotinterface.hh>
+
+std::vector<double> linspace(const double begin, const double end, const double samples)
+{
+    const double delta = (end-begin)/static_cast<double>(samples-1);
+    std::vector<double> vec(samples);
+    for (int i = 0; i < samples; ++i)
+        vec[i] = begin + i*delta;
+    return vec;
+}
+
+template<class Function>
+std::vector<double> eval(const Function& f, const std::vector<double>& x)
+{
+    auto y = x;
+    std::transform(x.begin(), x.end(), y.begin(), [&](const double x) { return f(x); });
+    return y;
+}
+
+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 = 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 = linspace(0.0, 5.0, 10);
+    const auto y = eval(f, samplePoints);
+
+    // create the spline
+    Dumux::MonotoneCubicSpline 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!");
+
+    // plot with Gnuplot (plot a bit more so we can see the linear extension)
+    const auto plotPoints = 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.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");
+
+    return 0;
+}