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

DenseAD: access the value and derivatives of Evaluations via accessor functions

Browse Source 
this allows some flexibility of how an Evaluation object is
represented internally. the main motivation for this is to allow using
SIMD instruction which expect special data types as their operants.
This commit is contained in:
Andreas Lauser
2016-10-20 17:11:12 +02:00
parent e047e85f54
commit 66ec787438
5 changed files with 221 additions and 204 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -190,9 +190,9 @@ public:
// the linear system of equations.
for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
for (unsigned pvIdx = 0; pvIdx < numEq; ++ pvIdx)
J[eqIdx][pvIdx] = defect[eqIdx].derivatives[pvIdx];
J[eqIdx][pvIdx] = defect[eqIdx].derivative(pvIdx);
b[eqIdx] = defect[eqIdx].value;
b[eqIdx] = defect[eqIdx].value();
}
Valgrind::CheckDefined(J);
Valgrind::CheckDefined(b);
@@ -274,7 +274,7 @@ protected:
FlashEval Slast = InputToolbox::createConstant(1.0);
for (unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx) {
FlashEval S = inputFluidState.saturation(phaseIdx);
S.derivatives[S0PvIdx + phaseIdx] = 1.0;
S.setDerivative(S0PvIdx + phaseIdx, 1.0);
Slast -= S;
@@ -285,7 +285,7 @@ protected:
// copy the pressures: the first pressure is the first primary variable, the
// remaining ones are given as p_beta = p_alpha + p_calpha,beta
FlashEval p0 = inputFluidState.pressure(0);
p0.derivatives[p0PvIdx] = 1.0;
p0.setDerivative(p0PvIdx, 1.0);
std::array<FlashEval, numPhases> pc;
MaterialLaw::capillaryPressures(pc, matParams, flashFluidState);
@@ -305,17 +305,17 @@ protected:
static void assignOutputFluidState_(const FlashFluidState& flashFluidState,
OutputFluidState& outputFluidState)
{
outputFluidState.setTemperature(flashFluidState.temperature(/*phaseIdx=*/0).value);
outputFluidState.setTemperature(flashFluidState.temperature(/*phaseIdx=*/0).value());
// copy the saturations, pressures and densities
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
const auto& S = flashFluidState.saturation(phaseIdx).value;
const auto& S = flashFluidState.saturation(phaseIdx).value();
outputFluidState.setSaturation(phaseIdx, S);
const auto& p = flashFluidState.pressure(phaseIdx).value;
const auto& p = flashFluidState.pressure(phaseIdx).value();
outputFluidState.setPressure(phaseIdx, p);
const auto& rho = flashFluidState.density(phaseIdx).value;
const auto& rho = flashFluidState.density(phaseIdx).value();
outputFluidState.setDensity(phaseIdx, rho);
}
}
@@ -387,8 +387,8 @@ protected:
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 = InnerEvalToolbox::min(0.5*fluidState.pressure(0).value(),
InnerEvalToolbox::max(-0.50*fluidState.pressure(0).value(), delta));
}
tmp -= delta;

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

@@ -209,9 +209,9 @@ public:
// the linear system of equations.
for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
for (unsigned pvIdx = 0; pvIdx < numEq; ++ pvIdx)
J[eqIdx][pvIdx] = defect[eqIdx].derivatives[pvIdx];
J[eqIdx][pvIdx] = defect[eqIdx].derivative(pvIdx);
b[eqIdx] = defect[eqIdx].value;
b[eqIdx] = defect[eqIdx].value();
}
Valgrind::CheckDefined(J);
Valgrind::CheckDefined(b);
@@ -330,7 +330,7 @@ protected:
FlashEval Slast = InputToolbox::createConstant(1.0);
for (unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx) {
FlashEval S = inputFluidState.saturation(phaseIdx);
S.derivatives[S0PvIdx + phaseIdx] = 1.0;
S.setDerivative(S0PvIdx + phaseIdx, 1.0);
Slast -= S;
@@ -341,7 +341,7 @@ protected:
// copy the pressures: the first pressure is the first primary variable, the
// remaining ones are given as p_beta = p_alpha + p_calpha,beta
FlashEval p0 = inputFluidState.pressure(0);
p0.derivatives[p0PvIdx] = 1.0;
p0.setDerivative(p0PvIdx, 1.0);
std::array<FlashEval, numPhases> pc;
MaterialLaw::capillaryPressures(pc, matParams, flashFluidState);
@@ -352,7 +352,7 @@ protected:
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
FlashEval x = inputFluidState.moleFraction(phaseIdx, compIdx);
x.derivatives[x00PvIdx + phaseIdx*numComponents + compIdx] = 1.0;
x.setDerivative(x00PvIdx + phaseIdx*numComponents + compIdx, 1.0);
flashFluidState.setMoleFraction(phaseIdx, compIdx, x);
}
}
@@ -376,17 +376,17 @@ protected:
static void assignOutputFluidState_(const FlashFluidState& flashFluidState,
OutputFluidState& outputFluidState)
{
outputFluidState.setTemperature(flashFluidState.temperature(/*phaseIdx=*/0).value);
outputFluidState.setTemperature(flashFluidState.temperature(/*phaseIdx=*/0).value());
// copy the saturations, pressures and densities
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
const auto& S = flashFluidState.saturation(phaseIdx).value;
const auto& S = flashFluidState.saturation(phaseIdx).value();
outputFluidState.setSaturation(phaseIdx, S);
const auto& p = flashFluidState.pressure(phaseIdx).value;
const auto& p = flashFluidState.pressure(phaseIdx).value();
outputFluidState.setPressure(phaseIdx, p);
const auto& rho = flashFluidState.density(phaseIdx).value;
const auto& rho = flashFluidState.density(phaseIdx).value();
outputFluidState.setDensity(phaseIdx, rho);
}
@@ -394,11 +394,11 @@ protected:
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
const auto& moleFrac =
flashFluidState.moleFraction(phaseIdx, compIdx).value;
flashFluidState.moleFraction(phaseIdx, compIdx).value();
outputFluidState.setMoleFraction(phaseIdx, compIdx, moleFrac);
const auto& fugCoeff =
flashFluidState.fugacityCoefficient(phaseIdx, compIdx).value;
flashFluidState.fugacityCoefficient(phaseIdx, compIdx).value();
outputFluidState.setFugacityCoefficient(phaseIdx, compIdx, fugCoeff);
}
}
@@ -493,8 +493,8 @@ protected:
}
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 = InnerEvalToolbox::min(0.5*fluidState.pressure(0).value(),
InnerEvalToolbox::max(-0.5*fluidState.pressure(0).value(),
delta));
}

150
opm/material/densead/Evaluation.hpp
View File

@@ -55,15 +55,15 @@ public:
enum { size = numVars };
Evaluation() : value(), derivatives()
Evaluation() : value_(), derivatives_()
{}
// copy other function evaluation
Evaluation(const Evaluation& other)
{
// copy evaluated function value and derivatives
value = other.value;
std::copy(other.derivatives.begin(), other.derivatives.end(), derivatives.begin());
value_ = other.value_;
std::copy(other.derivatives_.begin(), other.derivatives_.end(), derivatives_.begin());
}
// create an evaluation which represents a constant function
@@ -73,8 +73,8 @@ public:
template <class RhsValueType>
Evaluation(const RhsValueType& c)
{
value = c;
std::fill(derivatives.begin(), derivatives.end(), 0.0);
value_ = c;
std::fill(derivatives_.begin(), derivatives_.end(), 0.0);
}
// create a function evaluation for a "naked" depending variable (i.e., f(x) = x)
@@ -88,9 +88,9 @@ public:
// copy function value and set all derivatives to 0, except for the variable
// which is represented by the value (which is set to 1.0)
result.value = value;
std::fill(result.derivatives.begin(), result.derivatives.end(), 0.0);
result.derivatives[varPos] = 1.0;
result.value_ = value;
std::fill(result.derivatives_.begin(), result.derivatives_.end(), 0.0);
result.derivatives_[varPos] = 1.0;
return result;
}
@@ -101,19 +101,19 @@ public:
static Evaluation createConstant(const RhsValueType& value)
{
Evaluation result;
result.value = value;
std::fill(result.derivatives.begin(), result.derivatives.end(), 0.0);
Valgrind::CheckDefined(result.value);
Valgrind::CheckDefined(result.derivatives);
result.value_ = value;
std::fill(result.derivatives_.begin(), result.derivatives_.end(), 0.0);
Valgrind::CheckDefined(result.value_);
Valgrind::CheckDefined(result.derivatives_);
return result;
}
// print the value and the derivatives of the function evaluation
void print(std::ostream& os = std::cout) const
{
os << "v: " << value << " / d:";
for (int varIdx = 0; varIdx < derivatives.size(); ++varIdx)
os << " " << derivatives[varIdx];
os << "v: " << value_ << " / d:";
for (int varIdx = 0; varIdx < size; ++varIdx)
os << " " << derivatives_[varIdx];
}
void clearDerivatives()
@@ -129,9 +129,9 @@ public:
Evaluation& operator+=(const Evaluation& other)
{
// value and derivatives are added
this->value += other.value;
this->value_ += other.value_;
for (unsigned varIdx = 0; varIdx < size; ++varIdx)
this->derivatives[varIdx] += other.derivatives[varIdx];
this->derivatives_[varIdx] += other.derivatives_[varIdx];
return *this;
}
@@ -140,7 +140,7 @@ public:
Evaluation& operator+=(const RhsValueType& other)
{
// value is added, derivatives stay the same
this->value += other;
this->value_ += other;
return *this;
}
@@ -148,9 +148,9 @@ public:
Evaluation& operator-=(const Evaluation& other)
{
// value and derivatives are subtracted
this->value -= other.value;
this->value_ -= other.value_;
for (unsigned varIdx = 0; varIdx < size; ++varIdx)
this->derivatives[varIdx] -= other.derivatives[varIdx];
this->derivatives_[varIdx] -= other.derivatives_[varIdx];
return *this;
}
@@ -159,7 +159,7 @@ public:
Evaluation& operator-=(const RhsValueType& other)
{
// for constants, values are subtracted, derivatives stay the same
this->value -= other;
this->value_ -= other;
return *this;
}
@@ -168,15 +168,15 @@ public:
{
// while the values are multiplied, the derivatives follow the product rule,
// i.e., (u*v)' = (v'u + u'v).
const ValueType& u = this->value;
const ValueType& v = other.value;
const ValueType& u = this->value_;
const ValueType& v = other.value_;
for (unsigned varIdx = 0; varIdx < size; ++varIdx) {
const ValueType& uPrime = this->derivatives[varIdx];
const ValueType& vPrime = other.derivatives[varIdx];
const ValueType& uPrime = this->derivatives_[varIdx];
const ValueType& vPrime = other.derivatives_[varIdx];
this->derivatives[varIdx] = (v*uPrime + u*vPrime);
this->derivatives_[varIdx] = (v*uPrime + u*vPrime);
}
this->value *= v;
this->value_ *= v;
return *this;
}
@@ -185,9 +185,9 @@ public:
Evaluation& operator*=(RhsValueType other)
{
// values and derivatives are multiplied
this->value *= other;
this->value_ *= other;
for (unsigned varIdx = 0; varIdx < size; ++varIdx)
this->derivatives[varIdx] *= other;
this->derivatives_[varIdx] *= other;
return *this;
}
@@ -196,15 +196,15 @@ public:
{
// values are divided, derivatives follow the rule for division, i.e., (u/v)' = (v'u -
// u'v)/v^2.
const ValueType& u = this->value;
const ValueType& v = other.value;
const ValueType& u = this->value_;
const ValueType& v = other.value_;
for (unsigned varIdx = 0; varIdx < size; ++varIdx) {
const ValueType& uPrime = this->derivatives[varIdx];
const ValueType& vPrime = other.derivatives[varIdx];
const ValueType& uPrime = this->derivatives_[varIdx];
const ValueType& vPrime = other.derivatives_[varIdx];
this->derivatives[varIdx] = (v*uPrime - u*vPrime)/(v*v);
this->derivatives_[varIdx] = (v*uPrime - u*vPrime)/(v*v);
}
this->value /= v;
this->value_ /= v;
return *this;
}
@@ -214,9 +214,9 @@ public:
{
// values and derivatives are divided
auto tmp = 1.0/other;
this->value *= tmp;
this->value_ *= tmp;
for (unsigned varIdx = 0; varIdx < size; ++varIdx)
this->derivatives[varIdx] *= tmp;
this->derivatives_[varIdx] *= tmp;
return *this;
}
@@ -255,9 +255,9 @@ public:
Evaluation operator-() const
{
Evaluation result;
result.value = -this->value;
result.value_ = -this->value_;
for (unsigned varIdx = 0; varIdx < size; ++varIdx)
result.derivatives[varIdx] = - this->derivatives[varIdx];
result.derivatives_[varIdx] = - this->derivatives_[varIdx];
return result;
}
@@ -295,30 +295,30 @@ public:
template <class RhsValueType>
Evaluation& operator=(const RhsValueType& other)
{
this->value = other;
std::fill(this->derivatives.begin(), this->derivatives.end(), 0.0);
this->value_ = other;
std::fill(this->derivatives_.begin(), this->derivatives_.end(), 0.0);
return *this;
}
// copy assignment from evaluation
Evaluation& operator=(const Evaluation& other)
{
this->value = other.value;
std::copy(other.derivatives.begin(), other.derivatives.end(), this->derivatives.begin());
this->value_ = other.value_;
std::copy(other.derivatives_.begin(), other.derivatives_.end(), this->derivatives_.begin());
return *this;
}
template <class RhsValueType>
bool operator==(const RhsValueType& other) const
{ return this->value == other; }
{ return this->value_ == other; }
bool operator==(const Evaluation& other) const
{
if (this->value != other.value)
if (this->value_ != other.value_)
return false;
for (unsigned varIdx = 0; varIdx < size; ++varIdx)
if (this->derivatives[varIdx] != other.derivatives[varIdx])
if (this->derivatives_[varIdx] != other.derivatives_[varIdx])
return false;
return true;
@@ -329,35 +329,53 @@ public:
template <class RhsValueType>
bool operator>(RhsValueType other) const
{ return this->value > other; }
{ return this->value_ > other; }
bool operator>(const Evaluation& other) const
{ return this->value > other.value; }
{ return this->value_ > other.value_; }
template <class RhsValueType>
bool operator<(RhsValueType other) const
{ return this->value < other; }
{ return this->value_ < other; }
bool operator<(const Evaluation& other) const
{ return this->value < other.value; }
{ return this->value_ < other.value_; }
template <class RhsValueType>
bool operator>=(RhsValueType other) const
{ return this->value >= other; }
{ return this->value_ >= other; }
bool operator>=(const Evaluation& other) const
{ return this->value >= other.value; }
{ return this->value_ >= other.value_; }
template <class RhsValueType>
bool operator<=(RhsValueType other) const
{ return this->value <= other; }
{ return this->value_ <= other; }
bool operator<=(const Evaluation& other) const
{ return this->value <= other.value; }
{ return this->value_ <= other.value_; }
// maybe this should be made 'private'...
ValueType value;
std::array<ValueType, size> derivatives;
ValueType value() const
{ return value_; }
void setValue(const ValueType& val)
{ value_ = val; }
ValueType derivative(unsigned varIdx) const
{
assert(varIdx < numVars);
return derivatives_[varIdx];
}
void setDerivative(unsigned varIdx, const ValueType& derVal)
{
assert(varIdx < numVars);
derivatives_[varIdx] = derVal;
}
private:
ValueType value_;
std::array<ValueType, size> derivatives_;
};
template <class RhsValueType, class ValueType, int numVars>
@@ -378,7 +396,7 @@ bool operator>=(const RhsValueType& a, const Evaluation<ValueType, numVars> &b)
template <class RhsValueType, class ValueType, int numVars>
bool operator!=(const RhsValueType& a, const Evaluation<ValueType, numVars> &b)
{ return a != b.value; }
{ return a != b.value(); }
template <class RhsValueType, class ValueType, int numVars>
Evaluation<ValueType, numVars> operator+(const RhsValueType& a, const Evaluation<ValueType, numVars> &b)
@@ -395,9 +413,9 @@ Evaluation<ValueType, numVars> operator-(const RhsValueType& a, const Evaluation
{
Evaluation<ValueType, numVars> result;
result.value = a - b.value;
result.setValue(a - b.value());
for (unsigned varIdx = 0; varIdx < numVars; ++varIdx)
result.derivatives[varIdx] = - b.derivatives[varIdx];
result.setDerivative(varIdx, - b.derivative(varIdx));
return result;
}
@@ -407,12 +425,12 @@ Evaluation<ValueType, numVars> operator/(const RhsValueType& a, const Evaluation
{
Evaluation<ValueType, numVars> result;
result.value = a/b.value;
result.setValue(a/b.value());
// outer derivative
const ValueType& df_dg = - a/(b.value*b.value);
const ValueType& df_dg = - a/(b.value()*b.value());
for (unsigned varIdx = 0; varIdx < numVars; ++varIdx)
result.derivatives[varIdx] = df_dg*b.derivatives[varIdx];
result.setDerivative(varIdx, df_dg*b.derivative(varIdx));
return result;
}
@@ -422,9 +440,9 @@ Evaluation<ValueType, numVars> operator*(const RhsValueType& a, const Evaluation
{
Evaluation<ValueType, numVars> result;
result.value = a*b.value;
result.setValue(a*b.value());
for (unsigned varIdx = 0; varIdx < numVars; ++varIdx)
result.derivatives[varIdx] = a*b.derivatives[varIdx];
result.setDerivative(varIdx, a*b.derivative(varIdx));
return result;
}
@@ -432,7 +450,7 @@ Evaluation<ValueType, numVars> operator*(const RhsValueType& a, const Evaluation
template <class ValueType, int numVars>
std::ostream& operator<<(std::ostream& os, const Evaluation<ValueType, numVars>& eval)
{
os << eval.value;
os << eval.value();
return os;
}

129
opm/material/densead/Math.hpp
View File

@@ -84,13 +84,13 @@ Evaluation<ValueType, numVars> tan(const Evaluation<ValueType, numVars>& x)
Evaluation<ValueType, numVars> result;
const ValueType& tmp = ValueTypeToolbox::tan(x.value);
result.value = tmp;
const ValueType& tmp = ValueTypeToolbox::tan(x.value());
result.setValue(tmp);
// derivatives use the chain rule
const ValueType& df_dx = 1 + tmp*tmp;
for (unsigned curVarIdx = 0; curVarIdx < result.derivatives.size(); ++curVarIdx)
result.derivatives[curVarIdx] = df_dx*x.derivatives[curVarIdx];
for (unsigned curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
return result;
}
@@ -102,12 +102,12 @@ Evaluation<ValueType, numVars> atan(const Evaluation<ValueType, numVars>& x)
Evaluation<ValueType, numVars> result;
result.value = ValueTypeToolbox::atan(x.value);
result.setValue(ValueTypeToolbox::atan(x.value()));
// derivatives use the chain rule
const ValueType& df_dx = 1/(1 + x.value*x.value);
for (unsigned curVarIdx = 0; curVarIdx < result.derivatives.size(); ++curVarIdx)
result.derivatives[curVarIdx] = df_dx*x.derivatives[curVarIdx];
const ValueType& df_dx = 1/(1 + x.value()*x.value());
for (unsigned curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
return result;
}
@@ -120,15 +120,14 @@ Evaluation<ValueType, numVars> atan2(const Evaluation<ValueType, numVars>& x,
Evaluation<ValueType, numVars> result;
result.value = ValueTypeToolbox::atan2(x.value, y.value);
result.setValue(ValueTypeToolbox::atan2(x.value(), y.value()));
// derivatives use the chain rule
const ValueType& alpha = 1/(1 + (x.value*x.value)/(y.value*y.value));
for (unsigned curVarIdx = 0; curVarIdx < result.derivatives.size(); ++curVarIdx) {
result.derivatives[curVarIdx] =
alpha
/(y.value*y.value)
*(x.derivatives[curVarIdx]*y.value - x.value*y.derivatives[curVarIdx]);
const ValueType& alpha = 1/(1 + (x.value()*x.value())/(y.value()*y.value()));
for (unsigned curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) {
result.setDerivative(curVarIdx,
alpha/(y.value()*y.value())
*(x.derivative(curVarIdx)*y.value() - x.value()*y.derivative(curVarIdx)));
}
return result;
@@ -141,12 +140,12 @@ Evaluation<ValueType, numVars> sin(const Evaluation<ValueType, numVars>& x)
Evaluation<ValueType, numVars> result;
result.value = ValueTypeToolbox::sin(x.value);
result.setValue(ValueTypeToolbox::sin(x.value()));
// derivatives use the chain rule
const ValueType& df_dx = ValueTypeToolbox::cos(x.value);
for (unsigned curVarIdx = 0; curVarIdx < result.derivatives.size(); ++curVarIdx)
result.derivatives[curVarIdx] = df_dx*x.derivatives[curVarIdx];
const ValueType& df_dx = ValueTypeToolbox::cos(x.value());
for (unsigned curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
return result;
}
@@ -158,12 +157,12 @@ Evaluation<ValueType, numVars> asin(const Evaluation<ValueType, numVars>& x)
Evaluation<ValueType, numVars> result;
result.value = ValueTypeToolbox::asin(x.value);
result.setValue(ValueTypeToolbox::asin(x.value()));
// derivatives use the chain rule
const ValueType& df_dx = 1.0/ValueTypeToolbox::sqrt(1 - x.value*x.value);
for (unsigned curVarIdx = 0; curVarIdx < result.derivatives.size(); ++curVarIdx)
result.derivatives[curVarIdx] = df_dx*x.derivatives[curVarIdx];
const ValueType& df_dx = 1.0/ValueTypeToolbox::sqrt(1 - x.value()*x.value());
for (unsigned curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
return result;
}
@@ -175,12 +174,12 @@ Evaluation<ValueType, numVars> cos(const Evaluation<ValueType, numVars>& x)
Evaluation<ValueType, numVars> result;
result.value = ValueTypeToolbox::cos(x.value);
result.setValue(ValueTypeToolbox::cos(x.value()));
// derivatives use the chain rule
const ValueType& df_dx = -ValueTypeToolbox::sin(x.value);
for (unsigned curVarIdx = 0; curVarIdx < result.derivatives.size(); ++curVarIdx)
result.derivatives[curVarIdx] = df_dx*x.derivatives[curVarIdx];
const ValueType& df_dx = -ValueTypeToolbox::sin(x.value());
for (unsigned curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
return result;
}
@@ -192,12 +191,12 @@ Evaluation<ValueType, numVars> acos(const Evaluation<ValueType, numVars>& x)
Evaluation<ValueType, numVars> result;
result.value = ValueTypeToolbox::acos(x.value);
result.setValue(ValueTypeToolbox::acos(x.value()));
// derivatives use the chain rule
const ValueType& df_dx = - 1.0/ValueTypeToolbox::sqrt(1 - x.value*x.value);
for (unsigned curVarIdx = 0; curVarIdx < result.derivatives.size(); ++curVarIdx)
result.derivatives[curVarIdx] = df_dx*x.derivatives[curVarIdx];
const ValueType& df_dx = - 1.0/ValueTypeToolbox::sqrt(1 - x.value()*x.value());
for (unsigned curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
return result;
}
@@ -209,13 +208,13 @@ Evaluation<ValueType, numVars> sqrt(const Evaluation<ValueType, numVars>& x)
Evaluation<ValueType, numVars> result;
const ValueType& sqrt_x = ValueTypeToolbox::sqrt(x.value);
result.value = sqrt_x;
const ValueType& sqrt_x = ValueTypeToolbox::sqrt(x.value());
result.setValue(sqrt_x);
// derivatives use the chain rule
ValueType df_dx = 0.5/sqrt_x;
for (unsigned curVarIdx = 0; curVarIdx < result.derivatives.size(); ++curVarIdx) {
result.derivatives[curVarIdx] = df_dx*x.derivatives[curVarIdx];
for (unsigned curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) {
result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
}
return result;
@@ -227,13 +226,13 @@ Evaluation<ValueType, numVars> exp(const Evaluation<ValueType, numVars>& x)
typedef MathToolbox<ValueType> ValueTypeToolbox;
Evaluation<ValueType, numVars> result;
const ValueType& exp_x = ValueTypeToolbox::exp(x.value);
result.value = exp_x;
const ValueType& exp_x = ValueTypeToolbox::exp(x.value());
result.setValue(exp_x);
// derivatives use the chain rule
const ValueType& df_dx = exp_x;
for (unsigned curVarIdx = 0; curVarIdx < result.derivatives.size(); ++curVarIdx)
result.derivatives[curVarIdx] = df_dx*x.derivatives[curVarIdx];
for (unsigned curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
return result;
}
@@ -246,8 +245,8 @@ Evaluation<ValueType, numVars> pow(const Evaluation<ValueType, numVars>& base,
typedef MathToolbox<ValueType> ValueTypeToolbox;
Evaluation<ValueType, numVars> result;
const ValueType& pow_x = ValueTypeToolbox::pow(base.value, exp);
result.value = pow_x;
const ValueType& pow_x = ValueTypeToolbox::pow(base.value(), exp);
result.setValue(pow_x);
if (base == 0.0) {
// we special case the base 0 case because 0.0 is in the valid range of the
@@ -256,9 +255,9 @@ Evaluation<ValueType, numVars> pow(const Evaluation<ValueType, numVars>& base,
}
else {
// derivatives use the chain rule
const ValueType& df_dx = pow_x/base.value*exp;
for (unsigned curVarIdx = 0; curVarIdx < result.derivatives.size(); ++curVarIdx)
result.derivatives[curVarIdx] = df_dx*base.derivatives[curVarIdx];
const ValueType& df_dx = pow_x/base.value()*exp;
for (unsigned curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
result.setDerivative(curVarIdx, df_dx*base.derivative(curVarIdx));
}
return result;
@@ -280,12 +279,12 @@ Evaluation<ValueType, numVars> pow(const BaseType& base,
}
else {
const ValueType& lnBase = ValueTypeToolbox::log(base);
result.value = ValueTypeToolbox::exp(lnBase*exp.value);
result.setValue(ValueTypeToolbox::exp(lnBase*exp.value()));
// derivatives use the chain rule
const ValueType& df_dx = lnBase*result.value;
for (unsigned curVarIdx = 0; curVarIdx < result.derivatives.size(); ++curVarIdx)
result.derivatives[curVarIdx] = df_dx*exp.derivatives[curVarIdx];
const ValueType& df_dx = lnBase*result.value();
for (unsigned curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
result.setDerivative(curVarIdx, df_dx*exp.derivative(curVarIdx));
}
return result;
@@ -307,18 +306,18 @@ Evaluation<ValueType, numVars> pow(const Evaluation<ValueType, numVars>& base,
result = 0.0;
}
else {
ValueType valuePow = ValueTypeToolbox::pow(base.value, exp.value);
result.value = valuePow;
ValueType valuePow = ValueTypeToolbox::pow(base.value(), exp.value());
result.setValue(valuePow);
// use the chain rule for the derivatives. since both, the base and the exponent can
// potentially depend on the variable set, calculating these is quite elaborate...
const ValueType& f = base.value;
const ValueType& g = exp.value;
const ValueType& f = base.value();
const ValueType& g = exp.value();
const ValueType& logF = ValueTypeToolbox::log(f);
for (unsigned curVarIdx = 0; curVarIdx < result.derivatives.size(); ++curVarIdx) {
const ValueType& fPrime = base.derivatives[curVarIdx];
const ValueType& gPrime = exp.derivatives[curVarIdx];
result.derivatives[curVarIdx] = (g*fPrime/f + logF*gPrime) * valuePow;
for (unsigned curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) {
const ValueType& fPrime = base.derivative(curVarIdx);
const ValueType& gPrime = exp.derivative(curVarIdx);
result.setDerivative(curVarIdx, (g*fPrime/f + logF*gPrime) * valuePow);
}
}
@@ -332,12 +331,12 @@ Evaluation<ValueType, numVars> log(const Evaluation<ValueType, numVars>& x)
Evaluation<ValueType, numVars> result;
result.value = ValueTypeToolbox::log(x.value);
result.setValue(ValueTypeToolbox::log(x.value()));
// derivatives use the chain rule
const ValueType& df_dx = 1/x.value;
for (unsigned curVarIdx = 0; curVarIdx < result.derivatives.size(); ++curVarIdx)
result.derivatives[curVarIdx] = df_dx*x.derivatives[curVarIdx];
const ValueType& df_dx = 1/x.value();
for (unsigned curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
return result;
}
@@ -357,10 +356,10 @@ public:
typedef Opm::DenseAd::Evaluation<ValueType, numVars> Evaluation;
static ValueType value(const Evaluation& eval)
{ return eval.value; }
{ return eval.value(); }
static decltype(InnerToolbox::scalarValue(0.0)) scalarValue(const Evaluation& eval)
{ return InnerToolbox::scalarValue(eval.value); }
{ return InnerToolbox::scalarValue(eval.value()); }
static Evaluation createConstant(ValueType value)
{ return Evaluation::createConstant(value); }
@@ -384,7 +383,7 @@ public:
static typename std::enable_if<std::is_floating_point<LhsEval>::value,
LhsEval>::type
decay(const Evaluation& eval)
{ return eval.value; }
{ return eval.value(); }
// comparison
static bool isSame(const Evaluation& a, const Evaluation& b, Scalar tolerance)
@@ -392,12 +391,12 @@ public:
typedef MathToolbox<ValueType> ValueTypeToolbox;
// make sure that the value of the evaluation is identical
if (!ValueTypeToolbox::isSame(a.value, b.value, tolerance))
if (!ValueTypeToolbox::isSame(a.value(), b.value(), tolerance))
return false;
// make sure that the derivatives are identical
for (unsigned curVarIdx = 0; curVarIdx < numVars; ++curVarIdx)
if (!ValueTypeToolbox::isSame(a.derivatives[curVarIdx], b.derivatives[curVarIdx], tolerance))
if (!ValueTypeToolbox::isSame(a.derivative(curVarIdx), b.derivative(curVarIdx), tolerance))
return false;
return true;

100
tests/test_densead.cpp
View File

@@ -97,57 +97,57 @@ void testOperators(const Scalar tolerance)
// test the non-inplace operators
{
Eval a = xEval + yEval;
if (std::abs(a.value - (x + y)) > tolerance)
if (std::abs(a.value() - (x + y)) > tolerance)
throw std::logic_error("oops: operator+");
Eval b = xEval + c;
if (std::abs(b.value - (x + c)) > tolerance)
if (std::abs(b.value() - (x + c)) > tolerance)
throw std::logic_error("oops: operator+");
Eval d = xEval + cEval;
if (std::abs(d.value - (x + c)) > tolerance)
if (std::abs(d.value() - (x + c)) > tolerance)
throw std::logic_error("oops: operator+");
}
{
Eval a = xEval - yEval;
if (std::abs(a.value - (x - y)) > tolerance)
if (std::abs(a.value() - (x - y)) > tolerance)
throw std::logic_error("oops: operator-");
Eval b = xEval - c;
if (std::abs(b.value - (x - c)) > tolerance)
if (std::abs(b.value() - (x - c)) > tolerance)
throw std::logic_error("oops: operator-");
Eval d = xEval - cEval;
if (std::abs(d.value - (x - c)) > tolerance)
if (std::abs(d.value() - (x - c)) > tolerance)
throw std::logic_error("oops: operator-");
}
{
Eval a = xEval*yEval;
if (std::abs(a.value - (x*y)) > tolerance)
if (std::abs(a.value() - (x*y)) > tolerance)
throw std::logic_error("oops: operator*");
Eval b = xEval*c;
if (std::abs(b.value - (x*c)) > tolerance)
if (std::abs(b.value() - (x*c)) > tolerance)
throw std::logic_error("oops: operator*");
Eval d = xEval*cEval;
if (std::abs(d.value - (x*c)) > tolerance)
if (std::abs(d.value() - (x*c)) > tolerance)
throw std::logic_error("oops: operator*");
}
{
Eval a = xEval/yEval;
if (std::abs(a.value - (x/y)) > tolerance)
if (std::abs(a.value() - (x/y)) > tolerance)
throw std::logic_error("oops: operator/");
Eval b = xEval/c;
if (std::abs(b.value - (x/c)) > tolerance)
if (std::abs(b.value() - (x/c)) > tolerance)
throw std::logic_error("oops: operator/");
Eval d = xEval/cEval;
if (std::abs(d.value - (x/c)) > tolerance)
if (std::abs(d.value() - (x/c)) > tolerance)
throw std::logic_error("oops: operator/");
}
@@ -155,68 +155,68 @@ void testOperators(const Scalar tolerance)
{
Eval a = xEval;
a += yEval;
if (std::abs(a.value - (x + y)) > tolerance)
if (std::abs(a.value() - (x + y)) > tolerance)
throw std::logic_error("oops: operator+");
Eval b = xEval;
b += c;
if (std::abs(b.value - (x + c)) > tolerance)
if (std::abs(b.value() - (x + c)) > tolerance)
throw std::logic_error("oops: operator+");
Eval d = xEval;
d += cEval;
if (std::abs(d.value - (x + c)) > tolerance)
if (std::abs(d.value() - (x + c)) > tolerance)
throw std::logic_error("oops: operator+");
}
{
Eval a = xEval;
a -= yEval;
if (std::abs(a.value - (x - y)) > tolerance)
if (std::abs(a.value() - (x - y)) > tolerance)
throw std::logic_error("oops: operator-");
Eval b = xEval;
b -= c;
if (std::abs(b.value - (x - c)) > tolerance)
if (std::abs(b.value() - (x - c)) > tolerance)
throw std::logic_error("oops: operator-");
Eval d = xEval;
d -= cEval;
if (std::abs(d.value - (x - c)) > tolerance)
if (std::abs(d.value() - (x - c)) > tolerance)
throw std::logic_error("oops: operator-");
}
{
Eval a = xEval;
a *= yEval;
if (std::abs(a.value - (x*y)) > tolerance)
if (std::abs(a.value() - (x*y)) > tolerance)
throw std::logic_error("oops: operator*");
Eval b = xEval;
b *= c;
if (std::abs(b.value - (x*c)) > tolerance)
if (std::abs(b.value() - (x*c)) > tolerance)
throw std::logic_error("oops: operator*");
Eval d = xEval;
d *= cEval;
if (std::abs(d.value - (x*c)) > tolerance)
if (std::abs(d.value() - (x*c)) > tolerance)
throw std::logic_error("oops: operator*");
}
{
Eval a = xEval;
a /= yEval;
if (std::abs(a.value - (x/y)) > tolerance)
if (std::abs(a.value() - (x/y)) > tolerance)
throw std::logic_error("oops: operator/");
Eval b = xEval;
b /= c;
if (std::abs(b.value - (x/c)) > tolerance)
if (std::abs(b.value() - (x/c)) > tolerance)
throw std::logic_error("oops: operator/");
Eval d = xEval;
d /= cEval;
if (std::abs(d.value - (x/c)) > tolerance)
if (std::abs(d.value() - (x/c)) > tolerance)
throw std::logic_error("oops: operator/");
}
@@ -276,16 +276,16 @@ void test1DFunction(AdFn* adFn, ClassicFn* classicFn, Scalar xMin = 1e-6, Scalar
Scalar yStar2 = classicFn(x + eps);
Scalar yPrime = (yStar2 - yStar1)/(2*eps);
if (std::abs(y-yEval.value) > 5e-14)
if (std::abs(y-yEval.value()) > 5e-14)
throw std::logic_error("oops: value");
Scalar deltaAbs = std::abs(yPrime - yEval.derivatives[0]);
Scalar deltaAbs = std::abs(yPrime - yEval.derivative(0));
Scalar deltaRel = std::abs(deltaAbs/yPrime);
if (deltaAbs > 1000*eps && deltaRel > 1000*eps)
throw std::logic_error("oops: derivative @"+std::to_string((long double) x)+": "
+ std::to_string((long double) yPrime) + " vs "
+ std::to_string((long double) yEval.derivatives[0])
+ " delta: " + std::to_string((long double) std::abs(yPrime - yEval.derivatives[0])));
+ std::to_string((long double) yEval.derivative(0))
+ " delta: " + std::to_string((long double) std::abs(yPrime - yEval.derivative(0))));
}
}
@@ -312,17 +312,17 @@ void test2DFunction1(AdFn* adFn, ClassicFn* classicFn, Scalar xMin, Scalar xMax,
Scalar zStar2 = classicFn(x + eps, y);
Scalar zPrime = (zStar2 - zStar1)/(2.*eps);
if (std::abs(z - zEval.value)/std::abs(z + zEval.value)
if (std::abs(z - zEval.value())/std::abs(z + zEval.value())
> std::numeric_limits<Scalar>::epsilon()*1e2)
throw std::logic_error("oops: value");
Scalar deltaAbs = std::abs(zPrime - zEval.derivatives[0]);
Scalar deltaAbs = std::abs(zPrime - zEval.derivative(0));
Scalar deltaRel = std::abs(deltaAbs/zPrime);
if (deltaAbs > 1000*eps && deltaRel > 1000*eps)
throw std::logic_error("oops: derivative @"+std::to_string((long double) x)+": "
+ std::to_string((long double) zPrime) + " vs "
+ std::to_string((long double) zEval.derivatives[0])
+ " delta: " + std::to_string((long double) std::abs(zPrime - zEval.derivatives[0])));
+ std::to_string((long double) zEval.derivative(0))
+ " delta: " + std::to_string((long double) std::abs(zPrime - zEval.derivative(0))));
}
}
@@ -349,17 +349,17 @@ void test2DFunction2(AdFn* adFn, ClassicFn* classicFn, Scalar x, Scalar yMin, Sc
Scalar zStar2 = classicFn(x, y + eps);
Scalar zPrime = (zStar2 - zStar1)/(2*eps);
if (std::abs(z - zEval.value)/std::abs(z + zEval.value)
if (std::abs(z - zEval.value())/std::abs(z + zEval.value())
> std::numeric_limits<Scalar>::epsilon()*1e2)
throw std::logic_error("oops: value");
Scalar deltaAbs = std::abs(zPrime - zEval.derivatives[1]);
Scalar deltaAbs = std::abs(zPrime - zEval.derivative(1));
Scalar deltaRel = std::abs(deltaAbs/zPrime);
if (deltaAbs > 1000*eps && deltaRel > 1000*eps)
throw std::logic_error("oops: derivative @"+std::to_string((long double) x)+": "
+ std::to_string((long double) zPrime) + " vs "
+ std::to_string((long double) zEval.derivatives[1])
+ " delta: " + std::to_string((long double) std::abs(zPrime - zEval.derivatives[1])));
+ std::to_string((long double) zEval.derivative(1))
+ " delta: " + std::to_string((long double) std::abs(zPrime - zEval.derivative(1))));
}
}
@@ -388,19 +388,19 @@ void testPowBase(Scalar baseMin = 1e-2, Scalar baseMax = 100)
Scalar zStar2 = pow(base + eps, exp);
Scalar zPrime = (zStar2 - zStar1)/(2*eps);
if (std::abs(z - zEval2.value)/std::abs(z + zEval2.value)
if (std::abs(z - zEval2.value())/std::abs(z + zEval2.value())
> std::numeric_limits<Scalar>::epsilon()*1e2)
throw std::logic_error("oops: value");
Scalar deltaAbs = std::abs(zPrime - zEval1.derivatives[0]);
Scalar deltaAbs = std::abs(zPrime - zEval1.derivative(0));
Scalar deltaRel = std::abs(deltaAbs/zPrime);
if (deltaAbs > 1000*eps && deltaRel > 1000*eps)
throw std::logic_error("oops: derivative @"+std::to_string((long double) base)+": "
+ std::to_string((long double) zPrime) + " vs "
+ std::to_string((long double) zEval1.derivatives[0])
+ " delta: " + std::to_string((long double) std::abs(zPrime - zEval1.derivatives[0])));
+ std::to_string((long double) zEval1.derivative(0))
+ " delta: " + std::to_string((long double) std::abs(zPrime - zEval1.derivative(0))));
if (!EvalToolbox::isSame(zEval1, zEval2, /*tolerance=*/std::numeric_limits<Scalar>::epsilon()*1e3*zEval1.value))
if (!EvalToolbox::isSame(zEval1, zEval2, /*tolerance=*/std::numeric_limits<Scalar>::epsilon()*1e3*zEval1.value()))
throw std::logic_error("oops: pow(Eval, Scalar) != pow(Eval, Eval)");
}
}
@@ -430,19 +430,19 @@ void testPowExp(Scalar expMin = -100, Scalar expMax = 100)
Scalar zStar2 = pow(base, exp + eps);
Scalar zPrime = (zStar2 - zStar1)/(2*eps);
if (std::abs(z - zEval2.value)/std::abs(z + zEval2.value)
if (std::abs(z - zEval2.value())/std::abs(z + zEval2.value())
> std::numeric_limits<Scalar>::epsilon()*1e2)
throw std::logic_error("oops: value");
Scalar deltaAbs = std::abs(zPrime - zEval1.derivatives[1]);
Scalar deltaAbs = std::abs(zPrime - zEval1.derivative(1));
Scalar deltaRel = std::abs(deltaAbs/zPrime);
if (deltaAbs > 1000*eps && deltaRel > 1000*eps)
throw std::logic_error("oops: derivative @"+std::to_string((long double) base)+": "
+ std::to_string((long double) zPrime) + " vs "
+ std::to_string((long double) zEval1.derivatives[1])
+ " delta: " + std::to_string((long double) std::abs(zPrime - zEval1.derivatives[1])));
+ std::to_string((long double) zEval1.derivative(1))
+ " delta: " + std::to_string((long double) std::abs(zPrime - zEval1.derivative(1))));
if (!EvalToolbox::isSame(zEval1, zEval2, /*tolerance=*/std::numeric_limits<Scalar>::epsilon()*1e3*zEval1.value))
if (!EvalToolbox::isSame(zEval1, zEval2, /*tolerance=*/std::numeric_limits<Scalar>::epsilon()*1e3*zEval1.value()))
throw std::logic_error("oops: pow(Eval, Scalar) != pow(Eval, Eval)");
}
}
@@ -480,17 +480,17 @@ void testAtan2()
Scalar zStar2 = atan2(x + eps, y + eps);
Scalar zPrime = (zStar2 - zStar1)/(2*eps);
if (std::abs(z - zEval.value)/std::abs(z + zEval.value)
if (std::abs(z - zEval.value())/std::abs(z + zEval.value())
> std::numeric_limits<Scalar>::epsilon()*1e2)
throw std::logic_error("oops: value");
Scalar deltaAbs = std::abs(zPrime - zEval.derivatives[0]);
Scalar deltaAbs = std::abs(zPrime - zEval.derivative(0));
Scalar deltaRel = std::abs(deltaAbs/zPrime);
if (deltaAbs > 1000*eps && deltaRel > 1000*eps)
throw std::logic_error("oops: derivative @("+std::to_string((long double) x)+","+std::to_string((long double) y)+"): "
+ std::to_string((long double) zPrime) + " vs "
+ std::to_string((long double) zEval.derivatives[0])
+ " delta: " + std::to_string((long double) std::abs(zPrime - zEval.derivatives[0])));
+ std::to_string((long double) zEval.derivative(0))
+ " delta: " + std::to_string((long double) std::abs(zPrime - zEval.derivative(0))));
}
}