commit
4d253a1a5a
@ -53,7 +53,7 @@ int invertLinearPolynomial(SolContainer &sol,
|
||||
Scalar b)
|
||||
{
|
||||
typedef MathToolbox<Scalar> Toolbox;
|
||||
if (std::abs(Toolbox::value(a)) < 1e-30)
|
||||
if (std::abs(Toolbox::scalarValue(a)) < 1e-30)
|
||||
return 0;
|
||||
|
||||
sol[0] = -b/a;
|
||||
@ -85,7 +85,7 @@ int invertQuadraticPolynomial(SolContainer &sol,
|
||||
typedef MathToolbox<Scalar> Toolbox;
|
||||
|
||||
// check for a line
|
||||
if (std::abs(Toolbox::value(a)) < 1e-30)
|
||||
if (std::abs(Toolbox::scalarValue(a)) < 1e-30)
|
||||
return invertLinearPolynomial(sol, b, c);
|
||||
|
||||
// discriminant
|
||||
@ -121,12 +121,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::value(fPrime)) < 1e-30)
|
||||
if (std::abs(Toolbox::scalarValue(fPrime)) < 1e-30)
|
||||
continue;
|
||||
x -= fOld/fPrime;
|
||||
|
||||
Scalar fNew = d + x*(c + x*(b + x*a));
|
||||
if (std::abs(Toolbox::value(fNew)) < std::abs(Toolbox::value(fOld)))
|
||||
if (std::abs(Toolbox::scalarValue(fNew)) < std::abs(Toolbox::scalarValue(fOld)))
|
||||
sol[i] = x;
|
||||
}
|
||||
}
|
||||
@ -159,7 +159,7 @@ int invertCubicPolynomial(SolContainer *sol,
|
||||
typedef MathToolbox<Scalar> Toolbox;
|
||||
|
||||
// reduces to a quadratic polynomial
|
||||
if (std::abs(Toolbox::value(a)) < 1e-30)
|
||||
if (std::abs(Toolbox::scalarValue(a)) < 1e-30)
|
||||
return invertQuadraticPolynomial(sol, b, c, d);
|
||||
|
||||
// normalize the polynomial
|
||||
@ -172,7 +172,7 @@ int invertCubicPolynomial(SolContainer *sol,
|
||||
Scalar p = c - b*b/3;
|
||||
Scalar q = d + (2*b*b*b - 9*b*c)/27;
|
||||
|
||||
if (std::abs(Toolbox::value(p)) > 1e-30 && std::abs(Toolbox::value(q)) > 1e-30) {
|
||||
if (std::abs(Toolbox::scalarValue(p)) > 1e-30 && std::abs(Toolbox::scalarValue(q)) > 1e-30) {
|
||||
// At this point
|
||||
//
|
||||
// t^3 + p*t + q = 0
|
||||
@ -282,12 +282,12 @@ int invertCubicPolynomial(SolContainer *sol,
|
||||
}
|
||||
// Handle some (pretty unlikely) special cases required to avoid
|
||||
// divisions by zero in the code above...
|
||||
else if (std::abs(Toolbox::value(p)) < 1e-30 && std::abs(Toolbox::value(q)) < 1e-30) {
|
||||
else if (std::abs(Toolbox::scalarValue(p)) < 1e-30 && std::abs(Toolbox::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::value(p)) < 1e-30 && std::abs(Toolbox::value(q)) > 1e-30) {
|
||||
else if (std::abs(Toolbox::scalarValue(p)) < 1e-30 && std::abs(Toolbox::scalarValue(q)) > 1e-30) {
|
||||
// t^3 + q = 0,
|
||||
//
|
||||
// i. e. single real root at t=curt(q)
|
||||
@ -299,7 +299,7 @@ int invertCubicPolynomial(SolContainer *sol,
|
||||
return 1;
|
||||
}
|
||||
|
||||
assert(std::abs(Toolbox::value(p)) > 1e-30 && std::abs(Toolbox::value(q)) > 1e-30);
|
||||
assert(std::abs(Toolbox::scalarValue(p)) > 1e-30 && std::abs(Toolbox::scalarValue(q)) > 1e-30);
|
||||
|
||||
// t^3 + p*t = 0 = t*(t^2 + p),
|
||||
//
|
||||
|
@ -1042,9 +1042,9 @@ protected:
|
||||
size_t n = numSamples();
|
||||
|
||||
// create a vector containing 0...n-1
|
||||
std::vector<unsigned> idxVector(static_cast<unsigned>(n));
|
||||
std::vector<unsigned> idxVector(n);
|
||||
for (unsigned i = 0; i < n; ++i)
|
||||
idxVector[static_cast<unsigned>(i)] = i;
|
||||
idxVector[i] = i;
|
||||
|
||||
// sort the indices according to the x values of the sample
|
||||
// points
|
||||
@ -1148,8 +1148,10 @@ protected:
|
||||
M.solve(moments, d);
|
||||
|
||||
moments.resize(numSamples());
|
||||
for (int i = static_cast<int>(numSamples()) - 2; i >= 0; --i)
|
||||
moments[static_cast<unsigned>(i+1)] = moments[static_cast<unsigned>(i)];
|
||||
for (int i = static_cast<int>(numSamples()) - 2; i >= 0; --i) {
|
||||
unsigned ui = static_cast<unsigned>(i);
|
||||
moments[ui+1] = moments[ui];
|
||||
}
|
||||
moments[0] = moments[numSamples() - 1];
|
||||
|
||||
// convert the moments to slopes at the sample points
|
||||
@ -1731,7 +1733,7 @@ protected:
|
||||
{
|
||||
typedef Opm::MathToolbox<Evaluation> Toolbox;
|
||||
|
||||
Scalar x = Toolbox::value(xEval);
|
||||
Scalar x = Toolbox::scalarValue(xEval);
|
||||
|
||||
// bisection
|
||||
size_t iLow = 0;
|
||||
|
@ -259,7 +259,7 @@ public:
|
||||
segIdx = numSamples() - 2;
|
||||
else {
|
||||
assert(xValues_.front() <= x && x <= xValues_.back());
|
||||
segIdx = findSegmentIndex_(Toolbox::value(x));
|
||||
segIdx = findSegmentIndex_(Toolbox::scalarValue(x));
|
||||
}
|
||||
|
||||
Scalar x0 = xValues_[segIdx];
|
||||
|
@ -198,10 +198,10 @@ public:
|
||||
Evaluation alpha = xToI(x);
|
||||
Evaluation beta = yToJ(y);
|
||||
|
||||
unsigned i = std::max(0U, std::min(static_cast<unsigned>(numX()) - 2,
|
||||
static_cast<unsigned>(Toolbox::scalarValue(alpha))));
|
||||
unsigned j = std::max(0U, std::min(static_cast<unsigned>(numY()) - 2,
|
||||
static_cast<unsigned>(Toolbox::scalarValue(beta))));
|
||||
int i = std::max<int>(0, std::min(static_cast<int>(numX()) - 2,
|
||||
static_cast<int>(Toolbox::scalarValue(alpha))));
|
||||
int j = std::max<int>(0, std::min(static_cast<int>(numY()) - 2,
|
||||
static_cast<int>(Toolbox::scalarValue(beta))));
|
||||
|
||||
alpha -= i;
|
||||
beta -= j;
|
||||
|
@ -228,7 +228,7 @@ public:
|
||||
Scalar i = xToI(x, /*extrapolate=*/false);
|
||||
const auto &col1SamplePoints = samples_.at(unsigned(i));
|
||||
const auto &col2SamplePoints = samples_.at(unsigned(i));
|
||||
Scalar alpha = i - int(i);
|
||||
Scalar alpha = i - static_cast<int>(i);
|
||||
|
||||
Scalar minY =
|
||||
alpha*std::get<1>(col1SamplePoints.front()) +
|
||||
@ -252,7 +252,7 @@ public:
|
||||
typedef Opm::MathToolbox<Evaluation> Toolbox;
|
||||
|
||||
#ifndef NDEBUG
|
||||
if (!extrapolate && !applies(Toolbox::value(x), Toolbox::value(y))) {
|
||||
if (!extrapolate && !applies(Toolbox::scalarValue(x), Toolbox::scalarValue(y))) {
|
||||
OPM_THROW(NumericalProblem,
|
||||
"Attempt to get undefined table value (" << x << ", " << y << ")");
|
||||
};
|
||||
@ -262,8 +262,8 @@ public:
|
||||
// table ...
|
||||
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::value(alpha)))));
|
||||
static_cast<size_t>(std::max(0, std::min(static_cast<int>(numX()) - 2,
|
||||
static_cast<int>(Toolbox::scalarValue(alpha)))));
|
||||
alpha -= i;
|
||||
|
||||
Evaluation beta1;
|
||||
@ -272,15 +272,14 @@ public:
|
||||
beta1 = yToJ(i, y, extrapolate);
|
||||
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::value(beta1)))));
|
||||
size_t j2 = static_cast<size_t>(std::max(0, std::min(static_cast<int>(numY(i + 1) - 2),
|
||||
static_cast<int>(Toolbox::value(beta2)))));
|
||||
size_t j1 = static_cast<size_t>(std::max(0, std::min(static_cast<int>(numY(i)) - 2,
|
||||
static_cast<int>(Toolbox::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)))));
|
||||
|
||||
beta1 -= j1;
|
||||
beta2 -= j2;
|
||||
|
||||
|
||||
// evaluate the two function values for the same y value ...
|
||||
Evaluation s1, s2;
|
||||
s1 = valueAt(i, j1)*(1.0 - beta1) + valueAt(i, j1 + 1)*beta1;
|
||||
|
@ -151,23 +151,23 @@ public:
|
||||
{ +0.17965e-2, +0.71924e-3, -0.4900e-4 }
|
||||
};
|
||||
|
||||
Evaluation theta = temperature - 273.15;
|
||||
const Evaluation& theta = temperature - 273.15;
|
||||
|
||||
Scalar S = salinity;
|
||||
Scalar S_lSAT =
|
||||
Evaluation S = salinity;
|
||||
const Evaluation& S_lSAT =
|
||||
f[0]
|
||||
+ f[1]*Toolbox::value(theta)
|
||||
+ f[2]*std::pow(Toolbox::value(theta), 2)
|
||||
+ f[3]*std::pow(Toolbox::value(theta), 3);
|
||||
+ f[1]*theta
|
||||
+ f[2]*Toolbox::pow(theta, 2)
|
||||
+ f[3]*Toolbox::pow(theta, 3);
|
||||
|
||||
// Regularization
|
||||
if (S > S_lSAT)
|
||||
S = S_lSAT;
|
||||
|
||||
Evaluation hw = H2O::liquidEnthalpy(temperature, pressure)/1e3; // [kJ/kg]
|
||||
const Evaluation& hw = H2O::liquidEnthalpy(temperature, pressure)/1e3; // [kJ/kg]
|
||||
|
||||
// From Daubert and Danner
|
||||
Evaluation h_NaCl =
|
||||
const Evaluation& h_NaCl =
|
||||
(3.6710e4*temperature
|
||||
+ (6.2770e1/2)*temperature*temperature
|
||||
- (6.6670e-2/3)*temperature*temperature*temperature
|
||||
@ -179,14 +179,14 @@ public:
|
||||
Evaluation d_h = 0;
|
||||
for (int i = 0; i<=3; ++i) {
|
||||
for (int j = 0; j <= 2; ++j) {
|
||||
d_h += a[i][j] * std::pow(theta, i) * std::pow(m, j);
|
||||
d_h += a[i][j] * Toolbox::pow(theta, i) * std::pow(m, j);
|
||||
}
|
||||
}
|
||||
|
||||
Evaluation delta_h = 4.184/(1e3 + (58.44 * m))*d_h;
|
||||
const Evaluation& delta_h = 4.184/(1e3 + (58.44 * m))*d_h;
|
||||
|
||||
// Enthalpy of brine
|
||||
Evaluation h_ls = (1-S)*hw + S*h_NaCl + S*delta_h; // [kJ/kg]
|
||||
const Evaluation& h_ls = (1-S)*hw + S*h_NaCl + S*delta_h; // [kJ/kg]
|
||||
return h_ls*1e3; // convert to [J/kg]
|
||||
}
|
||||
|
||||
@ -290,12 +290,12 @@ public:
|
||||
// assume the pressure to be 10% higher than the vapor
|
||||
// pressure
|
||||
Evaluation pressure = 1.1*vaporPressure(temperature);
|
||||
Scalar eps = Toolbox::value(pressure)*1e-7;
|
||||
Scalar eps = Toolbox::scalarValue(pressure)*1e-7;
|
||||
|
||||
Evaluation deltaP = pressure*2;
|
||||
for (int i = 0;
|
||||
i < 5
|
||||
&& std::abs(Toolbox::value(pressure)*1e-9) < std::abs(Toolbox::value(deltaP));
|
||||
&& std::abs(Toolbox::scalarValue(pressure)*1e-9) < std::abs(Toolbox::scalarValue(deltaP));
|
||||
++i)
|
||||
{
|
||||
const Evaluation& f = liquidDensity(temperature, pressure) - density;
|
||||
|
@ -230,8 +230,6 @@ public:
|
||||
static Evaluation liquidEnthalpy(const Evaluation& temperature,
|
||||
const Evaluation& pressure)
|
||||
{
|
||||
typedef MathToolbox<Evaluation> Toolbox;
|
||||
|
||||
if (!Region1::isValid(temperature, pressure))
|
||||
{
|
||||
OPM_THROW(NumericalProblem,
|
||||
@ -247,7 +245,7 @@ public:
|
||||
const Evaluation& dh_dp =
|
||||
Rs * temperature*
|
||||
Region1::tau(temperature)*
|
||||
Region1::dpi_dp(Toolbox::value(pv))*
|
||||
Region1::dpi_dp(pv)*
|
||||
Region1::ddgamma_dtaudpi(temperature, pv);
|
||||
|
||||
return
|
||||
@ -590,7 +588,7 @@ public:
|
||||
|
||||
// calculate the partial derivative of the specific volume
|
||||
// to the pressure at the vapor pressure.
|
||||
Scalar eps = Toolbox::value(pv)*1e-8;
|
||||
Scalar eps = Toolbox::scalarValue(pv)*1e-8;
|
||||
Evaluation v0 = volumeRegion2_(temperature, pv);
|
||||
Evaluation v1 = volumeRegion2_(temperature, pv + eps);
|
||||
Evaluation dv_dp = (v1 - v0)/eps;
|
||||
@ -652,7 +650,7 @@ public:
|
||||
Evaluation deltaP = pressure*2;
|
||||
Valgrind::CheckDefined(pressure);
|
||||
Valgrind::CheckDefined(deltaP);
|
||||
for (int i = 0; i < 5 && std::abs(Toolbox::value(pressure)*1e-9) < std::abs(Toolbox::value(deltaP)); ++i) {
|
||||
for (int i = 0; i < 5 && std::abs(Toolbox::scalarValue(pressure)*1e-9) < std::abs(Toolbox::scalarValue(deltaP)); ++i) {
|
||||
Evaluation f = gasDensity(temperature, pressure) - density;
|
||||
|
||||
Evaluation df_dp;
|
||||
|
@ -547,7 +547,7 @@ private:
|
||||
if (alphaT < 0 || alphaT >= nTemp_ - 1)
|
||||
return std::numeric_limits<Scalar>::quiet_NaN();
|
||||
|
||||
unsigned iT = (unsigned) Toolbox::scalarValue(alphaT);
|
||||
size_t iT = static_cast<size_t>(Toolbox::scalarValue(alphaT));
|
||||
alphaT -= iT;
|
||||
|
||||
return
|
||||
@ -567,20 +567,17 @@ private:
|
||||
return Toolbox::createConstant(std::numeric_limits<Scalar>::quiet_NaN());
|
||||
}
|
||||
|
||||
int iT =
|
||||
std::max<int>(0,
|
||||
std::min<int>(nTemp_ - 2,
|
||||
static_cast<int>(Toolbox::scalarValue(alphaT))));
|
||||
size_t iT = static_cast<size_t>(Toolbox::scalarValue(alphaT));
|
||||
alphaT -= iT;
|
||||
|
||||
Evaluation alphaP1 = pressLiquidIdx_(p, iT);
|
||||
Evaluation alphaP2 = pressLiquidIdx_(p, iT + 1);
|
||||
|
||||
int iP1 =
|
||||
size_t iP1 =
|
||||
std::max<int>(0,
|
||||
std::min<int>(nPress_ - 2,
|
||||
static_cast<int>(Toolbox::scalarValue(alphaP1))));
|
||||
int iP2 =
|
||||
std::min<int>(nPress_ - 2,
|
||||
static_cast<int>(Toolbox::scalarValue(alphaP1))));
|
||||
size_t iP2 =
|
||||
std::max<int>(0,
|
||||
std::min<int>(nPress_ - 2,
|
||||
static_cast<int>(Toolbox::scalarValue(alphaP2))));
|
||||
@ -618,7 +615,7 @@ private:
|
||||
return Toolbox::createConstant(std::numeric_limits<Scalar>::quiet_NaN());
|
||||
}
|
||||
|
||||
int iT =
|
||||
size_t iT =
|
||||
std::max<int>(0,
|
||||
std::min<int>(nTemp_ - 2,
|
||||
static_cast<int>(Toolbox::scalarValue(alphaT))));
|
||||
@ -626,10 +623,10 @@ private:
|
||||
|
||||
Evaluation alphaP1 = pressGasIdx_(p, iT);
|
||||
Evaluation alphaP2 = pressGasIdx_(p, iT + 1);
|
||||
int iP1 =
|
||||
size_t iP1 =
|
||||
std::max<int>(0, std::min<int>(nPress_ - 2,
|
||||
static_cast<int>(Toolbox::scalarValue(alphaP1))));
|
||||
int iP2 =
|
||||
size_t iP2 =
|
||||
std::max<int>(0,
|
||||
std::min<int>(nPress_ - 2,
|
||||
static_cast<int>(Toolbox::scalarValue(alphaP2))));
|
||||
@ -661,7 +658,7 @@ private:
|
||||
static Evaluation interpolateGasTRho_(const Scalar *values, const Evaluation& T, const Evaluation& rho)
|
||||
{
|
||||
Evaluation alphaT = tempIdx_(T);
|
||||
int iT = std::max<int>(0, std::min<int>(nTemp_ - 2, (int) alphaT));
|
||||
unsigned iT = std::max<int>(0, std::min<int>(nTemp_ - 2, (int) alphaT));
|
||||
alphaT -= iT;
|
||||
|
||||
Evaluation alphaP1 = densityGasIdx_(rho, iT);
|
||||
@ -684,13 +681,13 @@ private:
|
||||
static Evaluation interpolateLiquidTRho_(const Scalar *values, const Evaluation& T, const Evaluation& rho)
|
||||
{
|
||||
Evaluation alphaT = tempIdx_(T);
|
||||
int iT = std::max<int>(0, std::min<int>(nTemp_ - 2, (int) alphaT));
|
||||
unsigned iT = std::max<int>(0, std::min<int>(nTemp_ - 2, static_cast<int>(alphaT)));
|
||||
alphaT -= iT;
|
||||
|
||||
Evaluation alphaP1 = densityLiquidIdx_(rho, iT);
|
||||
Evaluation alphaP2 = densityLiquidIdx_(rho, iT + 1);
|
||||
unsigned iP1 = std::max<int>(0, std::min<int>(nDensity_ - 2, (int) alphaP1));
|
||||
unsigned iP2 = std::max<int>(0, std::min<int>(nDensity_ - 2, (int) alphaP2));
|
||||
unsigned iP1 = std::max<int>(0, std::min<int>(nDensity_ - 2, static_cast<int>(alphaP1)));
|
||||
unsigned iP2 = std::max<int>(0, std::min<int>(nDensity_ - 2, static_cast<int>(alphaP2)));
|
||||
alphaP1 -= iP1;
|
||||
alphaP2 -= iP2;
|
||||
|
||||
|
@ -63,8 +63,8 @@ public:
|
||||
typedef MathToolbox<Evaluation> Toolbox;
|
||||
|
||||
return
|
||||
Toolbox::value(temperature) <= 623.15 &&
|
||||
Toolbox::value(pressure) <= 100e6;
|
||||
Toolbox::scalarValue(temperature) <= 623.15 &&
|
||||
Toolbox::scalarValue(pressure) <= 100e6;
|
||||
|
||||
// actually this is:
|
||||
/*
|
||||
@ -112,7 +112,8 @@ public:
|
||||
*
|
||||
* \param pressure temperature of component in \f$\mathrm{[Pa]}\f$
|
||||
*/
|
||||
static Scalar dpi_dp(Scalar /*pressure*/)
|
||||
template <class Evaluation>
|
||||
static Scalar dpi_dp(const Evaluation& /*pressure*/)
|
||||
{ return 1.0 / 16.53e6; }
|
||||
|
||||
/*!
|
||||
@ -121,7 +122,8 @@ public:
|
||||
*
|
||||
* \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
|
||||
*/
|
||||
static Scalar dp_dpi(Scalar /*pressure*/)
|
||||
template <class Evaluation>
|
||||
static Scalar dp_dpi(const Evaluation& /*pressure*/)
|
||||
{ return 16.53e6; }
|
||||
|
||||
/*!
|
||||
|
@ -241,7 +241,7 @@ protected:
|
||||
fluidState.setFugacityCoefficient(phaseIdx, i, phi);
|
||||
|
||||
defect[i] = targetFug[i] - f;
|
||||
absError = std::max(absError, std::abs(Toolbox::value(defect[i])));
|
||||
absError = std::max(absError, std::abs(Toolbox::scalarValue(defect[i])));
|
||||
}
|
||||
|
||||
// assemble jacobian matrix of the constraints for the composition
|
||||
@ -307,7 +307,7 @@ protected:
|
||||
Evaluation sumx = Toolbox::createConstant(0.0);
|
||||
for (unsigned i = 0; i < numComponents; ++i) {
|
||||
origComp[i] = fluidState.moleFraction(phaseIdx, i);
|
||||
relError = std::max(relError, std::abs(Toolbox::value(x[i])));
|
||||
relError = std::max(relError, std::abs(Toolbox::scalarValue(x[i])));
|
||||
|
||||
sumx += Toolbox::abs(fluidState.moleFraction(phaseIdx, i));
|
||||
sumDelta += Toolbox::abs(x[i]);
|
||||
|
@ -306,9 +306,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::value(Sw_ow) - Swco < epsilon) {
|
||||
if (Toolbox::scalarValue(Sw_ow) - Swco < epsilon) {
|
||||
Evaluation kro2 = (kro_ow + kro_go)/2;;
|
||||
if (Toolbox::value(Sw_ow) - Swco > epsilon/2) {
|
||||
if (Toolbox::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);
|
||||
@ -332,9 +332,9 @@ public:
|
||||
{
|
||||
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
|
||||
|
||||
Scalar Sw = FsToolbox::value(fluidState.saturation(waterPhaseIdx));
|
||||
Scalar So = FsToolbox::value(fluidState.saturation(oilPhaseIdx));
|
||||
Scalar Sg = FsToolbox::value(fluidState.saturation(gasPhaseIdx));
|
||||
Scalar Sw = FsToolbox::scalarValue(fluidState.saturation(waterPhaseIdx));
|
||||
Scalar So = FsToolbox::scalarValue(fluidState.saturation(oilPhaseIdx));
|
||||
Scalar Sg = FsToolbox::scalarValue(fluidState.saturation(gasPhaseIdx));
|
||||
|
||||
if (params.inconsistentHysteresisUpdate()) {
|
||||
Sg = std::min(Scalar(1.0), std::max(Scalar(0.0), Sg));
|
||||
|
@ -117,7 +117,8 @@ public:
|
||||
const std::vector<int>& compressedToCartesianElemIdx)
|
||||
{
|
||||
// get the number of saturation regions and the number of cells in the deck
|
||||
unsigned numSatRegions = static_cast<unsigned>(deck->getKeyword("TABDIMS").getRecord(0).getItem("NTSFUN").get< int >(0));
|
||||
int ntsFun = deck->getKeyword("TABDIMS").getRecord(0).getItem("NTSFUN").get<int>(0);
|
||||
unsigned numSatRegions = static_cast<unsigned>(ntsFun);
|
||||
size_t numCompressedElems = compressedToCartesianElemIdx.size();
|
||||
|
||||
// copy the SATNUM grid property. in some cases this is not necessary, but it
|
||||
@ -355,7 +356,7 @@ private:
|
||||
const std::vector<int>& compressedToCartesianElemIdx,
|
||||
const std::vector<int>& satnumRegionArray)
|
||||
{
|
||||
unsigned numSatRegions = static_cast<unsigned>(deck->getKeyword("TABDIMS").getRecord(0).getItem("NTSFUN").get< int >(0));
|
||||
unsigned numSatRegions = static_cast<unsigned>(deck->getKeyword("TABDIMS").getRecord(0).getItem("NTSFUN").get<int>(0));
|
||||
unsigned numCompressedElems = static_cast<unsigned>(compressedToCartesianElemIdx.size());
|
||||
|
||||
// read the end point scaling configuration. this needs to be done only once per
|
||||
|
@ -337,8 +337,8 @@ public:
|
||||
{
|
||||
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
|
||||
|
||||
Scalar Sw = FsToolbox::value(fluidState.saturation(waterPhaseIdx));
|
||||
Scalar Sg = FsToolbox::value(fluidState.saturation(gasPhaseIdx));
|
||||
Scalar Sw = FsToolbox::scalarValue(fluidState.saturation(waterPhaseIdx));
|
||||
Scalar Sg = FsToolbox::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);
|
||||
|
@ -316,8 +316,8 @@ public:
|
||||
{
|
||||
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
|
||||
|
||||
Scalar Sw = FsToolbox::value(fluidState.saturation(waterPhaseIdx));
|
||||
Scalar Sg = FsToolbox::value(fluidState.saturation(gasPhaseIdx));
|
||||
Scalar Sw = FsToolbox::scalarValue(fluidState.saturation(waterPhaseIdx));
|
||||
Scalar Sg = FsToolbox::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);
|
||||
|
@ -334,7 +334,7 @@ public:
|
||||
|
||||
switch (params.approach()) {
|
||||
case EclTwoPhaseGasOil: {
|
||||
Scalar So = FsToolbox::value(fluidState.saturation(oilPhaseIdx));
|
||||
Scalar So = FsToolbox::scalarValue(fluidState.saturation(oilPhaseIdx));
|
||||
|
||||
params.oilWaterParams().update(/*pcSw=*/0.0, /*krwSw=*/0.0, /*krnSw=*/0.0);
|
||||
params.gasOilParams().update(/*pcSw=*/So, /*krwSw=*/So, /*krnSw=*/So);
|
||||
@ -342,7 +342,7 @@ public:
|
||||
}
|
||||
|
||||
case EclTwoPhaseOilWater: {
|
||||
Scalar Sw = FsToolbox::value(fluidState.saturation(waterPhaseIdx));
|
||||
Scalar Sw = FsToolbox::scalarValue(fluidState.saturation(waterPhaseIdx));
|
||||
|
||||
params.oilWaterParams().update(/*pcSw=*/Sw, /*krwSw=*/Sw, /*krnSw=*/Sw);
|
||||
params.gasOilParams().update(/*pcSw=*/0.0, /*krwSw=*/0.0, /*krnSw=*/0.0);
|
||||
@ -350,7 +350,7 @@ public:
|
||||
}
|
||||
|
||||
case EclTwoPhaseGasWater: {
|
||||
Scalar Sw = FsToolbox::value(fluidState.saturation(waterPhaseIdx));
|
||||
Scalar Sw = FsToolbox::scalarValue(fluidState.saturation(waterPhaseIdx));
|
||||
|
||||
params.oilWaterParams().update(/*pcSw=*/Sw, /*krwSw=*/Sw, /*krnSw=*/0);
|
||||
params.gasOilParams().update(/*pcSw=*/1.0, /*krwSw=*/0.0, /*krnSw=*/Sw);
|
||||
|
@ -304,7 +304,7 @@ public:
|
||||
{
|
||||
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
|
||||
|
||||
Scalar Sw = FsToolbox::value(fs.saturation(Traits::wettingPhaseIdx));
|
||||
Scalar Sw = FsToolbox::scalarValue(fs.saturation(Traits::wettingPhaseIdx));
|
||||
|
||||
if (Sw > 1 - 1e-5) {
|
||||
// if the absolute saturation is almost 1,
|
||||
@ -392,7 +392,7 @@ public:
|
||||
typedef MathToolbox<Evaluation> Toolbox;
|
||||
|
||||
// calculate the current apparent saturation
|
||||
ScanningCurve *sc = findScanningCurve_(params, Toolbox::value(Sw));
|
||||
ScanningCurve *sc = findScanningCurve_(params, Toolbox::scalarValue(Sw));
|
||||
|
||||
// calculate the apparant saturation
|
||||
const Evaluation& Sw_app = absoluteToApparentSw_(params, Sw);
|
||||
|
@ -239,7 +239,7 @@ private:
|
||||
if (x >= xValues.back())
|
||||
return yValues.back();
|
||||
|
||||
size_t segIdx = findSegmentIndex_(xValues, Toolbox::value(x));
|
||||
size_t segIdx = findSegmentIndex_(xValues, Toolbox::scalarValue(x));
|
||||
|
||||
Scalar x0 = xValues[segIdx];
|
||||
Scalar x1 = xValues[segIdx + 1];
|
||||
@ -264,7 +264,7 @@ private:
|
||||
if (x <= xValues.back())
|
||||
return yValues.back();
|
||||
|
||||
size_t segIdx = findSegmentIndexDescending_(xValues, Toolbox::value(x));
|
||||
size_t segIdx = findSegmentIndexDescending_(xValues, Toolbox::scalarValue(x));
|
||||
|
||||
Scalar x0 = xValues[segIdx];
|
||||
Scalar x1 = xValues[segIdx + 1];
|
||||
@ -284,12 +284,12 @@ private:
|
||||
{
|
||||
typedef MathToolbox<Evaluation> Toolbox;
|
||||
|
||||
if (Toolbox::value(x) <= xValues.front())
|
||||
if (Toolbox::scalarValue(x) <= xValues.front())
|
||||
return Toolbox::createConstant(0.0);
|
||||
if (Toolbox::value(x) >= xValues.back())
|
||||
if (Toolbox::scalarValue(x) >= xValues.back())
|
||||
return Toolbox::createConstant(0.0);
|
||||
|
||||
size_t segIdx = findSegmentIndex_(xValues, Toolbox::value(x));
|
||||
size_t segIdx = findSegmentIndex_(xValues, Toolbox::scalarValue(x));
|
||||
|
||||
Scalar x0 = xValues[segIdx];
|
||||
Scalar x1 = xValues[segIdx + 1];
|
||||
|
@ -132,8 +132,8 @@ public:
|
||||
#ifndef NDEBUG
|
||||
typedef Opm::MathToolbox<Scalar> Toolbox;
|
||||
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
|
||||
assert(std::abs(FsToolbox::value(fs.temperature(phaseIdx))
|
||||
- Toolbox::value(temperature_)) < 1e-30);
|
||||
assert(std::abs(FsToolbox::scalarValue(fs.temperature(phaseIdx))
|
||||
- Toolbox::scalarValue(temperature_)) < 1e-30);
|
||||
}
|
||||
#endif
|
||||
}
|
||||
|
@ -579,7 +579,7 @@ private:
|
||||
theta = T - 273.15;
|
||||
|
||||
// Regularization
|
||||
Scalar scalarTheta = LhsToolbox::value(theta);
|
||||
Scalar scalarTheta = LhsToolbox::scalarValue(theta);
|
||||
Scalar S_lSAT = f[0] + scalarTheta*(f[1] + scalarTheta*(f[2] + scalarTheta*f[3]));
|
||||
if (S > S_lSAT)
|
||||
S = S_lSAT;
|
||||
|
@ -492,8 +492,8 @@ public:
|
||||
const Evaluation& delta = f/fPrime;
|
||||
pSat -= delta;
|
||||
|
||||
Scalar absDelta = std::abs(Toolbox::value(delta));
|
||||
if (absDelta < Toolbox::value(pSat) * 1e-10 || absDelta < 1e-4)
|
||||
Scalar absDelta = std::abs(Toolbox::scalarValue(delta));
|
||||
if (absDelta < Toolbox::scalarValue(pSat) * 1e-10 || absDelta < 1e-4)
|
||||
return pSat;
|
||||
}
|
||||
|
||||
|
@ -518,7 +518,7 @@ public:
|
||||
const Evaluation& delta = f/fPrime;
|
||||
pSat -= delta;
|
||||
|
||||
if (std::abs(Toolbox::value(delta)) < std::abs(Toolbox::value(pSat)) * 1e-10)
|
||||
if (std::abs(Toolbox::scalarValue(delta)) < std::abs(Toolbox::scalarValue(pSat)) * 1e-10)
|
||||
return pSat;
|
||||
}
|
||||
|
||||
|
@ -171,7 +171,7 @@ public:
|
||||
}
|
||||
|
||||
template <class RhsValueType>
|
||||
Evaluation& operator*=(const RhsValueType& other)
|
||||
Evaluation& operator*=(RhsValueType other)
|
||||
{
|
||||
// values and derivatives are multiplied
|
||||
this->value *= other;
|
||||
@ -202,9 +202,10 @@ public:
|
||||
Evaluation& operator/=(const RhsValueType& other)
|
||||
{
|
||||
// values and derivatives are divided
|
||||
this->value /= other;
|
||||
auto tmp = 1.0/other;
|
||||
this->value *= tmp;
|
||||
for (unsigned varIdx = 0; varIdx < size; ++varIdx)
|
||||
this->derivatives[varIdx] /= other;
|
||||
this->derivatives[varIdx] *= tmp;
|
||||
|
||||
return *this;
|
||||
}
|
||||
|
Loading…
Reference in New Issue
Block a user