From 14f7dc7efb42cb0eec0f6b4738692aa4a187893c Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Fri, 3 Feb 2017 14:22:06 +0100
Subject: [PATCH] [material] Introduce endPointPc method to material laws

The end point capillary pressure is the capillary pressure at sw=1.
For Brooks-Corey laws it is also called entry pressure.
---
 .../fluidmatrixinteractions/2p/brookscorey.hh | 10 ++++++++
 .../fluidmatrixinteractions/2p/efftoabslaw.hh | 10 ++++++++
 .../2p/linearmaterial.hh                      | 10 ++++++++
 .../2p/regularizedbrookscorey.hh              | 10 ++++++++
 .../2p/regularizedlinearmaterial.hh           | 10 ++++++++
 .../2p/regularizedvangenuchten.hh             | 24 ++++++++++++++++++-
 .../2p/vangenuchten.hh                        | 10 ++++++++
 7 files changed, 83 insertions(+), 1 deletion(-)

diff --git a/dumux/material/fluidmatrixinteractions/2p/brookscorey.hh b/dumux/material/fluidmatrixinteractions/2p/brookscorey.hh
index 35ef237461..370746e1c1 100644
--- a/dumux/material/fluidmatrixinteractions/2p/brookscorey.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/brookscorey.hh
@@ -107,6 +107,16 @@ public:
         return pow(pc/params.pe(), -params.lambda());
     }
 
+    /*!
+     * \brief The capillary pressure at Swe = 1.0 also called end point capillary pressure
+     *
+     * \param params A container object that is populated with the appropriate coefficients for the respective law.
+     *                  Therefore, in the (problem specific) spatialParameters  first, the material law is chosen, and then the params container
+     *                  is constructed accordingly. Afterwards the values are set there, too.
+     */
+    static Scalar endPointPc(const Params &params)
+    { return params.pe(); }
+
     /*!
      * \brief The partial derivative of the capillary
      *        pressure w.r.t. the effective saturation according to Brooks & Corey.
diff --git a/dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh b/dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh
index 93469ba29f..bb5b4164a8 100644
--- a/dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh
@@ -101,6 +101,16 @@ public:
         return sweToSw_(params, EffLaw::sw(params, pc));
     }
 
+     /*!
+     * \brief The capillary pressure at Swe = 1.0 also called end point capillary pressure
+     *
+     * \param params A container object that is populated with the appropriate coefficients for the respective law.
+     *                  Therefore, in the (problem specific) spatialParameters  first, the material law is chosen, and then the params container
+     *                  is constructed accordingly. Afterwards the values are set there, too.
+     */
+    static Scalar endPointPc(const Params &params)
+    { return EffLaw::endPointPc(params); }
+
     /*!
      * \brief Returns the partial derivative of the capillary
      *        pressure w.r.t the absolute saturation.
diff --git a/dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh b/dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh
index 8c28cf4268..cee892592b 100644
--- a/dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh
@@ -89,6 +89,16 @@ public:
         return 1 - (pc - params.entryPc())/(params.maxPc() - params.entryPc());
     }
 
+    /*!
+     * \brief The capillary pressure at Swe = 1.0 also called end point capillary pressure
+     *
+     * \param params A container object that is populated with the appropriate coefficients for the respective law.
+     *                  Therefore, in the (problem specific) spatialParameters  first, the material law is chosen, and then the params container
+     *                  is constructed accordingly. Afterwards the values are set there, too.
+     */
+    static Scalar endPointPc(const Params &params)
+    { return params.entryPc(); }
+
     /*!
      * \brief Returns the partial derivative of the capillary
      *        pressure w.r.t. the effective saturation.
diff --git a/dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh b/dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh
index cfc2b6f90d..22bef00571 100644
--- a/dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh
@@ -151,6 +151,16 @@ public:
         return BrooksCorey::sw(params, pc);
     }
 
+    /*!
+     * \brief The capillary pressure at Swe = 1.0 also called end point capillary pressure
+     *
+     * \param params A container object that is populated with the appropriate coefficients for the respective law.
+     *                  Therefore, in the (problem specific) spatialParameters  first, the material law is chosen, and then the params container
+     *                  is constructed accordingly. Afterwards the values are set there, too.
+     */
+    static Scalar endPointPc(const Params &params)
+    { return params.pe(); }
+
     /*!
      * \brief A regularized version of the partial derivative
      *        of the \f$\mathrm{p_c(\overline{S}_w)}\f$ w.r.t. effective saturation
diff --git a/dumux/material/fluidmatrixinteractions/2p/regularizedlinearmaterial.hh b/dumux/material/fluidmatrixinteractions/2p/regularizedlinearmaterial.hh
index b1382fe381..5947169d7d 100644
--- a/dumux/material/fluidmatrixinteractions/2p/regularizedlinearmaterial.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/regularizedlinearmaterial.hh
@@ -102,6 +102,16 @@ public:
         return LinearMaterial::sw(params, pc);
     }
 
+     /*!
+     * \brief The capillary pressure at Swe = 1.0 also called end point capillary pressure
+     *
+     * \param params A container object that is populated with the appropriate coefficients for the respective law.
+     *                  Therefore, in the (problem specific) spatialParameters  first, the material law is chosen, and then the params container
+     *                  is constructed accordingly. Afterwards the values are set there, too.
+     */
+    static Scalar endPointPc(const Params &params)
+    { return params.entryPc(); }
+
     /*!
      * \brief Returns the partial derivative of the capillary
      *        pressure to the effective saturation.
diff --git a/dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh b/dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh
index 30bf0d7c79..426ca41f98 100644
--- a/dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh
@@ -156,6 +156,12 @@ public:
         // Genuchten's law
         Scalar sw;
         if (pc <= 0) {
+            // for swThHigh = 1.0 the slope would get infinity
+            // swThHigh > 1.0 are not sensible threshold values
+            // setting swThHigh = 1.0 is a way to disable regularization
+            if (swThHigh > 1.0 - std::numeric_limits<Scalar>::epsilon())
+                return 1.0;
+
             // invert straight line for swe > 1.0
             Scalar yTh = VanGenuchten::pc(params, swThHigh);
             Scalar m1 = (0.0 - yTh)/(1.0 - swThHigh)*2;
@@ -187,6 +193,16 @@ public:
         return sw;
     }
 
+    /*!
+     * \brief The capillary pressure at Swe = 1.0 also called end point capillary pressure
+     *
+     * \param params A container object that is populated with the appropriate coefficients for the respective law.
+     *                  Therefore, in the (problem specific) spatialParameters  first, the material law is chosen, and then the params container
+     *                  is constructed accordingly. Afterwards the values are set there, too.
+     */
+    static Scalar endPointPc(const Params &params)
+    { return 0.0; }
+
     /*!
     * \brief A regularized version of the partial derivative
     *        of the \f$\mathrm{p_c(\overline{S}_w)}\f$ w.r.t. effective saturation
@@ -277,7 +293,7 @@ public:
 
         if (swe < 0)
             return 0;
-        else if (swe > 1)
+        else if (swe > 1 - std::numeric_limits<Scalar>::epsilon())
             return 1;
         else if (swe > swThHigh) {
             typedef Dumux::Spline<Scalar> Spline;
@@ -357,6 +373,12 @@ private:
     {
         const Scalar swThHigh = params.pcHighSw();
 
+        // for swThHigh = 1.0 the slope would get infinity
+        // swThHigh > 1.0 are not sensible threshold values
+        // setting swThHigh = 1.0 is a way to disable regularization
+        if (swThHigh > 1.0 - std::numeric_limits<Scalar>::epsilon())
+            return 0.0;
+
         Scalar pcswHigh = VanGenuchten::pc(params, swThHigh);
         return (0 - pcswHigh)/(1.0 - swThHigh);
     }
diff --git a/dumux/material/fluidmatrixinteractions/2p/vangenuchten.hh b/dumux/material/fluidmatrixinteractions/2p/vangenuchten.hh
index 75f1c768db..ffe43ef4b6 100644
--- a/dumux/material/fluidmatrixinteractions/2p/vangenuchten.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/vangenuchten.hh
@@ -106,6 +106,16 @@ public:
         return sw;
     }
 
+    /*!
+     * \brief The capillary pressure at Swe = 1.0 also called end point capillary pressure
+     *
+     * \param params A container object that is populated with the appropriate coefficients for the respective law.
+     *                  Therefore, in the (problem specific) spatialParameters  first, the material law is chosen, and then the params container
+     *                  is constructed accordingly. Afterwards the values are set there, too.
+     */
+    static Scalar endPointPc(const Params &params)
+    { return 0.0; }
+
     /*!
      * \brief The partial derivative of the capillary
      *        pressure w.r.t. the effective saturation according to van Genuchten.
-- 
GitLab