diff --git a/dumux/material/fluidmatrixinteractions/2p/brookscorey.hh b/dumux/material/fluidmatrixinteractions/2p/brookscorey.hh index a69963c302d342c3b284904148381461a2261abe..9b7a19caeab3cd1726bc01ae8de0a5c784dfbfc4 100644 --- a/dumux/material/fluidmatrixinteractions/2p/brookscorey.hh +++ b/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 { diff --git a/dumux/material/fluidmatrixinteractions/2p/vangenuchten.hh b/dumux/material/fluidmatrixinteractions/2p/vangenuchten.hh index 69aeb71099088c2905810378ac3fa211e63f0288..07e04277623d5556833c7010eb9b863175dc262c 100644 --- a/dumux/material/fluidmatrixinteractions/2p/vangenuchten.hh +++ b/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_;