OilfieldToolsNet/opm-common
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

DenseAD: make less fuzz about it

Browse Source 
this patch converts to code to use the convenience functions instead
of the math toolboxes whereever possible. the main advantage is that
Opm::foo(x) will work regardless of the type of `x`, but it also
reduces visual clutter.

also, constant Evaluations are now directly created by assigning
Scalars, which removes further visual noise.

while I hope it improves the readability of the code,
functionality-wise this patch should not change anything.
This commit is contained in:
Andreas Lauser 2017-06-12 16:36:14 +02:00
parent acfa7c43e5
commit 0f6540bdad
69 changed files with 738 additions and 1278 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
opm/material/binarycoefficients/Air_Mesitylene.hpp
View File

@ -58,14 +58,13 @@ public:
template <class Evaluation>
static Evaluation gasDiffCoeff(Evaluation temperature, Evaluation pressure)
{
typedef Opm::MathToolbox<Evaluation> Toolbox;
typedef Opm::Air<double> Air;
typedef Opm::Mesitylene<double> Mesitylene;
temperature = Toolbox::max(temperature, 1e-9); // regularization
temperature = Toolbox::min(temperature, 500.0); // regularization
pressure = Toolbox::max(pressure, 0.0); // regularization
pressure = Toolbox::min(pressure, 1e8); // regularization
temperature = Opm::max(temperature, 1e-9); // regularization
temperature = Opm::min(temperature, 500.0); // regularization
pressure = Opm::max(pressure, 0.0); // regularization
pressure = Opm::min(pressure, 1e8); // regularization
const double M_m = 1e3*Mesitylene::molarMass(); // [g/mol] molecular weight of mesitylene
const double M_a = 1e3*Air::molarMass(); // [g/mol] molecular weight of air
@ -79,13 +78,13 @@ public:
const double T_scal_am = std::sqrt(T_scal_a*T_scal_m);
Evaluation T_star = temperature/T_scal_am;
T_star = Toolbox::max(T_star, 1e-5); // regularization
T_star = Opm::max(T_star, 1e-5); // regularization
const Evaluation Omega = 1.06036/Toolbox::pow(T_star, 0.1561) + 0.193/Toolbox::exp(T_star*0.47635)
+ 1.03587/Toolbox::exp(T_star*1.52996) + 1.76474/Toolbox::exp(T_star*3.89411);
const Evaluation Omega = 1.06036/Opm::pow(T_star, 0.1561) + 0.193/Opm::exp(T_star*0.47635)
+ 1.03587/Opm::exp(T_star*1.52996) + 1.76474/Opm::exp(T_star*3.89411);
const double B_ = 0.00217 - 0.0005*std::sqrt(1.0/M_a + 1.0/M_m);
const double Mr = (M_a + M_m)/(M_a*M_m);
const Evaluation D_am = (B_*Toolbox::pow(temperature, 1.5) * std::sqrt(Mr))
const Evaluation D_am = (B_*Opm::pow(temperature, 1.5) * std::sqrt(Mr))
/(1e-5*pressure*std::pow(sigma_am, 2.0) * Omega); // [cm^2/s]
return 1e-4*D_am; // [m^2/s]

17
opm/material/binarycoefficients/Air_Xylene.hpp
View File

@ -56,14 +56,13 @@ public:
template <class Evaluation>
static Evaluation gasDiffCoeff(Evaluation temperature, Evaluation pressure)
{
typedef Opm::MathToolbox<Evaluation> Toolbox;
typedef Opm::Air<double> Air;
typedef Opm::Xylene<double> Xylene;
temperature = Toolbox::max(temperature, 1e-9); // regularization
temperature = Toolbox::min(temperature, 500.0); // regularization
pressure = Toolbox::max(pressure, 0.0); // regularization
pressure = Toolbox::min(pressure, 1e8); // regularization
temperature = Opm::max(temperature, 1e-9); // regularization
temperature = Opm::min(temperature, 500.0); // regularization
pressure = Opm::max(pressure, 0.0); // regularization
pressure = Opm::min(pressure, 1e8); // regularization
const double M_x = 1e3*Xylene::molarMass(); // [g/mol] molecular weight of xylene
const double M_a = 1e3*Air::molarMass(); // [g/mol] molecular weight of air
@ -76,14 +75,14 @@ public:
const double T_scal_x = 1.15*Tb_x;
const double T_scal_ax = std::sqrt(T_scal_a*T_scal_x);
const Evaluation& T_star = Toolbox::max(temperature/T_scal_ax, 1e-5);
const Evaluation& T_star = Opm::max(temperature/T_scal_ax, 1e-5);
const Evaluation& Omega = 1.06036/Toolbox::pow(T_star, 0.1561) + 0.193/Toolbox::exp(T_star*0.47635)
+ 1.03587/Toolbox::exp(T_star*1.52996) + 1.76474/Toolbox::exp(T_star*3.89411);
const Evaluation& Omega = 1.06036/Opm::pow(T_star, 0.1561) + 0.193/Opm::exp(T_star*0.47635)
+ 1.03587/Opm::exp(T_star*1.52996) + 1.76474/Opm::exp(T_star*3.89411);
const double B_ = 0.00217 - 0.0005*std::sqrt(1.0/M_a + 1.0/M_x);
const double Mr = (M_a + M_x)/(M_a*M_x);
return 1e-4
*(B_*Toolbox::pow(temperature,1.5)*std::sqrt(Mr))
*(B_*Opm::pow(temperature,1.5)*std::sqrt(Mr))
/(1e-5*pressure*std::pow(sigma_ax, 2.0)*Omega);
}

52
opm/material/binarycoefficients/Brine_CO2.hpp
View File

@ -153,8 +153,6 @@ public:
template <class Evaluation>
static Evaluation fugacityCoefficientCO2(const Evaluation& temperature, const Evaluation& pg)
{
typedef MathToolbox<Evaluation> Toolbox;
Valgrind::CheckDefined(temperature);
Valgrind::CheckDefined(pg);
@ -165,20 +163,20 @@ public:
Scalar R = IdealGas::R * 10.; // ideal gas constant with unit bar cm^3 /(K mol)
Evaluation lnPhiCO2;
lnPhiCO2 = Toolbox::log(V / (V - b_CO2));
lnPhiCO2 = Opm::log(V / (V - b_CO2));
lnPhiCO2 += b_CO2 / (V - b_CO2);
lnPhiCO2 -= 2 * a_CO2 / (R * Toolbox::pow(temperature, 1.5) * b_CO2) * log((V + b_CO2) / V);
lnPhiCO2 -= 2 * a_CO2 / (R * Opm::pow(temperature, 1.5) * b_CO2) * log((V + b_CO2) / V);
lnPhiCO2 +=
a_CO2 * b_CO2
/ (R
* Toolbox::pow(temperature, 1.5)
* Opm::pow(temperature, 1.5)
* b_CO2
* b_CO2)
* (Toolbox::log((V + b_CO2) / V)
* (Opm::log((V + b_CO2) / V)
- b_CO2 / (V + b_CO2));
lnPhiCO2 -= Toolbox::log(pg_bar * V / (R * temperature));
lnPhiCO2 -= Opm::log(pg_bar * V / (R * temperature));
return Toolbox::exp(lnPhiCO2); // fugacity coefficient of CO2
return Opm::exp(lnPhiCO2); // fugacity coefficient of CO2
}
/*!
@ -192,8 +190,6 @@ public:
template <class Evaluation>
static Evaluation fugacityCoefficientH2O(const Evaluation& temperature, const Evaluation& pg)
{
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& V = 1 / (CO2::gasDensity(temperature, pg) / CO2::molarMass()) * 1.e6; // molar volume in cm^3/mol
const Evaluation& pg_bar = pg / 1.e5; // gas phase pressure in bar
const Evaluation& a_CO2 = (7.54e7 - 4.13e4 * temperature);// mixture parameter of Redlich-Kwong equation
@ -204,13 +200,13 @@ public:
Evaluation lnPhiH2O;
lnPhiH2O =
Toolbox::log(V/(V - b_CO2))
Opm::log(V/(V - b_CO2))
+ b_H2O/(V - b_CO2) - 2*a_CO2_H2O
/ (R*Toolbox::pow(temperature, 1.5)*b_CO2)*Toolbox::log((V + b_CO2)/V)
+ a_CO2*b_H2O/(R*Toolbox::pow(temperature, 1.5)*b_CO2*b_CO2)
*(Toolbox::log((V + b_CO2)/V) - b_CO2/(V + b_CO2))
- Toolbox::log(pg_bar*V/(R*temperature));
return Toolbox::exp(lnPhiH2O); // fugacity coefficient of H2O
/ (R*Opm::pow(temperature, 1.5)*b_CO2)*Opm::log((V + b_CO2)/V)
+ a_CO2*b_H2O/(R*Opm::pow(temperature, 1.5)*b_CO2*b_CO2)
*(Opm::log((V + b_CO2)/V) - b_CO2/(V + b_CO2))
- Opm::log(pg_bar*V/(R*temperature));
return Opm::exp(lnPhiH2O); // fugacity coefficient of H2O
}
private:
@ -272,13 +268,11 @@ private:
const Evaluation& pg,
Scalar molalityNaCl)
{
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& lambda = computeLambda_(temperature, pg); // lambda_{CO2-Na+}
const Evaluation& xi = computeXi_(temperature, pg); // Xi_{CO2-Na+-Cl-}
const Evaluation& lnGammaStar =
2*molalityNaCl*lambda + xi*molalityNaCl*molalityNaCl;
return Toolbox::exp(lnGammaStar);
return Opm::exp(lnGammaStar);
}
/*!
@ -292,15 +286,13 @@ private:
template <class Evaluation>
static Evaluation computeA_(const Evaluation& temperature, const Evaluation& pg)
{
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& deltaP = pg / 1e5 - 1; // pressure range [bar] from p0 = 1bar to pg[bar]
Scalar v_av_H2O = 18.1; // average partial molar volume of H2O [cm^3/mol]
Scalar R = IdealGas::R * 10;
const Evaluation& k0_H2O = equilibriumConstantH2O_(temperature); // equilibrium constant for H2O at 1 bar
const Evaluation& phi_H2O = fugacityCoefficientH2O(temperature, pg); // fugacity coefficient of H2O for the water-CO2 system
const Evaluation& pg_bar = pg / 1.e5;
return k0_H2O/(phi_H2O*pg_bar)*Toolbox::exp(deltaP*v_av_H2O/(R*temperature));
return k0_H2O/(phi_H2O*pg_bar)*Opm::exp(deltaP*v_av_H2O/(R*temperature));
}
/*!
@ -314,15 +306,13 @@ private:
template <class Evaluation>
static Evaluation computeB_(const Evaluation& temperature, const Evaluation& pg)
{
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& deltaP = pg / 1e5 - 1; // pressure range [bar] from p0 = 1bar to pg[bar]
const Scalar v_av_CO2 = 32.6; // average partial molar volume of CO2 [cm^3/mol]
const Scalar R = IdealGas::R * 10;
const Evaluation& k0_CO2 = equilibriumConstantCO2_(temperature); // equilibrium constant for CO2 at 1 bar
const Evaluation& phi_CO2 = fugacityCoefficientCO2(temperature, pg); // fugacity coefficient of CO2 for the water-CO2 system
const Evaluation& pg_bar = pg / 1.e5;
return phi_CO2*pg_bar/(55.508*k0_CO2)*Toolbox::exp(-(deltaP*v_av_CO2)/(R*temperature));
return phi_CO2*pg_bar/(55.508*k0_CO2)*Opm::exp(-(deltaP*v_av_CO2)/(R*temperature));
}
/*!
@ -335,8 +325,6 @@ private:
template <class Evaluation>
static Evaluation computeLambda_(const Evaluation& temperature, const Evaluation& pg)
{
typedef MathToolbox<Evaluation> Toolbox;
static const Scalar c[6] =
{ -0.411370585, 6.07632013E-4, 97.5347708, -0.0237622469, 0.0170656236, 1.41335834E-5 };
@ -347,7 +335,7 @@ private:
+ c[2]/temperature
+ c[3]*pg_bar/temperature
+ c[4]*pg_bar/(630.0 - temperature)
+ c[5]*temperature*Toolbox::log(pg_bar);
+ c[5]*temperature*Opm::log(pg_bar);
}
/*!
@ -376,12 +364,10 @@ private:
template <class Evaluation>
static Evaluation equilibriumConstantCO2_(const Evaluation& temperature)
{
typedef MathToolbox<Evaluation> Toolbox;
Evaluation temperatureCelcius = temperature - 273.15;
static const Scalar c[3] = { 1.189, 1.304e-2, -5.446e-5 };
Evaluation logk0_CO2 = c[0] + temperatureCelcius*(c[1] + temperatureCelcius*c[2]);
Evaluation k0_CO2 = Toolbox::pow(10.0, logk0_CO2);
Evaluation k0_CO2 = Opm::pow(10.0, logk0_CO2);
return k0_CO2;
}
@ -394,13 +380,11 @@ private:
template <class Evaluation>
static Evaluation equilibriumConstantH2O_(const Evaluation& temperature)
{
typedef MathToolbox<Evaluation> Toolbox;
Evaluation temperatureCelcius = temperature - 273.15;
static const Scalar c[4] = { -2.209, 3.097e-2, -1.098e-4, 2.048e-7 };
Evaluation logk0_H2O =
c[0] + temperatureCelcius*(c[1] + temperatureCelcius*(c[2] + temperatureCelcius*c[3]));
return Toolbox::pow(10.0, logk0_H2O);
return Opm::pow(10.0, logk0_H2O);
}
};

4
opm/material/binarycoefficients/FullerMethod.hpp
View File

@ -58,14 +58,12 @@ inline Evaluation fullerMethod(const Scalar* M, // molar masses [g/mol]
const Evaluation& temperature, // [K]
const Evaluation& pressure) // [Pa]
{
typedef Opm::MathToolbox<Evaluation> Toolbox;
// "effective" molar mass in [g/m^3]
Scalar Mab = Opm::harmonicMean(M[0], M[1]);
// Fuller's method
const Evaluation& tmp = std::pow(SigmaNu[0], 1./3) + std::pow(SigmaNu[1], 1./3);
return 1e-4 * (143.0*Toolbox::pow(temperature, 1.75))/(pressure*std::sqrt(Mab)*tmp*tmp);
return 1e-4 * (143.0*Opm::pow(temperature, 1.75))/(pressure*std::sqrt(Mab)*tmp*tmp);
}
} // namespace BinaryCoeff

10
opm/material/binarycoefficients/H2O_Air.hpp
View File

@ -53,11 +53,7 @@ public:
*/
template <class Evaluation>
static Evaluation henry(const Evaluation& temperature)
{
typedef Opm::MathToolbox<Evaluation> Toolbox;
return 1.0/((0.8942+1.47*Toolbox::exp(-0.04394*(temperature-273.15)))*1.E-10);
}
{ return 1.0/((0.8942+1.47*Opm::exp(-0.04394*(temperature-273.15)))*1e-10); }
/*!
* \brief Binary diffusion coefficent \f$\mathrm{[m^2/s]}\f$ for molecular water and air
@ -73,14 +69,12 @@ public:
template <class Evaluation>
static Evaluation gasDiffCoeff(const Evaluation& temperature, const Evaluation& pressure)
{
typedef Opm::MathToolbox<Evaluation> Toolbox;
double Theta=1.8;
double Daw=2.13e-5; /* reference value */
double pg0=1.e5; /* reference pressure */
double T0=273.15; /* reference temperature */
return Daw*(pg0/pressure)*Toolbox::pow((temperature/T0),Theta);
return Daw*(pg0/pressure)*Opm::pow((temperature/T0),Theta);
}
/*!

18
opm/material/binarycoefficients/H2O_Mesitylene.hpp
View File

@ -66,15 +66,13 @@ public:
template <class Evaluation>
static Evaluation gasDiffCoeff(Evaluation temperature, Evaluation pressure)
{
typedef Opm::MathToolbox<Evaluation> Toolbox;
typedef Opm::H2O<double> H2O;
typedef Opm::Mesitylene<double> Mesitylene;
temperature = Toolbox::max(temperature, 1e-9); // regularization
temperature = Toolbox::min(temperature, 500.0); // regularization
pressure = Toolbox::max(pressure, 0.0); // regularization
pressure = Toolbox::min(pressure, 1e8); // regularization
temperature = Opm::max(temperature, 1e-9); // regularization
temperature = Opm::min(temperature, 500.0); // regularization
pressure = Opm::max(pressure, 0.0); // regularization
pressure = Opm::min(pressure, 1e8); // regularization
const double M_m = 1e3*Mesitylene::molarMass(); // [g/mol] molecular weight of mesitylene
const double M_w = 1e3*H2O::molarMass(); // [g/mol] molecular weight of water
@ -89,13 +87,13 @@ public:
const double T_scal_m = 1.15*Tb_m;
const double T_scal_wm = std::sqrt(T_scal_w*T_scal_m);
const Evaluation T_star = Toolbox::max(temperature/T_scal_wm, 1e-5);
const Evaluation T_star = Opm::max(temperature/T_scal_wm, 1e-5);
const Evaluation& Omega = 1.06036/Toolbox::pow(T_star,0.1561) + 0.193/Toolbox::exp(T_star*0.47635)
+ 1.03587/Toolbox::exp(T_star*1.52996) + 1.76474/Toolbox::exp(T_star*3.89411);
const Evaluation& Omega = 1.06036/Opm::pow(T_star,0.1561) + 0.193/Opm::exp(T_star*0.47635)
+ 1.03587/Opm::exp(T_star*1.52996) + 1.76474/Opm::exp(T_star*3.89411);
const double B_ = 0.00217 - 0.0005*std::sqrt(1.0/M_w + 1.0/M_m);
const double Mr = (M_w + M_m)/(M_w*M_m);
const Evaluation D_wm = (B_*Toolbox::pow(temperature,1.6)*std::sqrt(Mr))
const Evaluation D_wm = (B_*Opm::pow(temperature,1.6)*std::sqrt(Mr))
/(1e-5*pressure*std::pow(sigma_wm, 2.0)*Omega); // [cm^2/s]
return D_wm*1e-4; // [m^2/s]

17
opm/material/binarycoefficients/H2O_Xylene.hpp
View File

@ -65,14 +65,13 @@ public:
template <class Evaluation>
static Evaluation gasDiffCoeff(Evaluation temperature, Evaluation pressure)
{
typedef Opm::MathToolbox<Evaluation> Toolbox;
typedef Opm::H2O<double> H2O;
typedef Opm::Xylene<double> Xylene;
temperature = Toolbox::max(temperature, 1e-9); // regularization
temperature = Toolbox::min(temperature, 500.0); // regularization
pressure = Toolbox::max(pressure, 0.0); // regularization
pressure = Toolbox::min(pressure, 1e8); // regularization
temperature = Opm::max(temperature, 1e-9); // regularization
temperature = Opm::min(temperature, 500.0); // regularization
pressure = Opm::max(pressure, 0.0); // regularization
pressure = Opm::min(pressure, 1e8); // regularization
const double M_x = 1e3*Xylene::molarMass(); // [g/mol] molecular weight of xylene
const double M_w = 1e3*H2O::molarMass(); // [g/mol] molecular weight of water
@ -87,14 +86,14 @@ public:
const double T_scal_x = 1.15*Tb_x;
const double T_scal_wx = std::sqrt(T_scal_w*T_scal_x);
const Evaluation& T_star = Toolbox::max(temperature/T_scal_wx, 1e-5);
const Evaluation& T_star = Opm::max(temperature/T_scal_wx, 1e-5);
const Evaluation& Omega = 1.06036/Toolbox::pow(T_star, 0.1561) + 0.193/Toolbox::exp(T_star*0.47635)
+ 1.03587/Toolbox::exp(T_star*1.52996) + 1.76474/Toolbox::exp(T_star*3.89411);
const Evaluation& Omega = 1.06036/Opm::pow(T_star, 0.1561) + 0.193/Opm::exp(T_star*0.47635)
+ 1.03587/Opm::exp(T_star*1.52996) + 1.76474/Opm::exp(T_star*3.89411);
const double B_ = 0.00217 - 0.0005*std::sqrt(1.0/M_w + 1.0/M_x);
const double Mr = (M_w + M_x)/(M_w*M_x);
return 1e-4
*(B_*Toolbox::pow(temperature,1.6)*std::sqrt(Mr))
*(B_*Opm::pow(temperature,1.6)*std::sqrt(Mr))
/(1e-5*pressure*std::pow(sigma_wx, 2.0)*Omega);
}

9
opm/material/binarycoefficients/HenryIapws.hpp
View File

@ -48,7 +48,6 @@ inline Evaluation henryIAPWS(Scalar E,
Scalar H,
const Evaluation& temperature)
{
typedef Opm::MathToolbox<Evaluation> Toolbox;
typedef Opm::H2O<Evaluation> H2O;
Evaluation Tr = temperature/H2O::criticalTemperature();
@ -66,20 +65,20 @@ inline Evaluation henryIAPWS(Scalar E,
Evaluation f = 0;
for (int i = 0; i < 6; ++i) {
f += c[i]*Toolbox::pow(tau, d[i]);
f += c[i]*Opm::pow(tau, d[i]);
}
const Evaluation& exponent =
q*F +
E/temperature*f +
(F +
G*Toolbox::pow(tau, 2.0/3) +
G*Opm::pow(tau, 2.0/3) +
H*tau)*
Toolbox::exp((H2O::tripleTemperature() - temperature)/100);
Opm::exp((H2O::tripleTemperature() - temperature)/100);
// CAUTION: K_D is formulated in mole fractions. We have to
// multiply it with the vapor pressure of water in order to get
// derivative of the partial pressure.
return Toolbox::exp(exponent)*H2O::vaporPressure(temperature);
return Opm::exp(exponent)*H2O::vaporPressure(temperature);
}
} // namespace Opm

49
opm/material/common/PolynomialUtils.hpp
View File

@ -52,8 +52,7 @@ unsigned invertLinearPolynomial(SolContainer& sol,
Scalar a,
Scalar b)
{
typedef MathToolbox<Scalar> Toolbox;
if (std::abs(Toolbox::scalarValue(a)) < 1e-30)
if (std::abs(Opm::scalarValue(a)) < 1e-30)
return 0;
sol[0] = -b/a;
@ -82,10 +81,8 @@ unsigned invertQuadraticPolynomial(SolContainer& sol,
Scalar b,
Scalar c)
{
typedef MathToolbox<Scalar> Toolbox;
// check for a line
if (std::abs(Toolbox::scalarValue(a)) < 1e-30)
if (std::abs(Opm::scalarValue(a)) < 1e-30)
return invertLinearPolynomial(sol, b, c);
// discriminant
@ -93,7 +90,7 @@ unsigned invertQuadraticPolynomial(SolContainer& sol,
if (Delta < 0)
return 0; // no real roots
Delta = Toolbox::sqrt(Delta);
Delta = Opm::sqrt(Delta);
sol[0] = (- b + Delta)/(2*a);
sol[1] = (- b - Delta)/(2*a);
@ -112,8 +109,6 @@ void invertCubicPolynomialPostProcess_(SolContainer& sol,
Scalar c,
Scalar d)
{
typedef MathToolbox<Scalar> Toolbox;
// do one Newton iteration on the analytic solution if the
// precision is increased
for (int i = 0; i < numSol; ++i) {
@ -121,12 +116,12 @@ void invertCubicPolynomialPostProcess_(SolContainer& sol,
Scalar fOld = d + x*(c + x*(b + x*a));
Scalar fPrime = c + x*(2*b + x*3*a);
if (std::abs(Toolbox::scalarValue(fPrime)) < 1e-30)
if (std::abs(Opm::scalarValue(fPrime)) < 1e-30)
continue;
x -= fOld/fPrime;
Scalar fNew = d + x*(c + x*(b + x*a));
if (std::abs(Toolbox::scalarValue(fNew)) < std::abs(Toolbox::scalarValue(fOld)))
if (std::abs(Opm::scalarValue(fNew)) < std::abs(Opm::scalarValue(fOld)))
sol[i] = x;
}
}
@ -156,10 +151,8 @@ unsigned invertCubicPolynomial(SolContainer* sol,
Scalar c,
Scalar d)
{
typedef MathToolbox<Scalar> Toolbox;
// reduces to a quadratic polynomial
if (std::abs(Toolbox::scalarValue(a)) < 1e-30)
if (std::abs(Opm::scalarValue(a)) < 1e-30)
return invertQuadraticPolynomial(sol, b, c, d);
// normalize the polynomial
@ -172,7 +165,7 @@ unsigned invertCubicPolynomial(SolContainer* sol,
Scalar p = c - b*b/3;
Scalar q = d + (2*b*b*b - 9*b*c)/27;
if (std::abs(Toolbox::scalarValue(p)) > 1e-30 && std::abs(Toolbox::scalarValue(q)) > 1e-30) {
if (std::abs(Opm::scalarValue(p)) > 1e-30 && std::abs(Opm::scalarValue(q)) > 1e-30) {
// At this point
//
// t^3 + p*t + q = 0
@ -206,9 +199,9 @@ unsigned invertCubicPolynomial(SolContainer* sol,
Scalar wDisc = q*q/4 + p*p*p/27;
if (wDisc >= 0) { // the positive discriminant case:
// calculate the cube root of - q/2 + sqrt(q^2/4 + p^3/27)
Scalar u = - q/2 + Toolbox::sqrt(wDisc);
if (u < 0) u = - Toolbox::pow(-u, 1.0/3);
else u = Toolbox::pow(u, 1.0/3);
Scalar u = - q/2 + Opm::sqrt(wDisc);
if (u < 0) u = - Opm::pow(-u, 1.0/3);
else u = Opm::pow(u, 1.0/3);
// at this point, u != 0 since p^3 = 0 is necessary in order
// for u = 0 to hold, so
@ -221,10 +214,10 @@ unsigned invertCubicPolynomial(SolContainer* sol,
}
else { // the negative discriminant case:
Scalar uCubedRe = - q/2;
Scalar uCubedIm = Toolbox::sqrt(-wDisc);
Scalar uCubedIm = Opm::sqrt(-wDisc);
// calculate the cube root of - q/2 + sqrt(q^2/4 + p^3/27)
Scalar uAbs = Toolbox::pow(Toolbox::sqrt(uCubedRe*uCubedRe + uCubedIm*uCubedIm), 1.0/3);
Scalar phi = Toolbox::atan2(uCubedIm, uCubedRe)/3;
Scalar uAbs = Opm::pow(Opm::sqrt(uCubedRe*uCubedRe + uCubedIm*uCubedIm), 1.0/3);
Scalar phi = Opm::atan2(uCubedIm, uCubedRe)/3;
// calculate the length and the angle of the primitive root
@ -266,7 +259,7 @@ unsigned invertCubicPolynomial(SolContainer* sol,
// values for phi which differ by 2/3*pi. This allows to
// calculate the three real roots of the polynomial:
for (int i = 0; i < 3; ++i) {
sol[i] = Toolbox::cos(phi)*(uAbs - p/(3*uAbs)) - b/3;
sol[i] = Opm::cos(phi)*(uAbs - p/(3*uAbs)) - b/3;
phi += 2*M_PI/3;
}
@ -282,24 +275,24 @@ unsigned invertCubicPolynomial(SolContainer* sol,
}
// Handle some (pretty unlikely) special cases required to avoid
// divisions by zero in the code above...
else if (std::abs(Toolbox::scalarValue(p)) < 1e-30 && std::abs(Toolbox::scalarValue(q)) < 1e-30) {
else if (std::abs(Opm::scalarValue(p)) < 1e-30 && std::abs(Opm::scalarValue(q)) < 1e-30) {
// t^3 = 0, i.e. triple root at t = 0
sol[0] = sol[1] = sol[2] = 0.0 - b/3;
return 3;
}
else if (std::abs(Toolbox::scalarValue(p)) < 1e-30 && std::abs(Toolbox::scalarValue(q)) > 1e-30) {
else if (std::abs(Opm::scalarValue(p)) < 1e-30 && std::abs(Opm::scalarValue(q)) > 1e-30) {
// t^3 + q = 0,
//
// i. e. single real root at t=curt(q)
Scalar t;
if (-q > 0) t = Toolbox::pow(-q, 1./3);
else t = - Toolbox::pow(q, 1./3);
if (-q > 0) t = Opm::pow(-q, 1./3);
else t = - Opm::pow(q, 1./3);
sol[0] = t - b/3;
return 1;
}
assert(std::abs(Toolbox::scalarValue(p)) > 1e-30 && std::abs(Toolbox::scalarValue(q)) > 1e-30);
assert(std::abs(Opm::scalarValue(p)) > 1e-30 && std::abs(Opm::scalarValue(q)) > 1e-30);
// t^3 + p*t = 0 = t*(t^2 + p),
//
@ -310,9 +303,9 @@ unsigned invertCubicPolynomial(SolContainer* sol,
}
// two additional real roots at t = sqrt(-p) and t = -sqrt(-p)
sol[0] = -Toolbox::sqrt(-p) - b/3;;
sol[0] = -Opm::sqrt(-p) - b/3;;
sol[1] = 0.0 - b/3;
sol[2] = Toolbox::sqrt(-p) - b/3;
sol[2] = Opm::sqrt(-p) - b/3;
return 3;
}

16
opm/material/common/Spline.hpp
View File

@ -799,8 +799,6 @@ public:
template <class Evaluation>
Evaluation eval(const Evaluation& x, bool extrapolate = false) const
{
typedef Opm::MathToolbox<Evaluation> Toolbox;
if (!extrapolate && !applies(x))
OPM_THROW(Opm::NumericalProblem,
"Tried to evaluate a spline outside of its range");
@ -820,7 +818,7 @@ public:
}
}
return eval_(x, segmentIdx_(Toolbox::scalarValue(x)));
return eval_(x, segmentIdx_(Opm::scalarValue(x)));
}
/*!
@ -838,8 +836,6 @@ public:
template <class Evaluation>
Evaluation evalDerivative(const Evaluation& x, bool extrapolate = false) const
{
typedef Opm::MathToolbox<Evaluation> Toolbox;
if (!extrapolate && !applies(x))
OPM_THROW(Opm::NumericalProblem,
"Tried to evaluate the derivative of a spline outside of its range");
@ -852,7 +848,7 @@ public:
return evalDerivative_(xAt(numSamples() - 1), /*segmentIdx=*/numSamples() - 2);
}
return evalDerivative_(x, segmentIdx_(Toolbox::scalarValue(x)));
return evalDerivative_(x, segmentIdx_(Opm::scalarValue(x)));
}
/*!
@ -870,15 +866,13 @@ public:
template <class Evaluation>
Evaluation evalSecondDerivative(const Evaluation& x, bool extrapolate = false) const
{
typedef Opm::MathToolbox<Evaluation> Toolbox;
if (!extrapolate && !applies(x))
OPM_THROW(Opm::NumericalProblem,
"Tried to evaluate the second derivative of a spline outside of its range");
else if (extrapolate)
return 0.0;
return evalDerivative2_(x, segmentIdx_(Toolbox::scalarValue(x)));
return evalDerivative2_(x, segmentIdx_(Opm::scalarValue(x)));
}
/*!
@ -896,15 +890,13 @@ public:
template <class Evaluation>
Evaluation evalThirdDerivative(const Evaluation& x, bool extrapolate = false) const
{
typedef Opm::MathToolbox<Evaluation> Toolbox;
if (!extrapolate && !applies(x))
OPM_THROW(Opm::NumericalProblem,
"Tried to evaluate the third derivative of a spline outside of its range");
else if (extrapolate)
return 0.0;
return evalDerivative3_(x, segmentIdx_(Toolbox::scalarValue(x)));
return evalDerivative3_(x, segmentIdx_(Opm::scalarValue(x)));
}
/*!

8
opm/material/common/Tabulated1DFunction.hpp
View File

@ -251,8 +251,6 @@ public:
template <class Evaluation>
Evaluation eval(const Evaluation& x, bool extrapolate = false) const
{
typedef Opm::MathToolbox<Evaluation> Toolbox;
if (!extrapolate && !applies(x))
OPM_THROW(Opm::NumericalProblem,
"Tried to evaluate a tabulated function outside of its range");
@ -264,7 +262,7 @@ public:
else if (extrapolate && x > xValues_.back())
segIdx = numSamples() - 2;
else
segIdx = findSegmentIndex_(Toolbox::scalarValue(x));
segIdx = findSegmentIndex_(Opm::scalarValue(x));
Scalar x0 = xValues_[segIdx];
Scalar x1 = xValues_[segIdx + 1];
@ -289,15 +287,13 @@ public:
template <class Evaluation>
Evaluation evalDerivative(const Evaluation& x, bool extrapolate = false) const
{
typedef Opm::MathToolbox<Evaluation> Toolbox;
if (!extrapolate && !applies(x)) {
OPM_THROW(Opm::NumericalProblem,
"Tried to evaluate a derivative of a tabulated"
" function outside of its range");
}
unsigned segIdx = findSegmentIndex_(Toolbox::scalarValue(x));
unsigned segIdx = findSegmentIndex_(Opm::scalarValue(x));
return evalDerivative_(x, segIdx);
}

6
opm/material/common/UniformTabulated2DFunction.hpp
View File

@ -181,8 +181,6 @@ public:
template <class Evaluation>
Evaluation eval(const Evaluation& x, const Evaluation& y) const
{
typedef MathToolbox<Evaluation> Toolbox;
#ifndef NDEBUG
if (!applies(x,y))
{
@ -201,11 +199,11 @@ public:
unsigned i =
static_cast<unsigned>(
std::max(0, std::min(static_cast<int>(numX()) - 2,
static_cast<int>(Toolbox::scalarValue(alpha)))));
static_cast<int>(Opm::scalarValue(alpha)))));
unsigned j =
static_cast<unsigned>(
std::max(0, std::min(static_cast<int>(numY()) - 2,
static_cast<int>(Toolbox::scalarValue(beta)))));
static_cast<int>(Opm::scalarValue(beta)))));
alpha -= i;
beta -= j;

10
opm/material/common/UniformXTabulated2DFunction.hpp
View File

@ -249,10 +249,8 @@ public:
template <class Evaluation>
Evaluation eval(const Evaluation& x, const Evaluation& y, bool extrapolate=false) const
{
typedef Opm::MathToolbox<Evaluation> Toolbox;
#ifndef NDEBUG
if (!extrapolate && !applies(Toolbox::scalarValue(x), Toolbox::scalarValue(y))) {
if (!extrapolate && !applies(Opm::scalarValue(x), Opm::scalarValue(y))) {
OPM_THROW(NumericalProblem,
"Attempt to get undefined table value (" << x << ", " << y << ")");
};
@ -263,7 +261,7 @@ public:
Evaluation alpha = xToI(x, extrapolate);
size_t i =
static_cast<size_t>(std::max(0, std::min(static_cast<int>(numX()) - 2,
static_cast<int>(Toolbox::scalarValue(alpha)))));
static_cast<int>(Opm::scalarValue(alpha)))));
alpha -= i;
Evaluation beta1;
@ -273,9 +271,9 @@ public:
beta2 = yToJ(i + 1, y, extrapolate);
size_t j1 = static_cast<size_t>(std::max(0, std::min(static_cast<int>(numY(i)) - 2,
static_cast<int>(Toolbox::scalarValue(beta1)))));
static_cast<int>(Opm::scalarValue(beta1)))));
size_t j2 = static_cast<size_t>(std::max(0, std::min(static_cast<int>(numY(i + 1)) - 2,
static_cast<int>(Toolbox::scalarValue(beta2)))));
static_cast<int>(Opm::scalarValue(beta2)))));
beta1 -= j1;
beta2 -= j2;

38
opm/material/components/Air.hpp
View File

@ -137,8 +137,6 @@ public:
template <class Evaluation>
static Evaluation gasViscosity(const Evaluation& temperature, const Evaluation& /*pressure*/)
{
typedef MathToolbox<Evaluation> Toolbox;
Scalar Tc = criticalTemperature();
Scalar Vc = 84.525138; // critical specific volume [cm^3/mol]
Scalar omega = 0.078; // accentric factor
@ -152,23 +150,21 @@ public:
Scalar Fc = 1 - 0.2756*omega + 0.059035*mu_r4;
Evaluation Tstar = 1.2593 * temperature/Tc;
Evaluation Omega_v =
1.16145*Toolbox::pow(Tstar, -0.14874) +
0.52487*Toolbox::exp(- 0.77320*Tstar) +
2.16178*Toolbox::exp(- 2.43787*Tstar);
return 40.7851e-7*Fc*Toolbox::sqrt(M*temperature)/(std::pow(Vc, 2./3)*Omega_v);
1.16145*Opm::pow(Tstar, -0.14874) +
0.52487*Opm::exp(- 0.77320*Tstar) +
2.16178*Opm::exp(- 2.43787*Tstar);
return 40.7851e-7*Fc*Opm::sqrt(M*temperature)/(std::pow(Vc, 2./3)*Omega_v);
}
// simpler method, from old constrelAir.hh
template <class Evaluation>
static Evaluation simpleGasViscosity(const Evaluation& temperature, const Evaluation& /*pressure*/)
{
typedef MathToolbox<Evaluation> Toolbox;
if(temperature < 273.15 || temperature > 660.) {
OPM_THROW(NumericalProblem,
"Air: Temperature (" << temperature << "K) out of range");
}
return 1.496e-6*Toolbox::pow(temperature, 1.5)/(temperature + 120);
return 1.496e-6*Opm::pow(temperature, 1.5)/(temperature + 120);
}
/*!
@ -252,26 +248,24 @@ public:
static Evaluation gasHeatCapacity(const Evaluation& temperature,
const Evaluation& /*pressure*/)
{
typedef MathToolbox<Evaluation> Toolbox;
// scale temperature by reference temp of 100K
Evaluation phi = temperature/100;
Evaluation c_p =
0.661738E+01
-0.105885E+01 * phi
+0.201650E+00 * Toolbox::pow(phi,2.)
-0.196930E-01 * Toolbox::pow(phi,3.)
+0.106460E-02 * Toolbox::pow(phi,4.)
-0.303284E-04 * Toolbox::pow(phi,5.)
+0.355861E-06 * Toolbox::pow(phi,6.);
+0.201650E+00 * Opm::pow(phi,2.)
-0.196930E-01 * Opm::pow(phi,3.)
+0.106460E-02 * Opm::pow(phi,4.)
-0.303284E-04 * Opm::pow(phi,5.)
+0.355861E-06 * Opm::pow(phi,6.);
c_p +=
-0.549169E+01 * Toolbox::pow(phi,-1.)
+0.585171E+01* Toolbox::pow(phi,-2.)
-0.372865E+01* Toolbox::pow(phi,-3.)
+0.133981E+01* Toolbox::pow(phi,-4.)
-0.233758E+00* Toolbox::pow(phi,-5.)
+0.125718E-01* Toolbox::pow(phi,-6.);
-0.549169E+01 * Opm::pow(phi,-1.)
+0.585171E+01* Opm::pow(phi,-2.)
-0.372865E+01* Opm::pow(phi,-3.)
+0.133981E+01* Opm::pow(phi,-4.)
-0.233758E+00* Opm::pow(phi,-5.)
+0.125718E-01* Opm::pow(phi,-6.);
c_p *= IdealGas::R / (molarMass() * 1000); // in J/mol/K * mol / kg / 1000 = kJ/kg/K
return c_p;

22
opm/material/components/Brine.hpp
View File

@ -136,8 +136,6 @@ public:
static Evaluation liquidEnthalpy(const Evaluation& temperature,
const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
// Numerical coefficents from Palliser and McKibbin
static const Scalar f[] = {
2.63500e-1, 7.48368e-6, 1.44611e-6, -3.80860e-10
@ -157,8 +155,8 @@ public:
const Evaluation& S_lSAT =
f[0]
+ f[1]*theta
+ f[2]*Toolbox::pow(theta, 2)
+ f[3]*Toolbox::pow(theta, 3);
+ f[2]*Opm::pow(theta, 2)
+ f[3]*Opm::pow(theta, 3);
// Regularization
if (S > S_lSAT)
@ -179,7 +177,7 @@ public:
Evaluation d_h = 0;
for (int i = 0; i<=3; ++i) {
for (int j = 0; j <= 2; ++j) {
d_h += a[i][j] * Toolbox::pow(theta, i) * Toolbox::pow(m, j);
d_h += a[i][j] * Opm::pow(theta, i) * Opm::pow(m, j);
}
}
@ -284,18 +282,16 @@ public:
template <class Evaluation>
static Evaluation liquidPressure(const Evaluation& temperature, const Evaluation& density)
{
typedef MathToolbox<Evaluation> Toolbox;
// We use the newton method for this. For the initial value we
// assume the pressure to be 10% higher than the vapor
// pressure
Evaluation pressure = 1.1*vaporPressure(temperature);
Scalar eps = Toolbox::scalarValue(pressure)*1e-7;
Scalar eps = Opm::scalarValue(pressure)*1e-7;
Evaluation deltaP = pressure*2;
for (int i = 0;
i < 5
&& std::abs(Toolbox::scalarValue(pressure)*1e-9) < std::abs(Toolbox::scalarValue(deltaP));
&& std::abs(Opm::scalarValue(pressure)*1e-9) < std::abs(Opm::scalarValue(deltaP));
++i)
{
const Evaluation& f = liquidDensity(temperature, pressure) - density;
@ -330,14 +326,12 @@ public:
template <class Evaluation>
static Evaluation liquidViscosity(const Evaluation& temperature, const Evaluation& /*pressure*/)
{
typedef MathToolbox<Evaluation> Toolbox;
Evaluation T_C = temperature - 273.15;
if(temperature <= 275.) // regularization
T_C = Toolbox::createConstant(275.0);
T_C = 275.0;
Evaluation A = (0.42*std::pow((std::pow(salinity, 0.8)-0.17), 2) + 0.045)*Toolbox::pow(T_C, 0.8);
Evaluation mu_brine = 0.1 + 0.333*salinity + (1.65+91.9*salinity*salinity*salinity)*Toolbox::exp(-A);
Evaluation A = (0.42*std::pow((std::pow(salinity, 0.8)-0.17), 2) + 0.045)*Opm::pow(T_C, 0.8);
Evaluation mu_brine = 0.1 + 0.333*salinity + (1.65+91.9*salinity*salinity*salinity)*Opm::exp(-A);
return mu_brine/1000.0; // convert to [Pa s] (todo: check if correct cP->Pa s is times 10...)
}

22
opm/material/components/CO2.hpp
View File

@ -132,8 +132,6 @@ public:
template <class Evaluation>
static Evaluation vaporPressure(const Evaluation& T)
{
typedef MathToolbox<Evaluation> Toolbox;
static const Scalar a[4] =
{ -7.0602087, 1.9391218, -1.6463597, -3.2995634 };
static const Scalar t[4] =
@ -143,10 +141,10 @@ public:
Evaluation exponent = 0;
Evaluation Tred = T/criticalTemperature();
for (int i = 0; i < 5; ++i)
exponent += a[i]*Toolbox::pow(1 - Tred, t[i]);
exponent += a[i]*Opm::pow(1 - Tred, t[i]);
exponent *= 1.0/Tred;
return Toolbox::exp(exponent)*criticalPressure();
return Opm::exp(exponent)*criticalPressure();
}
@ -203,8 +201,6 @@ public:
template <class Evaluation>
static Evaluation gasViscosity(Evaluation temperature, const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
const Scalar a0 = 0.235156;
const Scalar a1 = -0.491266;
const Scalar a2 = 5.211155e-2;
@ -220,14 +216,14 @@ public:
const Scalar ESP = 251.196;
if(temperature < 275.) // regularization
temperature = Toolbox::createConstant(275.0);
temperature = 275.0;
Evaluation TStar = temperature/ESP;
// mu0: viscosity in zero-density limit
const Evaluation& logTStar = Toolbox::log(TStar);
Evaluation SigmaStar = Toolbox::exp(a0 + logTStar*(a1 + logTStar*(a2 + logTStar*(a3 + logTStar*a4))));
const Evaluation& logTStar = Opm::log(TStar);
Evaluation SigmaStar = Opm::exp(a0 + logTStar*(a1 + logTStar*(a2 + logTStar*(a3 + logTStar*a4))));
Evaluation mu0 = 1.00697*Toolbox::sqrt(temperature) / SigmaStar;
Evaluation mu0 = 1.00697*Opm::sqrt(temperature) / SigmaStar;
const Evaluation& rho = gasDensity(temperature, pressure); // CO2 mass density [kg/m^3]
@ -235,9 +231,9 @@ public:
Evaluation dmu =
d11*rho
+ d21*rho*rho
+ d64*Toolbox::pow(rho, 6.0)/(TStar*TStar*TStar)
+ d81*Toolbox::pow(rho, 8.0)
+ d82*Toolbox::pow(rho, 8.0)/TStar;
+ d64*Opm::pow(rho, 6.0)/(TStar*TStar*TStar)
+ d81*Opm::pow(rho, 8.0)
+ d82*Opm::pow(rho, 8.0)/TStar;
return (mu0 + dmu)/1.0e6; // conversion to [Pa s]
}

29
opm/material/components/H2O.hpp
View File

@ -177,8 +177,6 @@ public:
static Evaluation gasEnthalpy(const Evaluation& temperature,
const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
if (!Region2::isValid(temperature, pressure))
{
OPM_THROW(NumericalProblem,
@ -194,7 +192,7 @@ public:
// dependence on pressure, so we can just return the
// specific enthalpy at the point of regularization, i.e.
// the triple pressure - 100Pa
return enthalpyRegion2_(temperature, Toolbox::createConstant(triplePressure() - 100));
return enthalpyRegion2_<Evaluation>(temperature, triplePressure() - 100);
}
Evaluation pv = vaporPressure(temperature);
if (pressure > pv) {
@ -554,8 +552,6 @@ public:
template <class Evaluation>
static Evaluation gasDensity(const Evaluation& temperature, const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
if (!Region2::isValid(temperature, pressure))
{
OPM_THROW(NumericalProblem,
@ -588,7 +584,7 @@ public:
// calculate the partial derivative of the specific volume
// to the pressure at the vapor pressure.
Scalar eps = Toolbox::scalarValue(pv)*1e-8;
Scalar eps = Opm::scalarValue(pv)*1e-8;
Evaluation v0 = volumeRegion2_(temperature, pv);
Evaluation v1 = volumeRegion2_(temperature, pv + eps);
Evaluation dv_dp = (v1 - v0)/eps;
@ -637,8 +633,6 @@ public:
template <class Evaluation>
static Evaluation gasPressure(const Evaluation& temperature, Scalar density)
{
typedef MathToolbox<Evaluation> Toolbox;
Valgrind::CheckDefined(temperature);
Valgrind::CheckDefined(density);
@ -650,7 +644,7 @@ public:
Evaluation deltaP = pressure*2;
Valgrind::CheckDefined(pressure);
Valgrind::CheckDefined(deltaP);
for (int i = 0; i < 5 && std::abs(Toolbox::scalarValue(pressure)*1e-9) < std::abs(Toolbox::scalarValue(deltaP)); ++i) {
for (int i = 0; i < 5 && std::abs(Opm::scalarValue(pressure)*1e-9) < std::abs(Opm::scalarValue(deltaP)); ++i) {
Evaluation f = gasDensity(temperature, pressure) - density;
Evaluation df_dp;
@ -683,8 +677,6 @@ public:
template <class Evaluation>
static Evaluation liquidDensity(const Evaluation& temperature, const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
if (!Region1::isValid(temperature, pressure))
{
OPM_THROW(NumericalProblem,
@ -700,7 +692,7 @@ public:
// calculate the partial derivative of the specific volume
// to the pressure at the vapor pressure.
Scalar eps = Toolbox::scalarValue(pv)*1e-8;
Scalar eps = Opm::scalarValue(pv)*1e-8;
Evaluation v0 = volumeRegion1_(temperature, pv);
Evaluation v1 = volumeRegion1_(temperature, pv + eps);
Evaluation dv_dp = (v1 - v0)/eps;
@ -746,16 +738,14 @@ public:
template <class Evaluation>
static Evaluation liquidPressure(const Evaluation& temperature, Scalar density)
{
typedef MathToolbox<Evaluation> Toolbox;
// We use the Newton method for this. For the initial value we
// assume the pressure to be 10% higher than the vapor
// pressure
Evaluation pressure = 1.1*vaporPressure(temperature);
Scalar eps = Toolbox::scalarValue(pressure)*1e-7;
Scalar eps = Opm::scalarValue(pressure)*1e-7;
Evaluation deltaP = pressure*2;
for (int i = 0; i < 5 && std::abs(Toolbox::scalarValue(pressure)*1e-9) < std::abs(Toolbox::scalarValue(deltaP)); ++i) {
for (int i = 0; i < 5 && std::abs(Opm::scalarValue(pressure)*1e-9) < std::abs(Opm::scalarValue(deltaP)); ++i) {
Evaluation f = liquidDensity(temperature, pressure) - density;
Evaluation df_dp;
@ -879,9 +869,8 @@ private:
template <class Evaluation>
static Evaluation heatCap_p_Region1_(const Evaluation& temperature, const Evaluation& pressure)
{
typedef Opm::MathToolbox<Evaluation> Toolbox;
return
- Toolbox::pow(Region1::tau(temperature), 2.0) *
- Opm::pow(Region1::tau(temperature), 2.0) *
Region1::ddgamma_ddtau(temperature, pressure) *
Rs;
}
@ -944,10 +933,8 @@ private:
template <class Evaluation>
static Evaluation heatCap_p_Region2_(const Evaluation& temperature, const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
return
- Toolbox::pow(Region2::tau(temperature), 2 ) *
- Opm::pow(Region2::tau(temperature), 2 ) *
Region2::ddgamma_ddtau(temperature, pressure) *
Rs;
}

42
opm/material/components/Mesitylene.hpp
View File

@ -98,15 +98,13 @@ public:
template <class Evaluation>
static Evaluation vaporPressure(const Evaluation& temperature)
{
typedef MathToolbox<Evaluation> Toolbox;
const Scalar A = 7.07638;
const Scalar B = 1571.005;
const Scalar C = 209.728;
const Evaluation& T = temperature - 273.15;
return 100 * 1.334 * Toolbox::pow(10.0, A - (B / (T + C)));
return 100 * 1.334 * Opm::pow(10.0, A - (B / (T + C)));
}
@ -145,10 +143,8 @@ public:
template <class Evaluation>
static Evaluation heatVap(const Evaluation& temperature, const Evaluation& /*pressure*/)
{
typedef MathToolbox<Evaluation> Toolbox;
Evaluation T = Toolbox::min(temperature, criticalTemperature()); // regularization
T = Toolbox::max(T, 0.0); // regularization
Evaluation T = Opm::min(temperature, criticalTemperature()); // regularization
T = Opm::max(T, 0.0); // regularization
const Scalar T_crit = criticalTemperature();
const Scalar Tr1 = boilingTemperature()/criticalTemperature();
@ -163,7 +159,7 @@ public:
/* Variation with temp according to Watson relation eq 7-12.1*/
const Evaluation& Tr2 = T/criticalTemperature();
const Scalar n = 0.375;
const Evaluation& DH_vap = DH_v_boil * Toolbox::pow(((1.0 - Tr2)/(1.0 - Tr1)), n);
const Evaluation& DH_vap = DH_v_boil * Opm::pow(((1.0 - Tr2)/(1.0 - Tr1)), n);
return (DH_vap/molarMass()); // we need [J/kg]
}
@ -232,10 +228,8 @@ public:
template <class Evaluation>
static Evaluation gasViscosity(Evaluation temperature, const Evaluation& /*pressure*/, bool /*regularize*/=true)
{
typedef MathToolbox<Evaluation> Toolbox;
temperature = Toolbox::min(temperature, 500.0); // regularization
temperature = Toolbox::max(temperature, 250.0);
temperature = Opm::min(temperature, 500.0); // regularization
temperature = Opm::max(temperature, 250.0);
// reduced temperature
const Evaluation& Tr = temperature/criticalTemperature();
@ -243,9 +237,9 @@ public:
Scalar Fp0 = 1.0;
Scalar xi = 0.00474;
const Evaluation& eta_xi =
Fp0*(0.807*Toolbox::pow(Tr,0.618)
- 0.357*Toolbox::exp(-0.449*Tr)
+ 0.34*Toolbox::exp(-4.058*Tr)
Fp0*(0.807*Opm::pow(Tr,0.618)
- 0.357*Opm::exp(-0.449*Tr)
+ 0.34*Opm::exp(-4.058*Tr)
+ 0.018);
return eta_xi/xi/1e7; // [Pa s]
@ -260,15 +254,13 @@ public:
template <class Evaluation>
static Evaluation liquidViscosity(Evaluation temperature, const Evaluation& /*pressure*/)
{
typedef MathToolbox<Evaluation> Toolbox;
temperature = Toolbox::min(temperature, 500.0); // regularization
temperature = Toolbox::max(temperature, 250.0);
temperature = Opm::min(temperature, 500.0); // regularization
temperature = Opm::max(temperature, 250.0);
const Scalar A = -6.749;
const Scalar B = 2010.0;
return Toolbox::exp(A + B/temperature)*1e-3; // [Pa s]
return Opm::exp(A + B/temperature)*1e-3; // [Pa s]
}
/*!
@ -328,14 +320,12 @@ protected:
template <class Evaluation>
static Evaluation molarLiquidDensity_(Evaluation temperature)
{
typedef MathToolbox<Evaluation> Toolbox;
temperature = Toolbox::min(temperature, 500.0); // regularization
temperature = Toolbox::max(temperature, 250.0);
temperature = Opm::min(temperature, 500.0); // regularization
temperature = Opm::max(temperature, 250.0);
const Scalar Z_RA = 0.2556; // from equation
const Evaluation& expo = 1.0 + Toolbox::pow(1.0 - temperature/criticalTemperature(), 2.0/7.0);
const Evaluation& V = Consts::R*criticalTemperature()/criticalPressure()*Toolbox::pow(Z_RA, expo); // liquid molar volume [cm^3/mol]
const Evaluation& expo = 1.0 + Opm::pow(1.0 - temperature/criticalTemperature(), 2.0/7.0);
const Evaluation& V = Consts::R*criticalTemperature()/criticalPressure()*Opm::pow(Z_RA, expo); // liquid molar volume [cm^3/mol]
return 1.0/V; // molar density [mol/m^3]
}

16
opm/material/components/N2.hpp
View File

@ -103,8 +103,6 @@ public:
template <class Evaluation>
static Evaluation vaporPressure(const Evaluation& temperature)
{
typedef MathToolbox<Evaluation> Toolbox;
if (temperature > criticalTemperature())
return criticalPressure();
if (temperature < tripleTemperature())
@ -113,14 +111,14 @@ public:
// note: this is the ancillary equation given on page 1368
const Evaluation& sigma = 1.0 - temperature/criticalTemperature();
const Evaluation& sqrtSigma = Toolbox::sqrt(sigma);
const Evaluation& sqrtSigma = Opm::sqrt(sigma);
const Scalar N1 = -6.12445284;
const Scalar N2 = 1.26327220;
const Scalar N3 = -0.765910082;
const Scalar N4 = -1.77570564;
return
criticalPressure() *
Toolbox::exp(criticalTemperature()/temperature*
Opm::exp(criticalTemperature()/temperature*
(sigma*(N1 +
sqrtSigma*N2 +
sigma*(sqrtSigma*N3 +
@ -267,8 +265,6 @@ public:
template <class Evaluation>
static Evaluation gasViscosity(const Evaluation& temperature, const Evaluation& /*pressure*/)
{
typedef MathToolbox<Evaluation> Toolbox;
const Scalar Tc = criticalTemperature();
const Scalar Vc = 90.1; // critical specific volume [cm^3/mol]
const Scalar omega = 0.037; // accentric factor
@ -282,10 +278,10 @@ public:
Scalar Fc = 1 - 0.2756*omega + 0.059035*mu_r4;
const Evaluation& Tstar = 1.2593 * temperature/Tc;
const Evaluation& Omega_v =
1.16145*Toolbox::pow(Tstar, -0.14874) +
0.52487*Toolbox::exp(- 0.77320*Tstar) +
2.16178*Toolbox::exp(- 2.43787*Tstar);
const Evaluation& mu = 40.785*Fc*Toolbox::sqrt(M*temperature)/(std::pow(Vc, 2./3)*Omega_v);
1.16145*Opm::pow(Tstar, -0.14874) +
0.52487*Opm::exp(- 0.77320*Tstar) +
2.16178*Opm::exp(- 2.43787*Tstar);
const Evaluation& mu = 40.785*Fc*Opm::sqrt(M*temperature)/(std::pow(Vc, 2./3)*Omega_v);
// convertion from micro poise to Pa s
return mu/1e6 / 10;

4
opm/material/components/SimpleH2O.hpp
View File

@ -127,8 +127,6 @@ public:
template <class Evaluation>
static Evaluation vaporPressure(const Evaluation& T)
{
typedef Opm::MathToolbox<Evaluation> Toolbox;
if (T > criticalTemperature())
return criticalPressure();
if (T < tripleTemperature())
@ -147,7 +145,7 @@ public:
Evaluation B = (n[2]*sigma + n[3])*sigma + n[4];
Evaluation C = (n[5]*sigma + n[6])*sigma + n[7];
Evaluation tmp = 2.0*C/(Toolbox::sqrt(B*B - 4.0*A*C) - B);
Evaluation tmp = 2.0*C/(Opm::sqrt(B*B - 4.0*A*C) - B);
tmp *= tmp;
tmp *= tmp;

76
opm/material/components/TabulatedComponent.hpp
View File

@ -257,10 +257,8 @@ public:
template <class Evaluation>
static Evaluation vaporPressure(const Evaluation& temperature)
{
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& result = interpolateT_(vaporPressure_, temperature);
if (std::isnan(Toolbox::scalarValue(result)))
if (std::isnan(Opm::scalarValue(result)))
return RawComponent::vaporPressure(temperature);
return result;
}
@ -274,12 +272,10 @@ public:
template <class Evaluation>
static Evaluation gasEnthalpy(const Evaluation& temperature, const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& result = interpolateGasTP_(gasEnthalpy_,
temperature,
pressure);
if (std::isnan(Toolbox::scalarValue(result)))
if (std::isnan(Opm::scalarValue(result)))
return RawComponent::gasEnthalpy(temperature, pressure);
return result;
}
@ -293,12 +289,10 @@ public:
template <class Evaluation>
static Evaluation liquidEnthalpy(const Evaluation& temperature, const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& result = interpolateLiquidTP_(liquidEnthalpy_,
temperature,
pressure);
if (std::isnan(Toolbox::scalarValue(result)))
if (std::isnan(Opm::scalarValue(result)))
return RawComponent::liquidEnthalpy(temperature, pressure);
return result;
}
@ -312,12 +306,10 @@ public:
template <class Evaluation>
static Evaluation gasHeatCapacity(const Evaluation& temperature, const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& result = interpolateGasTP_(gasHeatCapacity_,
temperature,
pressure);
if (std::isnan(Toolbox::scalarValue(result)))
if (std::isnan(Opm::scalarValue(result)))
return RawComponent::gasHeatCapacity(temperature, pressure);
return result;
}
@ -331,12 +323,10 @@ public:
template <class Evaluation>
static Evaluation liquidHeatCapacity(const Evaluation& temperature, const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& result = interpolateLiquidTP_(liquidHeatCapacity_,
temperature,
pressure);
if (std::isnan(Toolbox::scalarValue(result)))
if (std::isnan(Opm::scalarValue(result)))
return RawComponent::liquidHeatCapacity(temperature, pressure);
return result;
}
@ -370,12 +360,10 @@ public:
template <class Evaluation>
static Evaluation gasPressure(const Evaluation& temperature, Scalar density)
{
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& result = interpolateGasTRho_(gasPressure_,
temperature,
density);
if (std::isnan(Toolbox::scalarValue(result)))
if (std::isnan(Opm::scalarValue(result)))
return RawComponent::gasPressure(temperature,
density);
return result;
@ -390,12 +378,10 @@ public:
template <class Evaluation>
static Evaluation liquidPressure(const Evaluation& temperature, Scalar density)
{
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& result = interpolateLiquidTRho_(liquidPressure_,
temperature,
density);
if (std::isnan(Toolbox::scalarValue(result)))
if (std::isnan(Opm::scalarValue(result)))
return RawComponent::liquidPressure(temperature,
density);
return result;
@ -430,12 +416,10 @@ public:
template <class Evaluation>
static Evaluation gasDensity(const Evaluation& temperature, const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& result = interpolateGasTP_(gasDensity_,
temperature,
pressure);
if (std::isnan(Toolbox::scalarValue(result)))
if (std::isnan(Opm::scalarValue(result)))
return RawComponent::gasDensity(temperature, pressure);
return result;
}
@ -450,12 +434,10 @@ public:
template <class Evaluation>
static Evaluation liquidDensity(const Evaluation& temperature, const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& result = interpolateLiquidTP_(liquidDensity_,
temperature,
pressure);
if (std::isnan(Toolbox::scalarValue(result)))
if (std::isnan(Opm::scalarValue(result)))
return RawComponent::liquidDensity(temperature, pressure);
return result;
}
@ -469,12 +451,10 @@ public:
template <class Evaluation>
static Evaluation gasViscosity(const Evaluation& temperature, const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& result = interpolateGasTP_(gasViscosity_,
temperature,
pressure);
if (std::isnan(Toolbox::scalarValue(result)))
if (std::isnan(Opm::scalarValue(result)))
return RawComponent::gasViscosity(temperature, pressure);
return result;
}
@ -488,12 +468,10 @@ public:
template <class Evaluation>
static Evaluation liquidViscosity(const Evaluation& temperature, const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& result = interpolateLiquidTP_(liquidViscosity_,
temperature,
pressure);
if (std::isnan(Toolbox::scalarValue(result)))
if (std::isnan(Opm::scalarValue(result)))
return RawComponent::liquidViscosity(temperature, pressure);
return result;
}
@ -507,12 +485,10 @@ public:
template <class Evaluation>
static Evaluation gasThermalConductivity(const Evaluation& temperature, const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& result = interpolateGasTP_(gasThermalConductivity_,
temperature,
pressure);
if (std::isnan(Toolbox::scalarValue(result)))
if (std::isnan(Opm::scalarValue(result)))
return RawComponent::gasThermalConductivity(temperature, pressure);
return result;
}
@ -526,12 +502,10 @@ public:
template <class Evaluation>
static Evaluation liquidThermalConductivity(const Evaluation& temperature, const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& result = interpolateLiquidTP_(liquidThermalConductivity_,
temperature,
pressure);
if (std::isnan(Toolbox::scalarValue(result)))
if (std::isnan(Opm::scalarValue(result)))
return RawComponent::liquidThermalConductivity(temperature, pressure);
return result;
}
@ -541,13 +515,11 @@ private:
template <class Evaluation>
static Evaluation interpolateT_(const Scalar* values, const Evaluation& T)
{
typedef Opm::MathToolbox<Evaluation> Toolbox;
Evaluation alphaT = tempIdx_(T);
if (alphaT < 0 || alphaT >= nTemp_ - 1)
return std::numeric_limits<Scalar>::quiet_NaN();
size_t iT = static_cast<size_t>(Toolbox::scalarValue(alphaT));
size_t iT = static_cast<size_t>(Opm::scalarValue(alphaT));
alphaT -= iT;
return
@ -560,13 +532,11 @@ private:
template <class Evaluation>
static Evaluation interpolateLiquidTP_(const Scalar* values, const Evaluation& T, const Evaluation& p)
{
typedef MathToolbox<Evaluation> Toolbox;
Evaluation alphaT = tempIdx_(T);
if (alphaT < 0 || alphaT >= nTemp_ - 1)
return Toolbox::createConstant(std::numeric_limits<Scalar>::quiet_NaN());
return std::numeric_limits<Scalar>::quiet_NaN();
size_t iT = static_cast<size_t>(Toolbox::scalarValue(alphaT));
size_t iT = static_cast<size_t>(Opm::scalarValue(alphaT));
alphaT -= iT;
Evaluation alphaP1 = pressLiquidIdx_(p, iT);
@ -575,11 +545,11 @@ private:
size_t iP1 =
static_cast<size_t>(
std::max<int>(0, std::min(static_cast<int>(nPress_) - 2,
static_cast<int>(Toolbox::scalarValue(alphaP1)))));
static_cast<int>(Opm::scalarValue(alphaP1)))));
size_t iP2 =
static_cast<size_t>(
std::max(0, std::min(static_cast<int>(nPress_) - 2,
static_cast<int>(Toolbox::scalarValue(alphaP2)))));
static_cast<int>(Opm::scalarValue(alphaP2)))));
alphaP1 -= iP1;
alphaP2 -= iP2;
@ -595,16 +565,14 @@ private:
template <class Evaluation>
static Evaluation interpolateGasTP_(const Scalar* values, const Evaluation& T, const Evaluation& p)
{
typedef MathToolbox<Evaluation> Toolbox;
Evaluation alphaT = tempIdx_(T);
if (alphaT < 0 || alphaT >= nTemp_ - 1)
return Toolbox::createConstant(std::numeric_limits<Scalar>::quiet_NaN());
return std::numeric_limits<Scalar>::quiet_NaN();
size_t iT =
static_cast<size_t>(
std::max(0, std::min(static_cast<int>(nTemp_) - 2,
static_cast<int>(Toolbox::scalarValue(alphaT)))));
static_cast<int>(Opm::scalarValue(alphaT)))));
alphaT -= iT;
Evaluation alphaP1 = pressGasIdx_(p, iT);
@ -612,11 +580,11 @@ private:
size_t iP1 =
static_cast<size_t>(
std::max(0, std::min(static_cast<int>(nPress_) - 2,
static_cast<int>(Toolbox::scalarValue(alphaP1)))));
static_cast<int>(Opm::scalarValue(alphaP1)))));
size_t iP2 =
static_cast<size_t>(
std::max(0, std::min(static_cast<int>(nPress_) - 2,
static_cast<int>(Toolbox::scalarValue(alphaP2)))));
static_cast<int>(Opm::scalarValue(alphaP2)))));
alphaP1 -= iP1;
alphaP2 -= iP2;

44
opm/material/components/Xylene.hpp
View File

@ -94,13 +94,11 @@ public:
template <class Evaluation>
static Evaluation vaporPressure(const Evaluation& temperature)
{
typedef MathToolbox<Evaluation> Toolbox;
const Scalar A = 7.00909;;
const Scalar B = 1462.266;;
const Scalar C = 215.110;;
return 100*1.334*Toolbox::pow(10.0, (A - (B/(temperature - 273.15 + C))));
return 100*1.334*Opm::pow(10.0, (A - (B/(temperature - 273.15 + C))));
}
/*!
@ -183,10 +181,8 @@ public:
static Evaluation heatVap(Evaluation temperature,
const Evaluation& /*pressure*/)
{
typedef MathToolbox<Evaluation> Toolbox;
temperature = Toolbox::min(temperature, criticalTemperature()); // regularization
temperature = Toolbox::max(temperature, 0.0); // regularization
temperature = Opm::min(temperature, criticalTemperature()); // regularization
temperature = Opm::max(temperature, 0.0); // regularization
const Scalar T_crit = criticalTemperature();
const Scalar Tr1 = boilingTemperature()/criticalTemperature();
@ -200,7 +196,7 @@ public:
/* Variation with temperature according to Watson relation eq 7-12.1*/
const Evaluation& Tr2 = temperature/criticalTemperature();
const Scalar n = 0.375;
const Evaluation& DH_vap = DH_v_boil * Toolbox::pow(((1.0 - Tr2)/(1.0 - Tr1)), n);
const Evaluation& DH_vap = DH_v_boil * Opm::pow(((1.0 - Tr2)/(1.0 - Tr1)), n);
return (DH_vap/molarMass()); // we need [J/kg]
}
@ -250,18 +246,16 @@ public:
template <class Evaluation>
static Evaluation molarLiquidDensity(Evaluation temperature, const Evaluation& /*pressure*/)
{
typedef MathToolbox<Evaluation> Toolbox;
// saturated molar volume according to Lide, CRC Handbook of
// Thermophysical and Thermochemical Data, CRC Press, 1994
// valid for 245 < Temperature < 600
temperature = Toolbox::min(temperature, 500.0); // regularization
temperature = Toolbox::max(temperature, 250.0); // regularization
temperature = Opm::min(temperature, 500.0); // regularization
temperature = Opm::max(temperature, 250.0); // regularization
const Scalar A1 = 0.25919; // from table
const Scalar A2 = 0.0014569; // from table
const Evaluation& expo = 1.0 + Toolbox::pow((1.0 - temperature/criticalTemperature()), (2.0/7.0));
return 1.0/(A2*Toolbox::pow(A1, expo));
const Evaluation& expo = 1.0 + Opm::pow((1.0 - temperature/criticalTemperature()), (2.0/7.0));
return 1.0/(A2*Opm::pow(A1, expo));
}
/*!
@ -295,17 +289,15 @@ public:
template <class Evaluation>
static Evaluation gasViscosity(Evaluation temperature, const Evaluation& /*pressure*/)
{
typedef MathToolbox<Evaluation> Toolbox;
temperature = Opm::min(temperature, 500.0); // regularization
temperature = Opm::max(temperature, 250.0); // regularization
temperature = Toolbox::min(temperature, 500.0); // regularization
temperature = Toolbox::max(temperature, 250.0); // regularization
const Evaluation& Tr = Toolbox::max(temperature/criticalTemperature(), 1e-10);
const Evaluation& Tr = Opm::max(temperature/criticalTemperature(), 1e-10);
const Scalar Fp0 = 1.0;
const Scalar xi = 0.004623;
const Evaluation& eta_xi = Fp0*(0.807*Toolbox::pow(Tr, 0.618)
- 0.357*Toolbox::exp(-0.449*Tr)
+ 0.34*Toolbox::exp(-4.058*Tr)
const Evaluation& eta_xi = Fp0*(0.807*Opm::pow(Tr, 0.618)
- 0.357*Opm::exp(-0.449*Tr)
+ 0.34*Opm::exp(-4.058*Tr)
+ 0.018);
return eta_xi/xi / 1e7; // [Pa s]
}
@ -316,17 +308,15 @@ public:
template <class Evaluation>
static Evaluation liquidViscosity(Evaluation temperature, const Evaluation& /*pressure*/)
{
typedef MathToolbox<Evaluation> Toolbox;
temperature = Toolbox::min(temperature, 500.0); // regularization
temperature = Toolbox::max(temperature, 250.0); // regularization
temperature = Opm::min(temperature, 500.0); // regularization
temperature = Opm::max(temperature, 250.0); // regularization
const Scalar A = -3.82;
const Scalar B = 1027.0;
const Scalar C = -6.38e-4;
const Scalar D = 4.52e-7;
return 1e-3*Toolbox::exp(A
return 1e-3*Opm::exp(A
+ B/temperature
+ C*temperature
+ D*temperature*temperature); // in [cP]

28
opm/material/components/iapws/Common.hpp
View File

@ -98,8 +98,6 @@ public:
template <class Evaluation>
static Evaluation viscosity(const Evaluation& temperature, const Evaluation& rho)
{
typedef MathToolbox<Evaluation> Toolbox;
Evaluation rhoBar = rho/322.0;
Evaluation TBar = temperature/criticalTemperature;
@ -126,10 +124,10 @@ public:
tmp3 *= 1.0/TBar - 1;
};
muBar *= rhoBar;
muBar = Toolbox::exp(muBar);
muBar = Opm::exp(muBar);
// muBar *= muBar_0
muBar *= 100*Toolbox::sqrt(TBar);
muBar *= 100*Opm::sqrt(TBar);
const Scalar H[4] = {
1.67752, 2.20462, 0.6366564, -0.241605
};
@ -160,8 +158,6 @@ public:
template <class Evaluation>
static Evaluation thermalConductivityIAPWS(const Evaluation& T, const Evaluation& rho)
{
typedef MathToolbox<Evaluation> Toolbox;
static const Scalar thcond_tstar = 647.26 ;
static const Scalar thcond_rhostar = 317.7 ;
/*static const Scalar thcond_kstar = 1.0 ;*/
@ -195,7 +191,7 @@ public:
Evaluation rhobar = rho / thcond_rhostar;
/* fast implementation... minimised calls to 'pow' routine... */
Evaluation Troot = Toolbox::sqrt(Tbar);
Evaluation Troot = Opm::sqrt(Tbar);
Evaluation Tpow = Troot;
Evaluation lam = 0;
@ -207,10 +203,10 @@ public:
lam +=
thcond_b0 + thcond_b1
* rhobar + thcond_b2
* Toolbox::exp(thcond_B1 * ((rhobar + thcond_B2)*(rhobar + thcond_B2)));
* Opm::exp(thcond_B1 * ((rhobar + thcond_B2)*(rhobar + thcond_B2)));
Evaluation DTbar = Toolbox::abs(Tbar - 1) + thcond_c4;
Evaluation DTbarpow = Toolbox::pow(DTbar, 3./5);
Evaluation DTbar = Opm::abs(Tbar - 1) + thcond_c4;
Evaluation DTbarpow = Opm::pow(DTbar, 3./5);
Evaluation Q = 2. + thcond_c5 / DTbarpow;
Evaluation S;
@ -219,16 +215,16 @@ public:
else
S = thcond_c6 / DTbarpow;
Evaluation rhobar18 = Toolbox::pow(rhobar, 1.8);
Evaluation rhobarQ = Toolbox::pow(rhobar, Q);
Evaluation rhobar18 = Opm::pow(rhobar, 1.8);
Evaluation rhobarQ = Opm::pow(rhobar, Q);
lam +=
(thcond_d1 / Toolbox::pow(Tbar,10.0) + thcond_d2) * rhobar18 *
Toolbox::exp(thcond_c1 * (1 - rhobar * rhobar18))
(thcond_d1 / Opm::pow(Tbar,10.0) + thcond_d2) * rhobar18 *
Opm::exp(thcond_c1 * (1 - rhobar * rhobar18))
+ thcond_d3 * S * rhobarQ *
Toolbox::exp((Q/(1+Q))*(1 - rhobar*rhobarQ))
Opm::exp((Q/(1+Q))*(1 - rhobar*rhobarQ))
+ thcond_d4 *
Toolbox::exp(thcond_c2 * Toolbox::pow(Troot,3.0) + thcond_c3 / Toolbox::pow(rhobar,5.0));
Opm::exp(thcond_c2 * Opm::pow(Troot,3.0) + thcond_c3 / Opm::pow(rhobar,5.0));
return /*thcond_kstar * */ lam;
}
};

50
opm/material/components/iapws/Region1.hpp
View File

@ -60,11 +60,7 @@ public:
template <class Evaluation>
static bool isValid(const Evaluation& temperature, const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
return
Toolbox::scalarValue(temperature) <= 623.15 &&
Toolbox::scalarValue(pressure) <= 100e6;
return temperature <= 623.15 && pressure <= 100e6;
// actually this is:
/*
@ -139,14 +135,12 @@ public:
template <class Evaluation>
static Evaluation gamma(const Evaluation& temperature, const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation tau_ = tau(temperature); /* reduced temperature */
const Evaluation pi_ = pi(pressure); /* reduced pressure */
Evaluation result = 0;
for (int i = 0; i < 34; ++i) {
result += n(i)*Toolbox::pow(7.1 - pi_, I(i))*Toolbox::pow(tau_ - 1.222, J(i));
result += n(i)*Opm::pow(7.1 - pi_, I(i))*Opm::pow(tau_ - 1.222, J(i));
}
return result;
@ -167,17 +161,15 @@ public:
template <class Evaluation>
static Evaluation dgamma_dtau(const Evaluation& temperature, const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation tau_ = tau(temperature); /* reduced temperature */
const Evaluation pi_ = pi(pressure); /* reduced pressure */
Evaluation result = Toolbox::createConstant(0.0);
Evaluation result = 0.0;
for (int i = 0; i < 34; i++) {
result +=
n(i) *
Toolbox::pow(7.1 - pi_, static_cast<Scalar>(I(i))) *
Toolbox::pow(tau_ - 1.222, static_cast<Scalar>(J(i)-1)) *
Opm::pow(7.1 - pi_, static_cast<Scalar>(I(i))) *
Opm::pow(tau_ - 1.222, static_cast<Scalar>(J(i)-1)) *
J(i);
}
@ -198,18 +190,16 @@ public:
template <class Evaluation>
static Evaluation dgamma_dpi(const Evaluation& temperature, const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation tau_ = tau(temperature); /* reduced temperature */
const Evaluation pi_ = pi(pressure); /* reduced pressure */
Evaluation result = Toolbox::createConstant(0.0);
Evaluation result = 0.0;
for (int i = 0; i < 34; i++) {
result +=
-n(i) *
I(i) *
Toolbox::pow(7.1 - pi_, static_cast<Scalar>(I(i) - 1)) *
Toolbox::pow(tau_ - 1.222, static_cast<Scalar>(J(i)));
Opm::pow(7.1 - pi_, static_cast<Scalar>(I(i) - 1)) *
Opm::pow(tau_ - 1.222, static_cast<Scalar>(J(i)));
}
return result;
@ -230,19 +220,17 @@ public:
template <class Evaluation>
static Evaluation ddgamma_dtaudpi(const Evaluation& temperature, const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation tau_ = tau(temperature); /* reduced temperature */
const Evaluation pi_ = pi(pressure); /* reduced pressure */
Evaluation result = Toolbox::createConstant(0.0);
Evaluation result = 0.0;
for (int i = 0; i < 34; i++) {
result +=
-n(i) *
I(i) *
J(i) *
Toolbox::pow(7.1 - pi_, static_cast<Scalar>(I(i) - 1)) *
Toolbox::pow(tau_ - 1.222, static_cast<Scalar>(J(i) - 1));
Opm::pow(7.1 - pi_, static_cast<Scalar>(I(i) - 1)) *
Opm::pow(tau_ - 1.222, static_cast<Scalar>(J(i) - 1));
}
return result;
@ -263,19 +251,17 @@ public:
template <class Evaluation>
static Evaluation ddgamma_ddpi(const Evaluation& temperature, const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation tau_ = tau(temperature); /* reduced temperature */
const Evaluation pi_ = pi(pressure); /* reduced pressure */
Evaluation result = Toolbox::createConstant(0.0);
Evaluation result = 0.0;
for (int i = 0; i < 34; i++) {
result +=
n(i) *
I(i) *
(I(i) - 1) *
Toolbox::pow(7.1 - pi_, I(i) - 2) *
Toolbox::pow(tau_ - 1.222, J(i));
Opm::pow(7.1 - pi_, I(i) - 2) *
Opm::pow(tau_ - 1.222, J(i));
}
return result;
@ -295,19 +281,17 @@ public:
template <class Evaluation>
static Evaluation ddgamma_ddtau(const Evaluation& temperature, const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation tau_ = tau(temperature); /* reduced temperature */
const Evaluation pi_ = pi(pressure); /* reduced pressure */
Evaluation result = Toolbox::createConstant(0.0);
Evaluation result = 0.0;
for (int i = 0; i < 34; i++) {
result +=
n(i) *
Toolbox::pow(7.1 - pi_, I(i)) *
Opm::pow(7.1 - pi_, I(i)) *
J(i) *
(J(i) - 1) *
Toolbox::pow(tau_ - 1.222, J(i) - 2);
Opm::pow(tau_ - 1.222, J(i) - 2);
}
return result;

48
opm/material/components/iapws/Region2.hpp
View File

@ -135,8 +135,6 @@ public:
template <class Evaluation>
static Evaluation gamma(const Evaluation& temperature, const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& tau_ = tau(temperature); /* reduced temperature */
const Evaluation& pi_ = pi(pressure); /* reduced pressure */
@ -145,14 +143,14 @@ public:
// ideal gas part
result = Opm::log(pi_);
for (int i = 0; i < 9; ++i)
result += n_g(i)*Toolbox::pow(tau_, J_g(i));
result += n_g(i)*Opm::pow(tau_, J_g(i));
// residual part
for (int i = 0; i < 43; ++i)
result +=
n_r(i)*
Toolbox::pow(pi_, I_r(i))*
Toolbox::pow(tau_ - 0.5, J_r(i));
Opm::pow(pi_, I_r(i))*
Opm::pow(tau_ - 0.5, J_r(i));
return result;
}
@ -171,27 +169,25 @@ public:
template <class Evaluation>
static Evaluation dgamma_dtau(const Evaluation& temperature, const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& tau_ = tau(temperature); /* reduced temperature */
const Evaluation& pi_ = pi(pressure); /* reduced pressure */
// ideal gas part
Evaluation result = Toolbox::createConstant(0.0);
Evaluation result = 0.0;
for (int i = 0; i < 9; i++) {
result +=
n_g(i) *
J_g(i) *
Toolbox::pow(tau_, static_cast<Scalar>(J_g(i) - 1));
Opm::pow(tau_, static_cast<Scalar>(J_g(i) - 1));
}
// residual part
for (int i = 0; i < 43; i++) {
result +=
n_r(i) *
Toolbox::pow(pi_, static_cast<Scalar>(I_r(i))) *
Opm::pow(pi_, static_cast<Scalar>(I_r(i))) *
J_r(i) *
Toolbox::pow(tau_ - 0.5, static_cast<Scalar>(J_r(i) - 1));
Opm::pow(tau_ - 0.5, static_cast<Scalar>(J_r(i) - 1));
}
return result;
@ -212,8 +208,6 @@ public:
template <class Evaluation>
static Evaluation dgamma_dpi(const Evaluation& temperature, const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& tau_ = tau(temperature); /* reduced temperature */
const Evaluation& pi_ = pi(pressure); /* reduced pressure */
@ -225,8 +219,8 @@ public:
result +=
n_r(i) *
I_r(i) *
Toolbox::pow(pi_, static_cast<Scalar>(I_r(i) - 1)) *
Toolbox::pow(tau_ - 0.5, static_cast<Scalar>(J_r(i)));
Opm::pow(pi_, static_cast<Scalar>(I_r(i) - 1)) *
Opm::pow(tau_ - 0.5, static_cast<Scalar>(J_r(i)));
}
return result;
@ -247,13 +241,11 @@ public:
template <class Evaluation>
static Evaluation ddgamma_dtaudpi(const Evaluation& temperature, const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& tau_ = tau(temperature); /* reduced temperature */
const Evaluation& pi_ = pi(pressure); /* reduced pressure */
// ideal gas part
Evaluation result = Toolbox::createConstant(0.0);
Evaluation result = 0.0;
// residual part
for (int i = 0; i < 43; i++) {
@ -261,8 +253,8 @@ public:
n_r(i) *
I_r(i) *
J_r(i) *
Toolbox::pow(pi_, static_cast<Scalar>(I_r(i) - 1)) *
Toolbox::pow(tau_ - 0.5, static_cast<Scalar>(J_r(i) - 1));
Opm::pow(pi_, static_cast<Scalar>(I_r(i) - 1)) *
Opm::pow(tau_ - 0.5, static_cast<Scalar>(J_r(i) - 1));
}
return result;
@ -283,8 +275,6 @@ public:
template <class Evaluation>
static Evaluation ddgamma_ddpi(const Evaluation& temperature, const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& tau_ = tau(temperature); /* reduced temperature */
const Evaluation& pi_ = pi(pressure); /* reduced pressure */
@ -297,8 +287,8 @@ public:
n_r(i) *
I_r(i) *
(I_r(i) - 1) *
Toolbox::pow(pi_, static_cast<Scalar>(I_r(i) - 2)) *
Toolbox::pow(tau_ - 0.5, static_cast<Scalar>(J_r(i)));
Opm::pow(pi_, static_cast<Scalar>(I_r(i) - 2)) *
Opm::pow(tau_ - 0.5, static_cast<Scalar>(J_r(i)));
}
return result;
@ -319,29 +309,27 @@ public:
template <class Evaluation>
static Evaluation ddgamma_ddtau(const Evaluation& temperature, const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& tau_ = tau(temperature); /* reduced temperature */
const Evaluation& pi_ = pi(pressure); /* reduced pressure */
// ideal gas part
Evaluation result = Toolbox::createConstant(0.0);
Evaluation result = 0.0;
for (int i = 0; i < 9; i++) {
result +=
n_g(i) *
J_g(i) *
(J_g(i) - 1) *
Toolbox::pow(tau_, static_cast<Scalar>(J_g(i) - 2));
Opm::pow(tau_, static_cast<Scalar>(J_g(i) - 2));
}
// residual part
for (int i = 0; i < 43; i++) {
result +=
n_r(i) *
Toolbox::pow(pi_, I_r(i)) *
Opm::pow(pi_, I_r(i)) *
J_r(i) *
(J_r(i) - 1.) *
Toolbox::pow(tau_ - 0.5, static_cast<Scalar>(J_r(i) - 2));
Opm::pow(tau_ - 0.5, static_cast<Scalar>(J_r(i) - 2));
}
return result;

10
opm/material/components/iapws/Region4.hpp
View File

@ -62,8 +62,6 @@ public:
template <class Evaluation>
static Evaluation saturationPressure(const Evaluation& temperature)
{
typedef MathToolbox<Evaluation> Toolbox;
static const Scalar n[10] = {
0.11670521452767e4, -0.72421316703206e6, -0.17073846940092e2,
0.12020824702470e5, -0.32325550322333e7, 0.14915108613530e2,
@ -77,7 +75,7 @@ public:
const Evaluation& B = (n[2]*sigma + n[3])*sigma + n[4];
const Evaluation& C = (n[5]*sigma + n[6])*sigma + n[7];
Evaluation tmp = 2*C/(Toolbox::sqrt(B*B - 4*A*C) - B);
Evaluation tmp = 2*C/(Opm::sqrt(B*B - 4*A*C) - B);
tmp *= tmp;
tmp *= tmp;
@ -95,8 +93,6 @@ public:
template <class Evaluation>
static Evaluation vaporTemperature(const Evaluation& pressure)
{
typedef MathToolbox<Evaluation> Toolbox;
static const Scalar n[10] = {
0.11670521452767e4, -0.72421316703206e6, -0.17073846940092e2,
0.12020824702470e5, -0.32325550322333e7, 0.14915108613530e2,
@ -109,9 +105,9 @@ public:
const Evaluation& F = n[0]*beta2 + n[3]*beta + n[6];
const Evaluation& G = n[1]*beta2 + n[4]*beta + n[7];
const Evaluation& D = ( 2.*G)/(-F -Toolbox::sqrt(pow(F,2.) - 4.*E*G));
const Evaluation& D = ( 2.*G)/(-F -Opm::sqrt(pow(F,2.) - 4.*E*G));
const Evaluation& temperature = (n[9] + D - Toolbox::sqrt(pow(n[9]+D , 2.) - 4.* (n[8] + n[9]*D)) ) * 0.5;
const Evaluation& temperature = (n[9] + D - Opm::sqrt(pow(n[9]+D , 2.) - 4.* (n[8] + n[9]*D)) ) * 0.5;
return temperature;
}

24
opm/material/constraintsolvers/CompositionFromFugacities.hpp
View File

@ -90,8 +90,6 @@ public:
unsigned phaseIdx,
const ComponentVector& targetFug)
{
typedef MathToolbox<Evaluation> Toolbox;
// use a much more efficient method in case the phase is an
// ideal mixture
if (FluidSystem::isIdealMixture(phaseIdx)) {
@ -140,7 +138,7 @@ public:
*/
// Solve J*x = b
x = Toolbox::createConstant(0.0);
x = 0.0;
try { J.solve(x, b); }
catch (Dune::FMatrixError e)
{ throw Opm::NumericalProblem(e.what()); }
@ -224,8 +222,6 @@ protected:
unsigned phaseIdx,
const ComponentVector& targetFug)
{
typedef MathToolbox<Evaluation> Toolbox;
// reset jacobian
J = 0;
@ -241,7 +237,7 @@ protected:
fluidState.setFugacityCoefficient(phaseIdx, i, phi);
defect[i] = targetFug[i] - f;
absError = std::max(absError, std::abs(Toolbox::scalarValue(defect[i])));
absError = std::max(absError, std::abs(Opm::scalarValue(defect[i])));
}
// assemble jacobian matrix of the constraints for the composition
@ -298,19 +294,17 @@ protected:
unsigned phaseIdx,
const Dune::FieldVector<Evaluation, numComponents>& targetFug)
{
typedef MathToolbox<Evaluation> Toolbox;
// store original composition and calculate relative error
Dune::FieldVector<Evaluation, numComponents> origComp;
Scalar relError = 0;
Evaluation sumDelta = Toolbox::createConstant(0.0);
Evaluation sumx = Toolbox::createConstant(0.0);
Evaluation sumDelta = 0.0;
Evaluation sumx = 0.0;
for (unsigned i = 0; i < numComponents; ++i) {
origComp[i] = fluidState.moleFraction(phaseIdx, i);
relError = std::max(relError, std::abs(Toolbox::scalarValue(x[i])));
relError = std::max(relError, std::abs(Opm::scalarValue(x[i])));
sumx += Toolbox::abs(fluidState.moleFraction(phaseIdx, i));
sumDelta += Toolbox::abs(x[i]);
sumx += Opm::abs(fluidState.moleFraction(phaseIdx, i));
sumDelta += Opm::abs(x[i]);
}
// chop update to at most 20% change in composition
@ -323,10 +317,10 @@ protected:
Evaluation newx = origComp[i] - x[i];
// only allow negative mole fractions if the target fugacity is negative
if (targetFug[i] > 0)
newx = Toolbox::max(0.0, newx);
newx = Opm::max(0.0, newx);
// only allow positive mole fractions if the target fugacity is positive
else if (targetFug[i] < 0)
newx = Toolbox::min(0.0, newx);
newx = Opm::min(0.0, newx);
// if the target fugacity is zero, the mole fraction must also be zero
else
newx = 0;

4
opm/material/constraintsolvers/ComputeFromReferencePhase.hpp
View File

@ -112,8 +112,6 @@ public:
bool setViscosity,
bool setEnthalpy)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
// compute the density and enthalpy of the
// reference phase
paramCache.updatePhase(fluidState, refPhaseIdx);
@ -152,7 +150,7 @@ public:
ComponentVector fugVec;
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
const auto& fug = fluidState.fugacity(refPhaseIdx, compIdx);
fugVec[compIdx] = FsToolbox::template decay<Evaluation>(fug);
fugVec[compIdx] = Opm::decay<Evaluation>(fug);
}
CompositionFromFugacities::solve(fluidState, paramCache, phaseIdx, fugVec);

18
opm/material/constraintsolvers/ImmiscibleFlash.hpp
View File

@ -259,8 +259,6 @@ protected:
typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& flashParamCache)
{
typedef typename FlashFluidState::Scalar FlashEval;
typedef typename InputFluidState::Scalar Evaluation;
typedef MathToolbox<Evaluation> InputToolbox;
// copy the temperature: even though the model which uses the flash solver might
// be non-isothermal, the flash solver does not consider energy. (it could be
@ -271,7 +269,7 @@ protected:
// copy the saturations: the first N-1 phases are primary variables, the last one
// is one minus the sum of the former.
FlashEval Slast = InputToolbox::createConstant(1.0);
FlashEval Slast = 1.0;
for (unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx) {
FlashEval S = inputFluidState.saturation(phaseIdx);
S.setDerivative(S0PvIdx + phaseIdx, 1.0);
@ -361,13 +359,12 @@ protected:
typedef Opm::MathToolbox<FlashEval> FlashEvalToolbox;
typedef typename FlashEvalToolbox::ValueType InnerEval;
typedef Opm::MathToolbox<InnerEval> InnerEvalToolbox;
#ifndef NDEBUG
// make sure we don't swallow non-finite update vectors
assert(deltaX.dimension == numEq);
for (unsigned i = 0; i < numEq; ++i)
assert(std::isfinite(InnerEvalToolbox::scalarValue(deltaX[i])));
assert(std::isfinite(Opm::scalarValue(deltaX[i])));
#endif
Scalar relError = 0;
@ -376,19 +373,19 @@ protected:
InnerEval delta = deltaX[pvIdx];
relError = std::max(relError,
std::abs(InnerEvalToolbox::scalarValue(delta))
std::abs(Opm::scalarValue(delta))
* quantityWeight_(fluidState, pvIdx));
if (isSaturationIdx_(pvIdx)) {
// dampen to at most 20% change in saturation per
// iteration
delta = InnerEvalToolbox::min(0.25, InnerEvalToolbox::max(-0.25, delta));
delta = Opm::min(0.25, Opm::max(-0.25, delta));
}
else if (isPressureIdx_(pvIdx)) {
// dampen to at most 30% change in pressure per
// iteration
delta = InnerEvalToolbox::min(0.5*fluidState.pressure(0).value(),
InnerEvalToolbox::max(-0.50*fluidState.pressure(0).value(), delta));
delta = Opm::min(0.5*fluidState.pressure(0).value(),
Opm::max(-0.50*fluidState.pressure(0).value(), delta));
}
tmp -= delta;
@ -408,11 +405,10 @@ protected:
typedef typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar> ParamCache;
typedef typename FlashFluidState::Scalar FlashEval;
typedef Opm::MathToolbox<FlashEval> FlashToolbox;
// calculate the saturation of the last phase as a function of
// the other saturations
FlashEval sumSat = FlashToolbox::createConstant(0.0);
FlashEval sumSat = 0.0;
for (unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
sumSat += flashFluidState.saturation(phaseIdx);
flashFluidState.setSaturation(/*phaseIdx=*/numPhases - 1, 1.0 - sumSat);

13
opm/material/constraintsolvers/MiscibleMultiPhaseComposition.hpp
View File

@ -170,7 +170,6 @@ public:
bool setViscosity,
bool setInternalEnergy)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
static_assert(std::is_same<typename FluidState::Scalar, Evaluation>::value,
"The scalar type of the fluid state must be 'Evaluation'");
@ -192,7 +191,7 @@ public:
// coefficients of the components cannot depend on
// composition, i.e. the parameters in the cache are valid
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
Evaluation fugCoeff = FsToolbox::template decay<Evaluation>(
Evaluation fugCoeff = Opm::decay<Evaluation>(
FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx));
fluidState.setFugacityCoefficient(phaseIdx, compIdx, fugCoeff);
}
@ -201,9 +200,9 @@ public:
// create the linear system of equations which defines the
// mole fractions
static const int numEq = numComponents*numPhases;
Dune::FieldMatrix<Evaluation, numEq, numEq> M(Toolbox::createConstant(0.0));
Dune::FieldVector<Evaluation, numEq> x(Toolbox::createConstant(0.0));
Dune::FieldVector<Evaluation, numEq> b(Toolbox::createConstant(0.0));
Dune::FieldMatrix<Evaluation, numEq, numEq> M(0.0);
Dune::FieldVector<Evaluation, numEq> x(0.0);
Dune::FieldVector<Evaluation, numEq> b(0.0);
// assemble the equations expressing the fact that the
// fugacities of each component are equal in all phases
@ -237,11 +236,11 @@ public:
unsigned rowIdx = numComponents*(numPhases - 1) + presentPhases;
presentPhases += 1;
b[rowIdx] = Toolbox::createConstant(1.0);
b[rowIdx] = 1.0;
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
unsigned colIdx = phaseIdx*numComponents + compIdx;
M[rowIdx][colIdx] = Toolbox::createConstant(1.0);
M[rowIdx][colIdx] = 1.0;
}
}

25
opm/material/constraintsolvers/NcpFlash.hpp
View File

@ -315,8 +315,6 @@ protected:
typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& flashParamCache)
{
typedef typename FlashFluidState::Scalar FlashEval;
typedef typename InputFluidState::Scalar Evaluation;
typedef MathToolbox<Evaluation> InputToolbox;
// copy the temperature: even though the model which uses the flash solver might
// be non-isothermal, the flash solver does not consider energy. (it could be
@ -327,7 +325,7 @@ protected:
// copy the saturations: the first N-1 phases are primary variables, the last one
// is one minus the sum of the former.
FlashEval Slast = InputToolbox::createConstant(1.0);
FlashEval Slast = 1.0;
for (unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx) {
FlashEval S = inputFluidState.saturation(phaseIdx);
S.setDerivative(S0PvIdx + phaseIdx, 1.0);
@ -410,7 +408,6 @@ protected:
const FlashComponentVector& globalMolarities)
{
typedef typename FlashFluidState::Scalar FlashEval;
typedef Opm::MathToolbox<FlashEval> FlashToolbox;
unsigned eqIdx = 0;
@ -443,7 +440,7 @@ protected:
// model assumptions (-> non-linear complementarity functions) must be adhered
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
FlashEval oneMinusSumMoleFrac = FlashToolbox::createConstant(1.0);
FlashEval oneMinusSumMoleFrac = 1.0;
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
oneMinusSumMoleFrac -= fluidState.moleFraction(phaseIdx, compIdx);
@ -465,13 +462,12 @@ protected:
// note that it is possible that FlashEval::Scalar is an Evaluation itself
typedef typename FlashFluidState::Scalar FlashEval;
typedef typename FlashEval::ValueType InnerEval;
typedef Opm::MathToolbox<InnerEval> InnerEvalToolbox;
#ifndef NDEBUG
// make sure we don't swallow non-finite update vectors
assert(deltaX.dimension == numEq);
for (unsigned i = 0; i < numEq; ++i)
assert(std::isfinite(InnerEvalToolbox::scalarValue(deltaX[i])));
assert(std::isfinite(Opm::scalarValue(deltaX[i])));
#endif
Scalar relError = 0;
@ -480,22 +476,22 @@ protected:
InnerEval delta = deltaX[pvIdx];
relError = std::max(relError,
std::abs(InnerEvalToolbox::scalarValue(delta))
std::abs(Opm::scalarValue(delta))
* quantityWeight_(fluidState, pvIdx));
if (isSaturationIdx_(pvIdx)) {
// dampen to at most 25% change in saturation per iteration
delta = InnerEvalToolbox::min(0.25, InnerEvalToolbox::max(-0.25, delta));
delta = Opm::min(0.25, Opm::max(-0.25, delta));
}
else if (isMoleFracIdx_(pvIdx)) {
// dampen to at most 20% change in mole fraction per iteration
delta = InnerEvalToolbox::min(0.20, InnerEvalToolbox::max(-0.20, delta));
delta = Opm::min(0.20, Opm::max(-0.20, delta));
}
else if (isPressureIdx_(pvIdx)) {
// dampen to at most 50% change in pressure per iteration
delta = InnerEvalToolbox::min(0.5*fluidState.pressure(0).value(),
InnerEvalToolbox::max(-0.5*fluidState.pressure(0).value(),
delta));
delta = Opm::min(0.5*fluidState.pressure(0).value(),
Opm::max(-0.5*fluidState.pressure(0).value(),
delta));
}
tmp -= delta;
@ -515,11 +511,10 @@ protected:
typedef typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar> ParamCache;
typedef typename FlashFluidState::Scalar FlashEval;
typedef Opm::MathToolbox<FlashEval> FlashToolbox;
// calculate the saturation of the last phase as a function of
// the other saturations
FlashEval sumSat = FlashToolbox::createConstant(0.0);
FlashEval sumSat = 0.0;
for (unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
sumSat += flashFluidState.saturation(phaseIdx);
flashFluidState.setSaturation(/*phaseIdx=*/numPhases - 1, 1.0 - sumSat);

56
opm/material/eos/PengRobinson.hpp
View File

@ -101,7 +101,6 @@ public:
template <class Evaluation, class Params>
static Evaluation computeVaporPressure(const Params& params, const Evaluation& T)
{
typedef MathToolbox<Evaluation> Toolbox;
typedef typename Params::Component Component;
if (T >= Component::criticalTemperature())
return Component::criticalPressure();
@ -130,7 +129,7 @@ public:
const Evaluation& delta = f/df_dp;
pVap = pVap - delta;
if (std::abs(Toolbox::scalarValue(delta/pVap)) < 1e-10)
if (std::abs(Opm::scalarValue(delta/pVap)) < 1e-10)
break;
}
@ -153,7 +152,6 @@ public:
Valgrind::CheckDefined(fs.pressure(phaseIdx));
typedef typename FluidState::Scalar Evaluation;
typedef MathToolbox<Evaluation> Toolbox;
Evaluation Vm = 0;
Valgrind::SetUndefined(Vm);
@ -164,10 +162,10 @@ public:
const Evaluation& a = params.a(phaseIdx); // "attractive factor"
const Evaluation& b = params.b(phaseIdx); // "co-volume"
if (!std::isfinite(Toolbox::scalarValue(a))
|| std::abs(Toolbox::scalarValue(a)) < 1e-30)
if (!std::isfinite(Opm::scalarValue(a))
|| std::abs(Opm::scalarValue(a)) < 1e-30)
return std::numeric_limits<Scalar>::quiet_NaN();
if (!std::isfinite(Toolbox::scalarValue(b)) || b <= 0)
if (!std::isfinite(Opm::scalarValue(b)) || b <= 0)
return std::numeric_limits<Scalar>::quiet_NaN();
const Evaluation& RT= R*T;
@ -229,7 +227,7 @@ public:
}
Valgrind::CheckDefined(Vm);
assert(std::isfinite(Toolbox::scalarValue(Vm)));
assert(std::isfinite(Opm::scalarValue(Vm)));
assert(Vm > 0);
return Vm;
}
@ -247,8 +245,6 @@ public:
template <class Evaluation, class Params>
static Evaluation computeFugacityCoeffient(const Params& params)
{
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& T = params.temperature();
const Evaluation& p = params.pressure();
const Evaluation& Vm = params.molarVolume();
@ -262,8 +258,8 @@ public:
(Vm + params.b()*(1 - std::sqrt(2)));
const Evaluation& expo = - params.a()/(RT * 2 * params.b() * std::sqrt(2));
const Evaluation& fugCoeff =
Toolbox::exp(Z - 1) / (Z - Bstar) *
Toolbox::pow(tmp, expo);
Opm::exp(Z - 1) / (Z - Bstar) *
Opm::pow(tmp, expo);
return fugCoeff;
}
@ -290,7 +286,6 @@ protected:
unsigned phaseIdx,
bool isGasPhase)
{
typedef MathToolbox<Evaluation> Toolbox;
Evaluation Tcrit, pcrit, Vcrit;
findCriticalPoint_(Tcrit,
pcrit,
@ -302,9 +297,9 @@ protected:
//Evaluation Vcrit = criticalMolarVolume_.eval(params.a(phaseIdx), params.b(phaseIdx));
if (isGasPhase)
Vm = Toolbox::max(Vm, Vcrit);
Vm = Opm::max(Vm, Vcrit);
else
Vm = Toolbox::min(Vm, Vcrit);
Vm = Opm::min(Vm, Vcrit);
}
template <class Evaluation>
@ -314,8 +309,6 @@ protected:
const Evaluation& a,
const Evaluation& b)
{
typedef MathToolbox<Evaluation> Toolbox;
Evaluation minVm(0);
Evaluation maxVm(1e30);
@ -356,14 +349,14 @@ protected:
const Scalar eps = - 1e-11;
bool hasExtrema OPM_OPTIM_UNUSED = findExtrema_(minVm, maxVm, minP, maxP, a, b, T + eps);
assert(hasExtrema);
assert(std::isfinite(Toolbox::scalarValue(maxVm)));
assert(std::isfinite(Opm::scalarValue(maxVm)));
Evaluation fStar = maxVm - minVm;
// derivative of the difference between the maximum's
// molar volume and the minimum's molar volume regarding
// temperature
Evaluation fPrime = (fStar - f)/eps;
if (std::abs(Toolbox::scalarValue(fPrime)) < 1e-40) {
if (std::abs(Opm::scalarValue(fPrime)) < 1e-40) {
Tcrit = T;
pcrit = (minP + maxP)/2;
Vcrit = (maxVm + minVm)/2;
@ -372,7 +365,7 @@ protected:
// update value for the current iteration
Evaluation delta = f/fPrime;
assert(std::isfinite(Toolbox::scalarValue(delta)));
assert(std::isfinite(Opm::scalarValue(delta)));
if (delta > 0)
delta = -10;
@ -415,8 +408,6 @@ protected:
const Evaluation& b,
const Evaluation& T)
{
typedef MathToolbox<Evaluation> Toolbox;
Scalar u = 2;
Scalar w = -1;
@ -430,11 +421,11 @@ protected:
const Evaluation& a4 = 2*RT*u*w*b*b*b + 2*u*a*b*b - 2*a*b*b;
const Evaluation& a5 = RT*w*w*b*b*b*b - u*a*b*b*b;
assert(std::isfinite(Toolbox::scalarValue(a1)));
assert(std::isfinite(Toolbox::scalarValue(a2)));
assert(std::isfinite(Toolbox::scalarValue(a3)));
assert(std::isfinite(Toolbox::scalarValue(a4)));
assert(std::isfinite(Toolbox::scalarValue(a5)));
assert(std::isfinite(Opm::scalarValue(a1)));
assert(std::isfinite(Opm::scalarValue(a2)));
assert(std::isfinite(Opm::scalarValue(a3)));
assert(std::isfinite(Opm::scalarValue(a4)));
assert(std::isfinite(Opm::scalarValue(a5)));
// Newton method to find first root
@ -443,11 +434,11 @@ protected:
// above the covolume
Evaluation V = b*1.1;
Evaluation delta = 1.0;
for (unsigned i = 0; std::abs(Toolbox::scalarValue(delta)) > 1e-12; ++i) {
for (unsigned i = 0; std::abs(Opm::scalarValue(delta)) > 1e-12; ++i) {
const Evaluation& f = a5 + V*(a4 + V*(a3 + V*(a2 + V*a1)));
const Evaluation& fPrime = a4 + V*(2*a3 + V*(3*a2 + V*4*a1));
if (std::abs(Toolbox::scalarValue(fPrime)) < 1e-20) {
if (std::abs(Opm::scalarValue(fPrime)) < 1e-20) {
// give up if the derivative is zero
return false;
}
@ -461,7 +452,7 @@ protected:
return false;
}
}
assert(std::isfinite(Toolbox::scalarValue(V)));
assert(std::isfinite(Opm::scalarValue(V)));
// polynomial division
Evaluation b1 = a1;
@ -506,16 +497,15 @@ protected:
template <class Evaluation, class Params>
static Evaluation ambroseWalton_(const Params& /*params*/, const Evaluation& T)
{
typedef MathToolbox<Evaluation> Toolbox;
typedef typename Params::Component Component;
const Evaluation& Tr = T / Component::criticalTemperature();
const Evaluation& tau = 1 - Tr;
const Evaluation& omega = Component::acentricFactor();
const Evaluation& f0 = (tau*(-5.97616 + Toolbox::sqrt(tau)*(1.29874 - tau*0.60394)) - 1.06841*Toolbox::pow(tau, 5))/Tr;
const Evaluation& f1 = (tau*(-5.03365 + Toolbox::sqrt(tau)*(1.11505 - tau*5.41217)) - 7.46628*Toolbox::pow(tau, 5))/Tr;
const Evaluation& f2 = (tau*(-0.64771 + Toolbox::sqrt(tau)*(2.41539 - tau*4.26979)) + 3.25259*Toolbox::pow(tau, 5))/Tr;
const Evaluation& f0 = (tau*(-5.97616 + Opm::sqrt(tau)*(1.29874 - tau*0.60394)) - 1.06841*Opm::pow(tau, 5))/Tr;
const Evaluation& f1 = (tau*(-5.03365 + Opm::sqrt(tau)*(1.11505 - tau*5.41217)) - 7.46628*Opm::pow(tau, 5))/Tr;
const Evaluation& f2 = (tau*(-0.64771 + Opm::sqrt(tau)*(2.41539 - tau*4.26979)) + 3.25259*Opm::pow(tau, 5))/Tr;
return Component::criticalPressure()*std::exp(f0 + omega * (f1 + omega*f2));
}

14
opm/material/eos/PengRobinsonMixture.hpp
View File

@ -93,8 +93,6 @@ public:
unsigned phaseIdx,
unsigned compIdx)
{
typedef MathToolbox<LhsEval> LhsToolbox;
// note that we normalize the component mole fractions, so
// that their sum is 100%. This increases numerical stability
// considerably if the fluid state is not physical.
@ -116,13 +114,13 @@ public:
LhsEval sumMoleFractions = 0.0;
for (unsigned compJIdx = 0; compJIdx < numComponents; ++compJIdx)
sumMoleFractions += fs.moleFraction(phaseIdx, compJIdx);
LhsEval deltai = 2*LhsToolbox::sqrt(params.aPure(phaseIdx, compIdx))/params.a(phaseIdx);
LhsEval deltai = 2*Opm::sqrt(params.aPure(phaseIdx, compIdx))/params.a(phaseIdx);
LhsEval tmp = 0;
for (unsigned compJIdx = 0; compJIdx < numComponents; ++compJIdx) {
tmp +=
fs.moleFraction(phaseIdx, compJIdx)
/ sumMoleFractions
* LhsToolbox::sqrt(params.aPure(phaseIdx, compJIdx))
* Opm::sqrt(params.aPure(phaseIdx, compJIdx))
* (1.0 - StaticParameters::interactionCoefficient(compIdx, compJIdx));
};
deltai *= tmp;
@ -133,8 +131,8 @@ public:
LhsEval expo = Astar/(Bstar*std::sqrt(u*u - 4*w))*(bi_b - deltai);
LhsEval fugCoeff =
LhsToolbox::exp(bi_b*(Z - 1))/LhsToolbox::max(1e-9, Z - Bstar) *
LhsToolbox::pow(base, expo);
Opm::exp(bi_b*(Z - 1))/Opm::max(1e-9, Z - Bstar) *
Opm::pow(base, expo);
////////
// limit the fugacity coefficient to a reasonable range:
@ -142,12 +140,12 @@ public:
// on one side, we want the mole fraction to be at
// least 10^-3 if the fugacity is at the current pressure
//
fugCoeff = LhsToolbox::min(1e10, fugCoeff);
fugCoeff = Opm::min(1e10, fugCoeff);
//
// on the other hand, if the mole fraction of the component is 100%, we want the
// fugacity to be at least 10^-3 Pa
//
fugCoeff = LhsToolbox::max(1e-10, fugCoeff);
fugCoeff = Opm::max(1e-10, fugCoeff);
///////////
return fugCoeff;

18
opm/material/eos/PengRobinsonParamsMixture.hpp
View File

@ -111,13 +111,13 @@ public:
Valgrind::CheckDefined(f_omega);
Scalar tmp = 1 + f_omega*(1 - Toolbox::sqrt(Tr));
Scalar tmp = 1 + f_omega*(1 - Opm::sqrt(Tr));
tmp = tmp*tmp;
Scalar newA = 0.4572355*RTc*RTc/pc * tmp;
Scalar newB = 0.0777961 * RTc / pc;
assert(std::isfinite(Toolbox::scalarValue(newA)));
assert(std::isfinite(Toolbox::scalarValue(newB)));
assert(std::isfinite(Opm::scalarValue(newA)));
assert(std::isfinite(Opm::scalarValue(newB)));
this->pureParams_[i].setA(newA);
this->pureParams_[i].setB(newB);
@ -151,23 +151,23 @@ public:
Scalar newB = 0;
for (unsigned compIIdx = 0; compIIdx < numComponents; ++compIIdx) {
const Scalar moleFracI = fs.moleFraction(phaseIdx, compIIdx);
Scalar xi = Toolbox::max(0.0, Toolbox::min(1.0, moleFracI));
Scalar xi = Opm::max(0.0, Opm::min(1.0, moleFracI));
Valgrind::CheckDefined(xi);
for (unsigned compJIdx = 0; compJIdx < numComponents; ++compJIdx) {
const Scalar moleFracJ = fs.moleFraction(phaseIdx, compJIdx );
Scalar xj = Toolbox::max(0.0, Toolbox::min(1.0, moleFracJ));
Scalar xj = Opm::max(0.0, Opm::min(1.0, moleFracJ));
Valgrind::CheckDefined(xj);
// mixing rule from Reid, page 82
newA += xi * xj * aCache_[compIIdx][compJIdx];
assert(std::isfinite(Toolbox::scalarValue(newA)));
assert(std::isfinite(Opm::scalarValue(newA)));
}
// mixing rule from Reid, page 82
newB += Toolbox::max(0.0, xi) * this->pureParams_[compIIdx].b();
assert(std::isfinite(Toolbox::scalarValue(newB)));
newB += Opm::max(0.0, xi) * this->pureParams_[compIIdx].b();
assert(std::isfinite(Opm::scalarValue(newB)));
}
// assert(newB > 0);
@ -236,7 +236,7 @@ private:
Scalar Psi = FluidSystem::interactionCoefficient(compIIdx, compJIdx);
aCache_[compIIdx][compJIdx] =
Toolbox::sqrt(this->pureParams_[compIIdx].a()
Opm::sqrt(this->pureParams_[compIIdx].a()
* this->pureParams_[compJIdx].a())
* (1 - Psi);
}

70
opm/material/fluidmatrixinteractions/BrooksCorey.hpp
View File

@ -114,7 +114,7 @@ public:
typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
values[Traits::wettingPhaseIdx] = Sw<FluidState, Evaluation>(params, fs);
values[Traits::nonWettingPhaseIdx] = 1 - values[Traits::wettingPhaseIdx];
values[Traits::nonWettingPhaseIdx] = 1.0 - values[Traits::wettingPhaseIdx];
}
/*!
@ -152,12 +152,10 @@ public:
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation pcnw(const Params& params, const FluidState& fs)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const Evaluation& Sw =
FsToolbox::template decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
Opm::decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
assert(0 <= Sw && Sw <= 1);
assert(0.0 <= Sw && Sw <= 1.0);
return twoPhaseSatPcnw(params, Sw);
}
@ -165,21 +163,17 @@ public:
template <class Evaluation>
static Evaluation twoPhaseSatPcnw(const Params& params, const Evaluation& Sw)
{
typedef MathToolbox<Evaluation> Toolbox;
assert(0.0 <= Sw && Sw <= 1.0);
assert(0 <= Sw && Sw <= 1);
return params.entryPressure()*Toolbox::pow(Sw, -1/params.lambda());
return params.entryPressure()*Opm::pow(Sw, -1/params.lambda());
}
template <class Evaluation>
static Evaluation twoPhaseSatPcnwInv(const Params& params, const Evaluation& pcnw)
{
typedef MathToolbox<Evaluation> Toolbox;
assert(pcnw > 0.0);
return Toolbox::pow(params.entryPressure()/pcnw, -params.lambda());
return Opm::pow(params.entryPressure()/pcnw, -params.lambda());
}
/*!
@ -197,22 +191,18 @@ public:
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation Sw(const Params& params, const FluidState& fs)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
Evaluation pC =
FsToolbox::template decay<Evaluation>(fs.pressure(Traits::nonWettingPhaseIdx))
- FsToolbox::template decay<Evaluation>(fs.pressure(Traits::wettingPhaseIdx));
Opm::decay<Evaluation>(fs.pressure(Traits::nonWettingPhaseIdx))
- Opm::decay<Evaluation>(fs.pressure(Traits::wettingPhaseIdx));
return twoPhaseSatSw(params, pC);
}
template <class Evaluation>
static Evaluation twoPhaseSatSw(const Params& params, const Evaluation& pc)
{
typedef MathToolbox<Evaluation> Toolbox;
assert(pc > 0.0); // if we don't assume that, std::pow will screw up!
assert(pc > 0); // if we don't assume that, std::pow will screw up!
return Toolbox::pow(pc/params.entryPressure(), -params.lambda());
return Opm::pow(pc/params.entryPressure(), -params.lambda());
}
/*!
@ -221,11 +211,11 @@ public:
*/
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation Sn(const Params& params, const FluidState& fs)
{ return 1 - Sw<FluidState, Evaluation>(params, fs); }
{ return 1.0 - Sw<FluidState, Evaluation>(params, fs); }
template <class Evaluation>
static Evaluation twoPhaseSatSn(const Params& params, const Evaluation& pc)
{ return 1 - twoPhaseSatSw(params, pc); }
{ return 1.0 - twoPhaseSatSw(params, pc); }
/*!
* \brief The relative permeability for the wetting phase of
* the medium implied by the Brooks-Corey
@ -237,10 +227,8 @@ public:
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krw(const Params& params, const FluidState& fs)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& Sw =
FsToolbox::template decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
Opm::decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
return twoPhaseSatKrw(params, Sw);
}
@ -248,19 +236,15 @@ public:
template <class Evaluation>
static Evaluation twoPhaseSatKrw(const Params& params, const Evaluation& Sw)
{
typedef MathToolbox<Evaluation> Toolbox;
assert(0.0 <= Sw && Sw <= 1.0);
assert(0 <= Sw && Sw <= 1);
return Toolbox::pow(Sw, 2.0/params.lambda() + 3);
return Opm::pow(Sw, 2.0/params.lambda() + 3.0);
}
template <class Evaluation>
static Evaluation twoPhaseSatKrwInv(const Params& params, const Evaluation& krw)
{
typedef MathToolbox<Evaluation> Toolbox;
return Toolbox::pow(krw, 1.0/(2.0/params.lambda() + 3));
return Opm::pow(krw, 1.0/(2.0/params.lambda() + 3.0));
}
/*!
@ -274,10 +258,8 @@ public:
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krn(const Params& params, const FluidState& fs)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const Evaluation& Sw =
1.0 - FsToolbox::template decay<Evaluation>(fs.saturation(Traits::nonWettingPhaseIdx));
1.0 - Opm::decay<Evaluation>(fs.saturation(Traits::nonWettingPhaseIdx));
return twoPhaseSatKrn(params, Sw);
}
@ -285,23 +267,19 @@ public:
template <class Evaluation>
static Evaluation twoPhaseSatKrn(const Params& params, const Evaluation& Sw)
{
typedef MathToolbox<Evaluation> Toolbox;
assert(0.0 <= Sw && Sw <= 1.0);
assert(0 <= Sw && Sw <= 1);
Scalar exponent = 2.0/params.lambda() + 1;
const Evaluation Sn = 1. - Sw;
return Sn*Sn*(1. - Toolbox::pow(Sw, exponent));
Scalar exponent = 2.0/params.lambda() + 1.0;
const Evaluation Sn = 1.0 - Sw;
return Sn*Sn*(1. - Opm::pow(Sw, exponent));
}
template <class Evaluation>
static Evaluation twoPhaseSatKrnInv(const Params& params, const Evaluation& krn)
{
typedef MathToolbox<Evaluation> Toolbox;
// since inverting the formula for krn is hard to do analytically, we use the
// Newton-Raphson method
Evaluation Sw = Toolbox::createConstant(0.5);
Evaluation Sw = 0.5;
Scalar eps = 1e-10;
for (int i = 0; i < 20; ++i) {
Evaluation f = krn - twoPhaseSatKrn(params, Sw);
@ -311,8 +289,8 @@ public:
Evaluation delta = f/fPrime;
Sw -= delta;
if (Sw < 0)
Sw = Toolbox::createConstant(0.0);
if (Toolbox::abs(delta) < 1e-10)
Sw = 0.0;
if (Opm::abs(delta) < 1e-10)
return Sw;
}

37
opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp
View File

@ -223,9 +223,7 @@ public:
static Evaluation pcgn(const Params& params,
const FluidState& fs)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& Sw = 1.0 - FsToolbox::template decay<Evaluation>(fs.saturation(gasPhaseIdx));
const auto& Sw = 1.0 - Opm::decay<Evaluation>(fs.saturation(gasPhaseIdx));
return GasOilMaterialLaw::twoPhaseSatPcnw(params.gasOilParams(), Sw);
}
@ -242,9 +240,7 @@ public:
static Evaluation pcnw(const Params& params,
const FluidState& fs)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& Sw = FsToolbox::template decay<Evaluation>(fs.saturation(waterPhaseIdx));
const auto& Sw = Opm::decay<Evaluation>(fs.saturation(waterPhaseIdx));
return OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw);
}
@ -323,9 +319,7 @@ public:
static Evaluation krg(const Params& params,
const FluidState& fluidState)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const Evaluation& Sw = 1 - FsToolbox::template decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
const Evaluation& Sw = 1 - Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
return GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), Sw);
}
@ -336,9 +330,7 @@ public:
static Evaluation krw(const Params& params,
const FluidState& fluidState)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const Evaluation& Sw = FsToolbox::template decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
const Evaluation& Sw = Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
return OilWaterMaterialLaw::twoPhaseSatKrw(params.oilWaterParams(), Sw);
}
@ -349,15 +341,12 @@ public:
static Evaluation krn(const Params& params,
const FluidState& fluidState)
{
typedef MathToolbox<Evaluation> Toolbox;
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
Scalar Swco = params.Swl();
Evaluation Sw =
Toolbox::max(Evaluation(Swco),
FsToolbox::template decay<Evaluation>(fluidState.saturation(waterPhaseIdx)));
Evaluation Sg = FsToolbox::template decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
Opm::max(Evaluation(Swco),
Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx)));
Evaluation Sg = Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
Evaluation Sw_ow = Sg + Sw;
Evaluation So_go = 1.0 - Sw_ow;
@ -368,9 +357,9 @@ public:
// < epsilon/2 and interpolate between the oridinary and the regularized kro between
// epsilon and epsilon/2
const Scalar epsilon = 1e-5;
if (Toolbox::scalarValue(Sw_ow) - Swco < epsilon) {
if (Opm::scalarValue(Sw_ow) - Swco < epsilon) {
Evaluation kro2 = (kro_ow + kro_go)/2;;
if (Toolbox::scalarValue(Sw_ow) - Swco > epsilon/2) {
if (Opm::scalarValue(Sw_ow) - Swco > epsilon/2) {
Evaluation kro1 = (Sg*kro_go + (Sw - Swco)*kro_ow)/(Sw_ow - Swco);
Evaluation alpha = (epsilon - (Sw_ow - Swco))/(epsilon/2);
return kro2*alpha + kro1*(1 - alpha);
@ -392,11 +381,9 @@ public:
template <class FluidState>
static void updateHysteresis(Params& params, const FluidState& fluidState)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
Scalar Sw = FsToolbox::scalarValue(fluidState.saturation(waterPhaseIdx));
Scalar So = FsToolbox::scalarValue(fluidState.saturation(oilPhaseIdx));
Scalar Sg = FsToolbox::scalarValue(fluidState.saturation(gasPhaseIdx));
Scalar Sw = Opm::scalarValue(fluidState.saturation(waterPhaseIdx));
Scalar So = Opm::scalarValue(fluidState.saturation(oilPhaseIdx));
Scalar Sg = Opm::scalarValue(fluidState.saturation(gasPhaseIdx));
if (params.inconsistentHysteresisUpdate()) {
Sg = std::min(Scalar(1.0), std::max(Scalar(0.0), Sg));

33
opm/material/fluidmatrixinteractions/EclStone1Material.hpp
View File

@ -221,9 +221,7 @@ public:
static Evaluation pcgn(const Params& params,
const FluidState& fs)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& Sw = 1.0 - FsToolbox::template decay<Evaluation>(fs.saturation(gasPhaseIdx));
const auto& Sw = 1.0 - Opm::decay<Evaluation>(fs.saturation(gasPhaseIdx));
return GasOilMaterialLaw::twoPhaseSatPcnw(params.gasOilParams(), Sw);
}
@ -240,9 +238,7 @@ public:
static Evaluation pcnw(const Params& params,
const FluidState& fs)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& Sw = FsToolbox::template decay<Evaluation>(fs.saturation(waterPhaseIdx));
const auto& Sw = Opm::decay<Evaluation>(fs.saturation(waterPhaseIdx));
Valgrind::CheckDefined(Sw);
const auto& result = OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw);
Valgrind::CheckDefined(result);
@ -324,9 +320,7 @@ public:
static Evaluation krg(const Params& params,
const FluidState& fluidState)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const Evaluation& Sw = 1 - FsToolbox::template decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
const Evaluation& Sw = 1 - Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
return GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), Sw);
}
@ -337,9 +331,7 @@ public:
static Evaluation krw(const Params& params,
const FluidState& fluidState)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const Evaluation& Sw = FsToolbox::template decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
const Evaluation& Sw = Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
return OilWaterMaterialLaw::twoPhaseSatKrw(params.oilWaterParams(), Sw);
}
@ -350,9 +342,6 @@ public:
static Evaluation krn(const Params& params,
const FluidState& fluidState)
{
typedef MathToolbox<Evaluation> Toolbox;
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
// the Eclipse docu is inconsistent with naming the variable of connate water: In
// some places the connate water saturation is represented by "Swl", in others
// "Swco" is used.
@ -361,8 +350,8 @@ public:
// oil relperm at connate water saturations (with Sg=0)
Scalar krocw = params.krocw();
const Evaluation& Sw = FsToolbox::template decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
const Evaluation& Sg = FsToolbox::template decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
const Evaluation& Sw = Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
const Evaluation& Sg = Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
Evaluation kro_ow = OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Sw);
Evaluation kro_go = GasOilMaterialLaw::twoPhaseSatKrw(params.gasOilParams(), 1 - Sg - Swco);
@ -381,10 +370,10 @@ public:
if (SSw >= 1.0 || SSg >= 1.0)
beta = 1.0;
else
beta = Toolbox::pow( SSo/((1 - SSw)*(1 - SSg)), params.eta());
beta = Opm::pow( SSo/((1 - SSw)*(1 - SSg)), params.eta());
}
return Toolbox::max(0.0, Toolbox::min(1.0, beta*kro_ow*kro_go/krocw));
return Opm::max(0.0, Opm::min(1.0, beta*kro_ow*kro_go/krocw));
}
/*!
@ -397,10 +386,8 @@ public:
template <class FluidState>
static void updateHysteresis(Params& params, const FluidState& fluidState)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
Scalar Sw = FsToolbox::scalarValue(fluidState.saturation(waterPhaseIdx));
Scalar Sg = FsToolbox::scalarValue(fluidState.saturation(gasPhaseIdx));
Scalar Sw = Opm::scalarValue(fluidState.saturation(waterPhaseIdx));
Scalar Sg = Opm::scalarValue(fluidState.saturation(gasPhaseIdx));
params.oilWaterParams().update(/*pcSw=*/Sw, /*krwSw=*/Sw, /*krnSw=*/Sw);
params.gasOilParams().update(/*pcSw=*/1 - Sg, /*krwSw=*/1 - Sg, /*krnSw=*/1 - Sg);

28
opm/material/fluidmatrixinteractions/EclStone2Material.hpp
View File

@ -222,9 +222,7 @@ public:
static Evaluation pcgn(const Params& params,
const FluidState& fs)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& Sw = 1.0 - FsToolbox::template decay<Evaluation>(fs.saturation(gasPhaseIdx));
const auto& Sw = 1.0 - Opm::decay<Evaluation>(fs.saturation(gasPhaseIdx));
return GasOilMaterialLaw::twoPhaseSatPcnw(params.gasOilParams(), Sw);
}
@ -241,9 +239,7 @@ public:
static Evaluation pcnw(const Params& params,
const FluidState& fs)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& Sw = FsToolbox::template decay<Evaluation>(fs.saturation(waterPhaseIdx));
const auto& Sw = Opm::decay<Evaluation>(fs.saturation(waterPhaseIdx));
Valgrind::CheckDefined(Sw);
const auto& result = OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw);
Valgrind::CheckDefined(result);
@ -325,9 +321,7 @@ public:
static Evaluation krg(const Params& params,
const FluidState& fluidState)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const Evaluation& Sw = 1 - FsToolbox::template decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
const Evaluation& Sw = 1 - Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
return GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), Sw);
}
@ -338,9 +332,7 @@ public:
static Evaluation krw(const Params& params,
const FluidState& fluidState)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const Evaluation& Sw = FsToolbox::template decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
const Evaluation& Sw = Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
return OilWaterMaterialLaw::twoPhaseSatKrw(params.oilWaterParams(), Sw);
}
@ -351,11 +343,9 @@ public:
static Evaluation krn(const Params& params,
const FluidState& fluidState)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
Scalar Swco = params.Swl();
const Evaluation& Sw = FsToolbox::template decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
const Evaluation& Sg = FsToolbox::template decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
const Evaluation& Sw = Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
const Evaluation& Sg = Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
Scalar krocw = OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Swco);
Evaluation krow = OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Sw);
@ -376,10 +366,8 @@ public:
template <class FluidState>
static void updateHysteresis(Params& params, const FluidState& fluidState)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
Scalar Sw = FsToolbox::scalarValue(fluidState.saturation(waterPhaseIdx));
Scalar Sg = FsToolbox::scalarValue(fluidState.saturation(gasPhaseIdx));
Scalar Sw = Opm::scalarValue(fluidState.saturation(waterPhaseIdx));
Scalar Sg = Opm::scalarValue(fluidState.saturation(gasPhaseIdx));
params.oilWaterParams().update(/*pcSw=*/Sw, /*krwSw=*/Sw, /*krnSw=*/Sw);
params.gasOilParams().update(/*pcSw=*/1 - Sg, /*krwSw=*/1 - Sg, /*krnSw=*/1 - Sg);

22
opm/material/fluidmatrixinteractions/EclTwoPhaseMaterial.hpp
View File

@ -128,12 +128,11 @@ public:
const FluidState& fluidState)
{
typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
switch (params.approach()) {
case EclTwoPhaseGasOil: {
const Evaluation& So =
FsToolbox::template decay<Evaluation>(fluidState.saturation(oilPhaseIdx));
Opm::decay<Evaluation>(fluidState.saturation(oilPhaseIdx));
values[oilPhaseIdx] = 0.0;
values[gasPhaseIdx] = GasOilMaterialLaw::twoPhaseSatPcnw(params.gasOilParams(), So);
@ -142,7 +141,7 @@ public:
case EclTwoPhaseOilWater: {
const Evaluation& Sw =
FsToolbox::template decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
values[waterPhaseIdx] = 0.0;
values[oilPhaseIdx] = OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw);
@ -151,7 +150,7 @@ public:
case EclTwoPhaseGasWater: {
const Evaluation& Sw =
FsToolbox::template decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
values[waterPhaseIdx] = 0.0;
values[gasPhaseIdx] =
@ -319,12 +318,11 @@ public:
const FluidState& fluidState)
{
typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
switch (params.approach()) {
case EclTwoPhaseGasOil: {
const Evaluation& So =
FsToolbox::template decay<Evaluation>(fluidState.saturation(oilPhaseIdx));
Opm::decay<Evaluation>(fluidState.saturation(oilPhaseIdx));
values[oilPhaseIdx] = GasOilMaterialLaw::twoPhaseSatKrw(params.gasOilParams(), So);
values[gasPhaseIdx] = GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), So);
@ -333,7 +331,7 @@ public:
case EclTwoPhaseOilWater: {
const Evaluation& Sw =
FsToolbox::template decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
values[waterPhaseIdx] = OilWaterMaterialLaw::twoPhaseSatKrw(params.oilWaterParams(), Sw);
values[oilPhaseIdx] = OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Sw);
@ -342,7 +340,7 @@ public:
case EclTwoPhaseGasWater: {
const Evaluation& Sw =
FsToolbox::template decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
values[waterPhaseIdx] = OilWaterMaterialLaw::twoPhaseSatKrw(params.oilWaterParams(), Sw);
values[gasPhaseIdx] = GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), Sw);
@ -392,25 +390,23 @@ public:
template <class FluidState>
static void updateHysteresis(Params& params, const FluidState& fluidState)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
switch (params.approach()) {
case EclTwoPhaseGasOil: {
Scalar So = FsToolbox::scalarValue(fluidState.saturation(oilPhaseIdx));
Scalar So = Opm::scalarValue(fluidState.saturation(oilPhaseIdx));
params.gasOilParams().update(/*pcSw=*/So, /*krwSw=*/So, /*krnSw=*/So);
break;
}
case EclTwoPhaseOilWater: {
Scalar Sw = FsToolbox::scalarValue(fluidState.saturation(waterPhaseIdx));
Scalar Sw = Opm::scalarValue(fluidState.saturation(waterPhaseIdx));
params.oilWaterParams().update(/*pcSw=*/Sw, /*krwSw=*/Sw, /*krnSw=*/Sw);
break;
}
case EclTwoPhaseGasWater: {
Scalar Sw = FsToolbox::scalarValue(fluidState.saturation(waterPhaseIdx));
Scalar Sw = Opm::scalarValue(fluidState.saturation(waterPhaseIdx));
params.oilWaterParams().update(/*pcSw=*/Sw, /*krwSw=*/Sw, /*krnSw=*/0);
params.gasOilParams().update(/*pcSw=*/1.0, /*krwSw=*/0.0, /*krnSw=*/Sw);

51
opm/material/fluidmatrixinteractions/LinearMaterial.hpp
View File

@ -102,11 +102,10 @@ public:
const FluidState& state)
{
typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
for (unsigned phaseIdx = 0; phaseIdx < Traits::numPhases; ++phaseIdx) {
const Evaluation& S =
FsToolbox::template decay<Evaluation>(state.saturation(phaseIdx));
Opm::decay<Evaluation>(state.saturation(phaseIdx));
Valgrind::CheckDefined(S);
values[phaseIdx] =
@ -135,15 +134,13 @@ public:
const FluidState& state)
{
typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
typedef MathToolbox<Evaluation> Toolbox;
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
for (unsigned phaseIdx = 0; phaseIdx < Traits::numPhases; ++phaseIdx) {
const Evaluation& S =
FsToolbox::template decay<Evaluation>(state.saturation(phaseIdx));
Opm::decay<Evaluation>(state.saturation(phaseIdx));
Valgrind::CheckDefined(S);
values[phaseIdx] = Toolbox::max(Toolbox::min(S,1.0),0.0);
values[phaseIdx] = Opm::max(Opm::min(S,1.0),0.0);
}
}
@ -153,9 +150,8 @@ public:
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation pcnw(const Params& params, const FluidState& fs)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const Evaluation& Sw =
FsToolbox::template decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
Opm::decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
Valgrind::CheckDefined(Sw);
const Evaluation& wPhasePressure =
@ -163,7 +159,7 @@ public:
(1.0 - Sw)*params.pcMinSat(Traits::wettingPhaseIdx);
const Evaluation& Sn =
FsToolbox::template decay<Evaluation>(fs.saturation(Traits::nonWettingPhaseIdx));
Opm::decay<Evaluation>(fs.saturation(Traits::nonWettingPhaseIdx));
Valgrind::CheckDefined(Sn);
const Evaluation& nPhasePressure =
@ -232,22 +228,15 @@ public:
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krw(const Params& /*params*/, const FluidState& fs)
{
typedef MathToolbox<Evaluation> Toolbox;
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const Evaluation& Sw =
FsToolbox::template decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
return Toolbox::max(0.0, Toolbox::min(1.0, Sw));
Opm::decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
return Opm::max(0.0, Opm::min(1.0, Sw));
}
template <class Evaluation = Scalar>
static typename std::enable_if<Traits::numPhases == 2, Evaluation>::type
twoPhaseSatKrw(const Params& /*params*/, const Evaluation& Sw)
{
typedef MathToolbox<Evaluation> Toolbox;
return Toolbox::max(0.0, Toolbox::min(1.0, Sw));
}
{ return Opm::max(0.0, Opm::min(1.0, Sw)); }
/*!
* \brief The relative permability of the liquid non-wetting phase
@ -255,21 +244,16 @@ public:
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krn(const Params& /*params*/, const FluidState& fs)
{
typedef MathToolbox<Evaluation> Toolbox;
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const Evaluation& Sn =
FsToolbox::template decay<Evaluation>(fs.saturation(Traits::nonWettingPhaseIdx));
return Toolbox::max(0.0, Toolbox::min(1.0, Sn));
Opm::decay<Evaluation>(fs.saturation(Traits::nonWettingPhaseIdx));
return Opm::max(0.0, Opm::min(1.0, Sn));
}
template <class Evaluation = Scalar>
static typename std::enable_if<Traits::numPhases == 2, Evaluation>::type
twoPhaseSatKrn(const Params& /*params*/, const Evaluation& Sw)
{
typedef MathToolbox<Evaluation> Toolbox;
return Toolbox::max(0.0, Toolbox::min(1.0, Sw));
return Opm::max(0.0, Opm::min(1.0, Sw));
}
/*!
@ -281,12 +265,9 @@ public:
static typename std::enable_if<Traits::numPhases == 3, Evaluation>::type
krg(const Params& /*params*/, const FluidState& fs)
{
typedef MathToolbox<Evaluation> Toolbox;
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const Evaluation& Sg =
FsToolbox::template decay<Evaluation>(fs.saturation(Traits::gasPhaseIdx));
return Toolbox::max(0.0, Toolbox::min(1.0, Sg));
Opm::decay<Evaluation>(fs.saturation(Traits::gasPhaseIdx));
return Opm::max(0.0, Opm::min(1.0, Sg));
}
/*!
@ -298,10 +279,8 @@ public:
static typename std::enable_if<Traits::numPhases == 3, Evaluation>::type
pcgn(const Params& params, const FluidState& fs)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const Evaluation& Sn =
FsToolbox::template decay<Evaluation>(fs.saturation(Traits::nonWettingPhaseIdx));
Opm::decay<Evaluation>(fs.saturation(Traits::nonWettingPhaseIdx));
Valgrind::CheckDefined(Sn);
const Evaluation& nPhasePressure =
@ -309,7 +288,7 @@ public:
(1.0 - Sn)*params.pcMinSat(Traits::nonWettingPhaseIdx);
const Evaluation& Sg =
FsToolbox::template decay<Evaluation>(fs.saturation(Traits::gasPhaseIdx));
Opm::decay<Evaluation>(fs.saturation(Traits::gasPhaseIdx));
Valgrind::CheckDefined(Sg);
const Evaluation& gPhasePressure =

39
opm/material/fluidmatrixinteractions/NullMaterial.hpp
View File

@ -115,13 +115,11 @@ public:
const FluidState& fluidState)
{
typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
typedef MathToolbox<Evaluation> Toolbox;
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
const Evaluation& S =
FsToolbox::template decay<Evaluation>(fluidState.saturation(phaseIdx));
values[phaseIdx] = Toolbox::max(Toolbox::min(S, 1.0), 0.0);
Opm::decay<Evaluation>(fluidState.saturation(phaseIdx));
values[phaseIdx] = Opm::max(Opm::min(S, 1.0), 0.0);
}
}
@ -182,23 +180,16 @@ public:
static typename std::enable_if<(numPhases > 1), Evaluation>::type
krw(const Params& /*params*/, const FluidState& fluidState)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& Sw =
FsToolbox::template decay<Evaluation>(fluidState.saturation(Traits::wettingPhaseIdx));
Opm::decay<Evaluation>(fluidState.saturation(Traits::wettingPhaseIdx));
return Toolbox::max(0.0, Toolbox::min(1.0, Sw));
return Opm::max(0.0, Opm::min(1.0, Sw));
}
template <class Evaluation>
static typename std::enable_if<numPhases == 2, Evaluation>::type
twoPhaseSatKrw(const Params& /*params*/, const Evaluation& Sw)
{
typedef MathToolbox<Evaluation> Toolbox;
return Toolbox::max(0.0, Toolbox::min(1.0, Sw));
}
{ return Opm::max(0.0, Opm::min(1.0, Sw)); }
/*!
* \brief The relative permability of the liquid non-wetting phase
@ -207,22 +198,17 @@ public:
static typename std::enable_if<(numPhases > 1), Evaluation>::type
krn(const Params& /*params*/, const FluidState& fluidState)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& Sn =
FsToolbox::template decay<Evaluation>(fluidState.saturation(Traits::nonWettingPhaseIdx));
Opm::decay<Evaluation>(fluidState.saturation(Traits::nonWettingPhaseIdx));
return Toolbox::max(0.0, Toolbox::min(1.0, Sn));
return Opm::max(0.0, Opm::min(1.0, Sn));
}
template <class Evaluation>
static typename std::enable_if<numPhases == 2, Evaluation>::type
twoPhaseSatKrn(const Params& /*params*/, const Evaluation& Sw)
{
typedef MathToolbox<Evaluation> Toolbox;
return Toolbox::max(0.0, Toolbox::min(1.0, 1.0 - Sw));
return Opm::max(0.0, Opm::min(1.0, 1.0 - Opm::decay<Evaluation>(Sw)));
}
/*!
@ -234,13 +220,10 @@ public:
static typename std::enable_if< (numPhases > 2), Evaluation>::type
krg(const Params& /*params*/, const FluidState& fluidState)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& Sg =
FsToolbox::template decay<Evaluation>(fluidState.saturation(Traits::gasPhaseIdx));
Opm::decay<Evaluation>(fluidState.saturation(Traits::gasPhaseIdx));
return Toolbox::max(0.0, Toolbox::min(1.0, Sg));
return Opm::max(0.0, Opm::min(1.0, Sg));
}
/*!
@ -251,7 +234,7 @@ public:
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static typename std::enable_if< (Traits::numPhases > 2), Evaluation>::type
pcgn(const Params& /*params*/, const FluidState& /*fluidState*/)
{ return 0; }
{ return 0.0; }
};
} // namespace Opm

20
opm/material/fluidmatrixinteractions/ParkerLenhard.hpp
View File

@ -302,9 +302,7 @@ public:
template <class FluidState>
static void update(Params& params, const FluidState& fs)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
Scalar Sw = FsToolbox::scalarValue(fs.saturation(Traits::wettingPhaseIdx));
Scalar Sw = Opm::scalarValue(fs.saturation(Traits::wettingPhaseIdx));
if (Sw > 1 - 1e-5) {
// if the absolute saturation is almost 1,
@ -378,10 +376,8 @@ public:
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation pcnw(const Params& params, const FluidState& fs)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const Evaluation& Sw =
FsToolbox::template decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
Opm::decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
return twoPhaseSatPcnw(params, Sw);
}
@ -389,10 +385,8 @@ public:
template <class Evaluation>
static Evaluation twoPhaseSatPcnw(const Params& params, const Evaluation& Sw)
{
typedef MathToolbox<Evaluation> Toolbox;
// calculate the current apparent saturation
ScanningCurve* sc = findScanningCurve_(params, Toolbox::scalarValue(Sw));
ScanningCurve* sc = findScanningCurve_(params, Opm::scalarValue(Sw));
// calculate the apparant saturation
const Evaluation& Sw_app = absoluteToApparentSw_(params, Sw);
@ -453,10 +447,8 @@ public:
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krw(const Params& params, const FluidState& fs)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const Evaluation& Sw =
FsToolbox::template decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
Opm::decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
return twoPhaseSatKrw(params, Sw);
}
@ -477,10 +469,8 @@ public:
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krn(const Params& params, const FluidState& fs)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const Evaluation& Sw =
FsToolbox::template decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
Opm::decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
return twoPhaseSatKrn(params, Sw);
}

31
opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterial.hpp
View File

@ -131,9 +131,8 @@ public:
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation pcnw(const Params& params, const FluidState& fs)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& Sw =
FsToolbox::template decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
Opm::decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
return twoPhaseSatPcnw(params, Sw);
}
@ -179,9 +178,8 @@ public:
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krw(const Params& params, const FluidState& fs)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& Sw =
FsToolbox::template decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
Opm::decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
return twoPhaseSatKrw(params, Sw);
}
@ -201,9 +199,8 @@ public:
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krn(const Params& params, const FluidState& fs)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& Sw =
FsToolbox::template decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
Opm::decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
return twoPhaseSatKrn(params, Sw);
}
@ -232,14 +229,12 @@ private:
const ValueVector& yValues,
const Evaluation& x)
{
typedef MathToolbox<Evaluation> Toolbox;
if (x <= xValues.front())
return yValues.front();
if (x >= xValues.back())
return yValues.back();
size_t segIdx = findSegmentIndex_(xValues, Toolbox::scalarValue(x));
size_t segIdx = findSegmentIndex_(xValues, Opm::scalarValue(x));
Scalar x0 = xValues[segIdx];
Scalar x1 = xValues[segIdx + 1];
@ -257,14 +252,12 @@ private:
const ValueVector& yValues,
const Evaluation& x)
{
typedef MathToolbox<Evaluation> Toolbox;
if (x >= xValues.front())
return yValues.front();
if (x <= xValues.back())
return yValues.back();
size_t segIdx = findSegmentIndexDescending_(xValues, Toolbox::scalarValue(x));
size_t segIdx = findSegmentIndexDescending_(xValues, Opm::scalarValue(x));
Scalar x0 = xValues[segIdx];
Scalar x1 = xValues[segIdx + 1];
@ -282,14 +275,12 @@ private:
const ValueVector& yValues,
const Evaluation& x)
{
typedef MathToolbox<Evaluation> Toolbox;
if (x <= xValues.front())
return 0.0;
if (x >= xValues.back())
return 0.0;
if (Toolbox::scalarValue(x) <= xValues.front())
return Toolbox::createConstant(0.0);
if (Toolbox::scalarValue(x) >= xValues.back())
return Toolbox::createConstant(0.0);
size_t segIdx = findSegmentIndex_(xValues, Toolbox::scalarValue(x));
size_t segIdx = findSegmentIndex_(xValues, Opm::scalarValue(x));
Scalar x0 = xValues[segIdx];
Scalar x1 = xValues[segIdx + 1];
@ -297,7 +288,7 @@ private:
Scalar y0 = yValues[segIdx];
Scalar y1 = yValues[segIdx + 1];
return Toolbox::createConstant((y1 - y0)/(x1 - x0));
return (y1 - y0)/(x1 - x0);
}
static size_t findSegmentIndex_(const ValueVector& xValues, Scalar x)

42
opm/material/fluidmatrixinteractions/RegularizedBrooksCorey.hpp
View File

@ -173,9 +173,7 @@ public:
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation pcnw(const Params& params, const FluidState& fs)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& Sw = FsToolbox::template decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
const auto& Sw = Opm::decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
return twoPhaseSatPcnw(params, Sw);
}
@ -209,11 +207,9 @@ public:
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation Sw(const Params& params, const FluidState& fs)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const Evaluation& pC =
FsToolbox::template decay<Evaluation>(fs.pressure(Traits::nonWettingPhaseIdx))
- FsToolbox::template decay<Evaluation>(fs.pressure(Traits::wettingPhaseIdx));
Opm::decay<Evaluation>(fs.pressure(Traits::nonWettingPhaseIdx))
- Opm::decay<Evaluation>(fs.pressure(Traits::wettingPhaseIdx));
return twoPhaseSatSw(params, pC);
}
@ -281,21 +277,17 @@ public:
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krw(const Params& params, const FluidState& fs)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& Sw = FsToolbox::template decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
const auto& Sw = Opm::decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
return twoPhaseSatKrw(params, Sw);
}
template <class Evaluation>
static Evaluation twoPhaseSatKrw(const Params& params, const Evaluation& Sw)
{
typedef MathToolbox<Evaluation> Toolbox;
if (Sw <= 0.0)
return Toolbox::createConstant(0.0);
return 0.0;
else if (Sw >= 1.0)
return Toolbox::createConstant(1.0);
return 1.0;
return BrooksCorey::twoPhaseSatKrw(params, Sw);
}
@ -303,12 +295,10 @@ public:
template <class Evaluation>
static Evaluation twoPhaseSatKrwInv(const Params& params, const Evaluation& krw)
{
typedef MathToolbox<Evaluation> Toolbox;
if (krw <= 0.0)
return Toolbox::createConstant(0.0);
return 0.0;
else if (krw >= 1.0)
return Toolbox::createConstant(1.0);
return 1.0;
return BrooksCorey::twoPhaseSatKrwInv(params, krw);
}
@ -330,22 +320,18 @@ public:
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krn(const Params& params, const FluidState& fs)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const Evaluation& Sw =
1.0 - FsToolbox::template decay<Evaluation>(fs.saturation(Traits::nonWettingPhaseIdx));
1.0 - Opm::decay<Evaluation>(fs.saturation(Traits::nonWettingPhaseIdx));
return twoPhaseSatKrn(params, Sw);
}
template <class Evaluation>
static Evaluation twoPhaseSatKrn(const Params& params, const Evaluation& Sw)
{
typedef MathToolbox<Evaluation> Toolbox;
if (Sw >= 1.0)
return Toolbox::createConstant(0.0);
return 0.0;
else if (Sw <= 0.0)
return Toolbox::createConstant(1.0);
return 1.0;
return BrooksCorey::twoPhaseSatKrn(params, Sw);
}
@ -353,12 +339,10 @@ public:
template <class Evaluation>
static Evaluation twoPhaseSatKrnInv(const Params& params, const Evaluation& krn)
{
typedef MathToolbox<Evaluation> Toolbox;
if (krn <= 0.0)
return Toolbox::createConstant(1.0);
return 1.0;
else if (krn >= 1.0)
return Toolbox::createConstant(0.0);
return 0.0;
return BrooksCorey::twoPhaseSatKrnInv(params, krn);
}

30
opm/material/fluidmatrixinteractions/RegularizedVanGenuchten.hpp
View File

@ -162,9 +162,7 @@ public:
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation pcnw(const Params& params, const FluidState& fs)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& Sw = FsToolbox::template decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
const auto& Sw = Opm::decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
return twoPhaseSatPcnw(params, Sw);
}
@ -223,11 +221,9 @@ public:
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation Sw(const Params& params, const FluidState& fs)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const Evaluation& pC =
FsToolbox::template decay<Evaluation>(fs.pressure(Traits::nonWettingPhaseIdx))
- FsToolbox::template decay<Evaluation>(fs.pressure(Traits::wettingPhaseIdx));
Opm::decay<Evaluation>(fs.pressure(Traits::nonWettingPhaseIdx))
- Opm::decay<Evaluation>(fs.pressure(Traits::wettingPhaseIdx));
return twoPhaseSatSw(params, pC);
}
@ -297,22 +293,18 @@ public:
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krw(const Params& params, const FluidState& fs)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& Sw = FsToolbox::template decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
const auto& Sw = Opm::decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
return twoPhaseSatKrw(params, Sw);
}
template <class Evaluation>
static Evaluation twoPhaseSatKrw(const Params& params, const Evaluation& Sw)
{
typedef MathToolbox<Evaluation> Toolbox;
// regularize
if (Sw <= 0)
return Toolbox::createConstant(0);
return 0.0;
else if (Sw >= 1)
return Toolbox::createConstant(1);
return 1.0;
return VanGenuchten::twoPhaseSatKrw(params, Sw);
}
@ -334,23 +326,19 @@ public:
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krn(const Params& params, const FluidState& fs)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const Evaluation& Sw =
1.0 - FsToolbox::template decay<Evaluation>(fs.saturation(Traits::nonWettingPhaseIdx));
1.0 - Opm::decay<Evaluation>(fs.saturation(Traits::nonWettingPhaseIdx));
return twoPhaseSatKrn(params, Sw);
}
template <class Evaluation>
static Evaluation twoPhaseSatKrn(const Params& params, const Evaluation& Sw)
{
typedef MathToolbox<Evaluation> Toolbox;
// regularize
if (Sw <= 0)
return Toolbox::createConstant(1);
return 1.0;
else if (Sw >= 1)
return Toolbox::createConstant(0);
return 0.0;
return VanGenuchten::twoPhaseSatKrn(params, Sw);
}

12
opm/material/fluidmatrixinteractions/SplineTwoPhaseMaterial.hpp
View File

@ -126,10 +126,8 @@ public:
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation pcnw(const Params& params, const FluidState& fluidState)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const Evaluation& Sw =
FsToolbox::template decay<Evaluation>(fluidState.saturation(Traits::wettingPhaseIdx));
Opm::decay<Evaluation>(fluidState.saturation(Traits::wettingPhaseIdx));
return twoPhaseSatPcnw(params, Sw);
}
@ -197,10 +195,8 @@ public:
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krw(const Params& params, const FluidState& fluidState)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const Evaluation& Sw =
FsToolbox::template decay<Evaluation>(fluidState.saturation(Traits::wettingPhaseIdx));
Opm::decay<Evaluation>(fluidState.saturation(Traits::wettingPhaseIdx));
return twoPhaseSatKrw(params, Sw);
}
@ -238,10 +234,8 @@ public:
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krn(const Params& params, const FluidState& fluidState)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const Evaluation& Sn =
FsToolbox::template decay<Evaluation>(fluidState.saturation(Traits::nonWettingPhaseIdx));
Opm::decay<Evaluation>(fluidState.saturation(Traits::nonWettingPhaseIdx));
return twoPhaseSatKrn(params, 1.0 - Sn);
}

69
opm/material/fluidmatrixinteractions/ThreePhaseParkerVanGenuchten.hpp
View File

@ -122,15 +122,12 @@ public:
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation pcgn(const Params& params, const FluidState& fluidState)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
typedef MathToolbox<Evaluation> Toolbox;
Scalar PC_VG_REG = 0.01;
// sum of liquid saturations
const auto& St =
FsToolbox::template decay<Evaluation>(fluidState.saturation(wettingPhaseIdx))
+ FsToolbox::template decay<Evaluation>(fluidState.saturation(nonWettingPhaseIdx));
Opm::decay<Evaluation>(fluidState.saturation(wettingPhaseIdx))
+ Opm::decay<Evaluation>(fluidState.saturation(nonWettingPhaseIdx));
Evaluation Se = (St - params.Swrx())/(1. - params.Swrx());
@ -142,8 +139,8 @@ public:
if (Se>PC_VG_REG && Se<1-PC_VG_REG)
{
const Evaluation& x = Toolbox::pow(Se,-1/params.vgM()) - 1;
return Toolbox::pow(x, 1.0 - params.vgM())/params.vgAlpha();
const Evaluation& x = Opm::pow(Se,-1/params.vgM()) - 1;
return Opm::pow(x, 1.0 - params.vgM())/params.vgAlpha();
}
// value and derivative at regularization point
@ -153,9 +150,9 @@ public:
else
Se_regu = 1-PC_VG_REG;
const Evaluation& x = std::pow(Se_regu,-1/params.vgM())-1;
const Evaluation& pc = Toolbox::pow(x, 1.0/params.vgN())/params.vgAlpha();
const Evaluation& pc = Opm::pow(x, 1.0/params.vgN())/params.vgAlpha();
const Evaluation& pc_prime =
Toolbox::pow(x, 1/params.vgN()-1)
Opm::pow(x, 1/params.vgN()-1)
* std::pow(Se_regu,-1/params.vgM()-1)
/ (-params.vgM())
/ params.vgAlpha()
@ -178,11 +175,8 @@ public:
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation pcnw(const Params& params, const FluidState& fluidState)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& Sw =
FsToolbox::template decay<Evaluation>(fluidState.saturation(wettingPhaseIdx));
Opm::decay<Evaluation>(fluidState.saturation(wettingPhaseIdx));
Evaluation Se = (Sw-params.Swr())/(1.-params.Snr());
Scalar PC_VG_REG = 0.01;
@ -194,8 +188,8 @@ public:
Se=1.0;
if (Se>PC_VG_REG && Se<1-PC_VG_REG) {
Evaluation x = Toolbox::pow(Se,-1/params.vgM()) - 1.0;
x = Toolbox::pow(x, 1 - params.vgM());
Evaluation x = Opm::pow(Se,-1/params.vgM()) - 1.0;
x = Opm::pow(x, 1 - params.vgM());
return x/params.vgAlpha();
}
@ -207,9 +201,9 @@ public:
Se_regu = 1.0 - PC_VG_REG;
const Evaluation& x = std::pow(Se_regu,-1/params.vgM())-1;
const Evaluation& pc = Toolbox::pow(x, 1/params.vgN())/params.vgAlpha();
const Evaluation& pc = Opm::pow(x, 1/params.vgN())/params.vgAlpha();
const Evaluation& pc_prime =
Toolbox::pow(x,1/params.vgN()-1)
Opm::pow(x,1/params.vgN()-1)
* std::pow(Se_regu, -1.0/params.vgM() - 1)
/ (-params.vgM())
/ params.vgAlpha()
@ -277,11 +271,8 @@ public:
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krw(const Params& params, const FluidState& fluidState)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& Sw =
FsToolbox::template decay<Evaluation>(fluidState.saturation(wettingPhaseIdx));
Opm::decay<Evaluation>(fluidState.saturation(wettingPhaseIdx));
// transformation to effective saturation
const Evaluation& Se = (Sw - params.Swr()) / (1-params.Swr());
@ -289,8 +280,8 @@ public:
if(Se > 1.0) return 1.;
if(Se < 0.0) return 0.;
const Evaluation& r = 1. - Toolbox::pow(1 - Toolbox::pow(Se, 1/params.vgM()), params.vgM());
return Toolbox::sqrt(Se)*r*r;
const Evaluation& r = 1. - Opm::pow(1 - Opm::pow(Se, 1/params.vgM()), params.vgM());
return Opm::sqrt(Se)*r*r;
}
/*!
@ -308,15 +299,12 @@ public:
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krn(const Params& params, const FluidState& fluidState)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& Sn =
FsToolbox::template decay<Evaluation>(fluidState.saturation(nonWettingPhaseIdx));
Opm::decay<Evaluation>(fluidState.saturation(nonWettingPhaseIdx));
const Evaluation& Sw =
FsToolbox::template decay<Evaluation>(fluidState.saturation(wettingPhaseIdx));
Evaluation Swe = Toolbox::min((Sw - params.Swr()) / (1 - params.Swr()), 1.);
Evaluation Ste = Toolbox::min((Sw + Sn - params.Swr()) / (1 - params.Swr()), 1.);
Opm::decay<Evaluation>(fluidState.saturation(wettingPhaseIdx));
Evaluation Swe = Opm::min((Sw - params.Swr()) / (1 - params.Swr()), 1.);
Evaluation Ste = Opm::min((Sw + Sn - params.Swr()) / (1 - params.Swr()), 1.);
// regularization
if(Swe <= 0.0) Swe = 0.;
@ -324,8 +312,8 @@ public:
if(Ste - Swe <= 0.0) return 0.;
Evaluation krn_;
krn_ = Toolbox::pow(1 - Toolbox::pow(Swe, 1/params.vgM()), params.vgM());
krn_ -= Toolbox::pow(1 - Toolbox::pow(Ste, 1/params.vgM()), params.vgM());
krn_ = Opm::pow(1 - Opm::pow(Swe, 1/params.vgM()), params.vgM());
krn_ -= Opm::pow(1 - Opm::pow(Ste, 1/params.vgM()), params.vgM());
krn_ *= krn_;
if (params.krRegardsSnr())
@ -333,11 +321,11 @@ public:
// regard Snr in the permeability of the non-wetting
// phase, see Helmig1997
const Evaluation& resIncluded =
Toolbox::max(Toolbox::min(Sw - params.Snr() / (1-params.Swr()), 1.0), 0.0);
krn_ *= Toolbox::sqrt(resIncluded );
Opm::max(Opm::min(Sw - params.Snr() / (1-params.Swr()), 1.0), 0.0);
krn_ *= Opm::sqrt(resIncluded );
}
else
krn_ *= Toolbox::sqrt(Sn / (1 - params.Swr()));
krn_ *= Opm::sqrt(Sn / (1 - params.Swr()));
return krn_;
}
@ -356,12 +344,9 @@ public:
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krg(const Params& params, const FluidState& fluidState)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& Sg =
FsToolbox::template decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
const Evaluation& Se = Toolbox::min(((1-Sg) - params.Sgr()) / (1 - params.Sgr()), 1.);
Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
const Evaluation& Se = Opm::min(((1-Sg) - params.Sgr()) / (1 - params.Sgr()), 1.);
// regularization
if(Se > 1.0)
@ -377,8 +362,8 @@ public:
}
return scaleFactor
* Toolbox::pow(1 - Se, 1.0/3.)
* Toolbox::pow(1 - Toolbox::pow(Se, 1/params.vgM()), 2*params.vgM());
* Opm::pow(1 - Se, 1.0/3.)
* Opm::pow(1 - Opm::pow(Se, 1/params.vgM()), 2*params.vgM());
}
};
} // namespace Opm

38
opm/material/fluidmatrixinteractions/VanGenuchten.hpp
View File

@ -168,10 +168,8 @@ public:
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation pcnw(const Params& params, const FluidState& fs)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const Evaluation& Sw =
FsToolbox::template decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
Opm::decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
assert(0 <= Sw && Sw <= 1);
@ -195,9 +193,7 @@ public:
template <class Evaluation>
static Evaluation twoPhaseSatPcnw(const Params& params, const Evaluation& Sw)
{
typedef MathToolbox<Evaluation> Toolbox;
return Toolbox::pow(Toolbox::pow(Sw, -1.0/params.vgM()) - 1, 1.0/params.vgN())/params.vgAlpha();
return Opm::pow(Opm::pow(Sw, -1.0/params.vgM()) - 1, 1.0/params.vgN())/params.vgAlpha();
}
/*!
@ -215,22 +211,18 @@ public:
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation Sw(const Params& params, const FluidState& fs)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
Evaluation pC =
FsToolbox::template decay<Evaluation>(fs.pressure(Traits::nonWettingPhaseIdx))
- FsToolbox::template decay<Evaluation>(fs.pressure(Traits::wettingPhaseIdx));
Opm::decay<Evaluation>(fs.pressure(Traits::nonWettingPhaseIdx))
- Opm::decay<Evaluation>(fs.pressure(Traits::wettingPhaseIdx));
return twoPhaseSatSw(params, pC);
}
template <class Evaluation>
static Evaluation twoPhaseSatSw(const Params& params, const Evaluation& pC)
{
typedef MathToolbox<Evaluation> Toolbox;
assert(pC >= 0);
return Toolbox::pow(Toolbox::pow(params.vgAlpha()*pC, params.vgN()) + 1, -params.vgM());
return Opm::pow(Opm::pow(params.vgAlpha()*pC, params.vgN()) + 1, -params.vgM());
}
/*!
@ -258,10 +250,8 @@ public:
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krw(const Params& params, const FluidState& fs)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const Evaluation& Sw =
FsToolbox::template decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
Opm::decay<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
return twoPhaseSatKrw(params, Sw);
}
@ -269,12 +259,10 @@ public:
template <class Evaluation>
static Evaluation twoPhaseSatKrw(const Params& params, const Evaluation& Sw)
{
typedef MathToolbox<Evaluation> Toolbox;
assert(0.0 <= Sw && Sw <= 1.0);
Evaluation r = 1.0 - Toolbox::pow(1.0 - Toolbox::pow(Sw, 1/params.vgM()), params.vgM());
return Toolbox::sqrt(Sw)*r*r;
Evaluation r = 1.0 - Opm::pow(1.0 - Opm::pow(Sw, 1/params.vgM()), params.vgM());
return Opm::sqrt(Sw)*r*r;
}
/*!
@ -289,10 +277,8 @@ public:
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krn(const Params& params, const FluidState& fs)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const Evaluation& Sw =
1.0 - FsToolbox::template decay<Evaluation>(fs.saturation(Traits::nonWettingPhaseIdx));
1.0 - Opm::decay<Evaluation>(fs.saturation(Traits::nonWettingPhaseIdx));
return twoPhaseSatKrn(params, Sw);
}
@ -300,13 +286,11 @@ public:
template <class Evaluation>
static Evaluation twoPhaseSatKrn(const Params& params, Evaluation Sw)
{
typedef MathToolbox<Evaluation> Toolbox;
assert(0 <= Sw && Sw <= 1);
return
Toolbox::pow(1 - Sw, 1.0/3) *
Toolbox::pow(1 - Toolbox::pow(Sw, 1/params.vgM()), 2*params.vgM());
Opm::pow(1 - Sw, 1.0/3) *
Opm::pow(1 - Opm::pow(Sw, 1/params.vgM()), 2*params.vgM());
}
};
} // namespace Opm

11
opm/material/fluidstates/FluidStateCompositionModules.hpp
View File

@ -74,13 +74,11 @@ public:
*/
Scalar massFraction(unsigned phaseIdx, unsigned compIdx) const
{
typedef Opm::MathToolbox<Scalar> Toolbox;
return
Toolbox::abs(sumMoleFractions_[phaseIdx])
Opm::abs(sumMoleFractions_[phaseIdx])
*moleFraction_[phaseIdx][compIdx]
*FluidSystem::molarMass(compIdx)
/ Toolbox::max(1e-40, Toolbox::abs(averageMolarMass_[phaseIdx]));
/ Opm::max(1e-40, Opm::abs(averageMolarMass_[phaseIdx]));
}
/*!
@ -136,15 +134,12 @@ public:
template <class FluidState>
void assign(const FluidState& fs)
{
typedef typename FluidState::Scalar FsScalar;
typedef Opm::MathToolbox<FsScalar> FsToolbox;
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
averageMolarMass_[phaseIdx] = 0;
sumMoleFractions_[phaseIdx] = 0;
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
moleFraction_[phaseIdx][compIdx] =
FsToolbox::template decay<Scalar>(fs.moleFraction(phaseIdx, compIdx));
Opm::decay<Scalar>(fs.moleFraction(phaseIdx, compIdx));
averageMolarMass_[phaseIdx] += moleFraction_[phaseIdx][compIdx]*FluidSystem::molarMass(compIdx);
sumMoleFractions_[phaseIdx] += moleFraction_[phaseIdx][compIdx];

4
opm/material/fluidstates/FluidStateDensityModules.hpp
View File

@ -82,10 +82,8 @@ public:
template <class FluidState>
void assign(const FluidState& fs)
{
typedef typename FluidState::Scalar FsScalar;
typedef Opm::MathToolbox<FsScalar> FsToolbox;
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
density_[phaseIdx] = FsToolbox::template decay<Scalar>(fs.density(phaseIdx));
density_[phaseIdx] = Opm::decay<Scalar>(fs.density(phaseIdx));
}
}

4
opm/material/fluidstates/FluidStateEnthalpyModules.hpp
View File

@ -75,10 +75,8 @@ public:
template <class FluidState>
void assign(const FluidState& fs)
{
typedef typename FluidState::Scalar FsScalar;
typedef Opm::MathToolbox<FsScalar> FsToolbox;
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
enthalpy_[phaseIdx] = FsToolbox::template decay<Scalar>(fs.enthalpy(phaseIdx));
enthalpy_[phaseIdx] = Opm::decay<Scalar>(fs.enthalpy(phaseIdx));
}
}

4
opm/material/fluidstates/FluidStatePressureModules.hpp
View File

@ -70,10 +70,8 @@ public:
template <class FluidState>
void assign(const FluidState& fs)
{
typedef typename FluidState::Scalar FsScalar;
typedef Opm::MathToolbox<FsScalar> FsToolbox;
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
pressure_[phaseIdx] = FsToolbox::template decay<Scalar>(fs.pressure(phaseIdx));
pressure_[phaseIdx] = Opm::decay<Scalar>(fs.pressure(phaseIdx));
}
}

4
opm/material/fluidstates/FluidStateSaturationModules.hpp
View File

@ -76,10 +76,8 @@ public:
template <class FluidState>
void assign(const FluidState& fs)
{
typedef typename FluidState::Scalar FsScalar;
typedef Opm::MathToolbox<FsScalar> FsToolbox;
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
saturation_[phaseIdx] = FsToolbox::template decay<Scalar>(fs.saturation(phaseIdx));
saturation_[phaseIdx] = Opm::decay<Scalar>(fs.saturation(phaseIdx));
}
}

8
opm/material/fluidstates/FluidStateTemperatureModules.hpp
View File

@ -125,15 +125,13 @@ public:
template <class FluidState>
void assign(const FluidState& fs)
{
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
temperature_ = FsToolbox::template decay<Scalar>(fs.temperature(/*phaseIdx=*/0));
temperature_ = Opm::decay<Scalar>(fs.temperature(/*phaseIdx=*/0));
#ifndef NDEBUG
typedef Opm::MathToolbox<Scalar> Toolbox;
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
assert(std::abs(FsToolbox::scalarValue(fs.temperature(phaseIdx))
- Toolbox::scalarValue(temperature_)) < 1e-30);
assert(std::abs(Opm::scalarValue(fs.temperature(phaseIdx))
- Opm::scalarValue(temperature_)) < 1e-30);
}
#endif
}

4
opm/material/fluidstates/FluidStateViscosityModules.hpp
View File

@ -69,10 +69,8 @@ public:
template <class FluidState>
void assign(const FluidState& fs)
{
typedef typename FluidState::Scalar FsScalar;
typedef Opm::MathToolbox<FsScalar> FsToolbox;
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
viscosity_[phaseIdx] = FsToolbox::template decay<Scalar>(fs.viscosity(phaseIdx));
viscosity_[phaseIdx] = Opm::decay<Scalar>(fs.viscosity(phaseIdx));
}
}

75
opm/material/fluidsystems/BlackOilFluidSystem.hpp
View File

@ -53,10 +53,8 @@ template <class FluidSystem, class LhsEval, class FluidState>
LhsEval getRs_(typename std::enable_if<!HasMember_Rs<FluidState>::value, const FluidState&>::type fluidState,
unsigned regionIdx)
{
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& XoG =
FsToolbox::template decay<LhsEval>(fluidState.massFraction(FluidSystem::oilPhaseIdx,
Opm::decay<LhsEval>(fluidState.massFraction(FluidSystem::oilPhaseIdx,
FluidSystem::gasCompIdx));
return FluidSystem::convertXoGToRs(XoG, regionIdx);
}
@ -67,18 +65,15 @@ auto getRs_(typename std::enable_if<HasMember_Rs<FluidState>::value, const Fluid
-> decltype(Opm::MathToolbox<typename FluidState::Scalar>
::template decay<LhsEval>(fluidState.Rs()))
{
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
return FsToolbox::template decay<LhsEval>(fluidState.Rs());
return Opm::decay<LhsEval>(fluidState.Rs());
}
template <class FluidSystem, class LhsEval, class FluidState>
LhsEval getRv_(typename std::enable_if<!HasMember_Rv<FluidState>::value, const FluidState&>::type fluidState,
unsigned regionIdx)
{
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& XgO =
FsToolbox::template decay<LhsEval>(fluidState.massFraction(FluidSystem::gasPhaseIdx,
Opm::decay<LhsEval>(fluidState.massFraction(FluidSystem::gasPhaseIdx,
FluidSystem::oilCompIdx));
return FluidSystem::convertXgOToRv(XgO, regionIdx);
}
@ -89,8 +84,7 @@ auto getRv_(typename std::enable_if<HasMember_Rv<FluidState>::value, const Fluid
-> decltype(Opm::MathToolbox<typename FluidState::Scalar>
::template decay<LhsEval>(fluidState.Rv()))
{
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
return FsToolbox::template decay<LhsEval>(fluidState.Rv());
return Opm::decay<LhsEval>(fluidState.Rv());
}
}
@ -135,9 +129,8 @@ public:
template <class OtherCache>
void assignPersistentData(const OtherCache& other)
{
typedef Opm::MathToolbox<typename OtherCache::Evaluation> OtherToolbox;
regionIdx_ = other.regionIndex();
maxOilSat_ = OtherToolbox::value(other.maxOilSat());
maxOilSat_ = other.maxOilSat();
}
/*!
@ -645,10 +638,8 @@ public:
assert(0 <= phaseIdx && phaseIdx <= numPhases);
assert(0 <= regionIdx && regionIdx <= numRegions());
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
switch (phaseIdx) {
case oilPhaseIdx: {
@ -659,7 +650,7 @@ public:
// interpolate between the saturated and undersaturated quantities to
// avoid a discontinuity
const auto& Rs = Opm::BlackOil::template getRs_<ThisType, LhsEval, FluidState>(fluidState, regionIdx);
const auto& alpha = FsToolbox::template decay<LhsEval>(fluidState.saturation(gasPhaseIdx))/1e-4;
const auto& alpha = Opm::decay<LhsEval>(fluidState.saturation(gasPhaseIdx))/1e-4;
const auto& bSat = oilPvt_->saturatedInverseFormationVolumeFactor(regionIdx, T, p);
const auto& bUndersat = oilPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rs);
return alpha*bSat + (1.0 - alpha)*bUndersat;
@ -683,7 +674,7 @@ public:
// interpolate between the saturated and undersaturated quantities to
// avoid a discontinuity
const auto& Rv = Opm::BlackOil::template getRv_<ThisType, LhsEval, FluidState>(fluidState, regionIdx);
const auto& alpha = FsToolbox::template decay<LhsEval>(fluidState.saturation(oilPhaseIdx))/1e-4;
const auto& alpha = Opm::decay<LhsEval>(fluidState.saturation(oilPhaseIdx))/1e-4;
const auto& bSat = gasPvt_->saturatedInverseFormationVolumeFactor(regionIdx, T, p);
const auto& bUndersat = gasPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rv);
return alpha*bSat + (1.0 - alpha)*bUndersat;
@ -720,10 +711,8 @@ public:
assert(0 <= phaseIdx && phaseIdx <= numPhases);
assert(0 <= regionIdx && regionIdx <= numRegions());
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
switch (phaseIdx) {
case oilPhaseIdx: return oilPvt_->saturatedInverseFormationVolumeFactor(regionIdx, T, p);
@ -744,10 +733,8 @@ public:
assert(0 <= compIdx && compIdx <= numComponents);
assert(0 <= regionIdx && regionIdx <= numRegions());
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
// for the fugacity coefficient of the oil component in the oil phase, we use
// some pseudo-realistic value for the vapor pressure to ease physical
@ -786,8 +773,8 @@ public:
const auto& x_oGSat = convertXoGToxoG(X_oGSat, regionIdx);
const auto& x_oOSat = 1.0 - x_oGSat;
const auto& p_o = FsToolbox::template decay<LhsEval>(fluidState.pressure(oilPhaseIdx));
const auto& p_g = FsToolbox::template decay<LhsEval>(fluidState.pressure(gasPhaseIdx));
const auto& p_o = Opm::decay<LhsEval>(fluidState.pressure(oilPhaseIdx));
const auto& p_g = Opm::decay<LhsEval>(fluidState.pressure(gasPhaseIdx));
return phi_oO*p_o*x_oOSat / (p_g*x_gOSat);
}
@ -823,8 +810,8 @@ public:
const auto& X_oGSat = convertRsToXoG(R_sSat, regionIdx);
const auto& x_oGSat = convertXoGToxoG(X_oGSat, regionIdx);
const auto& p_o = FsToolbox::template decay<LhsEval>(fluidState.pressure(oilPhaseIdx));
const auto& p_g = FsToolbox::template decay<LhsEval>(fluidState.pressure(gasPhaseIdx));
const auto& p_o = Opm::decay<LhsEval>(fluidState.pressure(oilPhaseIdx));
const auto& p_g = Opm::decay<LhsEval>(fluidState.pressure(gasPhaseIdx));
return phi_gG*p_g*x_gGSat / (p_o*x_oGSat);
}
@ -870,10 +857,8 @@ public:
assert(0 <= phaseIdx && phaseIdx <= numPhases);
assert(0 <= regionIdx && regionIdx <= numRegions());
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
switch (phaseIdx) {
case oilPhaseIdx: {
@ -884,7 +869,7 @@ public:
// interpolate between the saturated and undersaturated quantities to
// avoid a discontinuity
const auto& Rs = Opm::BlackOil::template getRs_<ThisType, LhsEval, FluidState>(fluidState, regionIdx);
const auto& alpha = FsToolbox::template decay<LhsEval>(fluidState.saturation(gasPhaseIdx))/1e-4;
const auto& alpha = Opm::decay<LhsEval>(fluidState.saturation(gasPhaseIdx))/1e-4;
const auto& muSat = oilPvt_->saturatedViscosity(regionIdx, T, p);
const auto& muUndersat = oilPvt_->viscosity(regionIdx, T, p, Rs);
return alpha*muSat + (1.0 - alpha)*muUndersat;
@ -909,7 +894,7 @@ public:
// interpolate between the saturated and undersaturated quantities to
// avoid a discontinuity
const auto& Rv = Opm::BlackOil::template getRv_<ThisType, LhsEval, FluidState>(fluidState, regionIdx);
const auto& alpha = FsToolbox::template decay<LhsEval>(fluidState.saturation(oilPhaseIdx))/1e-4;
const auto& alpha = Opm::decay<LhsEval>(fluidState.saturation(oilPhaseIdx))/1e-4;
const auto& muSat = gasPvt_->saturatedViscosity(regionIdx, T, p);
const auto& muUndersat = gasPvt_->viscosity(regionIdx, T, p, Rv);
return alpha*muSat + (1.0 - alpha)*muUndersat;
@ -950,11 +935,9 @@ public:
assert(0 <= phaseIdx && phaseIdx <= numPhases);
assert(0 <= regionIdx && regionIdx <= numRegions());
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& So = FsToolbox::template decay<LhsEval>(fluidState.saturation(oilPhaseIdx));
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& So = Opm::decay<LhsEval>(fluidState.saturation(oilPhaseIdx));
switch (phaseIdx) {
case oilPhaseIdx: return oilPvt_->saturatedGasDissolutionFactor(regionIdx, T, p, So, maxOilSaturation);
@ -980,10 +963,8 @@ public:
assert(0 <= phaseIdx && phaseIdx <= numPhases);
assert(0 <= regionIdx && regionIdx <= numRegions());
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
switch (phaseIdx) {
case oilPhaseIdx: return oilPvt_->saturatedGasDissolutionFactor(regionIdx, T, p);
@ -1032,9 +1013,7 @@ public:
assert(0 <= phaseIdx && phaseIdx <= numPhases);
assert(0 <= regionIdx && regionIdx <= numRegions());
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
switch (phaseIdx) {
case oilPhaseIdx: return oilPvt_->saturationPressure(regionIdx, T, Opm::BlackOil::template getRs_<ThisType, LhsEval, FluidState>(fluidState, regionIdx));

73
opm/material/fluidsystems/BrineCO2FluidSystem.hpp
View File

@ -234,20 +234,17 @@ public:
const ParameterCache<ParamCacheEval>& /*paramCache*/,
unsigned phaseIdx)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
typedef MathToolbox<LhsEval> LhsToolbox;
assert(0 <= phaseIdx && phaseIdx < numPhases);
const LhsEval& temperature = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const LhsEval& pressure = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const LhsEval& temperature = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const LhsEval& pressure = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
if (phaseIdx == liquidPhaseIdx) {
// use normalized composition for to calculate the density
// (the relations don't seem to take non-normalized
// compositions too well...)
LhsEval xlBrine = LhsToolbox::min(1.0, LhsToolbox::max(0.0, FsToolbox::template decay<LhsEval>(fluidState.moleFraction(liquidPhaseIdx, BrineIdx))));
LhsEval xlCO2 = LhsToolbox::min(1.0, LhsToolbox::max(0.0, FsToolbox::template decay<LhsEval>(fluidState.moleFraction(liquidPhaseIdx, CO2Idx))));
LhsEval xlBrine = Opm::min(1.0, Opm::max(0.0, Opm::decay<LhsEval>(fluidState.moleFraction(liquidPhaseIdx, BrineIdx))));
LhsEval xlCO2 = Opm::min(1.0, Opm::max(0.0, Opm::decay<LhsEval>(fluidState.moleFraction(liquidPhaseIdx, CO2Idx))));
LhsEval sumx = xlBrine + xlCO2;
xlBrine /= sumx;
xlCO2 /= sumx;
@ -266,8 +263,8 @@ public:
// use normalized composition for to calculate the density
// (the relations don't seem to take non-normalized
// compositions too well...)
LhsEval xgBrine = LhsToolbox::min(1.0, LhsToolbox::max(0.0, FsToolbox::template decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, BrineIdx))));
LhsEval xgCO2 = LhsToolbox::min(1.0, LhsToolbox::max(0.0, FsToolbox::template decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, CO2Idx))));
LhsEval xgBrine = Opm::min(1.0, Opm::max(0.0, Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, BrineIdx))));
LhsEval xgCO2 = Opm::min(1.0, Opm::max(0.0, Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, CO2Idx))));
LhsEval sumx = xgBrine + xgCO2;
xgBrine /= sumx;
xgCO2 /= sumx;
@ -288,12 +285,10 @@ public:
const ParameterCache<ParamCacheEval>& /*paramCache*/,
unsigned phaseIdx)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
assert(0 <= phaseIdx && phaseIdx < numPhases);
const LhsEval& temperature = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const LhsEval& pressure = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const LhsEval& temperature = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const LhsEval& pressure = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
if (phaseIdx == liquidPhaseIdx) {
// assume pure brine for the liquid phase. TODO: viscosity
@ -318,9 +313,6 @@ public:
unsigned phaseIdx,
unsigned compIdx)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
typedef MathToolbox<LhsEval> LhsToolbox;
assert(0 <= phaseIdx && phaseIdx < numPhases);
assert(0 <= compIdx && compIdx < numComponents);
@ -328,10 +320,10 @@ public:
// use the fugacity coefficients of an ideal gas. the
// actual value of the fugacity is not relevant, as long
// as the relative fluid compositions are observed,
return LhsToolbox::createConstant(1.0);
return 1.0;
const LhsEval& temperature = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const LhsEval& pressure = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const LhsEval& temperature = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const LhsEval& pressure = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
assert(temperature > 0);
assert(pressure > 0);
@ -348,8 +340,8 @@ public:
xgH2O);
// normalize the phase compositions
xlCO2 = LhsToolbox::max(0.0, LhsToolbox::min(1.0, xlCO2));
xgH2O = LhsToolbox::max(0.0, LhsToolbox::min(1.0, xgH2O));
xlCO2 = Opm::max(0.0, Opm::min(1.0, xlCO2));
xgH2O = Opm::max(0.0, Opm::min(1.0, xgH2O));
xlH2O = 1.0 - xlCO2;
xgCO2 = 1.0 - xgH2O;
@ -375,10 +367,8 @@ public:
unsigned phaseIdx,
unsigned /*compIdx*/)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const LhsEval& temperature = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const LhsEval& pressure = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const LhsEval& temperature = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const LhsEval& pressure = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
if (phaseIdx == liquidPhaseIdx)
return BinaryCoeffBrineCO2::liquidDiffCoeff(temperature, pressure);
@ -394,16 +384,13 @@ public:
const ParameterCache<ParamCacheEval>& /*paramCache*/,
unsigned phaseIdx)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
typedef MathToolbox<LhsEval> LhsToolbox;
assert(0 <= phaseIdx && phaseIdx < numPhases);
const LhsEval& temperature = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const LhsEval& pressure = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const LhsEval& temperature = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const LhsEval& pressure = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
if (phaseIdx == liquidPhaseIdx) {
const LhsEval& XlCO2 = FsToolbox::template decay<LhsEval>(fluidState.massFraction(phaseIdx, CO2Idx));
const LhsEval& XlCO2 = Opm::decay<LhsEval>(fluidState.massFraction(phaseIdx, CO2Idx));
const LhsEval& result = liquidEnthalpyBrineCO2_(temperature,
pressure,
Brine_IAPWS::salinity,
@ -412,10 +399,10 @@ public:
return result;
}
else {
const LhsEval& XCO2 = FsToolbox::template decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, CO2Idx));
const LhsEval& XBrine = FsToolbox::template decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, BrineIdx));
const LhsEval& XCO2 = Opm::decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, CO2Idx));
const LhsEval& XBrine = Opm::decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, BrineIdx));
LhsEval result = LhsToolbox::createConstant(0);
LhsEval result = 0;
result += XBrine * Brine::gasEnthalpy(temperature, pressure);
result += XCO2 * CO2::gasEnthalpy(temperature, pressure);
Valgrind::CheckDefined(result);
@ -431,14 +418,12 @@ public:
const ParameterCache<ParamCacheEval>& /*paramCache*/,
unsigned phaseIdx)
{
typedef MathToolbox<LhsEval> LhsToolbox;
// TODO way too simple!
if (phaseIdx == liquidPhaseIdx)
return LhsToolbox::createConstant(0.6); // conductivity of water[W / (m K ) ]
return 0.6; // conductivity of water[W / (m K ) ]
// gas phase
return LhsToolbox::createConstant(0.025); // conductivity of air [W / (m K ) ]
return 0.025; // conductivity of air [W / (m K ) ]
}
/*!
@ -458,12 +443,10 @@ public:
const ParameterCache<ParamCacheEval>& /*paramCache*/,
unsigned phaseIdx)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
assert(0 <= phaseIdx && phaseIdx < numPhases);
const LhsEval& temperature = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const LhsEval& pressure = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const LhsEval& temperature = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const LhsEval& pressure = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
if(phaseIdx == liquidPhaseIdx)
return H2O::liquidHeatCapacity(temperature, pressure);
@ -552,8 +535,6 @@ private:
Scalar S, // salinity
const LhsEval& X_CO2_w)
{
typedef MathToolbox<LhsEval> LhsToolbox;
/* X_CO2_w : mass fraction of CO2 in brine */
/* same function as enthalpy_brine, only extended by CO2 content */
@ -579,7 +560,7 @@ private:
theta = T - 273.15;
// Regularization
Scalar scalarTheta = LhsToolbox::scalarValue(theta);
Scalar scalarTheta = Opm::scalarValue(theta);
Scalar S_lSAT = f[0] + scalarTheta*(f[1] + scalarTheta*(f[2] + scalarTheta*f[3]));
if (S > S_lSAT)
S = S_lSAT;
@ -597,7 +578,7 @@ private:
for (i = 0; i<=3; i++) {
for (j=0; j<=2; j++) {
d_h = d_h + a[i][j] * LhsToolbox::pow(theta, static_cast<Scalar>(i)) * std::pow(m, j);
d_h = d_h + a[i][j] * Opm::pow(theta, static_cast<Scalar>(i)) * std::pow(m, j);
}
}
/* heat of dissolution for halite according to Michaelides 1971 */

70
opm/material/fluidsystems/H2OAirFluidSystem.hpp
View File

@ -259,15 +259,12 @@ public:
const ParameterCache<ParamCacheEval>& /*paramCache*/,
unsigned phaseIdx)
{
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
typedef Opm::MathToolbox<LhsEval> LhsToolbox;
assert(0 <= phaseIdx && phaseIdx < numPhases);
const auto& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
LhsEval p;
if (isCompressible(phaseIdx))
p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
else {
// random value which will hopefully cause things to blow
// up if it is used in a calculation!
@ -278,7 +275,7 @@ public:
LhsEval sumMoleFrac = 0;
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
sumMoleFrac += FsToolbox::template decay<LhsEval>(fluidState.moleFraction(phaseIdx, compIdx));
sumMoleFrac += Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, compIdx));
if (phaseIdx == liquidPhaseIdx)
{
@ -291,8 +288,8 @@ public:
const LhsEval& rholH2O = H2O::liquidDensity(T, p);
const LhsEval& clH2O = rholH2O/H2O::molarMass();
const auto& xlH2O = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(liquidPhaseIdx, H2OIdx));
const auto& xlAir = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(liquidPhaseIdx, AirIdx));
const auto& xlH2O = Opm::decay<LhsEval>(fluidState.moleFraction(liquidPhaseIdx, H2OIdx));
const auto& xlAir = Opm::decay<LhsEval>(fluidState.moleFraction(liquidPhaseIdx, AirIdx));
return clH2O*(H2O::molarMass()*xlH2O + Air::molarMass()*xlAir)/sumMoleFrac;
}
@ -303,16 +300,16 @@ public:
// for the gas phase assume an ideal gas
return
IdealGas::molarDensity(T, p)
* FsToolbox::template decay<LhsEval>(fluidState.averageMolarMass(gasPhaseIdx))
/ LhsToolbox::max(1e-5, sumMoleFrac);
* Opm::decay<LhsEval>(fluidState.averageMolarMass(gasPhaseIdx))
/ Opm::max(1e-5, sumMoleFrac);
LhsEval partialPressureH2O =
FsToolbox::template decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, H2OIdx))
*FsToolbox::template decay<LhsEval>(fluidState.pressure(gasPhaseIdx));
Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, H2OIdx))
*Opm::decay<LhsEval>(fluidState.pressure(gasPhaseIdx));
LhsEval partialPressureAir =
FsToolbox::template decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, AirIdx))
*FsToolbox::template decay<LhsEval>(fluidState.pressure(gasPhaseIdx));
Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, AirIdx))
*Opm::decay<LhsEval>(fluidState.pressure(gasPhaseIdx));
return H2O::gasDensity(T, partialPressureH2O) + Air::gasDensity(T, partialPressureAir);
}
@ -325,13 +322,10 @@ public:
const ParameterCache<ParamCacheEval>& /*paramCache*/,
unsigned phaseIdx)
{
typedef Opm::MathToolbox<LhsEval> LhsToolbox;
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
assert(0 <= phaseIdx && phaseIdx < numPhases);
const auto& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
if (phaseIdx == liquidPhaseIdx)
{
@ -372,14 +366,14 @@ public:
for (unsigned j = 0; j < numComponents; ++j) {
LhsEval phiIJ =
1 +
LhsToolbox::sqrt(mu[i]/mu[j]) * // 1 + (mu[i]/mu[j]^1/2
Opm::sqrt(mu[i]/mu[j]) * // 1 + (mu[i]/mu[j]^1/2
std::pow(M[j]/M[i], 1./4.0); // (M[i]/M[j])^1/4
phiIJ *= phiIJ;
phiIJ /= std::sqrt(8*(1 + M[i]/M[j]));
divisor += FsToolbox::template decay<LhsEval>(fluidState.moleFraction(phaseIdx, j))*phiIJ;
divisor += Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, j))*phiIJ;
}
const auto& xAlphaI = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(phaseIdx, i));
const auto& xAlphaI = Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, i));
muResult += xAlphaI*mu[i]/divisor;
}
return muResult;
@ -395,13 +389,11 @@ public:
unsigned phaseIdx,
unsigned compIdx)
{
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
assert(0 <= phaseIdx && phaseIdx < numPhases);
assert(0 <= compIdx && compIdx < numComponents);
const auto& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
if (phaseIdx == liquidPhaseIdx) {
if (compIdx == H2OIdx)
@ -421,10 +413,8 @@ public:
unsigned phaseIdx,
unsigned /*compIdx*/)
{
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
if (phaseIdx == liquidPhaseIdx)
return BinaryCoeff::H2O_Air::liquidDiffCoeff(T, p);
@ -439,10 +429,8 @@ public:
const ParameterCache<ParamCacheEval>& /*paramCache*/,
unsigned phaseIdx)
{
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
Valgrind::CheckDefined(T);
Valgrind::CheckDefined(p);
@ -457,11 +445,11 @@ public:
LhsEval result = 0.0;
result +=
H2O::gasEnthalpy(T, p) *
FsToolbox::template decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, H2OIdx));
Opm::decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, H2OIdx));
result +=
Air::gasEnthalpy(T, p) *
FsToolbox::template decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, AirIdx));
Opm::decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, AirIdx));
return result;
}
OPM_THROW(std::logic_error, "Invalid phase index " << phaseIdx);
@ -473,14 +461,12 @@ public:
const ParameterCache<ParamCacheEval>& /*paramCache*/,
unsigned phaseIdx)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
assert(0 <= phaseIdx && phaseIdx < numPhases);
const LhsEval& temperature =
FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const LhsEval& pressure =
FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
if (phaseIdx == liquidPhaseIdx)
return H2O::liquidThermalConductivity(temperature, pressure);
@ -489,9 +475,9 @@ public:
if (useComplexRelations){
const LhsEval& xAir =
FsToolbox::template decay<LhsEval>(fluidState.moleFraction(phaseIdx, AirIdx));
Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, AirIdx));
const LhsEval& xH2O =
FsToolbox::template decay<LhsEval>(fluidState.moleFraction(phaseIdx, H2OIdx));
Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, H2OIdx));
LhsEval lambdaAir = xAir*lambdaDryAir;
// Assuming Raoult's, Daltons law and ideal gas

78
opm/material/fluidsystems/H2OAirMesityleneFluidSystem.hpp
View File

@ -203,15 +203,13 @@ public:
const ParameterCache<ParamCacheEval>& /*paramCache*/,
unsigned phaseIdx)
{
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
const LhsEval& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const LhsEval& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
if (phaseIdx == waterPhaseIdx) {
// See: Ochs 2008
const LhsEval& p =
H2O::liquidIsCompressible()
? FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx))
? Opm::decay<LhsEval>(fluidState.pressure(phaseIdx))
: 1e30;
const LhsEval& rholH2O = H2O::liquidDensity(T, p);
@ -220,24 +218,24 @@ public:
// this assumes each dissolved molecule displaces exactly one
// water molecule in the liquid
return
clH2O*(H2O::molarMass()*FsToolbox::template decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, H2OIdx)) +
Air::molarMass()*FsToolbox::template decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, airIdx)) +
NAPL::molarMass()*FsToolbox::template decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, NAPLIdx)));
clH2O*(H2O::molarMass()*Opm::decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, H2OIdx)) +
Air::molarMass()*Opm::decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, airIdx)) +
NAPL::molarMass()*Opm::decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, NAPLIdx)));
}
else if (phaseIdx == naplPhaseIdx) {
// assume pure NAPL for the NAPL phase
const LhsEval& p =
NAPL::liquidIsCompressible()
? FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx))
? Opm::decay<LhsEval>(fluidState.pressure(phaseIdx))
: 1e30;
return NAPL::liquidDensity(T, p);
}
assert (phaseIdx == gasPhaseIdx);
const LhsEval& pg = FsToolbox::template decay<LhsEval>(fluidState.pressure(gasPhaseIdx));
const LhsEval& pH2O = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, H2OIdx))*pg;
const LhsEval& pAir = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, airIdx))*pg;
const LhsEval& pNAPL = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, NAPLIdx))*pg;
const LhsEval& pg = Opm::decay<LhsEval>(fluidState.pressure(gasPhaseIdx));
const LhsEval& pH2O = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, H2OIdx))*pg;
const LhsEval& pAir = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, airIdx))*pg;
const LhsEval& pNAPL = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, NAPLIdx))*pg;
return
H2O::gasDensity(T, pH2O) +
Air::gasDensity(T, pAir) +
@ -250,10 +248,8 @@ public:
const ParameterCache<ParamCacheEval>& /*paramCache*/,
unsigned phaseIdx)
{
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
const LhsEval& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const LhsEval& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const LhsEval& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const LhsEval& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
if (phaseIdx == waterPhaseIdx) {
// assume pure water viscosity
@ -291,9 +287,9 @@ public:
NAPL::molarMass()
};
const LhsEval& xgAir = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, airIdx));
const LhsEval& xgH2O = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, H2OIdx));
const LhsEval& xgNapl = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, NAPLIdx));
const LhsEval& xgAir = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, airIdx));
const LhsEval& xgH2O = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, H2OIdx));
const LhsEval& xgNapl = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, NAPLIdx));
const LhsEval& xgAW = xgAir + xgH2O;
const LhsEval& muAW = (mu[airIdx]*xgAir + mu[H2OIdx]*xgH2O)/xgAW;
const LhsEval& MAW = (xgAir*Air::molarMass() + xgH2O*H2O::molarMass())/xgAW;
@ -319,8 +315,8 @@ public:
#if 0
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
const LhsEval& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const LhsEval& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const LhsEval& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const LhsEval& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
LhsEval diffCont;
if (phaseIdx==gasPhaseIdx) {
@ -328,9 +324,9 @@ public:
const LhsEval& diffWC = Opm::BinaryCoeff::H2O_Mesitylene::gasDiffCoeff(T, p);
const LhsEval& diffAW = Opm::BinaryCoeff::H2O_Air::gasDiffCoeff(T, p);
const LhsEval& xga = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, airIdx));
const LhsEval& xgw = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, H2OIdx));
const LhsEval& xgc = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, NAPLIdx));
const LhsEval& xga = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, airIdx));
const LhsEval& xgw = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, H2OIdx));
const LhsEval& xgc = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, NAPLIdx));
if (compIdx==NAPLIdx) return (1 - xgw)/(xga/diffAW + xgc/diffWC);
else if (compIdx==H2OIdx) return (1 - xgc)/(xgw/diffWC + xga/diffAC);
@ -343,9 +339,9 @@ public:
const LhsEval& diffWCl = 1.e-9; // BinaryCoeff::H2O_Mesitylene::liquidDiffCoeff(temperature, pressure);
const LhsEval& diffAWl = 1.e-9; // BinaryCoeff::H2O_Air::liquidDiffCoeff(temperature, pressure);
const LhsEval& xwa = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, airIdx));
const LhsEval& xww = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, H2OIdx));
const LhsEval& xwc = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, NAPLIdx));
const LhsEval& xwa = Opm::decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, airIdx));
const LhsEval& xww = Opm::decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, H2OIdx));
const LhsEval& xwc = Opm::decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, NAPLIdx));
switch (compIdx) {
case NAPLIdx:
@ -376,13 +372,11 @@ public:
unsigned phaseIdx,
unsigned compIdx)
{
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
assert(0 <= phaseIdx && phaseIdx < numPhases);
assert(0 <= compIdx && compIdx < numComponents);
const LhsEval& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const LhsEval& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const LhsEval& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const LhsEval& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
Valgrind::CheckDefined(T);
Valgrind::CheckDefined(p);
@ -424,10 +418,8 @@ public:
const ParameterCache<ParamCacheEval>& /*paramCache*/,
unsigned phaseIdx)
{
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
const LhsEval& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const LhsEval& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const LhsEval& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const LhsEval& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
if (phaseIdx == waterPhaseIdx) {
return H2O::liquidEnthalpy(T, p);
@ -438,9 +430,9 @@ public:
else if (phaseIdx == gasPhaseIdx) {
// gas phase enthalpy depends strongly on composition
LhsEval result = 0;
result += H2O::gasEnthalpy(T, p) * FsToolbox::template decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, H2OIdx));
result += NAPL::gasEnthalpy(T, p) * FsToolbox::template decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, airIdx));
result += Air::gasEnthalpy(T, p) * FsToolbox::template decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, NAPLIdx));
result += H2O::gasEnthalpy(T, p) * Opm::decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, H2OIdx));
result += NAPL::gasEnthalpy(T, p) * Opm::decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, airIdx));
result += Air::gasEnthalpy(T, p) * Opm::decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, NAPLIdx));
return result;
}
@ -453,19 +445,17 @@ public:
const ParameterCache<ParamCacheEval>& /*paramCache*/,
unsigned phaseIdx)
{
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
assert(0 <= phaseIdx && phaseIdx < numPhases);
if (phaseIdx == waterPhaseIdx){ // water phase
const LhsEval& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const LhsEval& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const LhsEval& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const LhsEval& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
return H2O::liquidThermalConductivity(T, p);
}
else if (phaseIdx == gasPhaseIdx) { // gas phase
const LhsEval& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const LhsEval& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const LhsEval& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const LhsEval& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
return Air::gasThermalConductivity(T, p);
}

72
opm/material/fluidsystems/H2OAirXyleneFluidSystem.hpp
View File

@ -170,37 +170,35 @@ public:
const ParameterCache<ParamCacheEval>& /*paramCache*/,
unsigned phaseIdx)
{
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
if (phaseIdx == waterPhaseIdx) {
const auto& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
// See: Ochs 2008
// \todo: proper citation
const LhsEval& rholH2O = H2O::liquidDensity(T, p);
const LhsEval& clH2O = rholH2O/H2O::molarMass();
const auto& xwH2O = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, H2OIdx));
const auto& xwAir = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, airIdx));
const auto& xwNapl = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, NAPLIdx));
const auto& xwH2O = Opm::decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, H2OIdx));
const auto& xwAir = Opm::decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, airIdx));
const auto& xwNapl = Opm::decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, NAPLIdx));
// this assumes each dissolved molecule displaces exactly one
// water molecule in the liquid
return clH2O*(H2O::molarMass()*xwH2O + Air::molarMass()*xwAir + NAPL::molarMass()*xwNapl);
}
else if (phaseIdx == naplPhaseIdx) {
// assume pure NAPL for the NAPL phase
const auto& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
return NAPL::liquidDensity(T, LhsEval(1e30));
}
assert (phaseIdx == gasPhaseIdx);
const auto& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
const LhsEval& pH2O = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, H2OIdx))*p;
const LhsEval& pAir = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, airIdx))*p;
const LhsEval& pNAPL = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, NAPLIdx))*p;
const LhsEval& pH2O = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, H2OIdx))*p;
const LhsEval& pAir = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, airIdx))*p;
const LhsEval& pNAPL = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, NAPLIdx))*p;
return
H2O::gasDensity(T, pH2O) +
Air::gasDensity(T, pAir) +
@ -213,10 +211,8 @@ public:
const ParameterCache<ParamCacheEval>& /*paramCache*/,
unsigned phaseIdx)
{
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
if (phaseIdx == waterPhaseIdx) {
// assume pure water viscosity
@ -252,9 +248,9 @@ public:
NAPL::molarMass()
};
const auto& xgAir = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, airIdx));
const auto& xgH2O = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, H2OIdx));
const auto& xgNapl = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, NAPLIdx));
const auto& xgAir = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, airIdx));
const auto& xgH2O = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, H2OIdx));
const auto& xgNapl = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, NAPLIdx));
const LhsEval& xgAW = xgAir + xgH2O;
const LhsEval& muAW = (mu[airIdx]*xgAir + mu[H2OIdx]*xgH2O)/ xgAW;
@ -278,19 +274,17 @@ public:
unsigned phaseIdx,
unsigned compIdx)
{
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
if (phaseIdx==gasPhaseIdx) {
const auto& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
const LhsEval& diffAC = Opm::BinaryCoeff::Air_Xylene::gasDiffCoeff(T, p);
const LhsEval& diffWC = Opm::BinaryCoeff::H2O_Xylene::gasDiffCoeff(T, p);
const LhsEval& diffAW = Opm::BinaryCoeff::H2O_Air::gasDiffCoeff(T, p);
const LhsEval& xga = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, airIdx));
const LhsEval& xgw = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, H2OIdx));
const LhsEval& xgc = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, NAPLIdx));
const LhsEval& xga = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, airIdx));
const LhsEval& xgw = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, H2OIdx));
const LhsEval& xgc = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, NAPLIdx));
if (compIdx==NAPLIdx) return (1.- xgw)/(xga/diffAW + xgc/diffWC);
else if (compIdx==H2OIdx) return (1.- xgc)/(xgw/diffWC + xga/diffAC);
@ -302,9 +296,9 @@ public:
Scalar diffWCl = 1.e-9; // BinaryCoeff::H2O_Xylene::liquidDiffCoeff(temperature, pressure);
Scalar diffAWl = 1.e-9; // BinaryCoeff::H2O_Air::liquidDiffCoeff(temperature, pressure);
const LhsEval& xwa = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, airIdx));
const LhsEval& xww = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, H2OIdx));
const LhsEval& xwc = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, NAPLIdx));
const LhsEval& xwa = Opm::decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, airIdx));
const LhsEval& xww = Opm::decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, H2OIdx));
const LhsEval& xwc = Opm::decay<LhsEval>(fluidState.moleFraction(waterPhaseIdx, NAPLIdx));
switch (compIdx) {
case NAPLIdx:
@ -332,13 +326,11 @@ public:
unsigned phaseIdx,
unsigned compIdx)
{
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
assert(0 <= phaseIdx && phaseIdx < numPhases);
assert(0 <= compIdx && compIdx < numComponents);
const auto& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
if (phaseIdx == waterPhaseIdx) {
if (compIdx == H2OIdx)
@ -376,10 +368,8 @@ public:
const ParameterCache<ParamCacheEval>& /*paramCache*/,
unsigned phaseIdx)
{
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
if (phaseIdx == waterPhaseIdx) {
return H2O::liquidEnthalpy(T, p);
@ -393,9 +383,9 @@ public:
const LhsEval& hga = Air::gasEnthalpy(T, p);
LhsEval result = 0;
result += hgw * FsToolbox::template decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, H2OIdx));
result += hga * FsToolbox::template decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, airIdx));
result += hgc * FsToolbox::template decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, NAPLIdx));
result += hgw * Opm::decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, H2OIdx));
result += hga * Opm::decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, airIdx));
result += hgc * Opm::decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, NAPLIdx));
return result;
}

86
opm/material/fluidsystems/H2ON2FluidSystem.hpp
View File

@ -272,15 +272,12 @@ public:
{
assert(0 <= phaseIdx && phaseIdx < numPhases);
typedef Opm::MathToolbox<LhsEval> LhsToolbox;
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
LhsEval sumMoleFrac = 0;
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
sumMoleFrac += FsToolbox::template decay<LhsEval>(fluidState.moleFraction(phaseIdx, compIdx));
sumMoleFrac += Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, compIdx));
// liquid phase
if (phaseIdx == liquidPhaseIdx) {
@ -293,8 +290,8 @@ public:
const auto& rholH2O = H2O::liquidDensity(T, p);
const auto& clH2O = rholH2O/H2O::molarMass();
const auto& xlH2O = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(liquidPhaseIdx, H2OIdx));
const auto& xlN2 = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(liquidPhaseIdx, N2Idx));
const auto& xlH2O = Opm::decay<LhsEval>(fluidState.moleFraction(liquidPhaseIdx, H2OIdx));
const auto& xlN2 = Opm::decay<LhsEval>(fluidState.moleFraction(liquidPhaseIdx, N2Idx));
// this assumes each nitrogen molecule displaces exactly one
// water molecule in the liquid
@ -309,16 +306,16 @@ public:
// for the gas phase assume an ideal gas
return
IdealGas::molarDensity(T, p)
* FsToolbox::template decay<LhsEval>(fluidState.averageMolarMass(gasPhaseIdx))
/ LhsToolbox::max(1e-5, sumMoleFrac);
* Opm::decay<LhsEval>(fluidState.averageMolarMass(gasPhaseIdx))
/ Opm::max(1e-5, sumMoleFrac);
// assume ideal mixture: steam and nitrogen don't "see" each
// other
const auto& xgH2O = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, H2OIdx));
const auto& xgN2 = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, N2Idx));
const auto& xgH2O = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, H2OIdx));
const auto& xgN2 = Opm::decay<LhsEval>(fluidState.moleFraction(gasPhaseIdx, N2Idx));
const auto& rho_gH2O = H2O::gasDensity(T, p*xgH2O);
const auto& rho_gN2 = N2::gasDensity(T, p*xgN2);
return (rho_gH2O + rho_gN2)/LhsToolbox::max(1e-5, sumMoleFrac);
return (rho_gH2O + rho_gN2)/Opm::max(1e-5, sumMoleFrac);
}
//! \copydoc BaseFluidSystem::viscosity
@ -329,11 +326,8 @@ public:
{
assert(0 <= phaseIdx && phaseIdx < numPhases);
typedef Opm::MathToolbox<LhsEval> LhsToolbox;
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
// liquid phase
if (phaseIdx == liquidPhaseIdx)
@ -361,21 +355,21 @@ public:
LhsEval sumx = 0.0;
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
sumx += FsToolbox::template decay<LhsEval>(fluidState.moleFraction(phaseIdx, compIdx));
sumx = LhsToolbox::max(1e-10, sumx);
sumx += Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, compIdx));
sumx = Opm::max(1e-10, sumx);
for (unsigned i = 0; i < numComponents; ++i) {
LhsEval divisor = 0;
for (unsigned j = 0; j < numComponents; ++j) {
LhsEval phiIJ = 1 + LhsToolbox::sqrt(mu[i]/mu[j]) * std::pow(molarMass(j)/molarMass(i), 1/4.0);
LhsEval phiIJ = 1 + Opm::sqrt(mu[i]/mu[j]) * std::pow(molarMass(j)/molarMass(i), 1/4.0);
phiIJ *= phiIJ;
phiIJ /= std::sqrt(8*(1 + molarMass(i)/molarMass(j)));
divisor +=
FsToolbox::template decay<LhsEval>(fluidState.moleFraction(phaseIdx, j))
Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, j))
/sumx*phiIJ;
}
muResult +=
FsToolbox::template decay<LhsEval>(fluidState.moleFraction(phaseIdx, i))
Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, i))
/sumx*mu[i]/divisor;
}
return muResult;
@ -392,10 +386,8 @@ public:
assert(0 <= phaseIdx && phaseIdx < numPhases);
assert(0 <= compIdx && compIdx < numComponents);
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
// liquid phase
if (phaseIdx == liquidPhaseIdx) {
@ -419,10 +411,8 @@ public:
unsigned /*compIdx*/)
{
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
// liquid phase
if (phaseIdx == liquidPhaseIdx)
@ -439,10 +429,8 @@ public:
const ParameterCache<ParamCacheEval>& /*paramCache*/,
unsigned phaseIdx)
{
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
Valgrind::CheckDefined(T);
Valgrind::CheckDefined(p);
@ -459,8 +447,8 @@ public:
// "see" the molecules of the other component, which means
// that the total specific enthalpy is the sum of the
// "partial specific enthalpies" of the components.
const auto& XgH2O = FsToolbox::template decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, H2OIdx));
const auto& XgN2 = FsToolbox::template decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, N2Idx));
const auto& XgH2O = Opm::decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, H2OIdx));
const auto& XgN2 = Opm::decay<LhsEval>(fluidState.massFraction(gasPhaseIdx, N2Idx));
LhsEval hH2O = XgH2O*H2O::gasEnthalpy(T, p);
LhsEval hN2 = XgN2*N2::gasEnthalpy(T, p);
@ -475,10 +463,8 @@ public:
{
assert(0 <= phaseIdx && phaseIdx < numPhases);
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
if (phaseIdx == liquidPhaseIdx) // liquid phase
return H2O::liquidThermalConductivity(T, p);
@ -487,8 +473,8 @@ public:
if (useComplexRelations){
// return the sum of the partial conductivity of Nitrogen and Steam
const auto& xH2O = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(phaseIdx, H2OIdx));
const auto& xN2 = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(phaseIdx, N2Idx));
const auto& xH2O = Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, H2OIdx));
const auto& xN2 = Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, N2Idx));
// Assuming Raoult's, Daltons law and ideal gas in order to obtain the
// partial pressures in the gas phase
@ -508,14 +494,12 @@ public:
const ParameterCache<ParamCacheEval>& /*paramCache*/,
unsigned phaseIdx)
{
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& xAlphaH2O = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(phaseIdx, H2OIdx));
const auto& xAlphaN2 = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(phaseIdx, N2Idx));
const auto& XAlphaH2O = FsToolbox::template decay<LhsEval>(fluidState.massFraction(phaseIdx, H2OIdx));
const auto& XAlphaN2 = FsToolbox::template decay<LhsEval>(fluidState.massFraction(phaseIdx, N2Idx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& xAlphaH2O = Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, H2OIdx));
const auto& xAlphaN2 = Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, N2Idx));
const auto& XAlphaH2O = Opm::decay<LhsEval>(fluidState.massFraction(phaseIdx, H2OIdx));
const auto& XAlphaN2 = Opm::decay<LhsEval>(fluidState.massFraction(phaseIdx, N2Idx));
if (phaseIdx == liquidPhaseIdx)
return H2O::liquidHeatCapacity(T, p);

48
opm/material/fluidsystems/H2ON2LiquidPhaseFluidSystem.hpp
View File

@ -255,16 +255,14 @@ public:
const ParameterCache<ParamCacheEval>& /*paramCache*/,
unsigned phaseIdx)
{
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
assert(0 <= phaseIdx && phaseIdx < numPhases);
const auto& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
LhsEval sumMoleFrac = 0;
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
sumMoleFrac += FsToolbox::template decay<LhsEval>(fluidState.moleFraction(phaseIdx, compIdx));
sumMoleFrac += Opm::decay<LhsEval>(fluidState.moleFraction(phaseIdx, compIdx));
assert(phaseIdx == liquidPhaseIdx);
@ -277,8 +275,8 @@ public:
const LhsEval& rholH2O = H2O::liquidDensity(T, p);
const LhsEval& clH2O = rholH2O/H2O::molarMass();
const auto& xlH2O = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(liquidPhaseIdx, H2OIdx));
const auto& xlN2 = FsToolbox::template decay<LhsEval>(fluidState.moleFraction(liquidPhaseIdx, N2Idx));
const auto& xlH2O = Opm::decay<LhsEval>(fluidState.moleFraction(liquidPhaseIdx, H2OIdx));
const auto& xlN2 = Opm::decay<LhsEval>(fluidState.moleFraction(liquidPhaseIdx, N2Idx));
// this assumes each nitrogen molecule displaces exactly one
// water molecule in the liquid
@ -292,12 +290,10 @@ public:
const ParameterCache<ParamCacheEval>& /*paramCache*/,
unsigned phaseIdx)
{
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
assert(phaseIdx == liquidPhaseIdx);
const auto& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
// assume pure water for the liquid phase
return H2O::liquidViscosity(T, p);
@ -310,13 +306,11 @@ public:
unsigned phaseIdx,
unsigned compIdx)
{
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
assert(phaseIdx == liquidPhaseIdx);
assert(0 <= compIdx && compIdx < numComponents);
const auto& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
if (compIdx == H2OIdx)
return H2O::vaporPressure(T)/p;
@ -331,12 +325,10 @@ public:
unsigned /*compIdx*/)
{
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
assert(phaseIdx == liquidPhaseIdx);
const auto& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
return BinaryCoeff::H2O_N2::liquidDiffCoeff(T, p);
}
@ -347,12 +339,10 @@ public:
const ParameterCache<ParamCacheEval>& /*paramCache*/,
unsigned phaseIdx)
{
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
assert (phaseIdx == liquidPhaseIdx);
const auto& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
Valgrind::CheckDefined(T);
Valgrind::CheckDefined(p);
@ -366,13 +356,11 @@ public:
const ParameterCache<ParamCacheEval>& /*paramCache*/,
const unsigned phaseIdx)
{
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
assert(phaseIdx == liquidPhaseIdx);
if(useComplexRelations){
const auto& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
return H2O::liquidThermalConductivity(T, p);
}
else
@ -385,12 +373,10 @@ public:
const ParameterCache<ParamCacheEval>& /*paramCache*/,
unsigned phaseIdx)
{
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
assert (phaseIdx == liquidPhaseIdx);
const auto& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
return H2O::liquidHeatCapacity(T, p);
}

30
opm/material/fluidsystems/SinglePhaseFluidSystem.hpp
View File

@ -190,12 +190,10 @@ public:
const ParameterCache<ParamCacheEval>& /*paramCache*/,
unsigned phaseIdx)
{
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
assert(0 <= phaseIdx && phaseIdx < numPhases);
const auto& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
return Fluid::density(T, p);
}
@ -205,12 +203,10 @@ public:
const ParameterCache<ParamCacheEval>& /*paramCache*/,
unsigned phaseIdx)
{
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
assert(0 <= phaseIdx && phaseIdx < numPhases);
const auto& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
return Fluid::viscosity(T, p);
}
@ -239,12 +235,10 @@ public:
const ParameterCache<ParamCacheEval>& /*paramCache*/,
unsigned phaseIdx)
{
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
assert(0 <= phaseIdx && phaseIdx < numPhases);
const auto& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
return Fluid::enthalpy(T, p);
}
@ -254,12 +248,10 @@ public:
const ParameterCache<ParamCacheEval>& /*paramCache*/,
unsigned phaseIdx)
{
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
assert(0 <= phaseIdx && phaseIdx < numPhases);
const auto& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
return Fluid::thermalConductivity(T, p);
}
@ -269,12 +261,10 @@ public:
const ParameterCache<ParamCacheEval>& /*paramCache*/,
unsigned phaseIdx)
{
typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
assert(0 <= phaseIdx && phaseIdx < numPhases);
const auto& T = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
return Fluid::heatCapacity(T, p);
}
};

36
opm/material/fluidsystems/TwoPhaseImmiscibleFluidSystem.hpp
View File

@ -226,12 +226,10 @@ public:
const ParameterCache<ParamCacheEval>& /*paramCache*/,
unsigned phaseIdx)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
assert(0 <= phaseIdx && phaseIdx < numPhases);
const auto& temperature = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& pressure = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& temperature = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& pressure = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
if (phaseIdx == wettingPhaseIdx)
return WettingPhase::density(temperature, pressure);
return NonwettingPhase::density(temperature, pressure);
@ -243,12 +241,10 @@ public:
const ParameterCache<ParamCacheEval>& /*paramCache*/,
unsigned phaseIdx)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
assert(0 <= phaseIdx && phaseIdx < numPhases);
const auto& temperature = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& pressure = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& temperature = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& pressure = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
if (phaseIdx == wettingPhaseIdx)
return WettingPhase::viscosity(temperature, pressure);
return NonwettingPhase::viscosity(temperature, pressure);
@ -261,8 +257,6 @@ public:
unsigned phaseIdx,
unsigned compIdx)
{
typedef MathToolbox<LhsEval> LhsToolbox;
assert(0 <= phaseIdx && phaseIdx < numPhases);
assert(0 <= compIdx && compIdx < numComponents);
@ -271,8 +265,8 @@ public:
// the component in the fluid. Probably that's not worth
// the effort, since the fugacity coefficient of the other
// component is infinite anyway...
return LhsToolbox::createConstant(1.0);
return LhsToolbox::createConstant(std::numeric_limits<Scalar>::infinity());
return 1.0;
return std::numeric_limits<Scalar>::infinity();
}
//! \copydoc BaseFluidSystem::enthalpy
@ -281,12 +275,10 @@ public:
const ParameterCache<ParamCacheEval>& /*paramCache*/,
unsigned phaseIdx)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
assert(0 <= phaseIdx && phaseIdx < numPhases);
const auto& temperature = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& pressure = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& temperature = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& pressure = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
if (phaseIdx == wettingPhaseIdx)
return WettingPhase::enthalpy(temperature, pressure);
return NonwettingPhase::enthalpy(temperature, pressure);
@ -298,12 +290,10 @@ public:
const ParameterCache<ParamCacheEval>& /*paramCache*/,
unsigned phaseIdx)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
assert(0 <= phaseIdx && phaseIdx < numPhases);
const auto& temperature = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& pressure = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& temperature = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& pressure = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
if (phaseIdx == wettingPhaseIdx)
return WettingPhase::thermalConductivity(temperature, pressure);
return NonwettingPhase::thermalConductivity(temperature, pressure);
@ -315,12 +305,10 @@ public:
const ParameterCache<ParamCacheEval>& /*paramCache*/,
unsigned phaseIdx)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
assert(0 <= phaseIdx && phaseIdx < numPhases);
const auto& temperature = FsToolbox::template decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& pressure = FsToolbox::template decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& temperature = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
const auto& pressure = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
if (phaseIdx == wettingPhaseIdx)
return WettingPhase::heatCapacity(temperature, pressure);
return NonwettingPhase::heatCapacity(temperature, pressure);

9
opm/material/fluidsystems/blackoilpvt/LiveOilPvt.hpp
View File

@ -487,7 +487,6 @@ public:
const Evaluation& oilSaturation,
Scalar maxOilSaturation) const
{
typedef typename Opm::MathToolbox<Evaluation> Toolbox;
Evaluation tmp =
saturatedGasDissolutionFactorTable_[regionIdx].eval(pressure, /*extrapolate=*/true);
@ -495,8 +494,8 @@ public:
// keyword)
if (vapPar2_ > 0.0 && maxOilSaturation > 0.01 && oilSaturation < maxOilSaturation) {
static const Scalar sqrtEps = std::sqrt(std::numeric_limits<Scalar>::epsilon());
const Evaluation& So = Toolbox::max(oilSaturation, sqrtEps);
tmp *= Toolbox::pow(So/maxOilSaturation, vapPar2_);
const Evaluation& So = Opm::max(oilSaturation, sqrtEps);
tmp *= Opm::pow(So/maxOilSaturation, vapPar2_);
}
return tmp;
@ -531,7 +530,7 @@ public:
// If the derivative is "zero" Newton will not converge,
// so simply return our initial guess.
if (std::abs(Toolbox::scalarValue(fPrime)) < 1.0e-30) {
if (std::abs(Opm::scalarValue(fPrime)) < 1.0e-30) {
return pSat;
}
@ -549,7 +548,7 @@ public:
pSat = 0.0;
}
if (std::abs(Toolbox::scalarValue(delta)) < std::abs(Toolbox::scalarValue(pSat))*eps)
if (std::abs(Opm::scalarValue(delta)) < std::abs(Opm::scalarValue(pSat))*eps)
return pSat;
}

9
opm/material/fluidsystems/blackoilpvt/WetGasPvt.hpp
View File

@ -511,7 +511,6 @@ public:
const Evaluation& oilSaturation,
Scalar maxOilSaturation) const
{
typedef typename Opm::MathToolbox<Evaluation> Toolbox;
Evaluation tmp =
saturatedOilVaporizationFactorTable_[regionIdx].eval(pressure, /*extrapolate=*/true);
@ -519,8 +518,8 @@ public:
// keyword)
if (vapPar1_ > 0.0 && maxOilSaturation > 0.01 && oilSaturation < maxOilSaturation) {
static const Scalar sqrtEps = std::sqrt(std::numeric_limits<Scalar>::epsilon());
const Evaluation& So = Toolbox::max(oilSaturation, sqrtEps);
tmp *= Toolbox::pow(So/maxOilSaturation, vapPar1_);
const Evaluation& So = Opm::max(oilSaturation, sqrtEps);
tmp *= Opm::pow(So/maxOilSaturation, vapPar1_);
}
return tmp;
@ -559,7 +558,7 @@ public:
// If the derivative is "zero" Newton will not converge,
// so simply return our initial guess.
if (std::abs(Toolbox::scalarValue(fPrime)) < 1.0e-30) {
if (std::abs(Opm::scalarValue(fPrime)) < 1.0e-30) {
return pSat;
}
@ -577,7 +576,7 @@ public:
pSat = 0.0;
}
if (std::abs(Toolbox::scalarValue(delta)) < std::abs(Toolbox::scalarValue(pSat))*eps)
if (std::abs(Opm::scalarValue(delta)) < std::abs(Opm::scalarValue(pSat))*eps)
return pSat;
}

11
opm/material/heatconduction/Somerton.hpp
View File

@ -90,8 +90,6 @@ public:
static Evaluation heatConductivity(const Params& params,
const FluidState& fluidState)
{
typedef Opm::MathToolbox<Evaluation> Toolbox;
Valgrind::CheckDefined(params.vacuumLambda());
Evaluation lambda = 0;
@ -99,9 +97,9 @@ public:
Valgrind::CheckDefined(params.fullySaturatedLambda(phaseIdx));
if (FluidSystem::isLiquid(phaseIdx)) {
const auto& sat = Toolbox::template decay<Evaluation>(fluidState.saturation(phaseIdx));
const auto& sat = Opm::decay<Evaluation>(fluidState.saturation(phaseIdx));
lambda +=
regularizedSqrt_(Toolbox::max(0.0, Toolbox::min(1.0, sat)))
regularizedSqrt_(Opm::max(0.0, Opm::min(1.0, sat)))
* (params.fullySaturatedLambda(phaseIdx) - params.vacuumLambda());
}
else { // gas phase
@ -118,19 +116,18 @@ protected:
template <class Evaluation>
static Evaluation regularizedSqrt_(const Evaluation& x)
{
typedef Opm::MathToolbox<Evaluation> Toolbox;
typedef Opm::Spline<Scalar> Spline;
static const Scalar xMin = 1e-2;
static const Scalar sqrtXMin = std::sqrt(xMin);
static const Scalar fPrimeXMin = 1.0/(2*std::sqrt(xMin));
static const Scalar fPrime0 = 2*fPrimeXMin;
typedef Opm::Spline<Scalar> Spline;
static const Spline sqrtRegSpline(0, xMin, // x0, x1
0, sqrtXMin, // y0, y1
fPrime0, fPrimeXMin); // m0, m1
if (x > xMin)
return Toolbox::sqrt(x);
return Opm::sqrt(x);
else if (x <= 0)
return fPrime0 * x;
else

12
tests/test_eclmateriallawmanager.cpp
View File

@ -385,8 +385,8 @@ inline void testAll()
fs.setSaturation(oilPhaseIdx, So);
fs.setSaturation(gasPhaseIdx, Sg);
Scalar pcFam1[numPhases];
Scalar pcFam2[numPhases];
Scalar pcFam1[numPhases] = { 0.0, 0.0 };
Scalar pcFam2[numPhases] = { 0.0, 0.0 };
MaterialLaw::capillaryPressures(pcFam1,
materialLawManager.materialLawParams(elemIdx),
fs);
@ -394,8 +394,8 @@ inline void testAll()
fam2MaterialLawManager.materialLawParams(elemIdx),
fs);
Scalar krFam1[numPhases];
Scalar krFam2[numPhases];
Scalar krFam1[numPhases] = { 0.0, 0.0 };
Scalar krFam2[numPhases] = { 0.0, 0.0 };
MaterialLaw::relativePermeabilities(krFam1,
materialLawManager.materialLawParams(elemIdx),
fs);
@ -427,8 +427,8 @@ inline void testAll()
krnSwMdc_in[1],
elemIdx);
Scalar pcSwMdc_out[2];
Scalar krnSwMdc_out[2];
Scalar pcSwMdc_out[2] = { 0.0, 0.0 };
Scalar krnSwMdc_out[2] = { 0.0, 0.0 };
hysterMaterialLawManager.oilWaterHysteresisParams(
pcSwMdc_out[0],
krnSwMdc_out[0],

12
tests/test_fluidmatrixinteractions.cpp
View File

@ -156,14 +156,14 @@ void testTwoPhaseApi()
const FluidState fs;
const typename MaterialLaw::Params params;
Scalar v;
Scalar v OPM_UNUSED;
v = MaterialLaw::template pcnw<FluidState, Scalar>(params, fs);
v = MaterialLaw::template Sw<FluidState, Scalar>(params, fs);
v = MaterialLaw::template Sn<FluidState, Scalar>(params, fs);
v = MaterialLaw::template krw<FluidState, Scalar>(params, fs);
v = MaterialLaw::template krn<FluidState, Scalar>(params, fs);
typename FluidState::Scalar vEval;
typename FluidState::Scalar vEval OPM_UNUSED;
vEval = MaterialLaw::pcnw(params, fs);
vEval = MaterialLaw::Sw(params, fs);
vEval = MaterialLaw::Sn(params, fs);
@ -198,7 +198,7 @@ void testTwoPhaseSatApi()
const typename MaterialLaw::Params params;
Scalar Sw = 0;
Scalar v;
Scalar v OPM_UNUSED;
v = MaterialLaw::twoPhaseSatPcnw(params, Sw);
v = MaterialLaw::twoPhaseSatSw(params, Sw);
v = MaterialLaw::twoPhaseSatSn(params, Sw);
@ -206,7 +206,7 @@ void testTwoPhaseSatApi()
v = MaterialLaw::twoPhaseSatKrn(params, Sw);
typename FluidState::Scalar SwEval = 0;
typename FluidState::Scalar vEval;
typename FluidState::Scalar vEval OPM_UNUSED;
vEval = MaterialLaw::twoPhaseSatPcnw(params, SwEval);
vEval = MaterialLaw::twoPhaseSatSw(params, SwEval);
vEval = MaterialLaw::twoPhaseSatSn(params, SwEval);
@ -234,7 +234,7 @@ void testThreePhaseApi()
const FluidState fs;
const typename MaterialLaw::Params params;
Scalar v;
Scalar v OPM_UNUSED;
v = MaterialLaw::template pcnw<FluidState, Scalar>(params, fs);
v = MaterialLaw::template Sw<FluidState, Scalar>(params, fs);
v = MaterialLaw::template Sn<FluidState, Scalar>(params, fs);
@ -243,7 +243,7 @@ void testThreePhaseApi()
v = MaterialLaw::template krn<FluidState, Scalar>(params, fs);
v = MaterialLaw::template krg<FluidState, Scalar>(params, fs);
typename FluidState::Scalar vEval;
typename FluidState::Scalar vEval OPM_UNUSED;
vEval = MaterialLaw::pcnw(params, fs);
vEval = MaterialLaw::Sw(params, fs);
vEval = MaterialLaw::Sn(params, fs);