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

Evaluation: specialize evaluations for used number in Blackoil setting

Browse Source 
to increase performance. Also, change from unsigned to int because it's
better supported by compilers.
This commit is contained in:
Robert Kloefkorn
2017-03-10 15:49:04 +01:00
parent 2a9128a515
commit db46b907a2
8 changed files with 2402 additions and 150 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -26,6 +26,7 @@
* \brief Representation of an evaluation of a function and its derivatives w.r.t. a set
* of variables in the localized OPM automatic differentiation (AD) framework.
*/
#ifndef OPM_LOCAL_AD_EVALUATION_HPP
#define OPM_LOCAL_AD_EVALUATION_HPP
@@ -47,7 +48,7 @@ namespace DenseAd {
* \brief Represents a function evaluation and its derivatives w.r.t. a fixed set of
* variables.
*/
template <class ValueT, unsigned numVars>
template <class ValueT, int numVars>
class Evaluation
{
public:
@@ -55,18 +56,18 @@ public:
typedef ValueT ValueType;
//! number of derivatives
static constexpr unsigned size = numVars;
static constexpr int size = numVars;
protected:
//! length of internal data vector
static constexpr unsigned length_ = numVars + 1 ;
static constexpr int length_ = numVars + 1 ;
//! position index for value
static constexpr unsigned valuepos_ = 0;
static constexpr int valuepos_ = 0;
//! start index for derivatives
static constexpr unsigned dstart_ = 1;
static constexpr int dstart_ = 1;
//! end+1 index for derivatives
static constexpr unsigned dend_ = length_ ;
static constexpr int dend_ = length_ ;
public:
//! default constructor
@@ -95,7 +96,7 @@ public:
// i.e., f(x) = c. this implies an evaluation with the given value and all
// derivatives being zero.
template <class RhsValueType>
Evaluation(const RhsValueType& c, unsigned varPos)
Evaluation(const RhsValueType& c, int varPos)
{
setValue( c );
clearDerivatives();
@@ -109,19 +110,33 @@ public:
// set all derivatives to zero
void clearDerivatives()
{
for (unsigned i = dstart_; i < dend_; ++i)
for (int i = dstart_; i < dend_; ++i)
data_[ i ] = 0.0;
}
// create a function evaluation for a "naked" depending variable (i.e., f(x) = x)
template <class RhsValueType>
static Evaluation createVariable(const RhsValueType& value, unsigned varPos)
static Evaluation createVariable(const RhsValueType& value, int varPos)
{
// 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)
return Evaluation( value, varPos );
}
// create a function evaluation for a "naked" depending variable (i.e., f(x) = x)
template <class RhsValueType>
static Evaluation devide(const RhsValueType& a, const Evaluation& b )
{
Evaluation<ValueType, numVars> result;
result.setValue( a/b.value() );
const ValueType df_dg = - result.value()/b.value();
for (int idx = dstart_; idx < dend_; ++idx) {
result.data_[idx] = df_dg*b.data_[idx];
}
return result;
}
// "evaluate" a constant function (i.e. a function that does not depend on the set of
// relevant variables, f(x) = c).
template <class RhsValueType>
@@ -136,14 +151,14 @@ public:
// print value
os << "v: " << value() << " / d:";
// print derivatives
for (unsigned varIdx = 0; varIdx < numVars; ++varIdx)
for (int varIdx = 0; varIdx < numVars; ++varIdx)
os << " " << derivative(varIdx);
}
// copy all derivatives from other
void copyDerivatives(const Evaluation& other)
{
for (unsigned varIdx = dstart_; varIdx < dend_; ++varIdx)
for (int varIdx = dstart_; varIdx < dend_; ++varIdx)
data_[ varIdx ] = other.data_[ varIdx ];
}
@@ -152,7 +167,7 @@ public:
Evaluation& operator+=(const Evaluation& other)
{
// value and derivatives are added
for (unsigned varIdx = 0; varIdx < length_; ++ varIdx)
for (int varIdx = 0; varIdx < length_; ++ varIdx)
data_[ varIdx ] += other.data_[ varIdx ];
return *this;
@@ -171,7 +186,7 @@ public:
Evaluation& operator-=(const Evaluation& other)
{
// value and derivatives are subtracted
for (unsigned idx = 0 ; idx < length_ ; ++ idx)
for (int idx = 0 ; idx < length_ ; ++ idx)
data_[idx] -= other.data_[idx];
return *this;
@@ -192,25 +207,23 @@ public:
{
// while the values are multiplied, the derivatives follow the product rule,
// i.e., (u*v)' = (v'u + u'v).
ValueType& u = data_[ valuepos_ ];
const ValueType& v = other.value();
for (unsigned idx = dstart_; idx < dend_; ++idx) {
const ValueType& uPrime = data_[idx];
const ValueType& vPrime = other.data_[idx];
const ValueType u = value();
const ValueType v = other.value();
data_[idx] = (v*uPrime + u*vPrime);
data_[ valuepos_ ] *= v ;
for (int idx = dstart_; idx < dend_; ++idx) {
data_[idx] = data_[idx] * v + other.data_[idx] * u;
}
u *= v;
return *this;
}
// m(u*v)' = (v'u + u'v)
template <class RhsValueType>
Evaluation& operator*=(RhsValueType other)
Evaluation& operator*=(const RhsValueType& other)
{
// values and derivatives are multiplied
for (unsigned idx = 0 ; idx < length_ ; ++ idx)
for (int idx = 0 ; idx < length_ ; ++ idx)
data_[idx] *= other;
return *this;
@@ -219,17 +232,14 @@ public:
// m(u*v)' = (v'u + u'v)
Evaluation& operator/=(const Evaluation& other)
{
// values are divided, derivatives follow the rule for division, i.e., (u/v)' = (v'u -
// u'v)/v^2.
ValueType& u = data_[ valuepos_ ];
const ValueType& v = other.value();
for (unsigned idx = dstart_; idx < dend_; ++idx) {
const ValueType& uPrime = data_[idx];
const ValueType& vPrime = other.data_[idx];
// values are divided, derivatives follow the rule for division, i.e., (u/v)' = (v'u - u'v)/v^2.
const ValueType v_vv = 1.0 / other.value();
const ValueType u_vv = value() * v_vv * v_vv;
data_[idx] = (v*uPrime - u*vPrime)/(v*v);
data_[ valuepos_ ] *= v_vv;
for (int idx = dstart_; idx < dend_; ++idx) {
data_[idx] = data_[idx] * v_vv - other.data_[idx] * u_vv;
}
u /= v;
return *this;
}
@@ -239,11 +249,8 @@ public:
Evaluation& operator/=(const RhsValueType& other)
{
// values and derivatives are divided
ValueType factor = (1.0/other);
for (unsigned idx = 0; idx < length_; ++idx)
data_[idx] *= factor;
return *this;
const ValueType factor = (1.0/other);
return this->operator *=( factor );
}
// add two evaluation objects
@@ -285,7 +292,7 @@ public:
{
Evaluation result;
// set value and derivatives to negative
for (unsigned idx = 0; idx < length_; ++idx)
for (int idx = 0; idx < length_; ++idx)
result.data_[idx] = - data_[idx];
return result;
@@ -342,7 +349,7 @@ public:
bool operator==(const Evaluation& other) const
{
for (unsigned idx = 0; idx < length_; ++idx)
for (int idx = 0; idx < length_; ++idx)
if (data_[idx] != other.data_[idx])
return false;
@@ -389,14 +396,14 @@ public:
{ data_[valuepos_] = val; }
// return varIdx'th derivative
const ValueType& derivative(unsigned varIdx) const
const ValueType& derivative(int varIdx) const
{
assert(varIdx < numVars);
return data_[varIdx + dstart_];
}
// set derivative at position varIdx
void setDerivative(unsigned varIdx, const ValueType& derVal)
void setDerivative(int varIdx, const ValueType& derVal)
{
assert(varIdx < numVars);
data_[varIdx + dstart_] = derVal;
@@ -406,27 +413,27 @@ protected:
std::array<ValueType, length_> data_;
};
template <class RhsValueType, class ValueType, unsigned numVars>
template <class RhsValueType, class ValueType, int numVars>
bool operator<(const RhsValueType& a, const Evaluation<ValueType, numVars>& b)
{ return b > a; }
template <class RhsValueType, class ValueType, unsigned numVars>
template <class RhsValueType, class ValueType, int numVars>
bool operator>(const RhsValueType& a, const Evaluation<ValueType, numVars>& b)
{ return b < a; }
template <class RhsValueType, class ValueType, unsigned numVars>
template <class RhsValueType, class ValueType, int numVars>
bool operator<=(const RhsValueType& a, const Evaluation<ValueType, numVars>& b)
{ return b >= a; }
template <class RhsValueType, class ValueType, unsigned numVars>
template <class RhsValueType, class ValueType, int numVars>
bool operator>=(const RhsValueType& a, const Evaluation<ValueType, numVars>& b)
{ return b <= a; }
template <class RhsValueType, class ValueType, unsigned numVars>
template <class RhsValueType, class ValueType, int numVars>
bool operator!=(const RhsValueType& a, const Evaluation<ValueType, numVars>& b)
{ return a != b.value(); }
template <class RhsValueType, class ValueType, unsigned numVars>
template <class RhsValueType, class ValueType, int numVars>
Evaluation<ValueType, numVars> operator+(const RhsValueType& a, const Evaluation<ValueType, numVars>& b)
{
Evaluation<ValueType, numVars> result(b);
@@ -436,46 +443,29 @@ Evaluation<ValueType, numVars> operator+(const RhsValueType& a, const Evaluation
return result;
}
template <class RhsValueType, class ValueType, unsigned numVars>
template <class RhsValueType, class ValueType, int numVars>
Evaluation<ValueType, numVars> operator-(const RhsValueType& a, const Evaluation<ValueType, numVars>& b)
{
Evaluation<ValueType, numVars> result;
result.setValue(a - b.value());
for (unsigned varIdx = 0; varIdx < numVars; ++varIdx)
result.setDerivative(varIdx, - b.derivative(varIdx));
Evaluation<ValueType, numVars> result( a );
result -= b;
return result;
}
template <class RhsValueType, class ValueType, unsigned numVars>
template <class RhsValueType, class ValueType, int numVars>
Evaluation<ValueType, numVars> operator/(const RhsValueType& a, const Evaluation<ValueType, numVars>& b)
{
Evaluation<ValueType, numVars> result;
result.setValue(a/b.value());
// outer derivative
const ValueType& df_dg = - a/(b.value()*b.value());
for (unsigned varIdx = 0; varIdx < numVars; ++varIdx)
result.setDerivative(varIdx, df_dg*b.derivative(varIdx));
return result;
return Evaluation<ValueType, numVars>::devide( a, b );
}
template <class RhsValueType, class ValueType, unsigned numVars>
template <class RhsValueType, class ValueType, int numVars>
Evaluation<ValueType, numVars> operator*(const RhsValueType& a, const Evaluation<ValueType, numVars>& b)
{
Evaluation<ValueType, numVars> result;
result.setValue(a*b.value());
for (unsigned varIdx = 0; varIdx < numVars; ++varIdx)
result.setDerivative(varIdx, a*b.derivative(varIdx));
Evaluation<ValueType, numVars> result( b );
result *= a;
return result;
}
template <class ValueType, unsigned numVars>
template <class ValueType, int numVars>
std::ostream& operator<<(std::ostream& os, const Evaluation<ValueType, numVars>& eval)
{
os << eval.value();
@@ -511,12 +501,12 @@ std::ostream& operator<<(std::ostream& os, const Evaluation<ValueType, numVars>&
namespace Opm {
namespace DenseAd {
template <class ValueType, unsigned numVars>
template <class ValueType, int numVars>
Evaluation<ValueType, numVars> abs(const Evaluation<ValueType, numVars>&);
}}
namespace std {
template <class ValueType, unsigned numVars>
template <class ValueType, int numVars>
const Opm::DenseAd::Evaluation<ValueType, numVars> abs(const Opm::DenseAd::Evaluation<ValueType, numVars>& x)
{ return Opm::DenseAd::abs(x); }
@@ -535,7 +525,7 @@ const Opm::DenseAd::Evaluation<ValueType, numVars> abs(const Opm::DenseAd::Evalu
#include <dune/common/ftraits.hh>
namespace Dune {
template <class ValueType, unsigned numVars>
template <class ValueType, int numVars>
struct FieldTraits<Opm::DenseAd::Evaluation<ValueType, numVars> >
{
public:
@@ -547,4 +537,10 @@ public:
} // namespace Dune
#include <opm/material/densead/Evaluation1.hpp>
#include <opm/material/densead/Evaluation2.hpp>
#include <opm/material/densead/Evaluation3.hpp>
#include <opm/material/densead/Evaluation6.hpp>
#include <opm/material/densead/Evaluation12.hpp>
#endif

416
opm/material/densead/Evaluation1.hpp Normal file
View File

@@ -0,0 +1,416 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 2 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
*
* \brief Representation of an evaluation of a function and its derivatives w.r.t. a set
* of variables in the localized OPM automatic differentiation (AD) framework.
*/
#ifndef OPM_LOCAL_AD_EVALUATION_1_HPP
#define OPM_LOCAL_AD_EVALUATION_1_HPP
#include "Math.hpp"
#include <opm/common/Valgrind.hpp>
#include <dune/common/version.hh>
#include <array>
#include <cmath>
#include <cassert>
#include <iostream>
#include <algorithm>
namespace Opm {
namespace DenseAd {
/*!
* \brief Represents a function evaluation and its derivatives w.r.t. a fixed set of
* variables.
*/
template <class ValueT>
class Evaluation< ValueT, 1 >
{
static constexpr int numVars = 1;
public:
//! field type
typedef ValueT ValueType;
//! number of derivatives
static constexpr int size = numVars;
protected:
//! length of internal data vector
static constexpr int length_ = numVars + 1 ;
//! position index for value
static constexpr int valuepos_ = 0;
//! start index for derivatives
static constexpr int dstart_ = 1;
//! end+1 index for derivatives
static constexpr int dend_ = length_ ;
public:
//! default constructor
Evaluation() : data_()
{}
//! copy other function evaluation
Evaluation(const Evaluation& other) : data_( other.data_ )
{
}
// create an evaluation which represents a constant function
//
// i.e., f(x) = c. this implies an evaluation with the given value and all
// derivatives being zero.
template <class RhsValueType>
Evaluation(const RhsValueType& c)
{
setValue( c );
clearDerivatives();
Valgrind::CheckDefined( data_ );
}
// create an evaluation which represents a constant function
//
// i.e., f(x) = c. this implies an evaluation with the given value and all
// derivatives being zero.
template <class RhsValueType>
Evaluation(const RhsValueType& c, int varPos)
{
setValue( c );
clearDerivatives();
// The variable position must be in represented by the given variable descriptor
assert(0 <= varPos && varPos < numVars);
data_[varPos + dstart_] = 1.0;
Valgrind::CheckDefined(data_);
}
// set all derivatives to zero
void clearDerivatives()
{
data_[ 1 ] = 0;
}
// create a function evaluation for a "naked" depending variable (i.e., f(x) = x)
template <class RhsValueType>
static Evaluation devide(const RhsValueType& a, const Evaluation& b )
{
Evaluation<ValueType, numVars> result;
result.setValue( a/b.value() );
const ValueType df_dg = - result.value()/b.value();
for (int idx = dstart_; idx < dend_; ++idx) {
result.data_[idx] = df_dg*b.data_[idx];
}
return result;
}
// create a function evaluation for a "naked" depending variable (i.e., f(x) = x)
template <class RhsValueType>
static Evaluation createVariable(const RhsValueType& value, int varPos)
{
// 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)
return Evaluation( value, varPos );
}
// "evaluate" a constant function (i.e. a function that does not depend on the set of
// relevant variables, f(x) = c).
template <class RhsValueType>
static Evaluation createConstant(const RhsValueType& value)
{
return Evaluation( value );
}
// print the value and the derivatives of the function evaluation
void print(std::ostream& os = std::cout) const
{
// print value
os << "v: " << value() << " / d:";
// print derivatives
for (int varIdx = 0; varIdx < numVars; ++varIdx)
os << " " << derivative(varIdx);
}
// copy all derivatives from other
void copyDerivatives(const Evaluation& other)
{
data_[ 1 ] = other.data_[ 1 ];
}
// add value and derivatives from other to this values and derivatives
Evaluation& operator+=(const Evaluation& other)
{
data_[ 0 ] += other.data_[ 0 ];
data_[ 1 ] += other.data_[ 1 ];
return *this;
}
// add value from other to this values
template <class RhsValueType>
Evaluation& operator+=(const RhsValueType& other)
{
// value is added, derivatives stay the same
data_[valuepos_] += other;
return *this;
}
// subtract other's value and derivatives from this values
Evaluation& operator-=(const Evaluation& other)
{
// value and derivatives are subtracted
data_[ 0 ] -= other.data_[ 0 ];
data_[ 1 ] -= other.data_[ 1 ];
return *this;
}
// subtract other's value from this values
template <class RhsValueType>
Evaluation& operator-=(const RhsValueType& other)
{
// for constants, values are subtracted, derivatives stay the same
data_[ valuepos_ ] -= other;
return *this;
}
// multiply values and apply chain rule to derivatives: (u*v)' = (v'u + u'v)
Evaluation& operator*=(const Evaluation& other)
{
// while the values are multiplied, the derivatives follow the product rule,
// i.e., (u*v)' = (v'u + u'v).
const ValueType u = value();
const ValueType v = other.value();
data_[ 0 ] = u * v ;
data_[ 1 ] = data_[ 1 ] * v + other.data_[ 1 ] * u;
return *this;
}
// m(u*v)' = (v'u + u'v)
template <class RhsValueType>
Evaluation& operator*=(RhsValueType other)
{
data_[ 0 ] *= other;
data_[ 1 ] *= other;
return *this;
}
// m(u*v)' = (v'u + u'v)
Evaluation& operator/=(const Evaluation& other)
{
// values are divided, derivatives follow the rule for division, i.e., (u/v)' = (v'u -
// u'v)/v^2.
const ValueType v_vv = 1.0 / other.value();
const ValueType u_vv = value() * v_vv * v_vv;
data_[ 0 ] *= v_vv;
data_[ 1 ] = data_[ 1 ] * v_vv - other.data_[ 1 ] * u_vv ;
return *this;
}
// multiply value and derivatives by value of other
template <class RhsValueType>
Evaluation& operator/=(const RhsValueType& other)
{
// values and derivatives are divided
ValueType factor = (1.0/other);
data_[ 0 ] *= factor;
data_[ 1 ] *= factor;
return *this;
}
// add two evaluation objects
Evaluation operator+(const Evaluation& other) const
{
Evaluation result(*this);
result += other;
return result;
}
// add constant to this object
template <class RhsValueType>
Evaluation operator+(const RhsValueType& other) const
{
Evaluation result(*this);
result += other;
return result;
}
// subtract two evaluation objects
Evaluation operator-(const Evaluation& other) const
{
Evaluation result(*this);
result -= other;
return result;
}
// subtract constant from evaluation object
template <class RhsValueType>
Evaluation operator-(const RhsValueType& other) const
{
Evaluation result(*this);
result -= other;
return result;
}
// negation (unary minus) operator
Evaluation operator-() const
{
Evaluation result;
// set value and derivatives to negative
for (int idx = 0; idx < length_; ++idx)
result.data_[idx] = - data_[idx];
return result;
}
Evaluation operator*(const Evaluation& other) const
{
Evaluation result(*this);
result *= other;
return result;
}
template <class RhsValueType>
Evaluation operator*(const RhsValueType& other) const
{
Evaluation result(*this);
result *= other;
return result;
}
Evaluation operator/(const Evaluation& other) const
{
Evaluation result(*this);
result /= other;
return result;
}
template <class RhsValueType>
Evaluation operator/(const RhsValueType& other) const
{
Evaluation result(*this);
result /= other;
return result;
}
template <class RhsValueType>
Evaluation& operator=(const RhsValueType& other)
{
setValue( other );
clearDerivatives();
return *this;
}
// copy assignment from evaluation
Evaluation& operator=(const Evaluation& other)
{
data_ = other.data_;
return *this;
}
template <class RhsValueType>
bool operator==(const RhsValueType& other) const
{ return value() == other; }
bool operator==(const Evaluation& other) const
{
for (int idx = 0; idx < length_; ++idx)
if (data_[idx] != other.data_[idx])
return false;
return true;
}
bool operator!=(const Evaluation& other) const
{ return !operator==(other); }
template <class RhsValueType>
bool operator>(RhsValueType other) const
{ return value() > other; }
bool operator>(const Evaluation& other) const
{ return value() > other.value(); }
template <class RhsValueType>
bool operator<(RhsValueType other) const
{ return value() < other; }
bool operator<(const Evaluation& other) const
{ return value() < other.value(); }
template <class RhsValueType>
bool operator>=(RhsValueType other) const
{ return value() >= other; }
bool operator>=(const Evaluation& other) const
{ return value() >= other.value(); }
template <class RhsValueType>
bool operator<=(RhsValueType other) const
{ return value() <= other; }
bool operator<=(const Evaluation& other) const
{ return value() <= other.value(); }
// return value of variable
const ValueType& value() const
{ return data_[valuepos_]; }
// set value of variable
void setValue(const ValueType& val)
{ data_[valuepos_] = val; }
// return varIdx'th derivative
const ValueType& derivative(int varIdx) const
{
assert(varIdx < numVars);
return data_[varIdx + dstart_];
}
// set derivative at position varIdx
void setDerivative(int varIdx, const ValueType& derVal)
{
assert(varIdx < numVars);
data_[varIdx + dstart_] = derVal;
}
protected:
std::array<ValueType, length_> data_;
};
} // namespace DenseAD
} // namespace Dune
#endif

507
opm/material/densead/Evaluation12.hpp Normal file
View File

@@ -0,0 +1,507 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 2 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
*
* \brief Representation of an evaluation of a function and its derivatives w.r.t. a set
* of variables in the localized OPM automatic differentiation (AD) framework.
*/
#ifndef OPM_LOCAL_AD_EVALUATION_12_HPP
#define OPM_LOCAL_AD_EVALUATION_12_HPP
#include "Math.hpp"
#include <opm/common/Valgrind.hpp>
#include <dune/common/version.hh>
#include <array>
#include <cmath>
#include <cassert>
#include <iostream>
#include <algorithm>
namespace Opm {
namespace DenseAd {
/*!
* \brief Represents a function evaluation and its derivatives w.r.t. a fixed set of
* variables.
*/
template <class ValueT>
class Evaluation< ValueT, 12 >
{
static constexpr int numVars = 12;
public:
//! field type
typedef ValueT ValueType;
//! number of derivatives
static constexpr int size = numVars;
protected:
//! length of internal data vector
static constexpr int length_ = numVars + 1 ;
//! position index for value
static constexpr int valuepos_ = 0;
//! start index for derivatives
static constexpr int dstart_ = 1;
//! end+1 index for derivatives
static constexpr int dend_ = length_ ;
public:
//! default constructor
Evaluation() : data_()
{}
//! copy other function evaluation
Evaluation(const Evaluation& other) : data_( other.data_ )
{
}
// create an evaluation which represents a constant function
//
// i.e., f(x) = c. this implies an evaluation with the given value and all
// derivatives being zero.
template <class RhsValueType>
Evaluation(const RhsValueType& c)
{
setValue( c );
clearDerivatives();
Valgrind::CheckDefined( data_ );
}
// create a function evaluation for a "naked" depending variable (i.e., f(x) = x)
template <class RhsValueType>
static Evaluation devide(const RhsValueType& a, const Evaluation& b )
{
Evaluation<ValueType, numVars> result;
result.setValue( a/b.value() );
const ValueType df_dg = - result.value()/b.value();
for (int idx = dstart_; idx < dend_; ++idx) {
result.data_[idx] = df_dg*b.data_[idx];
}
return result;
}
// create an evaluation which represents a constant function
//
// i.e., f(x) = c. this implies an evaluation with the given value and all
// derivatives being zero.
template <class RhsValueType>
Evaluation(const RhsValueType& c, int varPos)
{
setValue( c );
clearDerivatives();
// The variable position must be in represented by the given variable descriptor
assert(0 <= varPos && varPos < numVars);
data_[varPos + dstart_] = 1.0;
Valgrind::CheckDefined(data_);
}
// set all derivatives to zero
void clearDerivatives()
{
data_[ 1 ] = 0;
data_[ 2 ] = 0;
data_[ 3 ] = 0;
data_[ 4 ] = 0;
data_[ 5 ] = 0;
data_[ 6 ] = 0;
data_[ 7 ] = 0;
data_[ 8 ] = 0;
data_[ 9 ] = 0;
data_[ 10 ] = 0;
data_[ 11 ] = 0;
data_[ 12 ] = 0;
}
// create a function evaluation for a "naked" depending variable (i.e., f(x) = x)
template <class RhsValueType>
static Evaluation createVariable(const RhsValueType& value, int varPos)
{
// 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)
return Evaluation( value, varPos );
}
// "evaluate" a constant function (i.e. a function that does not depend on the set of
// relevant variables, f(x) = c).
template <class RhsValueType>
static Evaluation createConstant(const RhsValueType& value)
{
return Evaluation( value );
}
// print the value and the derivatives of the function evaluation
void print(std::ostream& os = std::cout) const
{
// print value
os << "v: " << value() << " / d:";
// print derivatives
for (int varIdx = 0; varIdx < numVars; ++varIdx)
os << " " << derivative(varIdx);
}
// copy all derivatives from other
void copyDerivatives(const Evaluation& other)
{
data_[ 1 ] = other.data_[ 1 ];
data_[ 2 ] = other.data_[ 2 ];
data_[ 3 ] = other.data_[ 3 ];
data_[ 4 ] = other.data_[ 4 ];
data_[ 5 ] = other.data_[ 5 ];
data_[ 6 ] = other.data_[ 6 ];
data_[ 7 ] = other.data_[ 7 ];
data_[ 8 ] = other.data_[ 8 ];
data_[ 9 ] = other.data_[ 9 ];
data_[ 10 ] = other.data_[ 10 ];
data_[ 11 ] = other.data_[ 11 ];
data_[ 12 ] = other.data_[ 12 ];
}
// add value and derivatives from other to this values and derivatives
Evaluation& operator+=(const Evaluation& other)
{
data_[ 0 ] += other.data_[ 0 ];
data_[ 1 ] += other.data_[ 1 ];
data_[ 2 ] += other.data_[ 2 ];
data_[ 3 ] += other.data_[ 3 ];
data_[ 4 ] += other.data_[ 4 ];
data_[ 5 ] += other.data_[ 5 ];
data_[ 6 ] += other.data_[ 6 ];
data_[ 7 ] += other.data_[ 7 ];
data_[ 8 ] += other.data_[ 8 ];
data_[ 9 ] += other.data_[ 9 ];
data_[ 10 ] += other.data_[ 10 ];
data_[ 11 ] += other.data_[ 11 ];
data_[ 12 ] += other.data_[ 12 ];
return *this;
}
// add value from other to this values
template <class RhsValueType>
Evaluation& operator+=(const RhsValueType& other)
{
// value is added, derivatives stay the same
data_[valuepos_] += other;
return *this;
}
// subtract other's value and derivatives from this values
Evaluation& operator-=(const Evaluation& other)
{
// value and derivatives are subtracted
data_[ 0 ] -= other.data_[ 0 ];
data_[ 1 ] -= other.data_[ 1 ];
data_[ 2 ] -= other.data_[ 2 ];
data_[ 3 ] -= other.data_[ 3 ];
data_[ 4 ] -= other.data_[ 4 ];
data_[ 5 ] -= other.data_[ 5 ];
data_[ 6 ] -= other.data_[ 6 ];
data_[ 7 ] -= other.data_[ 7 ];
data_[ 8 ] -= other.data_[ 8 ];
data_[ 9 ] -= other.data_[ 9 ];
data_[ 10 ] -= other.data_[ 10 ];
data_[ 11 ] -= other.data_[ 11 ];
data_[ 12 ] -= other.data_[ 12 ];
return *this;
}
// subtract other's value from this values
template <class RhsValueType>
Evaluation& operator-=(const RhsValueType& other)
{
// for constants, values are subtracted, derivatives stay the same
data_[ valuepos_ ] -= other;
return *this;
}
// multiply values and apply chain rule to derivatives: (u*v)' = (v'u + u'v)
Evaluation& operator*=(const Evaluation& other)
{
// while the values are multiplied, the derivatives follow the product rule,
// i.e., (u*v)' = (v'u + u'v).
const ValueType u = value();
const ValueType v = other.value();
data_[ 0 ] = u * v ;
data_[ 1 ] = data_[ 1 ] * v + other.data_[ 1 ] * u;
data_[ 2 ] = data_[ 2 ] * v + other.data_[ 2 ] * u;
data_[ 3 ] = data_[ 3 ] * v + other.data_[ 3 ] * u;
data_[ 4 ] = data_[ 4 ] * v + other.data_[ 4 ] * u;
data_[ 5 ] = data_[ 5 ] * v + other.data_[ 5 ] * u;
data_[ 6 ] = data_[ 6 ] * v + other.data_[ 6 ] * u;
data_[ 7 ] = data_[ 7 ] * v + other.data_[ 7 ] * u;
data_[ 8 ] = data_[ 8 ] * v + other.data_[ 8 ] * u;
data_[ 9 ] = data_[ 9 ] * v + other.data_[ 9 ] * u;
data_[ 10 ] = data_[ 10 ] * v + other.data_[ 10 ] * u;
data_[ 11 ] = data_[ 11 ] * v + other.data_[ 11 ] * u;
data_[ 12 ] = data_[ 12 ] * v + other.data_[ 12 ] * u;
return *this;
}
// m(u*v)' = (v'u + u'v)
template <class RhsValueType>
Evaluation& operator*=(RhsValueType other)
{
data_[ 0 ] *= other;
data_[ 1 ] *= other;
data_[ 2 ] *= other;
data_[ 3 ] *= other;
data_[ 4 ] *= other;
data_[ 5 ] *= other;
data_[ 6 ] *= other;
data_[ 7 ] *= other;
data_[ 8 ] *= other;
data_[ 9 ] *= other;
data_[ 10 ] *= other;
data_[ 11 ] *= other;
data_[ 12 ] *= other;
return *this;
}
// m(u*v)' = (v'u + u'v)
Evaluation& operator/=(const Evaluation& other)
{
// values are divided, derivatives follow the rule for division, i.e., (u/v)' = (v'u -
// u'v)/v^2.
const ValueType v_vv = 1.0 / other.value();
const ValueType u_vv = value() * v_vv * v_vv;
data_[ 0 ] *= v_vv;
data_[ 1 ] = data_[ 1 ] * v_vv - other.data_[ 1 ] * u_vv ;
data_[ 2 ] = data_[ 2 ] * v_vv - other.data_[ 2 ] * u_vv ;
data_[ 3 ] = data_[ 3 ] * v_vv - other.data_[ 3 ] * u_vv ;
data_[ 4 ] = data_[ 4 ] * v_vv - other.data_[ 4 ] * u_vv ;
data_[ 5 ] = data_[ 5 ] * v_vv - other.data_[ 5 ] * u_vv ;
data_[ 6 ] = data_[ 6 ] * v_vv - other.data_[ 6 ] * u_vv ;
data_[ 7 ] = data_[ 7 ] * v_vv - other.data_[ 7 ] * u_vv ;
data_[ 8 ] = data_[ 8 ] * v_vv - other.data_[ 8 ] * u_vv ;
data_[ 9 ] = data_[ 9 ] * v_vv - other.data_[ 9 ] * u_vv ;
data_[ 10 ] = data_[ 10 ] * v_vv - other.data_[ 10 ] * u_vv ;
data_[ 11 ] = data_[ 11 ] * v_vv - other.data_[ 11 ] * u_vv ;
data_[ 12 ] = data_[ 12 ] * v_vv - other.data_[ 12 ] * u_vv ;
return *this;
}
// multiply value and derivatives by value of other
template <class RhsValueType>
Evaluation& operator/=(const RhsValueType& other)
{
// values and derivatives are divided
ValueType factor = (1.0/other);
data_[ 0 ] *= factor;
data_[ 1 ] *= factor;
data_[ 2 ] *= factor;
data_[ 3 ] *= factor;
data_[ 4 ] *= factor;
data_[ 5 ] *= factor;
data_[ 6 ] *= factor;
data_[ 7 ] *= factor;
data_[ 8 ] *= factor;
data_[ 9 ] *= factor;
data_[ 10 ] *= factor;
data_[ 11 ] *= factor;
data_[ 12 ] *= factor;
return *this;
}
// add two evaluation objects
Evaluation operator+(const Evaluation& other) const
{
Evaluation result(*this);
result += other;
return result;
}
// add constant to this object
template <class RhsValueType>
Evaluation operator+(const RhsValueType& other) const
{
Evaluation result(*this);
result += other;
return result;
}
// subtract two evaluation objects
Evaluation operator-(const Evaluation& other) const
{
Evaluation result(*this);
result -= other;
return result;
}
// subtract constant from evaluation object
template <class RhsValueType>
Evaluation operator-(const RhsValueType& other) const
{
Evaluation result(*this);
result -= other;
return result;
}
// negation (unary minus) operator
Evaluation operator-() const
{
Evaluation result;
// set value and derivatives to negative
for (int idx = 0; idx < length_; ++idx)
result.data_[idx] = - data_[idx];
return result;
}
Evaluation operator*(const Evaluation& other) const
{
Evaluation result(*this);
result *= other;
return result;
}
template <class RhsValueType>
Evaluation operator*(const RhsValueType& other) const
{
Evaluation result(*this);
result *= other;
return result;
}
Evaluation operator/(const Evaluation& other) const
{
Evaluation result(*this);
result /= other;
return result;
}
template <class RhsValueType>
Evaluation operator/(const RhsValueType& other) const
{
Evaluation result(*this);
result /= other;
return result;
}
template <class RhsValueType>
Evaluation& operator=(const RhsValueType& other)
{
setValue( other );
clearDerivatives();
return *this;
}
// copy assignment from evaluation
Evaluation& operator=(const Evaluation& other)
{
data_ = other.data_;
return *this;
}
template <class RhsValueType>
bool operator==(const RhsValueType& other) const
{ return value() == other; }
bool operator==(const Evaluation& other) const
{
for (int idx = 0; idx < length_; ++idx)
if (data_[idx] != other.data_[idx])
return false;
return true;
}
bool operator!=(const Evaluation& other) const
{ return !operator==(other); }
template <class RhsValueType>
bool operator>(RhsValueType other) const
{ return value() > other; }
bool operator>(const Evaluation& other) const
{ return value() > other.value(); }
template <class RhsValueType>
bool operator<(RhsValueType other) const
{ return value() < other; }
bool operator<(const Evaluation& other) const
{ return value() < other.value(); }
template <class RhsValueType>
bool operator>=(RhsValueType other) const
{ return value() >= other; }
bool operator>=(const Evaluation& other) const
{ return value() >= other.value(); }
template <class RhsValueType>
bool operator<=(RhsValueType other) const
{ return value() <= other; }
bool operator<=(const Evaluation& other) const
{ return value() <= other.value(); }
// return value of variable
const ValueType& value() const
{ return data_[valuepos_]; }
// set value of variable
void setValue(const ValueType& val)
{ data_[valuepos_] = val; }
// return varIdx'th derivative
const ValueType& derivative(int varIdx) const
{
assert(varIdx < numVars);
return data_[varIdx + dstart_];
}
// set derivative at position varIdx
void setDerivative(int varIdx, const ValueType& derVal)
{
assert(varIdx < numVars);
data_[varIdx + dstart_] = derVal;
}
protected:
std::array<ValueType, length_> data_;
};
} // namespace DenseAD
} // namespace Dune
// #include <opm/material/densead/EvaluationSIMD.hpp>
#endif

425
opm/material/densead/Evaluation2.hpp Normal file
View File

@@ -0,0 +1,425 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 2 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
*
* \brief Representation of an evaluation of a function and its derivatives w.r.t. a set
* of variables in the localized OPM automatic differentiation (AD) framework.
*/
#ifndef OPM_LOCAL_AD_EVALUATION_2_HPP
#define OPM_LOCAL_AD_EVALUATION_2_HPP
#include "Math.hpp"
#include <opm/common/Valgrind.hpp>
#include <dune/common/version.hh>
#include <array>
#include <cmath>
#include <cassert>
#include <iostream>
#include <algorithm>
namespace Opm {
namespace DenseAd {
/*!
* \brief Represents a function evaluation and its derivatives w.r.t. a fixed set of
* variables.
*/
template <class ValueT>
class Evaluation< ValueT, 2 >
{
static constexpr int numVars = 2;
public:
//! field type
typedef ValueT ValueType;
//! number of derivatives
static constexpr int size = numVars;
protected:
//! length of internal data vector
static constexpr int length_ = numVars + 1 ;
//! position index for value
static constexpr int valuepos_ = 0;
//! start index for derivatives
static constexpr int dstart_ = 1;
//! end+1 index for derivatives
static constexpr int dend_ = length_ ;
public:
//! default constructor
Evaluation() : data_()
{}
//! copy other function evaluation
Evaluation(const Evaluation& other) : data_( other.data_ )
{
}
// create an evaluation which represents a constant function
//
// i.e., f(x) = c. this implies an evaluation with the given value and all
// derivatives being zero.
template <class RhsValueType>
Evaluation(const RhsValueType& c)
{
setValue( c );
clearDerivatives();
Valgrind::CheckDefined( data_ );
}
// create an evaluation which represents a constant function
//
// i.e., f(x) = c. this implies an evaluation with the given value and all
// derivatives being zero.
template <class RhsValueType>
Evaluation(const RhsValueType& c, int varPos)
{
setValue( c );
clearDerivatives();
// The variable position must be in represented by the given variable descriptor
assert(0 <= varPos && varPos < numVars);
data_[varPos + dstart_] = 1.0;
Valgrind::CheckDefined(data_);
}
// set all derivatives to zero
void clearDerivatives()
{
data_[ 1 ] = 0;
data_[ 2 ] = 0;
}
// create a function evaluation for a "naked" depending variable (i.e., f(x) = x)
template <class RhsValueType>
static Evaluation devide(const RhsValueType& a, const Evaluation& b )
{
Evaluation<ValueType, numVars> result;
result.setValue( a/b.value() );
const ValueType df_dg = - result.value()/b.value();
for (int idx = dstart_; idx < dend_; ++idx) {
result.data_[idx] = df_dg*b.data_[idx];
}
return result;
}
// create a function evaluation for a "naked" depending variable (i.e., f(x) = x)
template <class RhsValueType>
static Evaluation createVariable(const RhsValueType& value, int varPos)
{
// 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)
return Evaluation( value, varPos );
}
// "evaluate" a constant function (i.e. a function that does not depend on the set of
// relevant variables, f(x) = c).
template <class RhsValueType>
static Evaluation createConstant(const RhsValueType& value)
{
return Evaluation( value );
}
// print the value and the derivatives of the function evaluation
void print(std::ostream& os = std::cout) const
{
// print value
os << "v: " << value() << " / d:";
// print derivatives
for (int varIdx = 0; varIdx < numVars; ++varIdx)
os << " " << derivative(varIdx);
}
// copy all derivatives from other
void copyDerivatives(const Evaluation& other)
{
data_[ 1 ] = other.data_[ 1 ];
data_[ 2 ] = other.data_[ 2 ];
}
// add value and derivatives from other to this values and derivatives
Evaluation& operator+=(const Evaluation& other)
{
data_[ 0 ] += other.data_[ 0 ];
data_[ 1 ] += other.data_[ 1 ];
data_[ 2 ] += other.data_[ 2 ];
return *this;
}
// add value from other to this values
template <class RhsValueType>
Evaluation& operator+=(const RhsValueType& other)
{
// value is added, derivatives stay the same
data_[valuepos_] += other;
return *this;
}
// subtract other's value and derivatives from this values
Evaluation& operator-=(const Evaluation& other)
{
// value and derivatives are subtracted
data_[ 0 ] -= other.data_[ 0 ];
data_[ 1 ] -= other.data_[ 1 ];
data_[ 2 ] -= other.data_[ 2 ];
return *this;
}
// subtract other's value from this values
template <class RhsValueType>
Evaluation& operator-=(const RhsValueType& other)
{
// for constants, values are subtracted, derivatives stay the same
data_[ valuepos_ ] -= other;
return *this;
}
// multiply values and apply chain rule to derivatives: (u*v)' = (v'u + u'v)
Evaluation& operator*=(const Evaluation& other)
{
// while the values are multiplied, the derivatives follow the product rule,
// i.e., (u*v)' = (v'u + u'v).
const ValueType u = value();
const ValueType v = other.value();
data_[ 0 ] = u * v ;
data_[ 1 ] = data_[ 1 ] * v + other.data_[ 1 ] * u;
data_[ 2 ] = data_[ 2 ] * v + other.data_[ 2 ] * u;
return *this;
}
// m(u*v)' = (v'u + u'v)
template <class RhsValueType>
Evaluation& operator*=(RhsValueType other)
{
data_[ 0 ] *= other;
data_[ 1 ] *= other;
data_[ 2 ] *= other;
return *this;
}
// m(u*v)' = (v'u + u'v)
Evaluation& operator/=(const Evaluation& other)
{
// values are divided, derivatives follow the rule for division, i.e., (u/v)' = (v'u -
// u'v)/v^2.
const ValueType v_vv = 1.0 / other.value();
const ValueType u_vv = value() * v_vv * v_vv;
data_[ 0 ] *= v_vv;
data_[ 1 ] = data_[ 1 ] * v_vv - other.data_[ 1 ] * u_vv ;
data_[ 2 ] = data_[ 2 ] * v_vv - other.data_[ 2 ] * u_vv ;
return *this;
}
// multiply value and derivatives by value of other
template <class RhsValueType>
Evaluation& operator/=(const RhsValueType& other)
{
// values and derivatives are divided
ValueType factor = (1.0/other);
data_[ 0 ] *= factor;
data_[ 1 ] *= factor;
data_[ 2 ] *= factor;
return *this;
}
// add two evaluation objects
Evaluation operator+(const Evaluation& other) const
{
Evaluation result(*this);
result += other;
return result;
}
// add constant to this object
template <class RhsValueType>
Evaluation operator+(const RhsValueType& other) const
{
Evaluation result(*this);
result += other;
return result;
}
// subtract two evaluation objects
Evaluation operator-(const Evaluation& other) const
{
Evaluation result(*this);
result -= other;
return result;
}
// subtract constant from evaluation object
template <class RhsValueType>
Evaluation operator-(const RhsValueType& other) const
{
Evaluation result(*this);
result -= other;
return result;
}
// negation (unary minus) operator
Evaluation operator-() const
{
Evaluation result;
// set value and derivatives to negative
for (int idx = 0; idx < length_; ++idx)
result.data_[idx] = - data_[idx];
return result;
}
Evaluation operator*(const Evaluation& other) const
{
Evaluation result(*this);
result *= other;
return result;
}
template <class RhsValueType>
Evaluation operator*(const RhsValueType& other) const
{
Evaluation result(*this);
result *= other;
return result;
}
Evaluation operator/(const Evaluation& other) const
{
Evaluation result(*this);
result /= other;
return result;
}
template <class RhsValueType>
Evaluation operator/(const RhsValueType& other) const
{
Evaluation result(*this);
result /= other;
return result;
}
template <class RhsValueType>
Evaluation& operator=(const RhsValueType& other)
{
setValue( other );
clearDerivatives();
return *this;
}
// copy assignment from evaluation
Evaluation& operator=(const Evaluation& other)
{
data_ = other.data_;
return *this;
}
template <class RhsValueType>
bool operator==(const RhsValueType& other) const
{ return value() == other; }
bool operator==(const Evaluation& other) const
{
for (int idx = 0; idx < length_; ++idx)
if (data_[idx] != other.data_[idx])
return false;
return true;
}
bool operator!=(const Evaluation& other) const
{ return !operator==(other); }
template <class RhsValueType>
bool operator>(RhsValueType other) const
{ return value() > other; }
bool operator>(const Evaluation& other) const
{ return value() > other.value(); }
template <class RhsValueType>
bool operator<(RhsValueType other) const
{ return value() < other; }
bool operator<(const Evaluation& other) const
{ return value() < other.value(); }
template <class RhsValueType>
bool operator>=(RhsValueType other) const
{ return value() >= other; }
bool operator>=(const Evaluation& other) const
{ return value() >= other.value(); }
template <class RhsValueType>
bool operator<=(RhsValueType other) const
{ return value() <= other; }
bool operator<=(const Evaluation& other) const
{ return value() <= other.value(); }
// return value of variable
const ValueType& value() const
{ return data_[valuepos_]; }
// set value of variable
void setValue(const ValueType& val)
{ data_[valuepos_] = val; }
// return varIdx'th derivative
const ValueType& derivative(int varIdx) const
{
assert(varIdx < numVars);
return data_[varIdx + dstart_];
}
// set derivative at position varIdx
void setDerivative(int varIdx, const ValueType& derVal)
{
assert(varIdx < numVars);
data_[varIdx + dstart_] = derVal;
}
protected:
std::array<ValueType, length_> data_;
};
} // namespace DenseAD
} // namespace Dune
// #include <opm/material/densead/EvaluationSIMD.hpp>
#endif

435
opm/material/densead/Evaluation3.hpp Normal file
View File

@@ -0,0 +1,435 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 2 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
*
* \brief Representation of an evaluation of a function and its derivatives w.r.t. a set
* of variables in the localized OPM automatic differentiation (AD) framework.
*/
#ifndef OPM_LOCAL_AD_EVALUATION_3_HPP
#define OPM_LOCAL_AD_EVALUATION_3_HPP
#include "Math.hpp"
#include <opm/common/Valgrind.hpp>
#include <dune/common/version.hh>
#include <array>
#include <cmath>
#include <cassert>
#include <iostream>
#include <algorithm>
namespace Opm {
namespace DenseAd {
/*!
* \brief Represents a function evaluation and its derivatives w.r.t. a fixed set of
* variables.
*/
template <class ValueT>
class Evaluation< ValueT, 3 >
{
static constexpr int numVars = 3;
public:
//! field type
typedef ValueT ValueType;
//! number of derivatives
static constexpr int size = numVars;
protected:
//! length of internal data vector
static constexpr int length_ = numVars + 1 ;
//! position index for value
static constexpr int valuepos_ = 0;
//! start index for derivatives
static constexpr int dstart_ = 1;
//! end+1 index for derivatives
static constexpr int dend_ = length_ ;
public:
//! default constructor
Evaluation() : data_()
{}
//! copy other function evaluation
Evaluation(const Evaluation& other) //: data_( other.data_ )
{
data_[ 0 ] = other.data_[ 0 ];
data_[ 1 ] = other.data_[ 1 ];
data_[ 2 ] = other.data_[ 2 ];
data_[ 3 ] = other.data_[ 3 ];
}
// create an evaluation which represents a constant function
//
// i.e., f(x) = c. this implies an evaluation with the given value and all
// derivatives being zero.
template <class RhsValueType>
Evaluation(const RhsValueType& c)
{
setValue( c );
clearDerivatives();
Valgrind::CheckDefined( data_ );
}
// create an evaluation which represents a constant function
//
// i.e., f(x) = c. this implies an evaluation with the given value and all
// derivatives being zero.
template <class RhsValueType>
Evaluation(const RhsValueType& c, int varPos)
{
setValue( c );
clearDerivatives();
// The variable position must be in represented by the given variable descriptor
assert(0 <= varPos && varPos < numVars);
data_[varPos + dstart_] = 1.0;
Valgrind::CheckDefined(data_);
}
// set all derivatives to zero
void clearDerivatives()
{
data_[ 1 ] = 0;
data_[ 2 ] = 0;
data_[ 3 ] = 0;
}
// create a function evaluation for a "naked" depending variable (i.e., f(x) = x)
template <class RhsValueType>
static Evaluation devide(const RhsValueType& a, const Evaluation& b )
{
Evaluation<ValueType, numVars> result;
result.setValue( a/b.value() );
const ValueType df_dg = - result.value()/b.value();
result.data_[ 1 ] = df_dg*b.data_[ 1 ];
result.data_[ 2 ] = df_dg*b.data_[ 2 ];
result.data_[ 3 ] = df_dg*b.data_[ 3 ];
return result;
}
// create a function evaluation for a "naked" depending variable (i.e., f(x) = x)
template <class RhsValueType>
static Evaluation createVariable(const RhsValueType& value, int varPos)
{
// 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)
return Evaluation( value, varPos );
}
// "evaluate" a constant function (i.e. a function that does not depend on the set of
// relevant variables, f(x) = c).
template <class RhsValueType>
static Evaluation createConstant(const RhsValueType& value)
{
return Evaluation( value );
}
// print the value and the derivatives of the function evaluation
void print(std::ostream& os = std::cout) const
{
// print value
os << "v: " << value() << " / d:";
// print derivatives
for (int varIdx = 0; varIdx < numVars; ++varIdx)
os << " " << derivative(varIdx);
}
// copy all derivatives from other
void copyDerivatives(const Evaluation& other)
{
data_[ 1 ] = other.data_[ 1 ];
data_[ 2 ] = other.data_[ 2 ];
data_[ 3 ] = other.data_[ 3 ];
}
// add value and derivatives from other to this values and derivatives
Evaluation& operator+=(const Evaluation& other)
{
data_[ 0 ] += other.data_[ 0 ];
data_[ 1 ] += other.data_[ 1 ];
data_[ 2 ] += other.data_[ 2 ];
data_[ 3 ] += other.data_[ 3 ];
return *this;
}
// add value from other to this values
template <class RhsValueType>
Evaluation& operator+=(const RhsValueType& other)
{
// value is added, derivatives stay the same
data_[valuepos_] += other;
return *this;
}
// subtract other's value and derivatives from this values
Evaluation& operator-=(const Evaluation& other)
{
// value and derivatives are subtracted
data_[ 0 ] -= other.data_[ 0 ];
data_[ 1 ] -= other.data_[ 1 ];
data_[ 2 ] -= other.data_[ 2 ];
data_[ 3 ] -= other.data_[ 3 ];
return *this;
}
// subtract other's value from this values
template <class RhsValueType>
Evaluation& operator-=(const RhsValueType& other)
{
// for constants, values are subtracted, derivatives stay the same
data_[ valuepos_ ] -= other;
return *this;
}
// multiply values and apply chain rule to derivatives: (u*v)' = (v'u + u'v)
Evaluation& operator*=(const Evaluation& other)
{
// while the values are multiplied, the derivatives follow the product rule,
// i.e., (u*v)' = (v'u + u'v).
const ValueType u = value();
const ValueType v = other.value();
data_[ 0 ] *= v ;
data_[ 1 ] = data_[ 1 ] * v + other.data_[ 1 ] * u;
data_[ 2 ] = data_[ 2 ] * v + other.data_[ 2 ] * u;
data_[ 3 ] = data_[ 3 ] * v + other.data_[ 3 ] * u;
return *this;
}
// m(u*v)' = (v'u + u'v)
template <class RhsValueType>
Evaluation& operator*=(RhsValueType other)
{
data_[ 0 ] *= other;
data_[ 1 ] *= other;
data_[ 2 ] *= other;
data_[ 3 ] *= other;
return *this;
}
// m(u*v)' = (v'u + u'v)
Evaluation& operator/=(const Evaluation& other)
{
// values are divided, derivatives follow the rule for division, i.e., (u/v)' = (v'u -
// u'v)/v^2.
const ValueType v_vv = 1.0 / other.value();
const ValueType u_vv = value() * v_vv * v_vv;
data_[ 0 ] *= v_vv;
data_[ 1 ] = data_[ 1 ] * v_vv - other.data_[ 1 ] * u_vv ;
data_[ 2 ] = data_[ 2 ] * v_vv - other.data_[ 2 ] * u_vv ;
data_[ 3 ] = data_[ 3 ] * v_vv - other.data_[ 3 ] * u_vv ;
return *this;
}
// multiply value and derivatives by value of other
template <class RhsValueType>
Evaluation& operator/=(const RhsValueType& other)
{
// values and derivatives are divided
const ValueType factor = (1.0/other);
return this->operator*=( factor );
}
// add two evaluation objects
Evaluation operator+(const Evaluation& other) const
{
Evaluation result(*this);
result += other;
return result;
}
// add constant to this object
template <class RhsValueType>
Evaluation operator+(const RhsValueType& other) const
{
Evaluation result(*this);
result += other;
return result;
}
// subtract two evaluation objects
Evaluation operator-(const Evaluation& other) const
{
Evaluation result(*this);
result -= other;
return result;
}
// subtract constant from evaluation object
template <class RhsValueType>
Evaluation operator-(const RhsValueType& other) const
{
Evaluation result(*this);
result -= other;
return result;
}
// negation (unary minus) operator
Evaluation operator-() const
{
Evaluation result;
// set value and derivatives to negative
for (int idx = 0; idx < length_; ++idx)
result.data_[idx] = - data_[idx];
return result;
}
Evaluation operator*(const Evaluation& other) const
{
Evaluation result(*this);
result *= other;
return result;
}
template <class RhsValueType>
Evaluation operator*(const RhsValueType& other) const
{
Evaluation result(*this);
result *= other;
return result;
}
Evaluation operator/(const Evaluation& other) const
{
Evaluation result(*this);
result /= other;
return result;
}
template <class RhsValueType>
Evaluation operator/(const RhsValueType& other) const
{
Evaluation result(*this);
result /= other;
return result;
}
template <class RhsValueType>
Evaluation& operator=(const RhsValueType& other)
{
setValue( other );
clearDerivatives();
return *this;
}
// copy assignment from evaluation
Evaluation& operator=(const Evaluation& other)
{
data_ = other.data_;
return *this;
}
template <class RhsValueType>
bool operator==(const RhsValueType& other) const
{ return value() == other; }
bool operator==(const Evaluation& other) const
{
for (int idx = 0; idx < length_; ++idx)
if (data_[idx] != other.data_[idx])
return false;
return true;
}
bool operator!=(const Evaluation& other) const
{ return !operator==(other); }
template <class RhsValueType>
bool operator>(const RhsValueType& other) const
{ return value() > other; }
bool operator>(const Evaluation& other) const
{ return value() > other.value(); }
template <class RhsValueType>
bool operator<(const RhsValueType& other) const
{ return value() < other; }
bool operator<(const Evaluation& other) const
{ return value() < other.value(); }
template <class RhsValueType>
bool operator>=(const RhsValueType& other) const
{ return value() >= other; }
bool operator>=(const Evaluation& other) const
{ return value() >= other.value(); }
template <class RhsValueType>
bool operator<=(const RhsValueType& other) const
{ return value() <= other; }
bool operator<=(const Evaluation& other) const
{ return value() <= other.value(); }
// return value of variable
const ValueType& value() const
{ return data_[valuepos_]; }
// set value of variable
void setValue(const ValueType& val)
{ data_[valuepos_] = val; }
// return varIdx'th derivative
const ValueType& derivative(int varIdx) const
{
assert(varIdx < numVars);
return data_[varIdx + dstart_];
}
// set derivative at position varIdx
void setDerivative(int varIdx, const ValueType& derVal)
{
assert(varIdx < numVars);
data_[varIdx + dstart_] = derVal;
}
protected:
std::array<ValueType, length_> data_;
};
} // namespace DenseAD
} // namespace Dune
// #include <opm/material/densead/EvaluationSIMD.hpp>
#endif

458
opm/material/densead/Evaluation6.hpp Normal file
View File

@@ -0,0 +1,458 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 2 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
*
* \brief Representation of an evaluation of a function and its derivatives w.r.t. a set
* of variables in the localized OPM automatic differentiation (AD) framework.
*/
#ifndef OPM_LOCAL_AD_EVALUATION_6_HPP
#define OPM_LOCAL_AD_EVALUATION_6_HPP
#include "Math.hpp"
#include <opm/common/Valgrind.hpp>
#include <dune/common/version.hh>
#include <array>
#include <cmath>
#include <cassert>
#include <iostream>
#include <algorithm>
namespace Opm {
namespace DenseAd {
/*!
* \brief Represents a function evaluation and its derivatives w.r.t. a fixed set of
* variables.
*/
template <class ValueT>
class Evaluation< ValueT, 6 >
{
static constexpr int numVars = 6;
public:
//! field type
typedef ValueT ValueType;
//! number of derivatives
static constexpr int size = numVars;
protected:
//! length of internal data vector
static constexpr int length_ = numVars + 1 ;
//! position index for value
static constexpr int valuepos_ = 0;
//! start index for derivatives
static constexpr int dstart_ = 1;
//! end+1 index for derivatives
static constexpr int dend_ = length_ ;
public:
//! default constructor
Evaluation() : data_()
{}
//! copy other function evaluation
Evaluation(const Evaluation& other) : data_( other.data_ )
{
}
// create an evaluation which represents a constant function
//
// i.e., f(x) = c. this implies an evaluation with the given value and all
// derivatives being zero.
template <class RhsValueType>
Evaluation(const RhsValueType& c)
{
setValue( c );
clearDerivatives();
Valgrind::CheckDefined( data_ );
}
// create an evaluation which represents a constant function
//
// i.e., f(x) = c. this implies an evaluation with the given value and all
// derivatives being zero.
template <class RhsValueType>
Evaluation(const RhsValueType& c, int varPos)
{
setValue( c );
clearDerivatives();
// The variable position must be in represented by the given variable descriptor
assert(0 <= varPos && varPos < numVars);
data_[varPos + dstart_] = 1.0;
Valgrind::CheckDefined(data_);
}
// create a function evaluation for a "naked" depending variable (i.e., f(x) = x)
template <class RhsValueType>
static Evaluation devide(const RhsValueType& a, const Evaluation& b )
{
Evaluation<ValueType, numVars> result;
result.setValue( a/b.value() );
const ValueType df_dg = - result.value()/b.value();
for (int idx = dstart_; idx < dend_; ++idx) {
result.data_[idx] = df_dg*b.data_[idx];
}
return result;
}
// set all derivatives to zero
void clearDerivatives()
{
data_[ 1 ] = 0;
data_[ 2 ] = 0;
data_[ 3 ] = 0;
data_[ 4 ] = 0;
data_[ 5 ] = 0;
data_[ 6 ] = 0;
}
// create a function evaluation for a "naked" depending variable (i.e., f(x) = x)
template <class RhsValueType>
static Evaluation createVariable(const RhsValueType& value, int varPos)
{
// 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)
return Evaluation( value, varPos );
}
// "evaluate" a constant function (i.e. a function that does not depend on the set of
// relevant variables, f(x) = c).
template <class RhsValueType>
static Evaluation createConstant(const RhsValueType& value)
{
return Evaluation( value );
}
// print the value and the derivatives of the function evaluation
void print(std::ostream& os = std::cout) const
{
// print value
os << "v: " << value() << " / d:";
// print derivatives
for (int varIdx = 0; varIdx < numVars; ++varIdx)
os << " " << derivative(varIdx);
}
// copy all derivatives from other
void copyDerivatives(const Evaluation& other)
{
data_[ 1 ] = other.data_[ 1 ];
data_[ 2 ] = other.data_[ 2 ];
data_[ 3 ] = other.data_[ 3 ];
data_[ 4 ] = other.data_[ 4 ];
data_[ 5 ] = other.data_[ 5 ];
data_[ 6 ] = other.data_[ 6 ];
}
// add value and derivatives from other to this values and derivatives
Evaluation& operator+=(const Evaluation& other)
{
data_[ 0 ] += other.data_[ 0 ];
data_[ 1 ] += other.data_[ 1 ];
data_[ 2 ] += other.data_[ 2 ];
data_[ 3 ] += other.data_[ 3 ];
data_[ 4 ] += other.data_[ 4 ];
data_[ 5 ] += other.data_[ 5 ];
data_[ 6 ] += other.data_[ 6 ];
return *this;
}
// add value from other to this values
template <class RhsValueType>
Evaluation& operator+=(const RhsValueType& other)
{
// value is added, derivatives stay the same
data_[valuepos_] += other;
return *this;
}
// subtract other's value and derivatives from this values
Evaluation& operator-=(const Evaluation& other)
{
// value and derivatives are subtracted
data_[ 0 ] -= other.data_[ 0 ];
data_[ 1 ] -= other.data_[ 1 ];
data_[ 2 ] -= other.data_[ 2 ];
data_[ 3 ] -= other.data_[ 3 ];
data_[ 4 ] -= other.data_[ 4 ];
data_[ 5 ] -= other.data_[ 5 ];
data_[ 6 ] -= other.data_[ 6 ];
return *this;
}
// subtract other's value from this values
template <class RhsValueType>
Evaluation& operator-=(const RhsValueType& other)
{
// for constants, values are subtracted, derivatives stay the same
data_[ valuepos_ ] -= other;
return *this;
}
// multiply values and apply chain rule to derivatives: (u*v)' = (v'u + u'v)
Evaluation& operator*=(const Evaluation& other)
{
// while the values are multiplied, the derivatives follow the product rule,
// i.e., (u*v)' = (v'u + u'v).
const ValueType u = value();
const ValueType v = other.value();
data_[ 0 ] = u * v ;
data_[ 1 ] = data_[ 1 ] * v + other.data_[ 1 ] * u;
data_[ 2 ] = data_[ 2 ] * v + other.data_[ 2 ] * u;
data_[ 3 ] = data_[ 3 ] * v + other.data_[ 3 ] * u;
data_[ 4 ] = data_[ 4 ] * v + other.data_[ 4 ] * u;
data_[ 5 ] = data_[ 5 ] * v + other.data_[ 5 ] * u;
data_[ 6 ] = data_[ 6 ] * v + other.data_[ 6 ] * u;
return *this;
}
// m(u*v)' = (v'u + u'v)
template <class RhsValueType>
Evaluation& operator*=(RhsValueType other)
{
data_[ 0 ] *= other;
data_[ 1 ] *= other;
data_[ 2 ] *= other;
data_[ 3 ] *= other;
data_[ 4 ] *= other;
data_[ 5 ] *= other;
data_[ 6 ] *= other;
return *this;
}
// m(u*v)' = (v'u + u'v)
Evaluation& operator/=(const Evaluation& other)
{
// values are divided, derivatives follow the rule for division, i.e., (u/v)' = (v'u -
// u'v)/v^2.
const ValueType v_vv = 1.0 / other.value();
const ValueType u_vv = value() * v_vv * v_vv;
data_[ 0 ] *= v_vv;
data_[ 1 ] = data_[ 1 ] * v_vv - other.data_[ 1 ] * u_vv ;
data_[ 2 ] = data_[ 2 ] * v_vv - other.data_[ 2 ] * u_vv ;
data_[ 3 ] = data_[ 3 ] * v_vv - other.data_[ 3 ] * u_vv ;
data_[ 4 ] = data_[ 4 ] * v_vv - other.data_[ 4 ] * u_vv ;
data_[ 5 ] = data_[ 5 ] * v_vv - other.data_[ 5 ] * u_vv ;
data_[ 6 ] = data_[ 6 ] * v_vv - other.data_[ 6 ] * u_vv ;
return *this;
}
// multiply value and derivatives by value of other
template <class RhsValueType>
Evaluation& operator/=(const RhsValueType& other)
{
// values and derivatives are divided
ValueType factor = (1.0/other);
data_[ 0 ] *= factor;
data_[ 1 ] *= factor;
data_[ 2 ] *= factor;
data_[ 3 ] *= factor;
data_[ 4 ] *= factor;
data_[ 5 ] *= factor;
data_[ 6 ] *= factor;
return *this;
}
// add two evaluation objects
Evaluation operator+(const Evaluation& other) const
{
Evaluation result(*this);
result += other;
return result;
}
// add constant to this object
template <class RhsValueType>
Evaluation operator+(const RhsValueType& other) const
{
Evaluation result(*this);
result += other;
return result;
}
// subtract two evaluation objects
Evaluation operator-(const Evaluation& other) const
{
Evaluation result(*this);
result -= other;
return result;
}
// subtract constant from evaluation object
template <class RhsValueType>
Evaluation operator-(const RhsValueType& other) const
{
Evaluation result(*this);
result -= other;
return result;
}
// negation (unary minus) operator
Evaluation operator-() const
{
Evaluation result;
// set value and derivatives to negative
for (int idx = 0; idx < length_; ++idx)
result.data_[idx] = - data_[idx];
return result;
}
Evaluation operator*(const Evaluation& other) const
{
Evaluation result(*this);
result *= other;
return result;
}
template <class RhsValueType>
Evaluation operator*(const RhsValueType& other) const
{
Evaluation result(*this);
result *= other;
return result;
}
Evaluation operator/(const Evaluation& other) const
{
Evaluation result(*this);
result /= other;
return result;
}
template <class RhsValueType>
Evaluation operator/(const RhsValueType& other) const
{
Evaluation result(*this);
result /= other;
return result;
}
template <class RhsValueType>
Evaluation& operator=(const RhsValueType& other)
{
setValue( other );
clearDerivatives();
return *this;
}
// copy assignment from evaluation
Evaluation& operator=(const Evaluation& other)
{
data_ = other.data_;
return *this;
}
template <class RhsValueType>
bool operator==(const RhsValueType& other) const
{ return value() == other; }
bool operator==(const Evaluation& other) const
{
for (int idx = 0; idx < length_; ++idx)
if (data_[idx] != other.data_[idx])
return false;
return true;
}
bool operator!=(const Evaluation& other) const
{ return !operator==(other); }
template <class RhsValueType>
bool operator>(RhsValueType other) const
{ return value() > other; }
bool operator>(const Evaluation& other) const
{ return value() > other.value(); }
template <class RhsValueType>
bool operator<(RhsValueType other) const
{ return value() < other; }
bool operator<(const Evaluation& other) const
{ return value() < other.value(); }
template <class RhsValueType>
bool operator>=(RhsValueType other) const
{ return value() >= other; }
bool operator>=(const Evaluation& other) const
{ return value() >= other.value(); }
template <class RhsValueType>
bool operator<=(RhsValueType other) const
{ return value() <= other; }
bool operator<=(const Evaluation& other) const
{ return value() <= other.value(); }
// return value of variable
const ValueType& value() const
{ return data_[valuepos_]; }
// set value of variable
void setValue(const ValueType& val)
{ data_[valuepos_] = val; }
// return varIdx'th derivative
const ValueType& derivative(int varIdx) const
{
assert(varIdx < numVars);
return data_[varIdx + dstart_];
}
// set derivative at position varIdx
void setDerivative(int varIdx, const ValueType& derVal)
{
assert(varIdx < numVars);
data_[varIdx + dstart_] = derVal;
}
protected:
std::array<ValueType, length_> data_;
};
} // namespace DenseAD
} // namespace Dune
// #include <opm/material/densead/EvaluationSIMD.hpp>
#endif

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

@@ -39,45 +39,45 @@
namespace Opm {
namespace DenseAd {
// forward declaration of the Evaluation template class
template <class ValueT, unsigned numVars>
template <class ValueT, int numVars>
class Evaluation;
// provide some algebraic functions
template <class ValueType, unsigned numVars>
template <class ValueType, int numVars>
Evaluation<ValueType, numVars> abs(const Evaluation<ValueType, numVars>& x)
{ return (x > 0.0)?x:-x; }
template <class ValueType, unsigned numVars>
template <class ValueType, int numVars>
Evaluation<ValueType, numVars> min(const Evaluation<ValueType, numVars>& x1,
const Evaluation<ValueType, numVars>& x2)
{ return (x1 < x2)?x1:x2; }
template <class Arg1ValueType, class ValueType, unsigned numVars>
template <class Arg1ValueType, class ValueType, int numVars>
Evaluation<ValueType, numVars> min(const Arg1ValueType& x1,
const Evaluation<ValueType, numVars>& x2)
{ return (x1 < x2)?x1:x2; }
template <class ValueType, unsigned numVars, class Arg2ValueType>
template <class ValueType, int numVars, class Arg2ValueType>
Evaluation<ValueType, numVars> min(const Evaluation<ValueType, numVars>& x1,
const Arg2ValueType& x2)
{ return min(x2, x1); }
template <class ValueType, unsigned numVars>
template <class ValueType, int numVars>
Evaluation<ValueType, numVars> max(const Evaluation<ValueType, numVars>& x1,
const Evaluation<ValueType, numVars>& x2)
{ return (x1 > x2)?x1:x2; }
template <class Arg1ValueType, class ValueType, unsigned numVars>
template <class Arg1ValueType, class ValueType, int numVars>
Evaluation<ValueType, numVars> max(const Arg1ValueType& x1,
const Evaluation<ValueType, numVars>& x2)
{ return (x1 > x2)?x1:x2; }
template <class ValueType, unsigned numVars, class Arg2ValueType>
template <class ValueType, int numVars, class Arg2ValueType>
Evaluation<ValueType, numVars> max(const Evaluation<ValueType, numVars>& x1,
const Arg2ValueType& x2)
{ return max(x2, x1); }
template <class ValueType, unsigned numVars>
template <class ValueType, int numVars>
Evaluation<ValueType, numVars> tan(const Evaluation<ValueType, numVars>& x)
{
typedef MathToolbox<ValueType> ValueTypeToolbox;
@@ -89,13 +89,13 @@ Evaluation<ValueType, numVars> tan(const Evaluation<ValueType, numVars>& x)
// derivatives use the chain rule
const ValueType& df_dx = 1 + tmp*tmp;
for (unsigned curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
return result;
}
template <class ValueType, unsigned numVars>
template <class ValueType, int numVars>
Evaluation<ValueType, numVars> atan(const Evaluation<ValueType, numVars>& x)
{
typedef MathToolbox<ValueType> ValueTypeToolbox;
@@ -106,13 +106,13 @@ Evaluation<ValueType, numVars> atan(const Evaluation<ValueType, numVars>& x)
// derivatives use the chain rule
const ValueType& df_dx = 1/(1 + x.value()*x.value());
for (unsigned curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
return result;
}
template <class ValueType, unsigned numVars>
template <class ValueType, int numVars>
Evaluation<ValueType, numVars> atan2(const Evaluation<ValueType, numVars>& x,
const Evaluation<ValueType, numVars>& y)
{
@@ -124,7 +124,7 @@ Evaluation<ValueType, numVars> atan2(const Evaluation<ValueType, numVars>& x,
// derivatives use the chain rule
const ValueType& alpha = 1/(1 + (x.value()*x.value())/(y.value()*y.value()));
for (unsigned curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) {
for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) {
result.setDerivative(curVarIdx,
alpha/(y.value()*y.value())
*(x.derivative(curVarIdx)*y.value() - x.value()*y.derivative(curVarIdx)));
@@ -133,7 +133,7 @@ Evaluation<ValueType, numVars> atan2(const Evaluation<ValueType, numVars>& x,
return result;
}
template <class ValueType, unsigned numVars>
template <class ValueType, int numVars>
Evaluation<ValueType, numVars> sin(const Evaluation<ValueType, numVars>& x)
{
typedef MathToolbox<ValueType> ValueTypeToolbox;
@@ -144,13 +144,13 @@ Evaluation<ValueType, numVars> sin(const Evaluation<ValueType, numVars>& x)
// derivatives use the chain rule
const ValueType& df_dx = ValueTypeToolbox::cos(x.value());
for (unsigned curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
return result;
}
template <class ValueType, unsigned numVars>
template <class ValueType, int numVars>
Evaluation<ValueType, numVars> asin(const Evaluation<ValueType, numVars>& x)
{
typedef MathToolbox<ValueType> ValueTypeToolbox;
@@ -161,13 +161,13 @@ Evaluation<ValueType, numVars> asin(const Evaluation<ValueType, numVars>& x)
// derivatives use the chain rule
const ValueType& df_dx = 1.0/ValueTypeToolbox::sqrt(1 - x.value()*x.value());
for (unsigned curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
return result;
}
template <class ValueType, unsigned numVars>
template <class ValueType, int numVars>
Evaluation<ValueType, numVars> cos(const Evaluation<ValueType, numVars>& x)
{
typedef MathToolbox<ValueType> ValueTypeToolbox;
@@ -178,13 +178,13 @@ Evaluation<ValueType, numVars> cos(const Evaluation<ValueType, numVars>& x)
// derivatives use the chain rule
const ValueType& df_dx = -ValueTypeToolbox::sin(x.value());
for (unsigned curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
return result;
}
template <class ValueType, unsigned numVars>
template <class ValueType, int numVars>
Evaluation<ValueType, numVars> acos(const Evaluation<ValueType, numVars>& x)
{
typedef MathToolbox<ValueType> ValueTypeToolbox;
@@ -195,13 +195,13 @@ Evaluation<ValueType, numVars> acos(const Evaluation<ValueType, numVars>& x)
// derivatives use the chain rule
const ValueType& df_dx = - 1.0/ValueTypeToolbox::sqrt(1 - x.value()*x.value());
for (unsigned curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
return result;
}
template <class ValueType, unsigned numVars>
template <class ValueType, int numVars>
Evaluation<ValueType, numVars> sqrt(const Evaluation<ValueType, numVars>& x)
{
typedef MathToolbox<ValueType> ValueTypeToolbox;
@@ -213,14 +213,14 @@ Evaluation<ValueType, numVars> sqrt(const Evaluation<ValueType, numVars>& x)
// derivatives use the chain rule
ValueType df_dx = 0.5/sqrt_x;
for (unsigned curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) {
for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) {
result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
}
return result;
}
template <class ValueType, unsigned numVars>
template <class ValueType, int numVars>
Evaluation<ValueType, numVars> exp(const Evaluation<ValueType, numVars>& x)
{
typedef MathToolbox<ValueType> ValueTypeToolbox;
@@ -231,14 +231,14 @@ Evaluation<ValueType, numVars> exp(const Evaluation<ValueType, numVars>& x)
// derivatives use the chain rule
const ValueType& df_dx = exp_x;
for (unsigned curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
return result;
}
// exponentiation of arbitrary base with a fixed constant
template <class ValueType, unsigned numVars, class ExpType>
template <class ValueType, int numVars, class ExpType>
Evaluation<ValueType, numVars> pow(const Evaluation<ValueType, numVars>& base,
const ExpType& exp)
{
@@ -256,7 +256,7 @@ 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.size; ++curVarIdx)
for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
result.setDerivative(curVarIdx, df_dx*base.derivative(curVarIdx));
}
@@ -264,7 +264,7 @@ Evaluation<ValueType, numVars> pow(const Evaluation<ValueType, numVars>& base,
}
// exponentiation of constant base with an arbitrary exponent
template <class BaseType, class ValueType, unsigned numVars>
template <class BaseType, class ValueType, int numVars>
Evaluation<ValueType, numVars> pow(const BaseType& base,
const Evaluation<ValueType, numVars>& exp)
{
@@ -283,7 +283,7 @@ Evaluation<ValueType, numVars> pow(const BaseType& base,
// derivatives use the chain rule
const ValueType& df_dx = lnBase*result.value();
for (unsigned curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
result.setDerivative(curVarIdx, df_dx*exp.derivative(curVarIdx));
}
@@ -292,7 +292,7 @@ Evaluation<ValueType, numVars> pow(const BaseType& base,
// this is the most expensive power function. Computationally it is pretty expensive, so
// one of the above two variants above should be preferred if possible.
template <class ValueType, unsigned numVars>
template <class ValueType, int numVars>
Evaluation<ValueType, numVars> pow(const Evaluation<ValueType, numVars>& base,
const Evaluation<ValueType, numVars>& exp)
{
@@ -314,7 +314,7 @@ Evaluation<ValueType, numVars> pow(const Evaluation<ValueType, numVars>& base,
const ValueType& f = base.value();
const ValueType& g = exp.value();
const ValueType& logF = ValueTypeToolbox::log(f);
for (unsigned curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) {
for (int 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);
@@ -324,7 +324,7 @@ Evaluation<ValueType, numVars> pow(const Evaluation<ValueType, numVars>& base,
return result;
}
template <class ValueType, unsigned numVars>
template <class ValueType, int numVars>
Evaluation<ValueType, numVars> log(const Evaluation<ValueType, numVars>& x)
{
typedef MathToolbox<ValueType> ValueTypeToolbox;
@@ -335,7 +335,7 @@ Evaluation<ValueType, numVars> log(const Evaluation<ValueType, numVars>& x)
// derivatives use the chain rule
const ValueType& df_dx = 1/x.value();
for (unsigned curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
return result;
@@ -345,7 +345,7 @@ Evaluation<ValueType, numVars> log(const Evaluation<ValueType, numVars>& x)
// a kind of traits class for the automatic differentiation case. (The toolbox for the
// scalar case is provided by the MathToolbox.hpp header file.)
template <class ValueT, unsigned numVars>
template <class ValueT, int numVars>
struct MathToolbox<Opm::DenseAd::Evaluation<ValueT, numVars> >
{
private:
@@ -364,7 +364,7 @@ public:
static Evaluation createConstant(ValueType value)
{ return Evaluation::createConstant(value); }
static Evaluation createVariable(ValueType value, unsigned varIdx)
static Evaluation createVariable(ValueType value, int varIdx)
{ return Evaluation::createVariable(value, varIdx); }
template <class LhsEval>
@@ -395,7 +395,7 @@ public:
return false;
// make sure that the derivatives are identical
for (unsigned curVarIdx = 0; curVarIdx < numVars; ++curVarIdx)
for (int curVarIdx = 0; curVarIdx < numVars; ++curVarIdx)
if (!ValueTypeToolbox::isSame(a.derivative(curVarIdx), b.derivative(curVarIdx), tolerance))
return false;
@@ -460,7 +460,7 @@ public:
if (!InnerToolbox::isfinite(arg.value()))
return false;
for (unsigned i = 0; i < numVars; ++i)
for (int i = 0; i < numVars; ++i)
if (!InnerToolbox::isfinite(arg.derivative(i)))
return false;
@@ -472,7 +472,7 @@ public:
if (InnerToolbox::isnan(arg.value()))
return true;
for (unsigned i = 0; i < numVars; ++i)
for (int i = 0; i < numVars; ++i)
if (InnerToolbox::isnan(arg.derivative(i)))
return true;

87
tests/test_densead.cpp
View File

@@ -50,9 +50,12 @@
#include <cassert>
#include <stdexcept>
static const int numVars = 3;
//static const int numVars = 3;
template <class Scalar, int numVars>
struct TestEnv
{
template <class Scalar>
void testOperators(const Scalar tolerance)
{
typedef Opm::DenseAd::Evaluation<Scalar, numVars> Eval;
@@ -256,7 +259,7 @@ void testOperators(const Scalar tolerance)
}
}
template <class Scalar, class AdFn, class ClassicFn>
template <class AdFn, class ClassicFn>
void test1DFunction(AdFn* adFn, ClassicFn* classicFn, Scalar xMin = 1e-6, Scalar xMax = 1000)
{
typedef Opm::DenseAd::Evaluation<Scalar, numVars> Eval;
@@ -289,8 +292,7 @@ void test1DFunction(AdFn* adFn, ClassicFn* classicFn, Scalar xMin = 1e-6, Scalar
}
}
template <class Scalar,
class AdFn,
template <class AdFn,
class ClassicFn>
void test2DFunction1(AdFn* adFn, ClassicFn* classicFn, Scalar xMin, Scalar xMax, Scalar y)
{
@@ -326,8 +328,7 @@ void test2DFunction1(AdFn* adFn, ClassicFn* classicFn, Scalar xMin, Scalar xMax,
}
}
template <class Scalar,
class AdFn,
template <class AdFn,
class ClassicFn>
void test2DFunction2(AdFn* adFn, ClassicFn* classicFn, Scalar x, Scalar yMin, Scalar yMax)
{
@@ -363,7 +364,6 @@ void test2DFunction2(AdFn* adFn, ClassicFn* classicFn, Scalar x, Scalar yMin, Sc
}
}
template <class Scalar>
void testPowBase(Scalar baseMin = 1e-2, Scalar baseMax = 100)
{
typedef Opm::DenseAd::Evaluation<Scalar, numVars> Eval;
@@ -405,7 +405,6 @@ void testPowBase(Scalar baseMin = 1e-2, Scalar baseMax = 100)
}
}
template <class Scalar>
void testPowExp(Scalar expMin = -100, Scalar expMax = 100)
{
typedef Opm::DenseAd::Evaluation<Scalar, numVars> Eval;
@@ -447,7 +446,6 @@ void testPowExp(Scalar expMin = -100, Scalar expMax = 100)
}
}
template <class Scalar>
void testAtan2()
{
typedef Opm::DenseAd::Evaluation<Scalar, numVars> Eval;
@@ -497,16 +495,12 @@ void testAtan2()
}
// prototypes
double myScalarMin(double a, double b);
double myScalarMax(double a, double b);
double myScalarMin(double a, double b)
static double myScalarMin(double a, double b)
{ return std::min(a, b); }
double myScalarMax(double a, double b)
static double myScalarMax(double a, double b)
{ return std::max(a, b); }
template <class Scalar>
inline void testAll()
{
// the following is commented out because it is supposed to produce a compiler
@@ -516,82 +510,82 @@ inline void testAll()
std::cout << "testing operators and constructors\n";
const Scalar eps = std::numeric_limits<Scalar>::epsilon()*1e3;
testOperators<Scalar>(eps);
testOperators(eps);
std::cout << "testing min()\n";
test2DFunction1<Scalar>(Opm::DenseAd::min<Scalar, numVars>,
test2DFunction1(Opm::DenseAd::min<Scalar, numVars>,
myScalarMin,
-1000, 1000,
/*p=*/1.234);
test2DFunction2<Scalar>(Opm::DenseAd::min<Scalar, numVars>,
test2DFunction2(Opm::DenseAd::min<Scalar, numVars>,
myScalarMin,
/*T=*/1.234,
-1000, 1000);
std::cout << "testing max()\n";
test2DFunction1<Scalar>(Opm::DenseAd::max<Scalar, numVars>,
test2DFunction1(Opm::DenseAd::max<Scalar, numVars>,
myScalarMax,
-1000, 1000,
/*p=*/1.234);
test2DFunction2<Scalar>(Opm::DenseAd::max<Scalar, numVars>,
test2DFunction2(Opm::DenseAd::max<Scalar, numVars>,
myScalarMax,
/*T=*/1.234,
-1000, 1000);
std::cout << "testing pow()\n";
testPowBase<Scalar>();
testPowExp<Scalar>();
testPowBase();
testPowExp();
std::cout << "testing abs()\n";
test1DFunction<Scalar>(Opm::DenseAd::abs<Scalar, numVars>,
test1DFunction(Opm::DenseAd::abs<Scalar, numVars>,
static_cast<Scalar (*)(Scalar)>(std::abs));
std::cout << "testing sqrt()\n";
test1DFunction<Scalar>(Opm::DenseAd::sqrt<Scalar, numVars>,
test1DFunction(Opm::DenseAd::sqrt<Scalar, numVars>,
static_cast<Scalar (*)(Scalar)>(std::sqrt));
std::cout << "testing sin()\n";
test1DFunction<Scalar>(Opm::DenseAd::sin<Scalar, numVars>,
test1DFunction(Opm::DenseAd::sin<Scalar, numVars>,
static_cast<Scalar (*)(Scalar)>(std::sin),
0, 2*M_PI);
std::cout << "testing asin()\n";
test1DFunction<Scalar>(Opm::DenseAd::asin<Scalar, numVars>,
test1DFunction(Opm::DenseAd::asin<Scalar, numVars>,
static_cast<Scalar (*)(Scalar)>(std::asin),
-1.0, 1.0);
std::cout << "testing cos()\n";
test1DFunction<Scalar>(Opm::DenseAd::cos<Scalar, numVars>,
test1DFunction(Opm::DenseAd::cos<Scalar, numVars>,
static_cast<Scalar (*)(Scalar)>(std::cos),
0, 2*M_PI);
std::cout << "testing acos()\n";
test1DFunction<Scalar>(Opm::DenseAd::acos<Scalar, numVars>,
test1DFunction(Opm::DenseAd::acos<Scalar, numVars>,
static_cast<Scalar (*)(Scalar)>(std::acos),
-1.0, 1.0);
std::cout << "testing tan()\n";
test1DFunction<Scalar>(Opm::DenseAd::tan<Scalar, numVars>,
test1DFunction(Opm::DenseAd::tan<Scalar, numVars>,
static_cast<Scalar (*)(Scalar)>(std::tan),
-M_PI / 2 * 0.95, M_PI / 2 * 0.95);
std::cout << "testing atan()\n";
test1DFunction<Scalar>(Opm::DenseAd::atan<Scalar, numVars>,
test1DFunction(Opm::DenseAd::atan<Scalar, numVars>,
static_cast<Scalar (*)(Scalar)>(std::atan),
-10*1000.0, 10*1000.0);
std::cout << "testing atan2()\n";
testAtan2<Scalar>();
testAtan2();
std::cout << "testing exp()\n";
test1DFunction<Scalar>(Opm::DenseAd::exp<Scalar, numVars>,
test1DFunction(Opm::DenseAd::exp<Scalar, numVars>,
static_cast<Scalar (*)(Scalar)>(std::exp),
-100, 100);
std::cout << "testing log()\n";
test1DFunction<Scalar>(Opm::DenseAd::log<Scalar, numVars>,
test1DFunction(Opm::DenseAd::log<Scalar, numVars>,
static_cast<Scalar (*)(Scalar)>(std::log),
1e-6, 1e9);
@@ -651,12 +645,33 @@ inline void testAll()
}
}
};//TestEnv
template <class Scalar, int numVars>
struct TestAll
{
static void test()
{
TestEnv< Scalar, numVars >().testAll();
TestAll< Scalar, numVars-1>::test();
}
};
template <class Scalar>
struct TestAll<Scalar, 2>
{
static void test()
{
TestEnv< Scalar, 2 >().testAll();
}
};
int main(int argc, char **argv)
{
Dune::MPIHelper::instance(argc, argv);
testAll<double>();
testAll<float>();
TestAll<double, 15>::test();
TestAll<float, 15>::test();
return 0;
}