diff --git a/dumux/material/fluidmatrixinteractions/2p/brookscorey.hh b/dumux/material/fluidmatrixinteractions/2p/brookscorey.hh
index e9c1b0a4e432a2e782d93e5f1ab46a250bf852d8..35ef237461161dfd0b1580abac1d8af113e58f1e 100644
--- a/dumux/material/fluidmatrixinteractions/2p/brookscorey.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/brookscorey.hh
@@ -66,10 +66,17 @@ public:
* 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.
* \return Capillary pressure calculated by Brooks & Corey constitutive relation.
+ *
+ * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number,
+ * by clamping the input.
*/
static Scalar pc(const Params ¶ms, Scalar swe)
{
- assert(0 <= swe && swe <= 1);
+ using std::pow;
+ using std::min;
+ using std::max;
+
+ swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
return params.pe()*pow(swe, -1.0/params.lambda());
}
@@ -86,13 +93,18 @@ public:
* 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.
* \return Effective wetting phase saturation calculated as inverse of BrooksCorey constitutive relation.
+ *
+ * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number,
+ * by clamping the input.
*/
static Scalar sw(const Params ¶ms, Scalar pc)
{
- assert(pc >= 0);
+ using std::pow;
+ using std::max;
- Scalar tmp = pow(pc/params.pe(), -params.lambda());
- return std::min(std::max(tmp, Scalar(0.0)), Scalar(1.0));
+ pc = max(pc, 0.0); // the equation below is undefined for negative pcs
+
+ return pow(pc/params.pe(), -params.lambda());
}
/*!
@@ -110,10 +122,17 @@ public:
* 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.
* \return Partial derivative of \f$\mathrm{[p_c]}\f$ w.r.t. effective saturation according to Brooks & Corey.
+ *
+ * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number,
+ * by clamping the input.
*/
static Scalar dpc_dswe(const Params ¶ms, Scalar swe)
{
- assert(0 <= swe && swe <= 1);
+ using std::pow;
+ using std::min;
+ using std::max;
+
+ swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
return - params.pe()/params.lambda() * pow(swe, -1/params.lambda() - 1);
}
@@ -127,10 +146,16 @@ public:
* 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.
* \return Partial derivative of effective saturation w.r.t. \f$\mathrm{[p_c]}\f$ according to Brooks & Corey.
+ *
+ * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number,
+ * by clamping the input.
*/
static Scalar dswe_dpc(const Params ¶ms, Scalar pc)
{
- assert(pc >= 0);
+ using std::pow;
+ using std::max;
+
+ pc = max(pc, 0.0); // the equation below is undefined for negative pcs
return -params.lambda()/params.pe() * pow(pc/params.pe(), - params.lambda() - 1);
}
@@ -145,10 +170,17 @@ public:
* 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.
* \return Relative permeability of the wetting phase calculated as implied by Brooks & Corey.
+ *
+ * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number,
+ * by clamping the input.
*/
static Scalar krw(const Params ¶ms, Scalar swe)
{
- assert(0 <= swe && swe <= 1);
+ using std::pow;
+ using std::min;
+ using std::max;
+
+ swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
return pow(swe, 2.0/params.lambda() + 3);
}
@@ -164,10 +196,17 @@ public:
* and then the params container is constructed accordingly. Afterwards the values are set there, too.
* \return Derivative of the relative permeability of the wetting phase w.r.t. effective wetting phase
* saturation calculated as implied by Brooks & Corey.
+ *
+ * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number,
+ * by clamping the input.
*/
static Scalar dkrw_dswe(const Params ¶ms, Scalar swe)
{
- assert(0 <= swe && swe <= 1);
+ using std::pow;
+ using std::min;
+ using std::max;
+
+ swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
return (2.0/params.lambda() + 3)*pow(swe, 2.0/params.lambda() + 2);
}
@@ -182,14 +221,21 @@ public:
* 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.
* \return Relative permeability of the non-wetting phase calculated as implied by Brooks & Corey.
+ *
+ * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number,
+ * by clamping the input.
*/
static Scalar krn(const Params ¶ms, Scalar swe)
{
- assert(0 <= swe && swe <= 1);
+ using std::pow;
+ using std::min;
+ using std::max;
+
+ swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
- Scalar exponent = 2.0/params.lambda() + 1;
- Scalar tmp = 1. - swe;
- return tmp*tmp*(1. - pow(swe, exponent));
+ const Scalar exponent = 2.0/params.lambda() + 1;
+ const Scalar tmp = 1.0 - swe;
+ return tmp*tmp*(1.0 - pow(swe, exponent));
}
/*!
@@ -204,22 +250,26 @@ public:
* and then the params container is constructed accordingly. Afterwards the values are set there, too.
* \return Derivative of the relative permeability of the non-wetting phase w.r.t. effective wetting phase
* saturation calculated as implied by Brooks & Corey.
+ *
+ * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number,
+ * by clamping the input.
*/
static Scalar dkrn_dswe(const Params ¶ms, Scalar swe)
{
- assert(0 <= swe && swe <= 1);
-
- return
- 2.0*(swe - 1)*(
- 1 +
- pow(swe, 2.0/params.lambda())*(
- 1.0/params.lambda() + 1.0/2 -
- swe*(1.0/params.lambda() + 1.0/2)
- )
- );
+ using std::pow;
+ using std::min;
+ using std::max;
+
+ swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
+
+ return 2.0*(swe - 1)*(1 + pow(swe, 2.0/params.lambda())*(1.0/params.lambda() + 1.0/2
+ - swe*(1.0/params.lambda() + 1.0/2)
+ )
+ );
}
};
-}
+
+} // end namespace Dumux
#endif // BROOKS_COREY_HH
diff --git a/dumux/material/fluidmatrixinteractions/2p/vangenuchten.hh b/dumux/material/fluidmatrixinteractions/2p/vangenuchten.hh
index 388e9aca0bce7225d7113dadb3cd730a28c67004..75f1c768db023c923efd6dcf088475a34c2aa224 100644
--- a/dumux/material/fluidmatrixinteractions/2p/vangenuchten.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/vangenuchten.hh
@@ -64,11 +64,19 @@ public:
* \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.
+ * \note Instead of undefined behaviour if swe is not in the valid range, we return a valid number,
+ * by clamping the input. Note that for pc(swe = 0.0) = inf, have a look at RegularizedVanGenuchten if this is a problem.
*/
static Scalar pc(const Params ¶ms, Scalar swe)
{
- assert(0 <= swe && swe <= 1);
- return pow(pow(swe, -1.0/params.vgm()) - 1, 1.0/params.vgn())/params.vgAlpha();
+ using std::pow;
+ using std::min;
+ using std::max;
+
+ swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
+
+ const Scalar pc = pow(pow(swe, -1.0/params.vgm()) - 1, 1.0/params.vgn())/params.vgAlpha();
+ return pc;
}
/*!
@@ -84,12 +92,18 @@ public:
* 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.
* \return The effective saturation of the wetting phase \f$\mathrm{\overline{S}_w}\f$
+ * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number,
+ * i.e. sw(pc < 0.0) = 0.0, by clamping the input to the physical bounds.
*/
static Scalar sw(const Params ¶ms, Scalar pc)
{
- assert(pc >= 0);
+ using std::pow;
+ using std::max;
+
+ pc = max(pc, 0.0); // the equation below is undefined for negative pcs
- return pow(pow(params.vgAlpha()*pc, params.vgn()) + 1, -params.vgm());
+ const Scalar sw = pow(pow(params.vgAlpha()*pc, params.vgn()) + 1, -params.vgm());
+ return sw;
}
/*!
@@ -107,14 +121,21 @@ public:
* \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.
- */
+ *
+ * \note Instead of undefined behaviour if swe is not in the valid range, we return a valid number,
+ * by clamping the input.
+ */
static Scalar dpc_dswe(const Params ¶ms, Scalar swe)
{
- assert(0 <= swe && swe <= 1);
+ using std::pow;
+ using std::min;
+ using std::max;
- Scalar powSwe = pow(swe, -1/params.vgm());
+ swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
+
+ const Scalar powSwe = pow(swe, -1/params.vgm());
return - 1.0/params.vgAlpha() * pow(powSwe - 1, 1.0/params.vgn() - 1)/params.vgn()
- * powSwe/swe/params.vgm();
+ * powSwe/swe/params.vgm();
}
/*!
@@ -125,14 +146,19 @@ public:
* \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.
+ *
+ * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number,
+ * by clamping the input.
*/
static Scalar dswe_dpc(const Params ¶ms, Scalar pc)
{
- assert(pc >= 0);
+ using std::pow;
+ using std::max;
+
+ pc = max(pc, 0.0); // the equation below is undefined for negative pcs
- Scalar powAlphaPc = pow(params.vgAlpha()*pc, params.vgn());
- return -pow(powAlphaPc + 1, -params.vgm()-1)*
- params.vgm()*powAlphaPc/pc*params.vgn();
+ const Scalar powAlphaPc = pow(params.vgAlpha()*pc, params.vgn());
+ return -pow(powAlphaPc + 1, -params.vgm()-1)*params.vgm()*powAlphaPc/pc*params.vgn();
}
/*!
@@ -143,12 +169,21 @@ public:
* \param swe The mobile saturation of the wetting phase.
* \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. */
+ * is constructed accordingly. Afterwards the values are set there, too.
+ *
+ * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number,
+ * by clamping the input.
+ */
static Scalar krw(const Params ¶ms, Scalar swe)
{
- assert(0 <= swe && swe <= 1);
+ using std::pow;
+ using std::sqrt;
+ using std::min;
+ using std::max;
+
+ swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
- Scalar r = 1.0 - pow(1.0 - pow(swe, 1.0/params.vgm()), params.vgm());
+ const Scalar r = 1.0 - pow(1.0 - pow(swe, 1.0/params.vgm()), params.vgm());
return sqrt(swe)*r*r;
}
@@ -161,14 +196,22 @@ public:
* \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.
+ *
+ * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number,
+ * by clamping the input.
*/
static Scalar dkrw_dswe(const Params ¶ms, Scalar swe)
{
- assert(0 <= swe && swe <= 1);
+ using std::pow;
+ using std::sqrt;
+ using std::min;
+ using std::max;
- const Scalar x = 1.0 - std::pow(swe, 1.0/params.vgm());
- const Scalar xToM = std::pow(x, params.vgm());
- return (1.0 - xToM)/std::sqrt(swe) * ( (1.0 - xToM)/2 + 2*xToM*(1.0-x)/x );
+ swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
+
+ const Scalar x = 1.0 - pow(swe, 1.0/params.vgm());
+ const Scalar xToM = pow(x, params.vgm());
+ return (1.0 - xToM)/sqrt(swe) * ( (1.0 - xToM)/2 + 2*xToM*(1.0-x)/x );
}
/*!
@@ -180,14 +223,19 @@ public:
* \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.
+ *
+ * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number,
+ * by clamping the input.
*/
static Scalar krn(const Params ¶ms, Scalar swe)
{
- assert(0 <= swe && swe <= 1);
+ using std::pow;
+ using std::min;
+ using std::max;
- return
- pow(1 - swe, 1.0/3) *
- pow(1 - pow(swe, 1.0/params.vgm()), 2*params.vgm());
+ swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
+
+ return pow(1 - swe, 1.0/3) * pow(1 - pow(swe, 1.0/params.vgm()), 2*params.vgm());
}
/*!
@@ -200,19 +248,24 @@ public:
* \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.
+ *
+ * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number,
+ * by clamping the input.
*/
static Scalar dkrn_dswe(const Params ¶ms, Scalar swe)
{
- assert(0 <= swe && swe <= 1);
+ using std::pow;
+ using std::min;
+ using std::max;
- const Scalar x = std::pow(swe, 1.0/params.vgm());
- return
- -std::pow(1.0 - x, 2*params.vgm())
- *std::pow(1.0 - swe, -2.0/3)
- *(1.0/3 + 2*x/swe);
+ swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= sw <= 1.0
+
+ const Scalar x = pow(swe, 1.0/params.vgm());
+ return -pow(1.0 - x, 2*params.vgm()) * pow(1.0 - swe, -2.0/3) * (1.0/3 + 2*x/swe);
}
};
-}
+
+} // end namespace Dumux
#endif