Merge branch 'fix/materiallaw-985' into 'master'

Fix/materiallaw vangenuchten spline assert

Closes #985

See merge request !2459

dumux/material/fluidmatrixinteractions/2p/brookscorey.hh

@@ -279,7 +279,6 @@ public:
 #include <algorithm>
 
 #include <dumux/common/parameters.hh>
-#include <dumux/common/spline.hh>
 #include <dumux/common/optionalscalar.hh>
 #include <dumux/material/fluidmatrixinteractions/2p/materiallaw.hh>
 
@@ -592,7 +591,6 @@ public:
         S pcLowSwe_ = 0.01;
     };
 
-    //! Initialize the spline
     template<class MaterialLaw>
     void init(const MaterialLaw* m, const std::string& paramGroup)
     {
@@ -624,8 +622,7 @@ public:
      * \brief The regularized capillary pressure-saturation curve
      * regularized part:
      *  - low saturation:  extend the \f$\mathrm{p_c(S_w)}\f$ curve with the slope at the regularization point (i.e. no kink).
-     *  - high saturation: connect the high regularization point with \f$\mathrm{\overline{S}_w =1}\f$
-     *                     with a spline and continue linearly for \f$\mathrm{\overline{S}_w > 1}\f$
+     *  - high saturation: continue linearly
      */
     OptionalScalar<Scalar> pc(const Scalar swe) const
     {

dumux/material/fluidmatrixinteractions/2p/vangenuchten.hh

@@ -848,11 +848,13 @@ private:
         pcLowSwePcValue_ = m->template pc<false>(lowSw);
         pcHighSwePcValue_ = m->template pc<false>(highSw);
 
-        pcSpline_ = Spline<Scalar>(highSwe, 1.0, // x0, x1
-                                   pcHighSwePcValue_, 0, // y0, y1
-                                   pcDerivativeHighSweThreshold_, pcDerivativeHighSweEnd_); // m0, m1
-
-
+        // Only initialize regularization spline if given parameters are in
+        // the admissible range. When constructing with non-sensible parameters
+        // the spline construction might fail (e.g. highSwe == 1.0)
+        if (highSwe < 1.0)
+            pcSpline_ = Spline<Scalar>(highSwe, 1.0, // x0, x1
+                                       pcHighSwePcValue_, 0, // y0, y1
+                                       pcDerivativeHighSweThreshold_, pcDerivativeHighSweEnd_); // m0, m1
     }
 
     template<class MaterialLaw>
@@ -868,14 +870,15 @@ private:
         const auto krnLowSw = m->template krn<false>(lowSw);
         const auto dkrnLowSw = m->template dkrn_dsw<false>(lowSw)*dsw_dswe;
 
-        krwSpline_ = Spline<Scalar>(highSwe, 1.0, // x0, x1
-                                    krwHighSw, 1.0, // y0, y1
-                                    dkrwHighSw, 0.0); // m0, m1
-
-        krnSpline_ = Spline<Scalar>(0.0, lowSwe, // x0, x1
-                                    1.0, krnLowSw, // y0, y1
-                                    0.0, dkrnLowSw); // m0, m1
+        if (highSwe < 1.0)
+            krwSpline_ = Spline<Scalar>(highSwe, 1.0, // x0, x1
+                                        krwHighSw, 1.0, // y0, y1
+                                        dkrwHighSw, 0.0); // m0, m1
 
+        if (lowSwe > 0.0)
+            krnSpline_ = Spline<Scalar>(0.0, lowSwe, // x0, x1
+                                        1.0, krnLowSw, // y0, y1
+                                        0.0, dkrnLowSw); // m0, m1
     }
 
     Scalar pcLowSwe_, pcHighSwe_;