add invertion methods to some of the two-phase fluid-matrix interactions

these are required for hysteresis

opm/material/fluidmatrixinteractions/BrooksCorey.hpp

@@ -30,6 +30,8 @@
 #include "BrooksCoreyParams.hpp"
 
 #include <opm/material/common/MathToolbox.hpp>
+#include <opm/material/common/ErrorMacros.hpp>
+#include <opm/material/common/Exceptions.hpp>
 
 #include <algorithm>
 #include <cassert>
@@ -167,7 +169,17 @@ public:
 
         assert(0 <= Sw && Sw <= 1);
 
-        return params.entryPressure()*Toolbox::pow(Sw, -1.0/params.lambda());
+        return params.entryPressure()*Toolbox::pow(Sw, -1/params.lambda());
+    }
+
+    template <class Evaluation>
+    static Evaluation twoPhaseSatPcnwInv(const Params &params, const Evaluation& pcnw)
+    {
+        typedef MathToolbox<Evaluation> Toolbox;
+
+        assert(pcnw > 0.0);
+
+        return Toolbox::pow(params.entryPressure()/pcnw, -params.lambda());
     }
 
     /*!
@@ -243,6 +255,14 @@ public:
         return Toolbox::pow(Sw, 2.0/params.lambda() + 3);
     }
 
+    template <class Evaluation>
+    static Evaluation twoPhaseSatKrwInv(const Params &params, const Evaluation& krw)
+    {
+        typedef MathToolbox<Evaluation> Toolbox;
+
+        return Toolbox::pow(krw, 1.0/(2.0/params.lambda() + 3));
+    }
+
     /*!
      * \brief The relative permeability for the non-wetting phase of
      *        the medium as implied by the Brooks-Corey
@@ -273,6 +293,33 @@ public:
         const Evaluation Sn = 1. - Sw;
         return Sn*Sn*(1. - Toolbox::pow(Sw, exponent));
     }
+
+    template <class Evaluation>
+    static Evaluation twoPhaseSatKrnInv(const Params &params, const Evaluation& krn)
+    {
+        typedef MathToolbox<Evaluation> Toolbox;
+
+        // since inverting the formula for krn is hard to do analytically, we use the
+        // Newton-Raphson method
+        Evaluation Sw = Toolbox::createConstant(0.5);
+        Scalar eps = 1e-10;
+        for (int i = 0; i < 20; ++i) {
+            Evaluation f = krn - twoPhaseSatKrn(params, Sw);
+            Evaluation fStar = krn - twoPhaseSatKrn(params, Sw + eps);
+            Evaluation fPrime = (fStar - f)/eps;
+
+            Evaluation delta = f/fPrime;
+            Sw -= delta;
+            if (Sw < 0)
+                Sw = Toolbox::createConstant(0.0);
+            if (Toolbox::abs(delta) < 1e-10)
+                return Sw;
+        }
+
+        OPM_THROW(NumericalIssue,
+                  "Couldn't invert the Brooks-Corey non-wetting phase"
+                  " relperm within 20 iterations");
+    }
 };
 } // namespace Opm
 

opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterial.hpp

@@ -143,6 +143,14 @@ public:
     static Evaluation twoPhaseSatPcnw(const Params &params, const Evaluation& Sw)
     { return eval_(params.SwSamples(), params.pcnwSamples(), Sw); }
 
+    template <class Evaluation>
+    static Evaluation twoPhaseSatPcnwInv(const Params &params, const Evaluation& pcnw)
+    {
+        return eval_(params.pcnwSamples(),
+                     params.SwSamples(),
+                     pcnw);
+    }
+
     /*!
      * \brief The saturation-capillary pressure curve
      */
@@ -181,7 +189,7 @@ public:
     }
 
     template <class Evaluation>
-    static Evaluation twoPhaseSatKrw(const Params &params, const Evaluation Sw)
+    static Evaluation twoPhaseSatKrw(const Params &params, const Evaluation& Sw)
     {
         typedef MathToolbox<Evaluation> Toolbox;
 
@@ -189,6 +197,14 @@ public:
         return Toolbox::max(0.0, Toolbox::min(1.0, res));
     }
 
+    template <class Evaluation>
+    static Evaluation twoPhaseSatKrwInv(const Params &params, const Evaluation& krw)
+    {
+        return eval_(params.krwSamples(),
+                     params.SwSamples(),
+                     krw);
+    }
+
     /*!
      * \brief The relative permeability for the non-wetting phase
      *        of the porous medium
@@ -213,18 +229,36 @@ public:
                                                          Sw)));
     }
 
+    template <class Evaluation>
+    static Evaluation twoPhaseSatKrnInv(const Params &params, const Evaluation& krn)
+    {
+        return eval_(params.krnSamples(),
+                     params.SwSamples(),
+                     krn);
+    }
+
 private:
     template <class Evaluation>
     static Evaluation eval_(const ValueVector &xValues,
                             const ValueVector &yValues,
                             const Evaluation& x)
+    {
+        if (xValues.front() < xValues.back())
+            return evalAscending_(xValues, yValues, x);
+        return evalDescending_(xValues, yValues, x);
+    }
+
+    template <class Evaluation>
+    static Evaluation evalAscending_(const ValueVector &xValues,
+                                     const ValueVector &yValues,
+                                     const Evaluation& x)
     {
         typedef MathToolbox<Evaluation> Toolbox;
 
-        if (Toolbox::value(x) < xValues.front())
+        if (x <= xValues.front())
-            return Toolbox::createConstant(yValues.front());
+            return yValues.front();
-        if (Toolbox::value(x) > xValues.back())
+        if (x >= xValues.back())
-            return Toolbox::createConstant(yValues.back());
+            return yValues.back();
 
         int segIdx = findSegmentIndex_(xValues, Toolbox::value(x));
 
@@ -239,6 +273,31 @@ private:
         return y0 + (x - x0)*m;
     }
 
+    template <class Evaluation>
+    static Evaluation evalDescending_(const ValueVector &xValues,
+                                      const ValueVector &yValues,
+                                      const Evaluation& x)
+    {
+        typedef MathToolbox<Evaluation> Toolbox;
+
+        if (x >= xValues.front())
+            return yValues.front();
+        if (x <= xValues.back())
+            return yValues.back();
+
+        int segIdx = findSegmentIndexDescending_(xValues, Toolbox::value(x));
+
+        Scalar x0 = xValues[segIdx];
+        Scalar x1 = xValues[segIdx + 1];
+
+        Scalar y0 = yValues[segIdx];
+        Scalar y1 = yValues[segIdx + 1];
+
+        Scalar m = (y1 - y0)/(x1 - x0);
+
+        return y0 + (x - x0)*m;
+    }
+
     template <class Evaluation>
     static Evaluation evalDeriv_(const ValueVector &xValues,
                                  const ValueVector &yValues,
@@ -283,6 +342,28 @@ private:
 
         return lowIdx;
     }
+
+    static int findSegmentIndexDescending_(const ValueVector &xValues, Scalar x)
+    {
+        int n = xValues.size() - 1;
+        assert(n >= 1); // we need at least two sampling points!
+        if (xValues[n] >= x)
+            return n - 1;
+        else if (xValues[0] <= x)
+            return 0;
+
+        // bisection
+        int lowIdx = 0, highIdx = n;
+        while (lowIdx + 1 < highIdx) {
+            int curIdx = (lowIdx + highIdx)/2;
+            if (xValues[curIdx] >= x)
+                lowIdx = curIdx;
+            else
+                highIdx = curIdx;
+        }
+
+        return lowIdx;
+    }
 };
 } // namespace Opm
 

opm/material/fluidmatrixinteractions/RegularizedBrooksCorey.hpp

@@ -299,6 +299,19 @@ public:
         return BrooksCorey::twoPhaseSatKrw(params, Sw);
     }
 
+    template <class Evaluation>
+    static Evaluation twoPhaseSatKrwInv(const Params &params, const Evaluation& krw)
+    {
+        typedef MathToolbox<Evaluation> Toolbox;
+
+        if (krw <= 0.0)
+            return Toolbox::createConstant(0.0);
+        else if (krw >= 1.0)
+            return Toolbox::createConstant(1.0);
+
+        return BrooksCorey::twoPhaseSatKrwInv(params, krw);
+    }
+
     /*!
      * \brief Regularized version of the relative permeability of the
      *        non-wetting phase of the Brooks-Corey curves.
@@ -335,6 +348,19 @@ public:
 
         return BrooksCorey::twoPhaseSatKrn(params, Sw);
     }
+
+    template <class Evaluation>
+    static Evaluation twoPhaseSatKrnInv(const Params &params, const Evaluation& krn)
+    {
+        typedef MathToolbox<Evaluation> Toolbox;
+
+        if (krn <= 0.0)
+            return Toolbox::createConstant(1.0);
+        else if (krn >= 1.0)
+            return Toolbox::createConstant(0.0);
+
+        return BrooksCorey::twoPhaseSatKrnInv(params, krn);
+    }
 };
 } // namespace Opm