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