[common] Add inverse of monotone cubic spline

dumux/common/monotonecubicspline.hh

@@ -35,6 +35,9 @@
 // Hermite basis functions
 #include <dumux/common/cubicsplinehermitebasis.hh>
 
+// for inversion
+#include <dumux/nonlinear/findscalarroot.hh>
+
 namespace Dumux {
 
 /*!
@@ -83,6 +86,9 @@ public:
         // the number of control points
         numPoints_ = x.size();
 
+        // whether we are increasing
+        increasing_ = y_.back() > y_.front();
+
         // the slope at every control point
         m_.resize(numPoints_);
 
@@ -150,11 +156,66 @@ public:
         }
     }
 
+    /*!
+     * \brief Evaluate the inverse function
+     * \param x the x-coordinate
+     * \note We extrapolate linearly if out of bounds
+     * \note Throws exception if inverse could not be found (e.g. not unique)
+     */
+    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);
+            }
+        }
+
+        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);
+            }
+        }
+    }
+
 private:
+    Scalar evalInverse_(const Scalar y, const std::size_t lookUpIndex) const
+    {
+        auto localPolynomial = [&](const auto x) {
+            // 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 - (y_[lookUpIndex-1]*Basis::h00(t) + h*m_[lookUpIndex-1]*Basis::h10(t)
+                   + y_[lookUpIndex]*Basis::h01(t) + h*m_[lookUpIndex]*Basis::h11(t));
+        };
+
+        // use an epsilon for the bracket
+        const auto eps = (x_[lookUpIndex]-x_[lookUpIndex-1])*1e-5;
+        return findScalarRootBrent(x_[lookUpIndex-1]-eps, x_[lookUpIndex]+eps, localPolynomial);
+    }
+
     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
 };
 
 } // end namespace Dumux

test/common/spline/test_monotonecubicspline.cc

@@ -79,6 +79,34 @@ int main(int argc, char** argv)
     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::max_element(testPoints.begin(), testPoints.end()));
+
+
+        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!");
+    }
+
     // 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);