Skip to content
Snippets Groups Projects
Commit 755f9a38 authored by Kilian Weishaupt's avatar Kilian Weishaupt
Browse files

Merge branch 'fix/clamp-input-material-laws' into 'master' 

Fix/clamp input material laws

See merge request !365
parent 77e83cff
No related branches found
No related tags found
Hide whitespace changes
Inline Side-by-side
dumux/material/fluidmatrixinteractions/2p/brookscorey.hh
+ 73
− 23
View file @ 755f9a38
......@@ -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 &params, 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 &params, 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 &params, 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 &params, 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 &params, 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 &params, 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 &params, 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 &params, 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
dumux/material/fluidmatrixinteractions/2p/vangenuchten.hh
+ 83
− 30
View file @ 755f9a38
......@@ -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 &params, 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 &params, 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 &params, 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 &params, 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 &params, 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 &params, 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 &params, 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 &params, 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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment