Merge pull request #3 from andlaus/generated_eval_specializations_v2

auto-generate and specialize the whole Evaluation class
This commit is contained in:
dr-robertk 2017-03-17 17:52:18 +01:00 committed by GitHub
commit efa158589f
14 changed files with 5731 additions and 1442 deletions

View File

@ -49,125 +49,609 @@ specializationTemplate = \
/*!
* \\file
*
{% if numDerivs < 0 %}
* \\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.
{% else %}
* \\brief This file specializes the dense-AD Evaluation class for {{ numDerivs }} derivatives.
{% endif %}
*
* \\attention THIS FILE GETS AUTOMATICALLY GENERATED BY THE "{{ scriptName }}"
* SCRIPT. DO NOT EDIT IT MANUALLY!
*/
{% if numDerivs < 0 %}
#ifndef OPM_DENSEAD_EVALUATION_HPP
#define OPM_DENSEAD_EVALUATION_HPP
{% else %}
#ifndef OPM_DENSEAD_EVALUATION{{numDerivs}}_HPP
#define OPM_DENSEAD_EVALUATION{{numDerivs}}_HPP
{% endif %}
#include "Math.hpp"
#include <opm/common/Valgrind.hpp>
#include <dune/common/version.hh>
#include <array>
#include <cmath>
#include <cassert>
#include <cstring>
#include <iostream>
#include <algorithm>
namespace Opm {
namespace DenseAd {
{% if numDerivs < 0 %}
/*!
* \\brief Represents a function evaluation and its derivatives w.r.t. a fixed set of
* variables.
*/
template <class ValueT, int numDerivs>
class Evaluation
{% else %}
template <class ValueT>
struct EvaluationOps<ValueT, {{numDerivs}}>
class Evaluation<ValueT, {{ numDerivs }}>
{% endif %}
{
private:
typedef Evaluation<ValueT, {{numDerivs}} > Eval;
public:
//! field type
typedef ValueT ValueType;
//! number of derivatives
{% if numDerivs < 0 %}
static constexpr int size = numDerivs;
{% else %}
static constexpr int size = {{ numDerivs }};
{% endif %}
protected:
//! length of internal data vector
static constexpr int length_ = size + 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:
typedef std::array<ValueT, {{numDerivs + 1}} > DataVector;
static inline void assign(Eval& a, const Eval& b)
//! default constructor
Evaluation() : data_()
{}
//! copy other function evaluation
Evaluation(const Evaluation& other)
{% if numDerivs < 0 %}
: data_(other.data_)
{ }
{% else %}
{
{% for i in range(0, numDerivs+1) %}
a.data_[{{i}}] = b.data_[{{i}}];{% endfor %}
data_[{{i}}] = other.data_[{{i}}];{% endfor %}
}
{% endif %}
// 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_ );
}
static inline void assignNegative(Eval& a, const Eval& b)
// 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)
{
{% for i in range(0, numDerivs+1) %}
a.data_[{{i}}] = - b.data_[{{i}}];{% endfor %}
setValue( c );
clearDerivatives();
// The variable position must be in represented by the given variable descriptor
assert(0 <= varPos && varPos < size);
data_[varPos + dstart_] = 1.0;
Valgrind::CheckDefined(data_);
}
static inline void clearDerivatives(Eval& a)
// set all derivatives to zero
void clearDerivatives()
{
{% if numDerivs < 0 %}
for (int i = dstart_; i < dend_; ++i)
data_[i] = 0.0;
{% else %}
{% for i in range(1, numDerivs+1) %}
a.data_[{{i}}] = 0.0;{% endfor %}
data_[{{i}}] = 0.0;{% endfor %}
{% endif %}
}
static inline Eval& addEq(Eval& a, const Eval& b)
// 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)
{
{% for i in range(0, numDerivs+1) %}
a.data_[{{i}}] += b.data_[{{i}}];{% endfor %}
return a;
// 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 );
}
static inline Eval& subEq(Eval& a, const Eval& b)
// "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)
{
{% for i in range(0, numDerivs+1) %}
a.data_[{{i}}] -= b.data_[{{i}}];{% endfor %}
return a;
return Evaluation( value );
}
static inline Eval& mulEq(Eval& a, const Eval& b)
// 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 < size; ++varIdx)
os << " " << derivative(varIdx);
}
// copy all derivatives from other
void copyDerivatives(const Evaluation& other)
{
for (int i = dstart_; 0 < dend_; ++i)
data_[i] = other.data_[i];
}
// add value and derivatives from other to this values and derivatives
Evaluation& operator+=(const Evaluation& other)
{
{% if numDerivs < 0 %}
for (int i = 0; i < length_; ++i)
data_[i] += other.data_[i];
{% else %}
{% for i in range(0, numDerivs+1) %}
data_[{{i}}] += other.data_[{{i}}];{% endfor %}
{% endif %}
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)
{
{% if numDerivs < 0 %}
for (int i = 0; i < length_; ++i)
data_[i] -= other.data_[i];
{% else %}
{% for i in range(0, numDerivs+1) %}
data_[{{i}}] -= other.data_[{{i}}];{% endfor %}
{% endif %}
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 ValueT u = a.value();
const ValueT v = b.value();
const ValueT u = this->value();
const ValueT v = other.value();
// value
a.data_[0] *= v ;
this->data_[valuepos_] *= v ;
// derivatives{% for i in range(1, numDerivs+1) %}
a.data_[{{i}}] = a.data_[{{i}}]*v + b.data_[{{i}}] * u;{% endfor %}
// derivatives
{% if numDerivs < 0 %}
for (int i = 0; i < length_; ++i)
this->data_[i] = this->data_[i]*v + other.data_[i] * u;
{% else %}
{% for i in range(1, numDerivs+1) %}
this->data_[{{i}}] = this->data_[{{i}}]*v + other.data_[{{i}}] * u;{% endfor %}
{% endif %}
return a;
return *this;
}
template <class RhsType>
static inline Eval& scalarMulEq(Eval& a, const RhsType& other)
// m(c*u)' = c*u'
template <class RhsValueType>
Evaluation& operator*=(const RhsValueType& other)
{
{% if numDerivs < 0 %}
for (int i = 0; i < length_; ++i)
data_[i] *= other;
{% else %}
{% for i in range(0, numDerivs+1) %}
a.data_[{{i}}] *= other;{% endfor %}
data_[{{i}}] *= other;{% endfor %}
{% endif %}
return a;
return *this;
}
static inline Eval& divEq(Eval& a, const Eval& b)
// 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 ValueT v_vv = 1.0 / b.value();
const ValueT u_vv = a.value() * v_vv * v_vv;
const ValueT v_vv = 1.0 / other.value();
const ValueT u_vv = value() * v_vv * v_vv;
// value
a.data_[0] *= v_vv;
data_[valuepos_] *= v_vv;
// derivatives{% for i in range(1, numDerivs+1) %}
a.data_[{{i}}] = a.data_[{{i}}]*v_vv - b.data_[{{i}}]*u_vv;{% endfor %}
// derivatives
{% if numDerivs < 0 %}
for (int i = 0; i < length_; ++i)
data_[i] = data_[i]*v_vv - other.data_[i]*u_vv;
{% else %}
{% for i in range(1, numDerivs+1) %}
data_[{{i}}] = data_[{{i}}]*v_vv - other.data_[{{i}}]*u_vv;{% endfor %}
{% endif %}
return a;
}
template <class RhsType>
static inline Eval& scalarDivEq(Eval& a, const RhsType& other)
{
{% for i in range(0, numDerivs+1) %}
a.data_[{{i}}] /= other;{% endfor %}
return a;
return *this;
}
// divide value and derivatives by value of other
template <class RhsValueType>
static inline Eval divide(const RhsValueType& a, const Eval& b)
Evaluation& operator/=(const RhsValueType& other)
{
Eval result;
result.setValue( a/b.value() );
const ValueT df_dg = - result.value()/b.value();
ValueType tmp = 1.0/other;
{% if numDerivs < 0 %}
for (int i = 0; i < length_; ++i)
data_[i] *= tmp;
{% else %}
{% for i in range(0, numDerivs+1) %}
data_[{{i}}] *= tmp;{% endfor %}
{% endif %}
return *this;
}
// division of a constant by an Evaluation
template <class RhsValueType>
static inline Evaluation divide(const RhsValueType& a, const Evaluation& b)
{
Evaluation result;
ValueType tmp = 1.0/b.value();
result.setValue( a*tmp );
const ValueT df_dg = - result.value()*tmp;
{% if numDerivs < 0 %}
for (int i = dstart_; i < dend_; ++i)
result.data_[i] = df_dg*b.data_[i];
{% else %}
{% for i in range(1, numDerivs+1) %}
result.data_[{{i}}] = df_dg*b.data_[{{i}}];{% endfor %}
{% endif %}
return result;
}
// 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);
return (result -= other);
}
// subtract constant from evaluation object
template <class RhsValueType>
Evaluation operator-(const RhsValueType& other) const
{
Evaluation result(*this);
return (result -= other);
}
// negation (unary minus) operator
Evaluation operator-() const
{
Evaluation result;
// set value and derivatives to negative
{% if numDerivs < 0 %}
for (int i = 0; i < length_; ++i)
result.data_[i] = - data_[i];
{% else %}
{% for i in range(0, numDerivs+1) %}
result.data_[{{i}}] = - data_[{{i}}];{% endfor %}
{% endif %}
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)
{
{% if numDerivs < 0 %}
for (int i = 0; i < length_; ++i)
data_[i] = other.data_[i];
{% else %}
{% for i in range(0, numDerivs+1) %}
data_[{{i}}] = other.data_[{{i}}];{% endfor %}
{% endif %}
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
template <class RhsValueType>
void setValue(const RhsValueType& val)
{ data_[valuepos_] = val; }
// return varIdx'th derivative
const ValueType& derivative(int varIdx) const
{
assert(0 <= varIdx && varIdx < size);
return data_[dstart_ + varIdx];
}
// set derivative at position varIdx
void setDerivative(int varIdx, const ValueType& derVal)
{
assert(0 <= varIdx && varIdx < size);
data_[dstart_ + varIdx] = derVal;
}
private:
std::array<ValueT, length_> data_;
};
{# the generic operators are only required for the unspecialized case #}
{% if numDerivs < 0 %}
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, int numVars>
bool operator>(const RhsValueType& a, const Evaluation<ValueType, numVars>& b)
{ return b < a; }
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, int numVars>
bool operator>=(const RhsValueType& a, const Evaluation<ValueType, numVars>& b)
{ return b <= a; }
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, int numVars>
Evaluation<ValueType, numVars> operator+(const RhsValueType& a, const Evaluation<ValueType, numVars>& b)
{
Evaluation<ValueType, numVars> result(b);
result += a;
return result;
}
template <class RhsValueType, class ValueType, int numVars>
Evaluation<ValueType, numVars> operator-(const RhsValueType& a, const Evaluation<ValueType, numVars>& b)
{
Evaluation<ValueType, numVars> result(a);
result -= b;
return result;
}
template <class RhsValueType, class ValueType, int numVars>
Evaluation<ValueType, numVars> operator/(const RhsValueType& a, const Evaluation<ValueType, numVars>& b)
{
return Evaluation<ValueType, numVars>::divide(a, b);
}
template <class RhsValueType, class ValueType, int numVars>
Evaluation<ValueType, numVars> operator*(const RhsValueType& a, const Evaluation<ValueType, numVars>& b)
{
Evaluation<ValueType, numVars> result(b);
result *= a;
return result;
}
template <class ValueType, int numVars>
std::ostream& operator<<(std::ostream& os, const Evaluation<ValueType, numVars>& eval)
{
os << eval.value();
return os;
}
{% endif %}
} } // namespace DenseAd, Opm
{% if numDerivs < 0 %}
// In Dune 2.3, the Evaluation.hpp header must be included before the fmatrix.hh
// header. Dune 2.4+ does not suffer from this because of some c++-foo.
//
// for those who are wondering: in C++ function templates cannot be partially
// specialized, and function argument overloads must be known _before_ they are used. The
// latter is what we do for the 'Dune::fvmeta::absreal()' function.
//
// consider the following test program:
//
// double foo(double i)
// { return i; }
//
// void bar()
// { std::cout << foo(0) << "\\n"; }
//
// int foo(int i)
// { return i + 1; }
//
// void foobar()
// { std::cout << foo(0) << "\\n"; }
//
// this will print '0' for bar() and '1' for foobar()...
#if !(DUNE_VERSION_NEWER(DUNE_COMMON, 2,4))
namespace Opm {
namespace DenseAd {
template <class ValueType, int numVars>
Evaluation<ValueType, numVars> abs(const Evaluation<ValueType, numVars>&);
}}
namespace std {
template <class ValueType, int numVars>
const Opm::DenseAd::Evaluation<ValueType, numVars> abs(const Opm::DenseAd::Evaluation<ValueType, numVars>& x)
{ return Opm::DenseAd::abs(x); }
} // namespace std
#if defined DUNE_DENSEMATRIX_HH
#warning \\
"Due to some C++ peculiarity regarding function overloads, the 'Evaluation.hpp'" \\
"header file must be included before Dune's 'densematrix.hh' for Dune < 2.4. " \\
"(If Evaluations are to be used in conjunction with a dense matrix.)"
#endif
#endif
// this makes the Dune matrix/vector classes happy...
#include <dune/common/ftraits.hh>
namespace Dune {
template <class ValueType, int numVars>
struct FieldTraits<Opm::DenseAd::Evaluation<ValueType, numVars> >
{
public:
typedef Opm::DenseAd::Evaluation<ValueType, numVars> field_type;
// setting real_type to field_type here potentially leads to slightly worse
// performance, but at least it makes things compile.
typedef field_type real_type;
};
} // namespace Dune
#include "EvaluationSpecializations.hpp"
#endif // OPM_DENSEAD_EVALUATION_HPP
{% else %}
#endif // OPM_DENSEAD_EVALUATION{{numDerivs}}_HPP
{% endif %}
"""
includeSpecializationsTemplate = \
@ -210,6 +694,15 @@ includeSpecializationsTemplate = \
#endif // OPM_DENSEAD_EVALUATION_SPECIALIZATIONS_HPP
"""
print "Generating generic template class"
fileName = "opm/material/densead/Evaluation.hpp"
template = jinja2.Template(specializationTemplate)
fileContents = template.render(numDerivs=-1, scriptName=os.path.basename(sys.argv[0]))
f = open(fileName, "w")
f.write(fileContents)
f.close()
for numDerivs in range(1, maxDerivs + 1):
print "Generating specialization for %d derivatives"%numDerivs

View File

@ -23,12 +23,18 @@
/*!
* \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.
*
* \attention THIS FILE GETS AUTOMATICALLY GENERATED BY THE "genEvalSpecializations.py"
* SCRIPT. DO NOT EDIT IT MANUALLY!
*/
#ifndef OPM_LOCAL_AD_EVALUATION_HPP
#define OPM_LOCAL_AD_EVALUATION_HPP
#ifndef OPM_DENSEAD_EVALUATION_HPP
#define OPM_DENSEAD_EVALUATION_HPP
#include "Math.hpp"
@ -45,154 +51,36 @@
namespace Opm {
namespace DenseAd {
template <class ValueT, int numVars>
class Evaluation;
/*!
* \brief Provides operations on Evaluations, so that these can be specialized without
* having to copy-and-paste the whole Evaluation class.
*
* \internal
*/
template <class ValueT, int numVars>
struct EvaluationOps
{
private:
typedef Evaluation<ValueT, numVars> Eval;
static constexpr int length_ = Eval::length_;
static constexpr int valuepos_ = Eval::valuepos_;
static constexpr int dstart_ = Eval::dstart_;
static constexpr int dend_ = Eval::dend_ ;
public:
typedef std::array<ValueT, length_> DataVector;
static inline void clearDerivatives(Eval& a)
{
for (int i = dstart_; i < dend_; ++i)
a.data_[i] = 0.0;
}
static inline void assign(Eval& a, const Eval& b)
{
a.data_ = b.data_;
}
static inline void assignNegative(Eval& a, const Eval& b)
{
for (int idx = 0; idx < length_; ++idx)
a.data_[idx] = - b.data_[idx];
}
static inline Eval& addEq(Eval& a, const Eval& b)
{
for (int i = 0; i < length_; ++i)
a.data_[i] += b.data_[i];
return a;
}
static inline Eval& subEq(Eval& a, const Eval& b)
{
for (int i = 0; i < length_; ++i)
a.data_[i] -= b.data_[i];
return a;
}
static inline Eval& mulEq(Eval& a, const Eval& b)
{
// while the values are multiplied, the derivatives follow the product rule,
// i.e., (u*v)' = (v'u + u'v).
const ValueT u = a.value();
const ValueT v = b.value();
a.data_[valuepos_] *= v ;
for (int idx = dstart_; idx < dend_; ++idx)
a.data_[idx] = a.data_[idx] * v + b.data_[idx] * u;
return a;
}
template <class RhsValueType>
static inline Eval& scalarMulEq(Eval& a, const RhsValueType& other)
{
// values and derivatives are multiplied
for (int idx = 0 ; idx < length_; ++ idx)
a.data_[idx] *= other;
return a;
}
static inline Eval& divEq(Eval& a, const Eval& b)
{
// values are divided, derivatives follow the rule for division, i.e., (u/v)' = (v'u - u'v)/v^2.
const ValueT v_vv = 1.0 / b.value();
const ValueT u_vv = a.value() * v_vv * v_vv;
a.data_[valuepos_] *= v_vv;
for (int idx = dstart_; idx < dend_; ++idx)
a.data_[idx] = a.data_[idx] * v_vv - b.data_[idx] * u_vv;
return a;
}
template <class RhsValueType>
static inline Eval& scalarDivEq(Eval& a, const RhsValueType& other)
{
// values and derivatives are divided
for (int idx = 0 ; idx <= length_; ++ idx)
a.data_[idx] /= other;
return a;
}
// a/b with 'a' being a scalar and 'b' an Evaluation
template <class RhsValueType>
static inline Eval divide(const RhsValueType& a, const Eval& b )
{
Eval result;
result.setValue( a/b.value() );
const ValueT df_dg = - result.value()/b.value();
for (int idx = dstart_; idx < dend_; ++idx)
result.data_[idx] = df_dg*b.data_[idx];
return result;
}
};
/*!
* \brief Represents a function evaluation and its derivatives w.r.t. a fixed set of
* variables.
*/
template <class ValueT, int numVars>
template <class ValueT, int numDerivs>
class Evaluation
{
typedef EvaluationOps<ValueT, numVars> Ops;
friend Ops;
public:
//! field type
typedef ValueT ValueType;
//! number of derivatives
static constexpr int size = numVars;
static constexpr int size = numDerivs;
protected:
//! length of internal data vector
static constexpr int length_ = numVars + 1 ;
static constexpr int length_ = size + 1;
//! position index for value
static constexpr int valuepos_ = 0;
//! start index for derivatives
static constexpr int dstart_ = 1;
static constexpr int dstart_ = 1;
//! end+1 index for derivatives
static constexpr int dend_ = length_ ;
static constexpr int dend_ = length_;
public:
//! default constructor
@ -201,9 +89,10 @@ public:
//! copy other function evaluation
Evaluation(const Evaluation& other)
{
Ops::assign(*this, other);
}
: data_(other.data_)
{ }
// create an evaluation which represents a constant function
//
@ -227,7 +116,7 @@ public:
setValue( c );
clearDerivatives();
// The variable position must be in represented by the given variable descriptor
assert(0 <= varPos && varPos < numVars);
assert(0 <= varPos && varPos < size);
data_[varPos + dstart_] = 1.0;
Valgrind::CheckDefined(data_);
@ -236,7 +125,10 @@ public:
// set all derivatives to zero
void clearDerivatives()
{
Ops::clearDerivatives(*this);
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)
@ -262,21 +154,28 @@ public:
// print value
os << "v: " << value() << " / d:";
// print derivatives
for (int varIdx = 0; varIdx < numVars; ++varIdx)
for (int varIdx = 0; varIdx < size; ++varIdx)
os << " " << derivative(varIdx);
}
// copy all derivatives from other
void copyDerivatives(const Evaluation& other)
{
for (int varIdx = dstart_; varIdx < dend_; ++varIdx)
data_[ varIdx ] = other.data_[ varIdx ];
for (int i = dstart_; 0 < dend_; ++i)
data_[i] = other.data_[i];
}
// add value and derivatives from other to this values and derivatives
Evaluation& operator+=(const Evaluation& other)
{ return Ops::addEq(*this, other); }
{
for (int i = 0; i < length_; ++i)
data_[i] += other.data_[i];
return *this;
}
// add value from other to this values
template <class RhsValueType>
@ -289,7 +188,13 @@ public:
// subtract other's value and derivatives from this values
Evaluation& operator-=(const Evaluation& other)
{ return Ops::subEq(*this, other); }
{
for (int i = 0; i < length_; ++i)
data_[i] -= other.data_[i];
return *this;
}
// subtract other's value from this values
template <class RhsValueType>
@ -303,21 +208,83 @@ public:
// multiply values and apply chain rule to derivatives: (u*v)' = (v'u + u'v)
Evaluation& operator*=(const Evaluation& other)
{ return Ops::mulEq(*this, other); }
{
// while the values are multiplied, the derivatives follow the product rule,
// i.e., (u*v)' = (v'u + u'v).
const ValueT u = this->value();
const ValueT v = other.value();
// m(u*v)' = (v'u + u'v)
// value
this->data_[valuepos_] *= v ;
// derivatives
for (int i = 0; i < length_; ++i)
this->data_[i] = this->data_[i]*v + other.data_[i] * u;
return *this;
}
// m(c*u)' = c*u'
template <class RhsValueType>
Evaluation& operator*=(const RhsValueType& other)
{ return Ops::scalarMulEq(*this, other); }
{
for (int i = 0; i < length_; ++i)
data_[i] *= other;
return *this;
}
// m(u*v)' = (v'u + u'v)
Evaluation& operator/=(const Evaluation& other)
{ return Ops::divEq(*this, other); }
{
// values are divided, derivatives follow the rule for division, i.e., (u/v)' = (v'u - u'v)/v^2.
const ValueT v_vv = 1.0 / other.value();
const ValueT u_vv = value() * v_vv * v_vv;
// multiply value and derivatives by value of other
// value
data_[valuepos_] *= v_vv;
// derivatives
for (int i = 0; i < length_; ++i)
data_[i] = data_[i]*v_vv - other.data_[i]*u_vv;
return *this;
}
// divide value and derivatives by value of other
template <class RhsValueType>
Evaluation& operator/=(const RhsValueType& other)
{ return Ops::scalarDivEq(*this, other); }
{
ValueType tmp = 1.0/other;
for (int i = 0; i < length_; ++i)
data_[i] *= tmp;
return *this;
}
// division of a constant by an Evaluation
template <class RhsValueType>
static inline Evaluation divide(const RhsValueType& a, const Evaluation& b)
{
Evaluation result;
ValueType tmp = 1.0/b.value();
result.setValue( a*tmp );
const ValueT df_dg = - result.value()*tmp;
for (int i = dstart_; i < dend_; ++i)
result.data_[i] = df_dg*b.data_[i];
return result;
}
// add two evaluation objects
Evaluation operator+(const Evaluation& other) const
@ -356,7 +323,10 @@ public:
{
Evaluation result;
// set value and derivatives to negative
Ops::assignNegative(result, *this);
for (int i = 0; i < length_; ++i)
result.data_[i] = - data_[i];
return result;
}
@ -402,7 +372,11 @@ public:
// copy assignment from evaluation
Evaluation& operator=(const Evaluation& other)
{
Ops::assign(*this, other);
for (int i = 0; i < length_; ++i)
data_[i] = other.data_[i];
return *this;
}
@ -462,21 +436,23 @@ public:
// return varIdx'th derivative
const ValueType& derivative(int varIdx) const
{
assert(varIdx < numVars);
return data_[varIdx + dstart_];
assert(0 <= varIdx && varIdx < size);
return data_[dstart_ + varIdx];
}
// set derivative at position varIdx
void setDerivative(int varIdx, const ValueType& derVal)
{
assert(varIdx < numVars);
data_[varIdx + dstart_] = derVal;
assert(0 <= varIdx && varIdx < size);
data_[dstart_ + varIdx] = derVal;
}
protected:
typename Ops::DataVector data_;
private:
std::array<ValueT, length_> data_;
};
template <class RhsValueType, class ValueType, int numVars>
bool operator<(const RhsValueType& a, const Evaluation<ValueType, numVars>& b)
{ return b > a; }
@ -516,9 +492,7 @@ Evaluation<ValueType, numVars> operator-(const RhsValueType& a, const Evaluation
template <class RhsValueType, class ValueType, int numVars>
Evaluation<ValueType, numVars> operator/(const RhsValueType& a, const Evaluation<ValueType, numVars>& b)
{
typedef EvaluationOps<ValueType, numVars> Ops;
return Ops::divide(a, b);
return Evaluation<ValueType, numVars>::divide(a, b);
}
template <class RhsValueType, class ValueType, int numVars>
@ -536,8 +510,8 @@ std::ostream& operator<<(std::ostream& os, const Evaluation<ValueType, numVars>&
return os;
}
} // namespace DenseAd
} // namespace Opm
} } // namespace DenseAd, Opm
// In Dune 2.3, the Evaluation.hpp header must be included before the fmatrix.hh
// header. Dune 2.4+ does not suffer from this because of some c++-foo.
@ -601,6 +575,6 @@ public:
} // namespace Dune
#include <opm/material/densead/EvaluationSpecializations.hpp>
#include "EvaluationSpecializations.hpp"
#endif
#endif // OPM_DENSEAD_EVALUATION_HPP

View File

@ -23,128 +23,441 @@
/*!
* \file
*
* \brief This file specializes the dense-AD Evaluation class for 1 derivatives.
*
* \attention THIS FILE GETS AUTOMATICALLY GENERATED BY THE "genEvalSpecializations.py"
* SCRIPT. DO NOT EDIT IT MANUALLY!
*/
#ifndef OPM_DENSEAD_EVALUATION1_HPP
#define OPM_DENSEAD_EVALUATION1_HPP
#include "Math.hpp"
#include <opm/common/Valgrind.hpp>
#include <dune/common/version.hh>
#include <array>
#include <cmath>
#include <cassert>
#include <cstring>
#include <iostream>
#include <algorithm>
namespace Opm {
namespace DenseAd {
template <class ValueT>
struct EvaluationOps<ValueT, 1>
class Evaluation<ValueT, 1>
{
private:
typedef Evaluation<ValueT, 1 > Eval;
public:
//! field type
typedef ValueT ValueType;
//! number of derivatives
static constexpr int size = 1;
protected:
//! length of internal data vector
static constexpr int length_ = size + 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:
typedef std::array<ValueT, 2 > DataVector;
static inline void assign(Eval& a, const Eval& b)
//! default constructor
Evaluation() : data_()
{}
//! copy other function evaluation
Evaluation(const Evaluation& other)
{
a.data_[0] = b.data_[0];
a.data_[1] = b.data_[1];
data_[0] = other.data_[0];
data_[1] = other.data_[1];
}
static inline void assignNegative(Eval& a, const Eval& b)
// 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 < size);
data_[varPos + dstart_] = 1.0;
Valgrind::CheckDefined(data_);
}
// set all derivatives to zero
void clearDerivatives()
{
a.data_[0] = - b.data_[0];
a.data_[1] = - b.data_[1];
data_[1] = 0.0;
}
static inline void clearDerivatives(Eval& a)
// 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 < size; ++varIdx)
os << " " << derivative(varIdx);
}
// copy all derivatives from other
void copyDerivatives(const Evaluation& other)
{
for (int i = dstart_; 0 < dend_; ++i)
data_[i] = other.data_[i];
}
// add value and derivatives from other to this values and derivatives
Evaluation& operator+=(const Evaluation& other)
{
a.data_[1] = 0.0;
data_[0] += other.data_[0];
data_[1] += other.data_[1];
return *this;
}
static inline Eval& addEq(Eval& a, const Eval& b)
// 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)
{
a.data_[0] += b.data_[0];
a.data_[1] += b.data_[1];
return a;
data_[0] -= other.data_[0];
data_[1] -= other.data_[1];
return *this;
}
static inline Eval& subEq(Eval& a, const Eval& b)
// 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;
a.data_[0] -= b.data_[0];
a.data_[1] -= b.data_[1];
return a;
return *this;
}
static inline Eval& mulEq(Eval& a, const Eval& b)
// 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 ValueT u = a.value();
const ValueT v = b.value();
const ValueT u = this->value();
const ValueT v = other.value();
// value
a.data_[0] *= v ;
this->data_[valuepos_] *= v ;
// derivatives
a.data_[1] = a.data_[1]*v + b.data_[1] * u;
return a;
this->data_[1] = this->data_[1]*v + other.data_[1] * u;
return *this;
}
template <class RhsType>
static inline Eval& scalarMulEq(Eval& a, const RhsType& other)
// m(c*u)' = c*u'
template <class RhsValueType>
Evaluation& operator*=(const RhsValueType& other)
{
a.data_[0] *= other;
a.data_[1] *= other;
return a;
data_[0] *= other;
data_[1] *= other;
return *this;
}
static inline Eval& divEq(Eval& a, const Eval& b)
// 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 ValueT v_vv = 1.0 / b.value();
const ValueT u_vv = a.value() * v_vv * v_vv;
const ValueT v_vv = 1.0 / other.value();
const ValueT u_vv = value() * v_vv * v_vv;
// value
a.data_[0] *= v_vv;
data_[valuepos_] *= v_vv;
// derivatives
a.data_[1] = a.data_[1]*v_vv - b.data_[1]*u_vv;
return a;
}
template <class RhsType>
static inline Eval& scalarDivEq(Eval& a, const RhsType& other)
{
a.data_[0] /= other;
a.data_[1] /= other;
return a;
data_[1] = data_[1]*v_vv - other.data_[1]*u_vv;
return *this;
}
// divide value and derivatives by value of other
template <class RhsValueType>
static inline Eval divide(const RhsValueType& a, const Eval& b)
Evaluation& operator/=(const RhsValueType& other)
{
Eval result;
result.setValue( a/b.value() );
const ValueT df_dg = - result.value()/b.value();
ValueType tmp = 1.0/other;
data_[0] *= tmp;
data_[1] *= tmp;
return *this;
}
// division of a constant by an Evaluation
template <class RhsValueType>
static inline Evaluation divide(const RhsValueType& a, const Evaluation& b)
{
Evaluation result;
ValueType tmp = 1.0/b.value();
result.setValue( a*tmp );
const ValueT df_dg = - result.value()*tmp;
result.data_[1] = df_dg*b.data_[1];
return result;
}
// 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);
return (result -= other);
}
// subtract constant from evaluation object
template <class RhsValueType>
Evaluation operator-(const RhsValueType& other) const
{
Evaluation result(*this);
return (result -= other);
}
// negation (unary minus) operator
Evaluation operator-() const
{
Evaluation result;
// set value and derivatives to negative
result.data_[0] = - data_[0];
result.data_[1] = - data_[1];
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_[0] = other.data_[0];
data_[1] = other.data_[1];
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
template <class RhsValueType>
void setValue(const RhsValueType& val)
{ data_[valuepos_] = val; }
// return varIdx'th derivative
const ValueType& derivative(int varIdx) const
{
assert(0 <= varIdx && varIdx < size);
return data_[dstart_ + varIdx];
}
// set derivative at position varIdx
void setDerivative(int varIdx, const ValueType& derVal)
{
assert(0 <= varIdx && varIdx < size);
data_[dstart_ + varIdx] = derVal;
}
private:
std::array<ValueT, length_> data_;
};
} } // namespace DenseAd, Opm
#endif // OPM_DENSEAD_EVALUATION1_HPP
#endif // OPM_DENSEAD_EVALUATION1_HPP

View File

@ -23,202 +23,337 @@
/*!
* \file
*
* \brief This file specializes the dense-AD Evaluation class for 10 derivatives.
*
* \attention THIS FILE GETS AUTOMATICALLY GENERATED BY THE "genEvalSpecializations.py"
* SCRIPT. DO NOT EDIT IT MANUALLY!
*/
#ifndef OPM_DENSEAD_EVALUATION10_HPP
#define OPM_DENSEAD_EVALUATION10_HPP
#include "Math.hpp"
#include <opm/common/Valgrind.hpp>
#include <dune/common/version.hh>
#include <array>
#include <cmath>
#include <cassert>
#include <cstring>
#include <iostream>
#include <algorithm>
namespace Opm {
namespace DenseAd {
template <class ValueT>
struct EvaluationOps<ValueT, 10>
class Evaluation<ValueT, 10>
{
private:
typedef Evaluation<ValueT, 10 > Eval;
public:
//! field type
typedef ValueT ValueType;
//! number of derivatives
static constexpr int size = 10;
protected:
//! length of internal data vector
static constexpr int length_ = size + 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:
typedef std::array<ValueT, 11 > DataVector;
static inline void assign(Eval& a, const Eval& b)
//! default constructor
Evaluation() : data_()
{}
//! copy other function evaluation
Evaluation(const Evaluation& other)
{
a.data_[0] = b.data_[0];
a.data_[1] = b.data_[1];
a.data_[2] = b.data_[2];
a.data_[3] = b.data_[3];
a.data_[4] = b.data_[4];
a.data_[5] = b.data_[5];
a.data_[6] = b.data_[6];
a.data_[7] = b.data_[7];
a.data_[8] = b.data_[8];
a.data_[9] = b.data_[9];
a.data_[10] = b.data_[10];
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];
}
static inline void assignNegative(Eval& a, const Eval& b)
// 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 < size);
data_[varPos + dstart_] = 1.0;
Valgrind::CheckDefined(data_);
}
// set all derivatives to zero
void clearDerivatives()
{
a.data_[0] = - b.data_[0];
a.data_[1] = - b.data_[1];
a.data_[2] = - b.data_[2];
a.data_[3] = - b.data_[3];
a.data_[4] = - b.data_[4];
a.data_[5] = - b.data_[5];
a.data_[6] = - b.data_[6];
a.data_[7] = - b.data_[7];
a.data_[8] = - b.data_[8];
a.data_[9] = - b.data_[9];
a.data_[10] = - b.data_[10];
data_[1] = 0.0;
data_[2] = 0.0;
data_[3] = 0.0;
data_[4] = 0.0;
data_[5] = 0.0;
data_[6] = 0.0;
data_[7] = 0.0;
data_[8] = 0.0;
data_[9] = 0.0;
data_[10] = 0.0;
}
static inline void clearDerivatives(Eval& a)
// 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 < size; ++varIdx)
os << " " << derivative(varIdx);
}
// copy all derivatives from other
void copyDerivatives(const Evaluation& other)
{
for (int i = dstart_; 0 < dend_; ++i)
data_[i] = other.data_[i];
}
// add value and derivatives from other to this values and derivatives
Evaluation& operator+=(const Evaluation& other)
{
a.data_[1] = 0.0;
a.data_[2] = 0.0;
a.data_[3] = 0.0;
a.data_[4] = 0.0;
a.data_[5] = 0.0;
a.data_[6] = 0.0;
a.data_[7] = 0.0;
a.data_[8] = 0.0;
a.data_[9] = 0.0;
a.data_[10] = 0.0;
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];
return *this;
}
static inline Eval& addEq(Eval& a, const Eval& b)
// 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)
{
a.data_[0] += b.data_[0];
a.data_[1] += b.data_[1];
a.data_[2] += b.data_[2];
a.data_[3] += b.data_[3];
a.data_[4] += b.data_[4];
a.data_[5] += b.data_[5];
a.data_[6] += b.data_[6];
a.data_[7] += b.data_[7];
a.data_[8] += b.data_[8];
a.data_[9] += b.data_[9];
a.data_[10] += b.data_[10];
return a;
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];
return *this;
}
static inline Eval& subEq(Eval& a, const Eval& b)
// 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;
a.data_[0] -= b.data_[0];
a.data_[1] -= b.data_[1];
a.data_[2] -= b.data_[2];
a.data_[3] -= b.data_[3];
a.data_[4] -= b.data_[4];
a.data_[5] -= b.data_[5];
a.data_[6] -= b.data_[6];
a.data_[7] -= b.data_[7];
a.data_[8] -= b.data_[8];
a.data_[9] -= b.data_[9];
a.data_[10] -= b.data_[10];
return a;
return *this;
}
static inline Eval& mulEq(Eval& a, const Eval& b)
// 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 ValueT u = a.value();
const ValueT v = b.value();
const ValueT u = this->value();
const ValueT v = other.value();
// value
a.data_[0] *= v ;
this->data_[valuepos_] *= v ;
// derivatives
a.data_[1] = a.data_[1]*v + b.data_[1] * u;
a.data_[2] = a.data_[2]*v + b.data_[2] * u;
a.data_[3] = a.data_[3]*v + b.data_[3] * u;
a.data_[4] = a.data_[4]*v + b.data_[4] * u;
a.data_[5] = a.data_[5]*v + b.data_[5] * u;
a.data_[6] = a.data_[6]*v + b.data_[6] * u;
a.data_[7] = a.data_[7]*v + b.data_[7] * u;
a.data_[8] = a.data_[8]*v + b.data_[8] * u;
a.data_[9] = a.data_[9]*v + b.data_[9] * u;
a.data_[10] = a.data_[10]*v + b.data_[10] * u;
return a;
this->data_[1] = this->data_[1]*v + other.data_[1] * u;
this->data_[2] = this->data_[2]*v + other.data_[2] * u;
this->data_[3] = this->data_[3]*v + other.data_[3] * u;
this->data_[4] = this->data_[4]*v + other.data_[4] * u;
this->data_[5] = this->data_[5]*v + other.data_[5] * u;
this->data_[6] = this->data_[6]*v + other.data_[6] * u;
this->data_[7] = this->data_[7]*v + other.data_[7] * u;
this->data_[8] = this->data_[8]*v + other.data_[8] * u;
this->data_[9] = this->data_[9]*v + other.data_[9] * u;
this->data_[10] = this->data_[10]*v + other.data_[10] * u;
return *this;
}
template <class RhsType>
static inline Eval& scalarMulEq(Eval& a, const RhsType& other)
// m(c*u)' = c*u'
template <class RhsValueType>
Evaluation& operator*=(const RhsValueType& other)
{
a.data_[0] *= other;
a.data_[1] *= other;
a.data_[2] *= other;
a.data_[3] *= other;
a.data_[4] *= other;
a.data_[5] *= other;
a.data_[6] *= other;
a.data_[7] *= other;
a.data_[8] *= other;
a.data_[9] *= other;
a.data_[10] *= other;
return a;
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;
return *this;
}
static inline Eval& divEq(Eval& a, const Eval& b)
// 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 ValueT v_vv = 1.0 / b.value();
const ValueT u_vv = a.value() * v_vv * v_vv;
const ValueT v_vv = 1.0 / other.value();
const ValueT u_vv = value() * v_vv * v_vv;
// value
a.data_[0] *= v_vv;
data_[valuepos_] *= v_vv;
// derivatives
a.data_[1] = a.data_[1]*v_vv - b.data_[1]*u_vv;
a.data_[2] = a.data_[2]*v_vv - b.data_[2]*u_vv;
a.data_[3] = a.data_[3]*v_vv - b.data_[3]*u_vv;
a.data_[4] = a.data_[4]*v_vv - b.data_[4]*u_vv;
a.data_[5] = a.data_[5]*v_vv - b.data_[5]*u_vv;
a.data_[6] = a.data_[6]*v_vv - b.data_[6]*u_vv;
a.data_[7] = a.data_[7]*v_vv - b.data_[7]*u_vv;
a.data_[8] = a.data_[8]*v_vv - b.data_[8]*u_vv;
a.data_[9] = a.data_[9]*v_vv - b.data_[9]*u_vv;
a.data_[10] = a.data_[10]*v_vv - b.data_[10]*u_vv;
return a;
}
template <class RhsType>
static inline Eval& scalarDivEq(Eval& a, const RhsType& other)
{
a.data_[0] /= other;
a.data_[1] /= other;
a.data_[2] /= other;
a.data_[3] /= other;
a.data_[4] /= other;
a.data_[5] /= other;
a.data_[6] /= other;
a.data_[7] /= other;
a.data_[8] /= other;
a.data_[9] /= other;
a.data_[10] /= other;
return a;
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;
return *this;
}
// divide value and derivatives by value of other
template <class RhsValueType>
static inline Eval divide(const RhsValueType& a, const Eval& b)
Evaluation& operator/=(const RhsValueType& other)
{
Eval result;
result.setValue( a/b.value() );
const ValueT df_dg = - result.value()/b.value();
ValueType tmp = 1.0/other;
data_[0] *= tmp;
data_[1] *= tmp;
data_[2] *= tmp;
data_[3] *= tmp;
data_[4] *= tmp;
data_[5] *= tmp;
data_[6] *= tmp;
data_[7] *= tmp;
data_[8] *= tmp;
data_[9] *= tmp;
data_[10] *= tmp;
return *this;
}
// division of a constant by an Evaluation
template <class RhsValueType>
static inline Evaluation divide(const RhsValueType& a, const Evaluation& b)
{
Evaluation result;
ValueType tmp = 1.0/b.value();
result.setValue( a*tmp );
const ValueT df_dg = - result.value()*tmp;
result.data_[1] = df_dg*b.data_[1];
result.data_[2] = df_dg*b.data_[2];
@ -233,8 +368,195 @@ public:
return result;
}
// 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);
return (result -= other);
}
// subtract constant from evaluation object
template <class RhsValueType>
Evaluation operator-(const RhsValueType& other) const
{
Evaluation result(*this);
return (result -= other);
}
// negation (unary minus) operator
Evaluation operator-() const
{
Evaluation result;
// set value and derivatives to negative
result.data_[0] = - data_[0];
result.data_[1] = - data_[1];
result.data_[2] = - data_[2];
result.data_[3] = - data_[3];
result.data_[4] = - data_[4];
result.data_[5] = - data_[5];
result.data_[6] = - data_[6];
result.data_[7] = - data_[7];
result.data_[8] = - data_[8];
result.data_[9] = - data_[9];
result.data_[10] = - data_[10];
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_[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];
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
template <class RhsValueType>
void setValue(const RhsValueType& val)
{ data_[valuepos_] = val; }
// return varIdx'th derivative
const ValueType& derivative(int varIdx) const
{
assert(0 <= varIdx && varIdx < size);
return data_[dstart_ + varIdx];
}
// set derivative at position varIdx
void setDerivative(int varIdx, const ValueType& derVal)
{
assert(0 <= varIdx && varIdx < size);
data_[dstart_ + varIdx] = derVal;
}
private:
std::array<ValueT, length_> data_;
};
} } // namespace DenseAd, Opm
#endif // OPM_DENSEAD_EVALUATION10_HPP
#endif // OPM_DENSEAD_EVALUATION10_HPP

View File

@ -23,211 +23,345 @@
/*!
* \file
*
* \brief This file specializes the dense-AD Evaluation class for 11 derivatives.
*
* \attention THIS FILE GETS AUTOMATICALLY GENERATED BY THE "genEvalSpecializations.py"
* SCRIPT. DO NOT EDIT IT MANUALLY!
*/
#ifndef OPM_DENSEAD_EVALUATION11_HPP
#define OPM_DENSEAD_EVALUATION11_HPP
#include "Math.hpp"
#include <opm/common/Valgrind.hpp>
#include <dune/common/version.hh>
#include <array>
#include <cmath>
#include <cassert>
#include <cstring>
#include <iostream>
#include <algorithm>
namespace Opm {
namespace DenseAd {
template <class ValueT>
struct EvaluationOps<ValueT, 11>
class Evaluation<ValueT, 11>
{
private:
typedef Evaluation<ValueT, 11 > Eval;
public:
//! field type
typedef ValueT ValueType;
//! number of derivatives
static constexpr int size = 11;
protected:
//! length of internal data vector
static constexpr int length_ = size + 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:
typedef std::array<ValueT, 12 > DataVector;
static inline void assign(Eval& a, const Eval& b)
//! default constructor
Evaluation() : data_()
{}
//! copy other function evaluation
Evaluation(const Evaluation& other)
{
a.data_[0] = b.data_[0];
a.data_[1] = b.data_[1];
a.data_[2] = b.data_[2];
a.data_[3] = b.data_[3];
a.data_[4] = b.data_[4];
a.data_[5] = b.data_[5];
a.data_[6] = b.data_[6];
a.data_[7] = b.data_[7];
a.data_[8] = b.data_[8];
a.data_[9] = b.data_[9];
a.data_[10] = b.data_[10];
a.data_[11] = b.data_[11];
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];
}
static inline void assignNegative(Eval& a, const Eval& b)
// 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 < size);
data_[varPos + dstart_] = 1.0;
Valgrind::CheckDefined(data_);
}
// set all derivatives to zero
void clearDerivatives()
{
a.data_[0] = - b.data_[0];
a.data_[1] = - b.data_[1];
a.data_[2] = - b.data_[2];
a.data_[3] = - b.data_[3];
a.data_[4] = - b.data_[4];
a.data_[5] = - b.data_[5];
a.data_[6] = - b.data_[6];
a.data_[7] = - b.data_[7];
a.data_[8] = - b.data_[8];
a.data_[9] = - b.data_[9];
a.data_[10] = - b.data_[10];
a.data_[11] = - b.data_[11];
data_[1] = 0.0;
data_[2] = 0.0;
data_[3] = 0.0;
data_[4] = 0.0;
data_[5] = 0.0;
data_[6] = 0.0;
data_[7] = 0.0;
data_[8] = 0.0;
data_[9] = 0.0;
data_[10] = 0.0;
data_[11] = 0.0;
}
static inline void clearDerivatives(Eval& a)
// 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 < size; ++varIdx)
os << " " << derivative(varIdx);
}
// copy all derivatives from other
void copyDerivatives(const Evaluation& other)
{
for (int i = dstart_; 0 < dend_; ++i)
data_[i] = other.data_[i];
}
// add value and derivatives from other to this values and derivatives
Evaluation& operator+=(const Evaluation& other)
{
a.data_[1] = 0.0;
a.data_[2] = 0.0;
a.data_[3] = 0.0;
a.data_[4] = 0.0;
a.data_[5] = 0.0;
a.data_[6] = 0.0;
a.data_[7] = 0.0;
a.data_[8] = 0.0;
a.data_[9] = 0.0;
a.data_[10] = 0.0;
a.data_[11] = 0.0;
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];
return *this;
}
static inline Eval& addEq(Eval& a, const Eval& b)
// 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)
{
a.data_[0] += b.data_[0];
a.data_[1] += b.data_[1];
a.data_[2] += b.data_[2];
a.data_[3] += b.data_[3];
a.data_[4] += b.data_[4];
a.data_[5] += b.data_[5];
a.data_[6] += b.data_[6];
a.data_[7] += b.data_[7];
a.data_[8] += b.data_[8];
a.data_[9] += b.data_[9];
a.data_[10] += b.data_[10];
a.data_[11] += b.data_[11];
return a;
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];
return *this;
}
static inline Eval& subEq(Eval& a, const Eval& b)
// 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;
a.data_[0] -= b.data_[0];
a.data_[1] -= b.data_[1];
a.data_[2] -= b.data_[2];
a.data_[3] -= b.data_[3];
a.data_[4] -= b.data_[4];
a.data_[5] -= b.data_[5];
a.data_[6] -= b.data_[6];
a.data_[7] -= b.data_[7];
a.data_[8] -= b.data_[8];
a.data_[9] -= b.data_[9];
a.data_[10] -= b.data_[10];
a.data_[11] -= b.data_[11];
return a;
return *this;
}
static inline Eval& mulEq(Eval& a, const Eval& b)
// 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 ValueT u = a.value();
const ValueT v = b.value();
const ValueT u = this->value();
const ValueT v = other.value();
// value
a.data_[0] *= v ;
this->data_[valuepos_] *= v ;
// derivatives
a.data_[1] = a.data_[1]*v + b.data_[1] * u;
a.data_[2] = a.data_[2]*v + b.data_[2] * u;
a.data_[3] = a.data_[3]*v + b.data_[3] * u;
a.data_[4] = a.data_[4]*v + b.data_[4] * u;
a.data_[5] = a.data_[5]*v + b.data_[5] * u;
a.data_[6] = a.data_[6]*v + b.data_[6] * u;
a.data_[7] = a.data_[7]*v + b.data_[7] * u;
a.data_[8] = a.data_[8]*v + b.data_[8] * u;
a.data_[9] = a.data_[9]*v + b.data_[9] * u;
a.data_[10] = a.data_[10]*v + b.data_[10] * u;
a.data_[11] = a.data_[11]*v + b.data_[11] * u;
return a;
this->data_[1] = this->data_[1]*v + other.data_[1] * u;
this->data_[2] = this->data_[2]*v + other.data_[2] * u;
this->data_[3] = this->data_[3]*v + other.data_[3] * u;
this->data_[4] = this->data_[4]*v + other.data_[4] * u;
this->data_[5] = this->data_[5]*v + other.data_[5] * u;
this->data_[6] = this->data_[6]*v + other.data_[6] * u;
this->data_[7] = this->data_[7]*v + other.data_[7] * u;
this->data_[8] = this->data_[8]*v + other.data_[8] * u;
this->data_[9] = this->data_[9]*v + other.data_[9] * u;
this->data_[10] = this->data_[10]*v + other.data_[10] * u;
this->data_[11] = this->data_[11]*v + other.data_[11] * u;
return *this;
}
template <class RhsType>
static inline Eval& scalarMulEq(Eval& a, const RhsType& other)
// m(c*u)' = c*u'
template <class RhsValueType>
Evaluation& operator*=(const RhsValueType& other)
{
a.data_[0] *= other;
a.data_[1] *= other;
a.data_[2] *= other;
a.data_[3] *= other;
a.data_[4] *= other;
a.data_[5] *= other;
a.data_[6] *= other;
a.data_[7] *= other;
a.data_[8] *= other;
a.data_[9] *= other;
a.data_[10] *= other;
a.data_[11] *= other;
return a;
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;
return *this;
}
static inline Eval& divEq(Eval& a, const Eval& b)
// 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 ValueT v_vv = 1.0 / b.value();
const ValueT u_vv = a.value() * v_vv * v_vv;
const ValueT v_vv = 1.0 / other.value();
const ValueT u_vv = value() * v_vv * v_vv;
// value
a.data_[0] *= v_vv;
data_[valuepos_] *= v_vv;
// derivatives
a.data_[1] = a.data_[1]*v_vv - b.data_[1]*u_vv;
a.data_[2] = a.data_[2]*v_vv - b.data_[2]*u_vv;
a.data_[3] = a.data_[3]*v_vv - b.data_[3]*u_vv;
a.data_[4] = a.data_[4]*v_vv - b.data_[4]*u_vv;
a.data_[5] = a.data_[5]*v_vv - b.data_[5]*u_vv;
a.data_[6] = a.data_[6]*v_vv - b.data_[6]*u_vv;
a.data_[7] = a.data_[7]*v_vv - b.data_[7]*u_vv;
a.data_[8] = a.data_[8]*v_vv - b.data_[8]*u_vv;
a.data_[9] = a.data_[9]*v_vv - b.data_[9]*u_vv;
a.data_[10] = a.data_[10]*v_vv - b.data_[10]*u_vv;
a.data_[11] = a.data_[11]*v_vv - b.data_[11]*u_vv;
return a;
}
template <class RhsType>
static inline Eval& scalarDivEq(Eval& a, const RhsType& other)
{
a.data_[0] /= other;
a.data_[1] /= other;
a.data_[2] /= other;
a.data_[3] /= other;
a.data_[4] /= other;
a.data_[5] /= other;
a.data_[6] /= other;
a.data_[7] /= other;
a.data_[8] /= other;
a.data_[9] /= other;
a.data_[10] /= other;
a.data_[11] /= other;
return a;
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;
return *this;
}
// divide value and derivatives by value of other
template <class RhsValueType>
static inline Eval divide(const RhsValueType& a, const Eval& b)
Evaluation& operator/=(const RhsValueType& other)
{
Eval result;
result.setValue( a/b.value() );
const ValueT df_dg = - result.value()/b.value();
ValueType tmp = 1.0/other;
data_[0] *= tmp;
data_[1] *= tmp;
data_[2] *= tmp;
data_[3] *= tmp;
data_[4] *= tmp;
data_[5] *= tmp;
data_[6] *= tmp;
data_[7] *= tmp;
data_[8] *= tmp;
data_[9] *= tmp;
data_[10] *= tmp;
data_[11] *= tmp;
return *this;
}
// division of a constant by an Evaluation
template <class RhsValueType>
static inline Evaluation divide(const RhsValueType& a, const Evaluation& b)
{
Evaluation result;
ValueType tmp = 1.0/b.value();
result.setValue( a*tmp );
const ValueT df_dg = - result.value()*tmp;
result.data_[1] = df_dg*b.data_[1];
result.data_[2] = df_dg*b.data_[2];
@ -243,8 +377,197 @@ public:
return result;
}
// 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);
return (result -= other);
}
// subtract constant from evaluation object
template <class RhsValueType>
Evaluation operator-(const RhsValueType& other) const
{
Evaluation result(*this);
return (result -= other);
}
// negation (unary minus) operator
Evaluation operator-() const
{
Evaluation result;
// set value and derivatives to negative
result.data_[0] = - data_[0];
result.data_[1] = - data_[1];
result.data_[2] = - data_[2];
result.data_[3] = - data_[3];
result.data_[4] = - data_[4];
result.data_[5] = - data_[5];
result.data_[6] = - data_[6];
result.data_[7] = - data_[7];
result.data_[8] = - data_[8];
result.data_[9] = - data_[9];
result.data_[10] = - data_[10];
result.data_[11] = - data_[11];
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_[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];
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
template <class RhsValueType>
void setValue(const RhsValueType& val)
{ data_[valuepos_] = val; }
// return varIdx'th derivative
const ValueType& derivative(int varIdx) const
{
assert(0 <= varIdx && varIdx < size);
return data_[dstart_ + varIdx];
}
// set derivative at position varIdx
void setDerivative(int varIdx, const ValueType& derVal)
{
assert(0 <= varIdx && varIdx < size);
data_[dstart_ + varIdx] = derVal;
}
private:
std::array<ValueT, length_> data_;
};
} } // namespace DenseAd, Opm
#endif // OPM_DENSEAD_EVALUATION11_HPP
#endif // OPM_DENSEAD_EVALUATION11_HPP

View File

@ -23,220 +23,353 @@
/*!
* \file
*
* \brief This file specializes the dense-AD Evaluation class for 12 derivatives.
*
* \attention THIS FILE GETS AUTOMATICALLY GENERATED BY THE "genEvalSpecializations.py"
* SCRIPT. DO NOT EDIT IT MANUALLY!
*/
#ifndef OPM_DENSEAD_EVALUATION12_HPP
#define OPM_DENSEAD_EVALUATION12_HPP
#include "Math.hpp"
#include <opm/common/Valgrind.hpp>
#include <dune/common/version.hh>
#include <array>
#include <cmath>
#include <cassert>
#include <cstring>
#include <iostream>
#include <algorithm>
namespace Opm {
namespace DenseAd {
template <class ValueT>
struct EvaluationOps<ValueT, 12>
class Evaluation<ValueT, 12>
{
private:
typedef Evaluation<ValueT, 12 > Eval;
public:
//! field type
typedef ValueT ValueType;
//! number of derivatives
static constexpr int size = 12;
protected:
//! length of internal data vector
static constexpr int length_ = size + 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:
typedef std::array<ValueT, 13 > DataVector;
static inline void assign(Eval& a, const Eval& b)
//! default constructor
Evaluation() : data_()
{}
//! copy other function evaluation
Evaluation(const Evaluation& other)
{
a.data_[0] = b.data_[0];
a.data_[1] = b.data_[1];
a.data_[2] = b.data_[2];
a.data_[3] = b.data_[3];
a.data_[4] = b.data_[4];
a.data_[5] = b.data_[5];
a.data_[6] = b.data_[6];
a.data_[7] = b.data_[7];
a.data_[8] = b.data_[8];
a.data_[9] = b.data_[9];
a.data_[10] = b.data_[10];
a.data_[11] = b.data_[11];
a.data_[12] = b.data_[12];
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];
}
static inline void assignNegative(Eval& a, const Eval& b)
// 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 < size);
data_[varPos + dstart_] = 1.0;
Valgrind::CheckDefined(data_);
}
// set all derivatives to zero
void clearDerivatives()
{
a.data_[0] = - b.data_[0];
a.data_[1] = - b.data_[1];
a.data_[2] = - b.data_[2];
a.data_[3] = - b.data_[3];
a.data_[4] = - b.data_[4];
a.data_[5] = - b.data_[5];
a.data_[6] = - b.data_[6];
a.data_[7] = - b.data_[7];
a.data_[8] = - b.data_[8];
a.data_[9] = - b.data_[9];
a.data_[10] = - b.data_[10];
a.data_[11] = - b.data_[11];
a.data_[12] = - b.data_[12];
data_[1] = 0.0;
data_[2] = 0.0;
data_[3] = 0.0;
data_[4] = 0.0;
data_[5] = 0.0;
data_[6] = 0.0;
data_[7] = 0.0;
data_[8] = 0.0;
data_[9] = 0.0;
data_[10] = 0.0;
data_[11] = 0.0;
data_[12] = 0.0;
}
static inline void clearDerivatives(Eval& a)
// 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 < size; ++varIdx)
os << " " << derivative(varIdx);
}
// copy all derivatives from other
void copyDerivatives(const Evaluation& other)
{
for (int i = dstart_; 0 < dend_; ++i)
data_[i] = other.data_[i];
}
// add value and derivatives from other to this values and derivatives
Evaluation& operator+=(const Evaluation& other)
{
a.data_[1] = 0.0;
a.data_[2] = 0.0;
a.data_[3] = 0.0;
a.data_[4] = 0.0;
a.data_[5] = 0.0;
a.data_[6] = 0.0;
a.data_[7] = 0.0;
a.data_[8] = 0.0;
a.data_[9] = 0.0;
a.data_[10] = 0.0;
a.data_[11] = 0.0;
a.data_[12] = 0.0;
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;
}
static inline Eval& addEq(Eval& a, const Eval& b)
// 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)
{
a.data_[0] += b.data_[0];
a.data_[1] += b.data_[1];
a.data_[2] += b.data_[2];
a.data_[3] += b.data_[3];
a.data_[4] += b.data_[4];
a.data_[5] += b.data_[5];
a.data_[6] += b.data_[6];
a.data_[7] += b.data_[7];
a.data_[8] += b.data_[8];
a.data_[9] += b.data_[9];
a.data_[10] += b.data_[10];
a.data_[11] += b.data_[11];
a.data_[12] += b.data_[12];
return a;
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;
}
static inline Eval& subEq(Eval& a, const Eval& b)
// 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;
a.data_[0] -= b.data_[0];
a.data_[1] -= b.data_[1];
a.data_[2] -= b.data_[2];
a.data_[3] -= b.data_[3];
a.data_[4] -= b.data_[4];
a.data_[5] -= b.data_[5];
a.data_[6] -= b.data_[6];
a.data_[7] -= b.data_[7];
a.data_[8] -= b.data_[8];
a.data_[9] -= b.data_[9];
a.data_[10] -= b.data_[10];
a.data_[11] -= b.data_[11];
a.data_[12] -= b.data_[12];
return a;
return *this;
}
static inline Eval& mulEq(Eval& a, const Eval& b)
// 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 ValueT u = a.value();
const ValueT v = b.value();
const ValueT u = this->value();
const ValueT v = other.value();
// value
a.data_[0] *= v ;
this->data_[valuepos_] *= v ;
// derivatives
a.data_[1] = a.data_[1]*v + b.data_[1] * u;
a.data_[2] = a.data_[2]*v + b.data_[2] * u;
a.data_[3] = a.data_[3]*v + b.data_[3] * u;
a.data_[4] = a.data_[4]*v + b.data_[4] * u;
a.data_[5] = a.data_[5]*v + b.data_[5] * u;
a.data_[6] = a.data_[6]*v + b.data_[6] * u;
a.data_[7] = a.data_[7]*v + b.data_[7] * u;
a.data_[8] = a.data_[8]*v + b.data_[8] * u;
a.data_[9] = a.data_[9]*v + b.data_[9] * u;
a.data_[10] = a.data_[10]*v + b.data_[10] * u;
a.data_[11] = a.data_[11]*v + b.data_[11] * u;
a.data_[12] = a.data_[12]*v + b.data_[12] * u;
return a;
this->data_[1] = this->data_[1]*v + other.data_[1] * u;
this->data_[2] = this->data_[2]*v + other.data_[2] * u;
this->data_[3] = this->data_[3]*v + other.data_[3] * u;
this->data_[4] = this->data_[4]*v + other.data_[4] * u;
this->data_[5] = this->data_[5]*v + other.data_[5] * u;
this->data_[6] = this->data_[6]*v + other.data_[6] * u;
this->data_[7] = this->data_[7]*v + other.data_[7] * u;
this->data_[8] = this->data_[8]*v + other.data_[8] * u;
this->data_[9] = this->data_[9]*v + other.data_[9] * u;
this->data_[10] = this->data_[10]*v + other.data_[10] * u;
this->data_[11] = this->data_[11]*v + other.data_[11] * u;
this->data_[12] = this->data_[12]*v + other.data_[12] * u;
return *this;
}
template <class RhsType>
static inline Eval& scalarMulEq(Eval& a, const RhsType& other)
// m(c*u)' = c*u'
template <class RhsValueType>
Evaluation& operator*=(const RhsValueType& other)
{
a.data_[0] *= other;
a.data_[1] *= other;
a.data_[2] *= other;
a.data_[3] *= other;
a.data_[4] *= other;
a.data_[5] *= other;
a.data_[6] *= other;
a.data_[7] *= other;
a.data_[8] *= other;
a.data_[9] *= other;
a.data_[10] *= other;
a.data_[11] *= other;
a.data_[12] *= other;
return a;
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;
}
static inline Eval& divEq(Eval& a, const Eval& b)
// 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 ValueT v_vv = 1.0 / b.value();
const ValueT u_vv = a.value() * v_vv * v_vv;
const ValueT v_vv = 1.0 / other.value();
const ValueT u_vv = value() * v_vv * v_vv;
// value
a.data_[0] *= v_vv;
data_[valuepos_] *= v_vv;
// derivatives
a.data_[1] = a.data_[1]*v_vv - b.data_[1]*u_vv;
a.data_[2] = a.data_[2]*v_vv - b.data_[2]*u_vv;
a.data_[3] = a.data_[3]*v_vv - b.data_[3]*u_vv;
a.data_[4] = a.data_[4]*v_vv - b.data_[4]*u_vv;
a.data_[5] = a.data_[5]*v_vv - b.data_[5]*u_vv;
a.data_[6] = a.data_[6]*v_vv - b.data_[6]*u_vv;
a.data_[7] = a.data_[7]*v_vv - b.data_[7]*u_vv;
a.data_[8] = a.data_[8]*v_vv - b.data_[8]*u_vv;
a.data_[9] = a.data_[9]*v_vv - b.data_[9]*u_vv;
a.data_[10] = a.data_[10]*v_vv - b.data_[10]*u_vv;
a.data_[11] = a.data_[11]*v_vv - b.data_[11]*u_vv;
a.data_[12] = a.data_[12]*v_vv - b.data_[12]*u_vv;
return a;
}
template <class RhsType>
static inline Eval& scalarDivEq(Eval& a, const RhsType& other)
{
a.data_[0] /= other;
a.data_[1] /= other;
a.data_[2] /= other;
a.data_[3] /= other;
a.data_[4] /= other;
a.data_[5] /= other;
a.data_[6] /= other;
a.data_[7] /= other;
a.data_[8] /= other;
a.data_[9] /= other;
a.data_[10] /= other;
a.data_[11] /= other;
a.data_[12] /= other;
return a;
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;
}
// divide value and derivatives by value of other
template <class RhsValueType>
static inline Eval divide(const RhsValueType& a, const Eval& b)
Evaluation& operator/=(const RhsValueType& other)
{
Eval result;
result.setValue( a/b.value() );
const ValueT df_dg = - result.value()/b.value();
ValueType tmp = 1.0/other;
data_[0] *= tmp;
data_[1] *= tmp;
data_[2] *= tmp;
data_[3] *= tmp;
data_[4] *= tmp;
data_[5] *= tmp;
data_[6] *= tmp;
data_[7] *= tmp;
data_[8] *= tmp;
data_[9] *= tmp;
data_[10] *= tmp;
data_[11] *= tmp;
data_[12] *= tmp;
return *this;
}
// division of a constant by an Evaluation
template <class RhsValueType>
static inline Evaluation divide(const RhsValueType& a, const Evaluation& b)
{
Evaluation result;
ValueType tmp = 1.0/b.value();
result.setValue( a*tmp );
const ValueT df_dg = - result.value()*tmp;
result.data_[1] = df_dg*b.data_[1];
result.data_[2] = df_dg*b.data_[2];
@ -253,8 +386,199 @@ public:
return result;
}
// 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);
return (result -= other);
}
// subtract constant from evaluation object
template <class RhsValueType>
Evaluation operator-(const RhsValueType& other) const
{
Evaluation result(*this);
return (result -= other);
}
// negation (unary minus) operator
Evaluation operator-() const
{
Evaluation result;
// set value and derivatives to negative
result.data_[0] = - data_[0];
result.data_[1] = - data_[1];
result.data_[2] = - data_[2];
result.data_[3] = - data_[3];
result.data_[4] = - data_[4];
result.data_[5] = - data_[5];
result.data_[6] = - data_[6];
result.data_[7] = - data_[7];
result.data_[8] = - data_[8];
result.data_[9] = - data_[9];
result.data_[10] = - data_[10];
result.data_[11] = - data_[11];
result.data_[12] = - data_[12];
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_[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;
}
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
template <class RhsValueType>
void setValue(const RhsValueType& val)
{ data_[valuepos_] = val; }
// return varIdx'th derivative
const ValueType& derivative(int varIdx) const
{
assert(0 <= varIdx && varIdx < size);
return data_[dstart_ + varIdx];
}
// set derivative at position varIdx
void setDerivative(int varIdx, const ValueType& derVal)
{
assert(0 <= varIdx && varIdx < size);
data_[dstart_ + varIdx] = derVal;
}
private:
std::array<ValueT, length_> data_;
};
} } // namespace DenseAd, Opm
#endif // OPM_DENSEAD_EVALUATION12_HPP
#endif // OPM_DENSEAD_EVALUATION12_HPP

View File

@ -23,138 +23,452 @@
/*!
* \file
*
* \brief This file specializes the dense-AD Evaluation class for 2 derivatives.
*
* \attention THIS FILE GETS AUTOMATICALLY GENERATED BY THE "genEvalSpecializations.py"
* SCRIPT. DO NOT EDIT IT MANUALLY!
*/
#ifndef OPM_DENSEAD_EVALUATION2_HPP
#define OPM_DENSEAD_EVALUATION2_HPP
#include "Math.hpp"
#include <opm/common/Valgrind.hpp>
#include <dune/common/version.hh>
#include <array>
#include <cmath>
#include <cassert>
#include <cstring>
#include <iostream>
#include <algorithm>
namespace Opm {
namespace DenseAd {
template <class ValueT>
struct EvaluationOps<ValueT, 2>
class Evaluation<ValueT, 2>
{
private:
typedef Evaluation<ValueT, 2 > Eval;
public:
//! field type
typedef ValueT ValueType;
//! number of derivatives
static constexpr int size = 2;
protected:
//! length of internal data vector
static constexpr int length_ = size + 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:
typedef std::array<ValueT, 3 > DataVector;
static inline void assign(Eval& a, const Eval& b)
//! default constructor
Evaluation() : data_()
{}
//! copy other function evaluation
Evaluation(const Evaluation& other)
{
a.data_[0] = b.data_[0];
a.data_[1] = b.data_[1];
a.data_[2] = b.data_[2];
data_[0] = other.data_[0];
data_[1] = other.data_[1];
data_[2] = other.data_[2];
}
static inline void assignNegative(Eval& a, const Eval& b)
// 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 < size);
data_[varPos + dstart_] = 1.0;
Valgrind::CheckDefined(data_);
}
// set all derivatives to zero
void clearDerivatives()
{
a.data_[0] = - b.data_[0];
a.data_[1] = - b.data_[1];
a.data_[2] = - b.data_[2];
data_[1] = 0.0;
data_[2] = 0.0;
}
static inline void clearDerivatives(Eval& a)
// 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 < size; ++varIdx)
os << " " << derivative(varIdx);
}
// copy all derivatives from other
void copyDerivatives(const Evaluation& other)
{
for (int i = dstart_; 0 < dend_; ++i)
data_[i] = other.data_[i];
}
// add value and derivatives from other to this values and derivatives
Evaluation& operator+=(const Evaluation& other)
{
a.data_[1] = 0.0;
a.data_[2] = 0.0;
data_[0] += other.data_[0];
data_[1] += other.data_[1];
data_[2] += other.data_[2];
return *this;
}
static inline Eval& addEq(Eval& a, const Eval& b)
// 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)
{
a.data_[0] += b.data_[0];
a.data_[1] += b.data_[1];
a.data_[2] += b.data_[2];
return a;
data_[0] -= other.data_[0];
data_[1] -= other.data_[1];
data_[2] -= other.data_[2];
return *this;
}
static inline Eval& subEq(Eval& a, const Eval& b)
// 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;
a.data_[0] -= b.data_[0];
a.data_[1] -= b.data_[1];
a.data_[2] -= b.data_[2];
return a;
return *this;
}
static inline Eval& mulEq(Eval& a, const Eval& b)
// 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 ValueT u = a.value();
const ValueT v = b.value();
const ValueT u = this->value();
const ValueT v = other.value();
// value
a.data_[0] *= v ;
this->data_[valuepos_] *= v ;
// derivatives
a.data_[1] = a.data_[1]*v + b.data_[1] * u;
a.data_[2] = a.data_[2]*v + b.data_[2] * u;
return a;
this->data_[1] = this->data_[1]*v + other.data_[1] * u;
this->data_[2] = this->data_[2]*v + other.data_[2] * u;
return *this;
}
template <class RhsType>
static inline Eval& scalarMulEq(Eval& a, const RhsType& other)
// m(c*u)' = c*u'
template <class RhsValueType>
Evaluation& operator*=(const RhsValueType& other)
{
a.data_[0] *= other;
a.data_[1] *= other;
a.data_[2] *= other;
return a;
data_[0] *= other;
data_[1] *= other;
data_[2] *= other;
return *this;
}
static inline Eval& divEq(Eval& a, const Eval& b)
// 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 ValueT v_vv = 1.0 / b.value();
const ValueT u_vv = a.value() * v_vv * v_vv;
const ValueT v_vv = 1.0 / other.value();
const ValueT u_vv = value() * v_vv * v_vv;
// value
a.data_[0] *= v_vv;
data_[valuepos_] *= v_vv;
// derivatives
a.data_[1] = a.data_[1]*v_vv - b.data_[1]*u_vv;
a.data_[2] = a.data_[2]*v_vv - b.data_[2]*u_vv;
return a;
}
template <class RhsType>
static inline Eval& scalarDivEq(Eval& a, const RhsType& other)
{
a.data_[0] /= other;
a.data_[1] /= other;
a.data_[2] /= other;
return a;
data_[1] = data_[1]*v_vv - other.data_[1]*u_vv;
data_[2] = data_[2]*v_vv - other.data_[2]*u_vv;
return *this;
}
// divide value and derivatives by value of other
template <class RhsValueType>
static inline Eval divide(const RhsValueType& a, const Eval& b)
Evaluation& operator/=(const RhsValueType& other)
{
Eval result;
result.setValue( a/b.value() );
const ValueT df_dg = - result.value()/b.value();
ValueType tmp = 1.0/other;
data_[0] *= tmp;
data_[1] *= tmp;
data_[2] *= tmp;
return *this;
}
// division of a constant by an Evaluation
template <class RhsValueType>
static inline Evaluation divide(const RhsValueType& a, const Evaluation& b)
{
Evaluation result;
ValueType tmp = 1.0/b.value();
result.setValue( a*tmp );
const ValueT df_dg = - result.value()*tmp;
result.data_[1] = df_dg*b.data_[1];
result.data_[2] = df_dg*b.data_[2];
return result;
}
// 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);
return (result -= other);
}
// subtract constant from evaluation object
template <class RhsValueType>
Evaluation operator-(const RhsValueType& other) const
{
Evaluation result(*this);
return (result -= other);
}
// negation (unary minus) operator
Evaluation operator-() const
{
Evaluation result;
// set value and derivatives to negative
result.data_[0] = - data_[0];
result.data_[1] = - data_[1];
result.data_[2] = - data_[2];
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_[0] = other.data_[0];
data_[1] = other.data_[1];
data_[2] = other.data_[2];
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
template <class RhsValueType>
void setValue(const RhsValueType& val)
{ data_[valuepos_] = val; }
// return varIdx'th derivative
const ValueType& derivative(int varIdx) const
{
assert(0 <= varIdx && varIdx < size);
return data_[dstart_ + varIdx];
}
// set derivative at position varIdx
void setDerivative(int varIdx, const ValueType& derVal)
{
assert(0 <= varIdx && varIdx < size);
data_[dstart_ + varIdx] = derVal;
}
private:
std::array<ValueT, length_> data_;
};
} } // namespace DenseAd, Opm
#endif // OPM_DENSEAD_EVALUATION2_HPP
#endif // OPM_DENSEAD_EVALUATION2_HPP

View File

@ -23,139 +23,281 @@
/*!
* \file
*
* \brief This file specializes the dense-AD Evaluation class for 3 derivatives.
*
* \attention THIS FILE GETS AUTOMATICALLY GENERATED BY THE "genEvalSpecializations.py"
* SCRIPT. DO NOT EDIT IT MANUALLY!
*/
#ifndef OPM_DENSEAD_EVALUATION3_HPP
#define OPM_DENSEAD_EVALUATION3_HPP
#include "Math.hpp"
#include <opm/common/Valgrind.hpp>
#include <dune/common/version.hh>
#include <array>
#include <cmath>
#include <cassert>
#include <cstring>
#include <iostream>
#include <algorithm>
namespace Opm {
namespace DenseAd {
template <class ValueT>
struct EvaluationOps<ValueT, 3>
class Evaluation<ValueT, 3>
{
private:
typedef Evaluation<ValueT, 3 > Eval;
public:
//! field type
typedef ValueT ValueType;
//! number of derivatives
static constexpr int size = 3;
protected:
//! length of internal data vector
static constexpr int length_ = size + 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:
typedef std::array<ValueT, 4 > DataVector;
static inline void assign(Eval& a, const Eval& b)
//! default constructor
Evaluation() : data_()
{}
//! copy other function evaluation
Evaluation(const Evaluation& other)
{
a.data_[0] = b.data_[0];
a.data_[1] = b.data_[1];
a.data_[2] = b.data_[2];
a.data_[3] = b.data_[3];
data_[0] = other.data_[0];
data_[1] = other.data_[1];
data_[2] = other.data_[2];
data_[3] = other.data_[3];
}
static inline void assignNegative(Eval& a, const Eval& b)
// 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 < size);
data_[varPos + dstart_] = 1.0;
Valgrind::CheckDefined(data_);
}
// set all derivatives to zero
void clearDerivatives()
{
a.data_[0] = - b.data_[0];
a.data_[1] = - b.data_[1];
a.data_[2] = - b.data_[2];
a.data_[3] = - b.data_[3];
data_[1] = 0.0;
data_[2] = 0.0;
data_[3] = 0.0;
}
static inline void clearDerivatives(Eval& a)
// 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 < size; ++varIdx)
os << " " << derivative(varIdx);
}
// copy all derivatives from other
void copyDerivatives(const Evaluation& other)
{
for (int i = dstart_; 0 < dend_; ++i)
data_[i] = other.data_[i];
}
// add value and derivatives from other to this values and derivatives
Evaluation& operator+=(const Evaluation& other)
{
a.data_[1] = 0.0;
a.data_[2] = 0.0;
a.data_[3] = 0.0;
data_[0] += other.data_[0];
data_[1] += other.data_[1];
data_[2] += other.data_[2];
data_[3] += other.data_[3];
return *this;
}
static inline Eval& addEq(Eval& a, const Eval& b)
// 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)
{
a.data_[0] += b.data_[0];
a.data_[1] += b.data_[1];
a.data_[2] += b.data_[2];
a.data_[3] += b.data_[3];
return a;
data_[0] -= other.data_[0];
data_[1] -= other.data_[1];
data_[2] -= other.data_[2];
data_[3] -= other.data_[3];
return *this;
}
static inline Eval& subEq(Eval& a, const Eval& b)
// 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;
a.data_[0] -= b.data_[0];
a.data_[1] -= b.data_[1];
a.data_[2] -= b.data_[2];
a.data_[3] -= b.data_[3];
return a;
return *this;
}
static inline Eval& mulEq(Eval& a, const Eval& b)
// 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 ValueT u = a.value();
const ValueT v = b.value();
const ValueT u = this->value();
const ValueT v = other.value();
// value
a.data_[0] *= v ;
this->data_[valuepos_] *= v ;
// derivatives
a.data_[1] = a.data_[1]*v + b.data_[1] * u;
a.data_[2] = a.data_[2]*v + b.data_[2] * u;
a.data_[3] = a.data_[3]*v + b.data_[3] * u;
return a;
this->data_[1] = this->data_[1]*v + other.data_[1] * u;
this->data_[2] = this->data_[2]*v + other.data_[2] * u;
this->data_[3] = this->data_[3]*v + other.data_[3] * u;
return *this;
}
template <class RhsType>
static inline Eval& scalarMulEq(Eval& a, const RhsType& other)
// m(c*u)' = c*u'
template <class RhsValueType>
Evaluation& operator*=(const RhsValueType& other)
{
a.data_[0] *= other;
a.data_[1] *= other;
a.data_[2] *= other;
a.data_[3] *= other;
return a;
data_[0] *= other;
data_[1] *= other;
data_[2] *= other;
data_[3] *= other;
return *this;
}
static inline Eval& divEq(Eval& a, const Eval& b)
// 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 ValueT v_vv = 1.0 / b.value();
const ValueT u_vv = a.value() * v_vv * v_vv;
const ValueT v_vv = 1.0 / other.value();
const ValueT u_vv = value() * v_vv * v_vv;
// value
a.data_[0] *= v_vv;
data_[valuepos_] *= v_vv;
// derivatives
a.data_[1] = a.data_[1]*v_vv - b.data_[1]*u_vv;
a.data_[2] = a.data_[2]*v_vv - b.data_[2]*u_vv;
a.data_[3] = a.data_[3]*v_vv - b.data_[3]*u_vv;
return a;
}
template <class RhsType>
static inline Eval& scalarDivEq(Eval& a, const RhsType& other)
{
a.data_[0] /= other;
a.data_[1] /= other;
a.data_[2] /= other;
a.data_[3] /= other;
return a;
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;
}
// divide value and derivatives by value of other
template <class RhsValueType>
static inline Eval divide(const RhsValueType& a, const Eval& b)
Evaluation& operator/=(const RhsValueType& other)
{
Eval result;
result.setValue( a/b.value() );
const ValueT df_dg = - result.value()/b.value();
ValueType tmp = 1.0/other;
data_[0] *= tmp;
data_[1] *= tmp;
data_[2] *= tmp;
data_[3] *= tmp;
return *this;
}
// division of a constant by an Evaluation
template <class RhsValueType>
static inline Evaluation divide(const RhsValueType& a, const Evaluation& b)
{
Evaluation result;
ValueType tmp = 1.0/b.value();
result.setValue( a*tmp );
const ValueT df_dg = - result.value()*tmp;
result.data_[1] = df_dg*b.data_[1];
result.data_[2] = df_dg*b.data_[2];
@ -163,8 +305,181 @@ public:
return result;
}
// 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);
return (result -= other);
}
// subtract constant from evaluation object
template <class RhsValueType>
Evaluation operator-(const RhsValueType& other) const
{
Evaluation result(*this);
return (result -= other);
}
// negation (unary minus) operator
Evaluation operator-() const
{
Evaluation result;
// set value and derivatives to negative
result.data_[0] = - data_[0];
result.data_[1] = - data_[1];
result.data_[2] = - data_[2];
result.data_[3] = - data_[3];
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_[0] = other.data_[0];
data_[1] = other.data_[1];
data_[2] = other.data_[2];
data_[3] = other.data_[3];
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
template <class RhsValueType>
void setValue(const RhsValueType& val)
{ data_[valuepos_] = val; }
// return varIdx'th derivative
const ValueType& derivative(int varIdx) const
{
assert(0 <= varIdx && varIdx < size);
return data_[dstart_ + varIdx];
}
// set derivative at position varIdx
void setDerivative(int varIdx, const ValueType& derVal)
{
assert(0 <= varIdx && varIdx < size);
data_[dstart_ + varIdx] = derVal;
}
private:
std::array<ValueT, length_> data_;
};
} } // namespace DenseAd, Opm
#endif // OPM_DENSEAD_EVALUATION3_HPP
#endif // OPM_DENSEAD_EVALUATION3_HPP

View File

@ -23,148 +23,289 @@
/*!
* \file
*
* \brief This file specializes the dense-AD Evaluation class for 4 derivatives.
*
* \attention THIS FILE GETS AUTOMATICALLY GENERATED BY THE "genEvalSpecializations.py"
* SCRIPT. DO NOT EDIT IT MANUALLY!
*/
#ifndef OPM_DENSEAD_EVALUATION4_HPP
#define OPM_DENSEAD_EVALUATION4_HPP
#include "Math.hpp"
#include <opm/common/Valgrind.hpp>
#include <dune/common/version.hh>
#include <array>
#include <cmath>
#include <cassert>
#include <cstring>
#include <iostream>
#include <algorithm>
namespace Opm {
namespace DenseAd {
template <class ValueT>
struct EvaluationOps<ValueT, 4>
class Evaluation<ValueT, 4>
{
private:
typedef Evaluation<ValueT, 4 > Eval;
public:
//! field type
typedef ValueT ValueType;
//! number of derivatives
static constexpr int size = 4;
protected:
//! length of internal data vector
static constexpr int length_ = size + 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:
typedef std::array<ValueT, 5 > DataVector;
static inline void assign(Eval& a, const Eval& b)
//! default constructor
Evaluation() : data_()
{}
//! copy other function evaluation
Evaluation(const Evaluation& other)
{
a.data_[0] = b.data_[0];
a.data_[1] = b.data_[1];
a.data_[2] = b.data_[2];
a.data_[3] = b.data_[3];
a.data_[4] = b.data_[4];
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];
}
static inline void assignNegative(Eval& a, const Eval& b)
// 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 < size);
data_[varPos + dstart_] = 1.0;
Valgrind::CheckDefined(data_);
}
// set all derivatives to zero
void clearDerivatives()
{
a.data_[0] = - b.data_[0];
a.data_[1] = - b.data_[1];
a.data_[2] = - b.data_[2];
a.data_[3] = - b.data_[3];
a.data_[4] = - b.data_[4];
data_[1] = 0.0;
data_[2] = 0.0;
data_[3] = 0.0;
data_[4] = 0.0;
}
static inline void clearDerivatives(Eval& a)
// 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 < size; ++varIdx)
os << " " << derivative(varIdx);
}
// copy all derivatives from other
void copyDerivatives(const Evaluation& other)
{
for (int i = dstart_; 0 < dend_; ++i)
data_[i] = other.data_[i];
}
// add value and derivatives from other to this values and derivatives
Evaluation& operator+=(const Evaluation& other)
{
a.data_[1] = 0.0;
a.data_[2] = 0.0;
a.data_[3] = 0.0;
a.data_[4] = 0.0;
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];
return *this;
}
static inline Eval& addEq(Eval& a, const Eval& b)
// 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)
{
a.data_[0] += b.data_[0];
a.data_[1] += b.data_[1];
a.data_[2] += b.data_[2];
a.data_[3] += b.data_[3];
a.data_[4] += b.data_[4];
return a;
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];
return *this;
}
static inline Eval& subEq(Eval& a, const Eval& b)
// 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;
a.data_[0] -= b.data_[0];
a.data_[1] -= b.data_[1];
a.data_[2] -= b.data_[2];
a.data_[3] -= b.data_[3];
a.data_[4] -= b.data_[4];
return a;
return *this;
}
static inline Eval& mulEq(Eval& a, const Eval& b)
// 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 ValueT u = a.value();
const ValueT v = b.value();
const ValueT u = this->value();
const ValueT v = other.value();
// value
a.data_[0] *= v ;
this->data_[valuepos_] *= v ;
// derivatives
a.data_[1] = a.data_[1]*v + b.data_[1] * u;
a.data_[2] = a.data_[2]*v + b.data_[2] * u;
a.data_[3] = a.data_[3]*v + b.data_[3] * u;
a.data_[4] = a.data_[4]*v + b.data_[4] * u;
return a;
this->data_[1] = this->data_[1]*v + other.data_[1] * u;
this->data_[2] = this->data_[2]*v + other.data_[2] * u;
this->data_[3] = this->data_[3]*v + other.data_[3] * u;
this->data_[4] = this->data_[4]*v + other.data_[4] * u;
return *this;
}
template <class RhsType>
static inline Eval& scalarMulEq(Eval& a, const RhsType& other)
// m(c*u)' = c*u'
template <class RhsValueType>
Evaluation& operator*=(const RhsValueType& other)
{
a.data_[0] *= other;
a.data_[1] *= other;
a.data_[2] *= other;
a.data_[3] *= other;
a.data_[4] *= other;
return a;
data_[0] *= other;
data_[1] *= other;
data_[2] *= other;
data_[3] *= other;
data_[4] *= other;
return *this;
}
static inline Eval& divEq(Eval& a, const Eval& b)
// 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 ValueT v_vv = 1.0 / b.value();
const ValueT u_vv = a.value() * v_vv * v_vv;
const ValueT v_vv = 1.0 / other.value();
const ValueT u_vv = value() * v_vv * v_vv;
// value
a.data_[0] *= v_vv;
data_[valuepos_] *= v_vv;
// derivatives
a.data_[1] = a.data_[1]*v_vv - b.data_[1]*u_vv;
a.data_[2] = a.data_[2]*v_vv - b.data_[2]*u_vv;
a.data_[3] = a.data_[3]*v_vv - b.data_[3]*u_vv;
a.data_[4] = a.data_[4]*v_vv - b.data_[4]*u_vv;
return a;
}
template <class RhsType>
static inline Eval& scalarDivEq(Eval& a, const RhsType& other)
{
a.data_[0] /= other;
a.data_[1] /= other;
a.data_[2] /= other;
a.data_[3] /= other;
a.data_[4] /= other;
return a;
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;
return *this;
}
// divide value and derivatives by value of other
template <class RhsValueType>
static inline Eval divide(const RhsValueType& a, const Eval& b)
Evaluation& operator/=(const RhsValueType& other)
{
Eval result;
result.setValue( a/b.value() );
const ValueT df_dg = - result.value()/b.value();
ValueType tmp = 1.0/other;
data_[0] *= tmp;
data_[1] *= tmp;
data_[2] *= tmp;
data_[3] *= tmp;
data_[4] *= tmp;
return *this;
}
// division of a constant by an Evaluation
template <class RhsValueType>
static inline Evaluation divide(const RhsValueType& a, const Evaluation& b)
{
Evaluation result;
ValueType tmp = 1.0/b.value();
result.setValue( a*tmp );
const ValueT df_dg = - result.value()*tmp;
result.data_[1] = df_dg*b.data_[1];
result.data_[2] = df_dg*b.data_[2];
@ -173,8 +314,183 @@ public:
return result;
}
// 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);
return (result -= other);
}
// subtract constant from evaluation object
template <class RhsValueType>
Evaluation operator-(const RhsValueType& other) const
{
Evaluation result(*this);
return (result -= other);
}
// negation (unary minus) operator
Evaluation operator-() const
{
Evaluation result;
// set value and derivatives to negative
result.data_[0] = - data_[0];
result.data_[1] = - data_[1];
result.data_[2] = - data_[2];
result.data_[3] = - data_[3];
result.data_[4] = - data_[4];
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_[0] = other.data_[0];
data_[1] = other.data_[1];
data_[2] = other.data_[2];
data_[3] = other.data_[3];
data_[4] = other.data_[4];
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
template <class RhsValueType>
void setValue(const RhsValueType& val)
{ data_[valuepos_] = val; }
// return varIdx'th derivative
const ValueType& derivative(int varIdx) const
{
assert(0 <= varIdx && varIdx < size);
return data_[dstart_ + varIdx];
}
// set derivative at position varIdx
void setDerivative(int varIdx, const ValueType& derVal)
{
assert(0 <= varIdx && varIdx < size);
data_[dstart_ + varIdx] = derVal;
}
private:
std::array<ValueT, length_> data_;
};
} } // namespace DenseAd, Opm
#endif // OPM_DENSEAD_EVALUATION4_HPP
#endif // OPM_DENSEAD_EVALUATION4_HPP

View File

@ -23,157 +23,297 @@
/*!
* \file
*
* \brief This file specializes the dense-AD Evaluation class for 5 derivatives.
*
* \attention THIS FILE GETS AUTOMATICALLY GENERATED BY THE "genEvalSpecializations.py"
* SCRIPT. DO NOT EDIT IT MANUALLY!
*/
#ifndef OPM_DENSEAD_EVALUATION5_HPP
#define OPM_DENSEAD_EVALUATION5_HPP
#include "Math.hpp"
#include <opm/common/Valgrind.hpp>
#include <dune/common/version.hh>
#include <array>
#include <cmath>
#include <cassert>
#include <cstring>
#include <iostream>
#include <algorithm>
namespace Opm {
namespace DenseAd {
template <class ValueT>
struct EvaluationOps<ValueT, 5>
class Evaluation<ValueT, 5>
{
private:
typedef Evaluation<ValueT, 5 > Eval;
public:
//! field type
typedef ValueT ValueType;
//! number of derivatives
static constexpr int size = 5;
protected:
//! length of internal data vector
static constexpr int length_ = size + 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:
typedef std::array<ValueT, 6 > DataVector;
static inline void assign(Eval& a, const Eval& b)
//! default constructor
Evaluation() : data_()
{}
//! copy other function evaluation
Evaluation(const Evaluation& other)
{
a.data_[0] = b.data_[0];
a.data_[1] = b.data_[1];
a.data_[2] = b.data_[2];
a.data_[3] = b.data_[3];
a.data_[4] = b.data_[4];
a.data_[5] = b.data_[5];
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];
}
static inline void assignNegative(Eval& a, const Eval& b)
// 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 < size);
data_[varPos + dstart_] = 1.0;
Valgrind::CheckDefined(data_);
}
// set all derivatives to zero
void clearDerivatives()
{
a.data_[0] = - b.data_[0];
a.data_[1] = - b.data_[1];
a.data_[2] = - b.data_[2];
a.data_[3] = - b.data_[3];
a.data_[4] = - b.data_[4];
a.data_[5] = - b.data_[5];
data_[1] = 0.0;
data_[2] = 0.0;
data_[3] = 0.0;
data_[4] = 0.0;
data_[5] = 0.0;
}
static inline void clearDerivatives(Eval& a)
// 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 < size; ++varIdx)
os << " " << derivative(varIdx);
}
// copy all derivatives from other
void copyDerivatives(const Evaluation& other)
{
for (int i = dstart_; 0 < dend_; ++i)
data_[i] = other.data_[i];
}
// add value and derivatives from other to this values and derivatives
Evaluation& operator+=(const Evaluation& other)
{
a.data_[1] = 0.0;
a.data_[2] = 0.0;
a.data_[3] = 0.0;
a.data_[4] = 0.0;
a.data_[5] = 0.0;
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];
return *this;
}
static inline Eval& addEq(Eval& a, const Eval& b)
// 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)
{
a.data_[0] += b.data_[0];
a.data_[1] += b.data_[1];
a.data_[2] += b.data_[2];
a.data_[3] += b.data_[3];
a.data_[4] += b.data_[4];
a.data_[5] += b.data_[5];
return a;
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];
return *this;
}
static inline Eval& subEq(Eval& a, const Eval& b)
// 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;
a.data_[0] -= b.data_[0];
a.data_[1] -= b.data_[1];
a.data_[2] -= b.data_[2];
a.data_[3] -= b.data_[3];
a.data_[4] -= b.data_[4];
a.data_[5] -= b.data_[5];
return a;
return *this;
}
static inline Eval& mulEq(Eval& a, const Eval& b)
// 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 ValueT u = a.value();
const ValueT v = b.value();
const ValueT u = this->value();
const ValueT v = other.value();
// value
a.data_[0] *= v ;
this->data_[valuepos_] *= v ;
// derivatives
a.data_[1] = a.data_[1]*v + b.data_[1] * u;
a.data_[2] = a.data_[2]*v + b.data_[2] * u;
a.data_[3] = a.data_[3]*v + b.data_[3] * u;
a.data_[4] = a.data_[4]*v + b.data_[4] * u;
a.data_[5] = a.data_[5]*v + b.data_[5] * u;
return a;
this->data_[1] = this->data_[1]*v + other.data_[1] * u;
this->data_[2] = this->data_[2]*v + other.data_[2] * u;
this->data_[3] = this->data_[3]*v + other.data_[3] * u;
this->data_[4] = this->data_[4]*v + other.data_[4] * u;
this->data_[5] = this->data_[5]*v + other.data_[5] * u;
return *this;
}
template <class RhsType>
static inline Eval& scalarMulEq(Eval& a, const RhsType& other)
// m(c*u)' = c*u'
template <class RhsValueType>
Evaluation& operator*=(const RhsValueType& other)
{
a.data_[0] *= other;
a.data_[1] *= other;
a.data_[2] *= other;
a.data_[3] *= other;
a.data_[4] *= other;
a.data_[5] *= other;
return a;
data_[0] *= other;
data_[1] *= other;
data_[2] *= other;
data_[3] *= other;
data_[4] *= other;
data_[5] *= other;
return *this;
}
static inline Eval& divEq(Eval& a, const Eval& b)
// 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 ValueT v_vv = 1.0 / b.value();
const ValueT u_vv = a.value() * v_vv * v_vv;
const ValueT v_vv = 1.0 / other.value();
const ValueT u_vv = value() * v_vv * v_vv;
// value
a.data_[0] *= v_vv;
data_[valuepos_] *= v_vv;
// derivatives
a.data_[1] = a.data_[1]*v_vv - b.data_[1]*u_vv;
a.data_[2] = a.data_[2]*v_vv - b.data_[2]*u_vv;
a.data_[3] = a.data_[3]*v_vv - b.data_[3]*u_vv;
a.data_[4] = a.data_[4]*v_vv - b.data_[4]*u_vv;
a.data_[5] = a.data_[5]*v_vv - b.data_[5]*u_vv;
return a;
}
template <class RhsType>
static inline Eval& scalarDivEq(Eval& a, const RhsType& other)
{
a.data_[0] /= other;
a.data_[1] /= other;
a.data_[2] /= other;
a.data_[3] /= other;
a.data_[4] /= other;
a.data_[5] /= other;
return a;
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;
return *this;
}
// divide value and derivatives by value of other
template <class RhsValueType>
static inline Eval divide(const RhsValueType& a, const Eval& b)
Evaluation& operator/=(const RhsValueType& other)
{
Eval result;
result.setValue( a/b.value() );
const ValueT df_dg = - result.value()/b.value();
ValueType tmp = 1.0/other;
data_[0] *= tmp;
data_[1] *= tmp;
data_[2] *= tmp;
data_[3] *= tmp;
data_[4] *= tmp;
data_[5] *= tmp;
return *this;
}
// division of a constant by an Evaluation
template <class RhsValueType>
static inline Evaluation divide(const RhsValueType& a, const Evaluation& b)
{
Evaluation result;
ValueType tmp = 1.0/b.value();
result.setValue( a*tmp );
const ValueT df_dg = - result.value()*tmp;
result.data_[1] = df_dg*b.data_[1];
result.data_[2] = df_dg*b.data_[2];
@ -183,8 +323,185 @@ public:
return result;
}
// 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);
return (result -= other);
}
// subtract constant from evaluation object
template <class RhsValueType>
Evaluation operator-(const RhsValueType& other) const
{
Evaluation result(*this);
return (result -= other);
}
// negation (unary minus) operator
Evaluation operator-() const
{
Evaluation result;
// set value and derivatives to negative
result.data_[0] = - data_[0];
result.data_[1] = - data_[1];
result.data_[2] = - data_[2];
result.data_[3] = - data_[3];
result.data_[4] = - data_[4];
result.data_[5] = - data_[5];
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_[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];
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
template <class RhsValueType>
void setValue(const RhsValueType& val)
{ data_[valuepos_] = val; }
// return varIdx'th derivative
const ValueType& derivative(int varIdx) const
{
assert(0 <= varIdx && varIdx < size);
return data_[dstart_ + varIdx];
}
// set derivative at position varIdx
void setDerivative(int varIdx, const ValueType& derVal)
{
assert(0 <= varIdx && varIdx < size);
data_[dstart_ + varIdx] = derVal;
}
private:
std::array<ValueT, length_> data_;
};
} } // namespace DenseAd, Opm
#endif // OPM_DENSEAD_EVALUATION5_HPP
#endif // OPM_DENSEAD_EVALUATION5_HPP

View File

@ -23,166 +23,305 @@
/*!
* \file
*
* \brief This file specializes the dense-AD Evaluation class for 6 derivatives.
*
* \attention THIS FILE GETS AUTOMATICALLY GENERATED BY THE "genEvalSpecializations.py"
* SCRIPT. DO NOT EDIT IT MANUALLY!
*/
#ifndef OPM_DENSEAD_EVALUATION6_HPP
#define OPM_DENSEAD_EVALUATION6_HPP
#include "Math.hpp"
#include <opm/common/Valgrind.hpp>
#include <dune/common/version.hh>
#include <array>
#include <cmath>
#include <cassert>
#include <cstring>
#include <iostream>
#include <algorithm>
namespace Opm {
namespace DenseAd {
template <class ValueT>
struct EvaluationOps<ValueT, 6>
class Evaluation<ValueT, 6>
{
private:
typedef Evaluation<ValueT, 6 > Eval;
public:
//! field type
typedef ValueT ValueType;
//! number of derivatives
static constexpr int size = 6;
protected:
//! length of internal data vector
static constexpr int length_ = size + 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:
typedef std::array<ValueT, 7 > DataVector;
static inline void assign(Eval& a, const Eval& b)
//! default constructor
Evaluation() : data_()
{}
//! copy other function evaluation
Evaluation(const Evaluation& other)
{
a.data_[0] = b.data_[0];
a.data_[1] = b.data_[1];
a.data_[2] = b.data_[2];
a.data_[3] = b.data_[3];
a.data_[4] = b.data_[4];
a.data_[5] = b.data_[5];
a.data_[6] = b.data_[6];
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];
}
static inline void assignNegative(Eval& a, const Eval& b)
// 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 < size);
data_[varPos + dstart_] = 1.0;
Valgrind::CheckDefined(data_);
}
// set all derivatives to zero
void clearDerivatives()
{
a.data_[0] = - b.data_[0];
a.data_[1] = - b.data_[1];
a.data_[2] = - b.data_[2];
a.data_[3] = - b.data_[3];
a.data_[4] = - b.data_[4];
a.data_[5] = - b.data_[5];
a.data_[6] = - b.data_[6];
data_[1] = 0.0;
data_[2] = 0.0;
data_[3] = 0.0;
data_[4] = 0.0;
data_[5] = 0.0;
data_[6] = 0.0;
}
static inline void clearDerivatives(Eval& a)
// 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 < size; ++varIdx)
os << " " << derivative(varIdx);
}
// copy all derivatives from other
void copyDerivatives(const Evaluation& other)
{
for (int i = dstart_; 0 < dend_; ++i)
data_[i] = other.data_[i];
}
// add value and derivatives from other to this values and derivatives
Evaluation& operator+=(const Evaluation& other)
{
a.data_[1] = 0.0;
a.data_[2] = 0.0;
a.data_[3] = 0.0;
a.data_[4] = 0.0;
a.data_[5] = 0.0;
a.data_[6] = 0.0;
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;
}
static inline Eval& addEq(Eval& a, const Eval& b)
// 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)
{
a.data_[0] += b.data_[0];
a.data_[1] += b.data_[1];
a.data_[2] += b.data_[2];
a.data_[3] += b.data_[3];
a.data_[4] += b.data_[4];
a.data_[5] += b.data_[5];
a.data_[6] += b.data_[6];
return a;
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;
}
static inline Eval& subEq(Eval& a, const Eval& b)
// 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;
a.data_[0] -= b.data_[0];
a.data_[1] -= b.data_[1];
a.data_[2] -= b.data_[2];
a.data_[3] -= b.data_[3];
a.data_[4] -= b.data_[4];
a.data_[5] -= b.data_[5];
a.data_[6] -= b.data_[6];
return a;
return *this;
}
static inline Eval& mulEq(Eval& a, const Eval& b)
// 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 ValueT u = a.value();
const ValueT v = b.value();
const ValueT u = this->value();
const ValueT v = other.value();
// value
a.data_[0] *= v ;
this->data_[valuepos_] *= v ;
// derivatives
a.data_[1] = a.data_[1]*v + b.data_[1] * u;
a.data_[2] = a.data_[2]*v + b.data_[2] * u;
a.data_[3] = a.data_[3]*v + b.data_[3] * u;
a.data_[4] = a.data_[4]*v + b.data_[4] * u;
a.data_[5] = a.data_[5]*v + b.data_[5] * u;
a.data_[6] = a.data_[6]*v + b.data_[6] * u;
return a;
this->data_[1] = this->data_[1]*v + other.data_[1] * u;
this->data_[2] = this->data_[2]*v + other.data_[2] * u;
this->data_[3] = this->data_[3]*v + other.data_[3] * u;
this->data_[4] = this->data_[4]*v + other.data_[4] * u;
this->data_[5] = this->data_[5]*v + other.data_[5] * u;
this->data_[6] = this->data_[6]*v + other.data_[6] * u;
return *this;
}
template <class RhsType>
static inline Eval& scalarMulEq(Eval& a, const RhsType& other)
// m(c*u)' = c*u'
template <class RhsValueType>
Evaluation& operator*=(const RhsValueType& other)
{
a.data_[0] *= other;
a.data_[1] *= other;
a.data_[2] *= other;
a.data_[3] *= other;
a.data_[4] *= other;
a.data_[5] *= other;
a.data_[6] *= other;
return a;
data_[0] *= other;
data_[1] *= other;
data_[2] *= other;
data_[3] *= other;
data_[4] *= other;
data_[5] *= other;
data_[6] *= other;
return *this;
}
static inline Eval& divEq(Eval& a, const Eval& b)
// 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 ValueT v_vv = 1.0 / b.value();
const ValueT u_vv = a.value() * v_vv * v_vv;
const ValueT v_vv = 1.0 / other.value();
const ValueT u_vv = value() * v_vv * v_vv;
// value
a.data_[0] *= v_vv;
data_[valuepos_] *= v_vv;
// derivatives
a.data_[1] = a.data_[1]*v_vv - b.data_[1]*u_vv;
a.data_[2] = a.data_[2]*v_vv - b.data_[2]*u_vv;
a.data_[3] = a.data_[3]*v_vv - b.data_[3]*u_vv;
a.data_[4] = a.data_[4]*v_vv - b.data_[4]*u_vv;
a.data_[5] = a.data_[5]*v_vv - b.data_[5]*u_vv;
a.data_[6] = a.data_[6]*v_vv - b.data_[6]*u_vv;
return a;
}
template <class RhsType>
static inline Eval& scalarDivEq(Eval& a, const RhsType& other)
{
a.data_[0] /= other;
a.data_[1] /= other;
a.data_[2] /= other;
a.data_[3] /= other;
a.data_[4] /= other;
a.data_[5] /= other;
a.data_[6] /= other;
return a;
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;
}
// divide value and derivatives by value of other
template <class RhsValueType>
static inline Eval divide(const RhsValueType& a, const Eval& b)
Evaluation& operator/=(const RhsValueType& other)
{
Eval result;
result.setValue( a/b.value() );
const ValueT df_dg = - result.value()/b.value();
ValueType tmp = 1.0/other;
data_[0] *= tmp;
data_[1] *= tmp;
data_[2] *= tmp;
data_[3] *= tmp;
data_[4] *= tmp;
data_[5] *= tmp;
data_[6] *= tmp;
return *this;
}
// division of a constant by an Evaluation
template <class RhsValueType>
static inline Evaluation divide(const RhsValueType& a, const Evaluation& b)
{
Evaluation result;
ValueType tmp = 1.0/b.value();
result.setValue( a*tmp );
const ValueT df_dg = - result.value()*tmp;
result.data_[1] = df_dg*b.data_[1];
result.data_[2] = df_dg*b.data_[2];
@ -193,8 +332,187 @@ public:
return result;
}
// 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);
return (result -= other);
}
// subtract constant from evaluation object
template <class RhsValueType>
Evaluation operator-(const RhsValueType& other) const
{
Evaluation result(*this);
return (result -= other);
}
// negation (unary minus) operator
Evaluation operator-() const
{
Evaluation result;
// set value and derivatives to negative
result.data_[0] = - data_[0];
result.data_[1] = - data_[1];
result.data_[2] = - data_[2];
result.data_[3] = - data_[3];
result.data_[4] = - data_[4];
result.data_[5] = - data_[5];
result.data_[6] = - data_[6];
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_[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;
}
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
template <class RhsValueType>
void setValue(const RhsValueType& val)
{ data_[valuepos_] = val; }
// return varIdx'th derivative
const ValueType& derivative(int varIdx) const
{
assert(0 <= varIdx && varIdx < size);
return data_[dstart_ + varIdx];
}
// set derivative at position varIdx
void setDerivative(int varIdx, const ValueType& derVal)
{
assert(0 <= varIdx && varIdx < size);
data_[dstart_ + varIdx] = derVal;
}
private:
std::array<ValueT, length_> data_;
};
} } // namespace DenseAd, Opm
#endif // OPM_DENSEAD_EVALUATION6_HPP
#endif // OPM_DENSEAD_EVALUATION6_HPP

View File

@ -23,175 +23,313 @@
/*!
* \file
*
* \brief This file specializes the dense-AD Evaluation class for 7 derivatives.
*
* \attention THIS FILE GETS AUTOMATICALLY GENERATED BY THE "genEvalSpecializations.py"
* SCRIPT. DO NOT EDIT IT MANUALLY!
*/
#ifndef OPM_DENSEAD_EVALUATION7_HPP
#define OPM_DENSEAD_EVALUATION7_HPP
#include "Math.hpp"
#include <opm/common/Valgrind.hpp>
#include <dune/common/version.hh>
#include <array>
#include <cmath>
#include <cassert>
#include <cstring>
#include <iostream>
#include <algorithm>
namespace Opm {
namespace DenseAd {
template <class ValueT>
struct EvaluationOps<ValueT, 7>
class Evaluation<ValueT, 7>
{
private:
typedef Evaluation<ValueT, 7 > Eval;
public:
//! field type
typedef ValueT ValueType;
//! number of derivatives
static constexpr int size = 7;
protected:
//! length of internal data vector
static constexpr int length_ = size + 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:
typedef std::array<ValueT, 8 > DataVector;
static inline void assign(Eval& a, const Eval& b)
//! default constructor
Evaluation() : data_()
{}
//! copy other function evaluation
Evaluation(const Evaluation& other)
{
a.data_[0] = b.data_[0];
a.data_[1] = b.data_[1];
a.data_[2] = b.data_[2];
a.data_[3] = b.data_[3];
a.data_[4] = b.data_[4];
a.data_[5] = b.data_[5];
a.data_[6] = b.data_[6];
a.data_[7] = b.data_[7];
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];
}
static inline void assignNegative(Eval& a, const Eval& b)
// 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 < size);
data_[varPos + dstart_] = 1.0;
Valgrind::CheckDefined(data_);
}
// set all derivatives to zero
void clearDerivatives()
{
a.data_[0] = - b.data_[0];
a.data_[1] = - b.data_[1];
a.data_[2] = - b.data_[2];
a.data_[3] = - b.data_[3];
a.data_[4] = - b.data_[4];
a.data_[5] = - b.data_[5];
a.data_[6] = - b.data_[6];
a.data_[7] = - b.data_[7];
data_[1] = 0.0;
data_[2] = 0.0;
data_[3] = 0.0;
data_[4] = 0.0;
data_[5] = 0.0;
data_[6] = 0.0;
data_[7] = 0.0;
}
static inline void clearDerivatives(Eval& a)
// 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 < size; ++varIdx)
os << " " << derivative(varIdx);
}
// copy all derivatives from other
void copyDerivatives(const Evaluation& other)
{
for (int i = dstart_; 0 < dend_; ++i)
data_[i] = other.data_[i];
}
// add value and derivatives from other to this values and derivatives
Evaluation& operator+=(const Evaluation& other)
{
a.data_[1] = 0.0;
a.data_[2] = 0.0;
a.data_[3] = 0.0;
a.data_[4] = 0.0;
a.data_[5] = 0.0;
a.data_[6] = 0.0;
a.data_[7] = 0.0;
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];
return *this;
}
static inline Eval& addEq(Eval& a, const Eval& b)
// 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)
{
a.data_[0] += b.data_[0];
a.data_[1] += b.data_[1];
a.data_[2] += b.data_[2];
a.data_[3] += b.data_[3];
a.data_[4] += b.data_[4];
a.data_[5] += b.data_[5];
a.data_[6] += b.data_[6];
a.data_[7] += b.data_[7];
return a;
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];
return *this;
}
static inline Eval& subEq(Eval& a, const Eval& b)
// 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;
a.data_[0] -= b.data_[0];
a.data_[1] -= b.data_[1];
a.data_[2] -= b.data_[2];
a.data_[3] -= b.data_[3];
a.data_[4] -= b.data_[4];
a.data_[5] -= b.data_[5];
a.data_[6] -= b.data_[6];
a.data_[7] -= b.data_[7];
return a;
return *this;
}
static inline Eval& mulEq(Eval& a, const Eval& b)
// 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 ValueT u = a.value();
const ValueT v = b.value();
const ValueT u = this->value();
const ValueT v = other.value();
// value
a.data_[0] *= v ;
this->data_[valuepos_] *= v ;
// derivatives
a.data_[1] = a.data_[1]*v + b.data_[1] * u;
a.data_[2] = a.data_[2]*v + b.data_[2] * u;
a.data_[3] = a.data_[3]*v + b.data_[3] * u;
a.data_[4] = a.data_[4]*v + b.data_[4] * u;
a.data_[5] = a.data_[5]*v + b.data_[5] * u;
a.data_[6] = a.data_[6]*v + b.data_[6] * u;
a.data_[7] = a.data_[7]*v + b.data_[7] * u;
return a;
this->data_[1] = this->data_[1]*v + other.data_[1] * u;
this->data_[2] = this->data_[2]*v + other.data_[2] * u;
this->data_[3] = this->data_[3]*v + other.data_[3] * u;
this->data_[4] = this->data_[4]*v + other.data_[4] * u;
this->data_[5] = this->data_[5]*v + other.data_[5] * u;
this->data_[6] = this->data_[6]*v + other.data_[6] * u;
this->data_[7] = this->data_[7]*v + other.data_[7] * u;
return *this;
}
template <class RhsType>
static inline Eval& scalarMulEq(Eval& a, const RhsType& other)
// m(c*u)' = c*u'
template <class RhsValueType>
Evaluation& operator*=(const RhsValueType& other)
{
a.data_[0] *= other;
a.data_[1] *= other;
a.data_[2] *= other;
a.data_[3] *= other;
a.data_[4] *= other;
a.data_[5] *= other;
a.data_[6] *= other;
a.data_[7] *= other;
return a;
data_[0] *= other;
data_[1] *= other;
data_[2] *= other;
data_[3] *= other;
data_[4] *= other;
data_[5] *= other;
data_[6] *= other;
data_[7] *= other;
return *this;
}
static inline Eval& divEq(Eval& a, const Eval& b)
// 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 ValueT v_vv = 1.0 / b.value();
const ValueT u_vv = a.value() * v_vv * v_vv;
const ValueT v_vv = 1.0 / other.value();
const ValueT u_vv = value() * v_vv * v_vv;
// value
a.data_[0] *= v_vv;
data_[valuepos_] *= v_vv;
// derivatives
a.data_[1] = a.data_[1]*v_vv - b.data_[1]*u_vv;
a.data_[2] = a.data_[2]*v_vv - b.data_[2]*u_vv;
a.data_[3] = a.data_[3]*v_vv - b.data_[3]*u_vv;
a.data_[4] = a.data_[4]*v_vv - b.data_[4]*u_vv;
a.data_[5] = a.data_[5]*v_vv - b.data_[5]*u_vv;
a.data_[6] = a.data_[6]*v_vv - b.data_[6]*u_vv;
a.data_[7] = a.data_[7]*v_vv - b.data_[7]*u_vv;
return a;
}
template <class RhsType>
static inline Eval& scalarDivEq(Eval& a, const RhsType& other)
{
a.data_[0] /= other;
a.data_[1] /= other;
a.data_[2] /= other;
a.data_[3] /= other;
a.data_[4] /= other;
a.data_[5] /= other;
a.data_[6] /= other;
a.data_[7] /= other;
return a;
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;
return *this;
}
// divide value and derivatives by value of other
template <class RhsValueType>
static inline Eval divide(const RhsValueType& a, const Eval& b)
Evaluation& operator/=(const RhsValueType& other)
{
Eval result;
result.setValue( a/b.value() );
const ValueT df_dg = - result.value()/b.value();
ValueType tmp = 1.0/other;
data_[0] *= tmp;
data_[1] *= tmp;
data_[2] *= tmp;
data_[3] *= tmp;
data_[4] *= tmp;
data_[5] *= tmp;
data_[6] *= tmp;
data_[7] *= tmp;
return *this;
}
// division of a constant by an Evaluation
template <class RhsValueType>
static inline Evaluation divide(const RhsValueType& a, const Evaluation& b)
{
Evaluation result;
ValueType tmp = 1.0/b.value();
result.setValue( a*tmp );
const ValueT df_dg = - result.value()*tmp;
result.data_[1] = df_dg*b.data_[1];
result.data_[2] = df_dg*b.data_[2];
@ -203,8 +341,189 @@ public:
return result;
}
// 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);
return (result -= other);
}
// subtract constant from evaluation object
template <class RhsValueType>
Evaluation operator-(const RhsValueType& other) const
{
Evaluation result(*this);
return (result -= other);
}
// negation (unary minus) operator
Evaluation operator-() const
{
Evaluation result;
// set value and derivatives to negative
result.data_[0] = - data_[0];
result.data_[1] = - data_[1];
result.data_[2] = - data_[2];
result.data_[3] = - data_[3];
result.data_[4] = - data_[4];
result.data_[5] = - data_[5];
result.data_[6] = - data_[6];
result.data_[7] = - data_[7];
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_[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];
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
template <class RhsValueType>
void setValue(const RhsValueType& val)
{ data_[valuepos_] = val; }
// return varIdx'th derivative
const ValueType& derivative(int varIdx) const
{
assert(0 <= varIdx && varIdx < size);
return data_[dstart_ + varIdx];
}
// set derivative at position varIdx
void setDerivative(int varIdx, const ValueType& derVal)
{
assert(0 <= varIdx && varIdx < size);
data_[dstart_ + varIdx] = derVal;
}
private:
std::array<ValueT, length_> data_;
};
} } // namespace DenseAd, Opm
#endif // OPM_DENSEAD_EVALUATION7_HPP
#endif // OPM_DENSEAD_EVALUATION7_HPP

View File

@ -23,184 +23,321 @@
/*!
* \file
*
* \brief This file specializes the dense-AD Evaluation class for 8 derivatives.
*
* \attention THIS FILE GETS AUTOMATICALLY GENERATED BY THE "genEvalSpecializations.py"
* SCRIPT. DO NOT EDIT IT MANUALLY!
*/
#ifndef OPM_DENSEAD_EVALUATION8_HPP
#define OPM_DENSEAD_EVALUATION8_HPP
#include "Math.hpp"
#include <opm/common/Valgrind.hpp>
#include <dune/common/version.hh>
#include <array>
#include <cmath>
#include <cassert>
#include <cstring>
#include <iostream>
#include <algorithm>
namespace Opm {
namespace DenseAd {
template <class ValueT>
struct EvaluationOps<ValueT, 8>
class Evaluation<ValueT, 8>
{
private:
typedef Evaluation<ValueT, 8 > Eval;
public:
//! field type
typedef ValueT ValueType;
//! number of derivatives
static constexpr int size = 8;
protected:
//! length of internal data vector
static constexpr int length_ = size + 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:
typedef std::array<ValueT, 9 > DataVector;
static inline void assign(Eval& a, const Eval& b)
//! default constructor
Evaluation() : data_()
{}
//! copy other function evaluation
Evaluation(const Evaluation& other)
{
a.data_[0] = b.data_[0];
a.data_[1] = b.data_[1];
a.data_[2] = b.data_[2];
a.data_[3] = b.data_[3];
a.data_[4] = b.data_[4];
a.data_[5] = b.data_[5];
a.data_[6] = b.data_[6];
a.data_[7] = b.data_[7];
a.data_[8] = b.data_[8];
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];
}
static inline void assignNegative(Eval& a, const Eval& b)
// 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 < size);
data_[varPos + dstart_] = 1.0;
Valgrind::CheckDefined(data_);
}
// set all derivatives to zero
void clearDerivatives()
{
a.data_[0] = - b.data_[0];
a.data_[1] = - b.data_[1];
a.data_[2] = - b.data_[2];
a.data_[3] = - b.data_[3];
a.data_[4] = - b.data_[4];
a.data_[5] = - b.data_[5];
a.data_[6] = - b.data_[6];
a.data_[7] = - b.data_[7];
a.data_[8] = - b.data_[8];
data_[1] = 0.0;
data_[2] = 0.0;
data_[3] = 0.0;
data_[4] = 0.0;
data_[5] = 0.0;
data_[6] = 0.0;
data_[7] = 0.0;
data_[8] = 0.0;
}
static inline void clearDerivatives(Eval& a)
// 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 < size; ++varIdx)
os << " " << derivative(varIdx);
}
// copy all derivatives from other
void copyDerivatives(const Evaluation& other)
{
for (int i = dstart_; 0 < dend_; ++i)
data_[i] = other.data_[i];
}
// add value and derivatives from other to this values and derivatives
Evaluation& operator+=(const Evaluation& other)
{
a.data_[1] = 0.0;
a.data_[2] = 0.0;
a.data_[3] = 0.0;
a.data_[4] = 0.0;
a.data_[5] = 0.0;
a.data_[6] = 0.0;
a.data_[7] = 0.0;
a.data_[8] = 0.0;
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];
return *this;
}
static inline Eval& addEq(Eval& a, const Eval& b)
// 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)
{
a.data_[0] += b.data_[0];
a.data_[1] += b.data_[1];
a.data_[2] += b.data_[2];
a.data_[3] += b.data_[3];
a.data_[4] += b.data_[4];
a.data_[5] += b.data_[5];
a.data_[6] += b.data_[6];
a.data_[7] += b.data_[7];
a.data_[8] += b.data_[8];
return a;
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];
return *this;
}
static inline Eval& subEq(Eval& a, const Eval& b)
// 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;
a.data_[0] -= b.data_[0];
a.data_[1] -= b.data_[1];
a.data_[2] -= b.data_[2];
a.data_[3] -= b.data_[3];
a.data_[4] -= b.data_[4];
a.data_[5] -= b.data_[5];
a.data_[6] -= b.data_[6];
a.data_[7] -= b.data_[7];
a.data_[8] -= b.data_[8];
return a;
return *this;
}
static inline Eval& mulEq(Eval& a, const Eval& b)
// 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 ValueT u = a.value();
const ValueT v = b.value();
const ValueT u = this->value();
const ValueT v = other.value();
// value
a.data_[0] *= v ;
this->data_[valuepos_] *= v ;
// derivatives
a.data_[1] = a.data_[1]*v + b.data_[1] * u;
a.data_[2] = a.data_[2]*v + b.data_[2] * u;
a.data_[3] = a.data_[3]*v + b.data_[3] * u;
a.data_[4] = a.data_[4]*v + b.data_[4] * u;
a.data_[5] = a.data_[5]*v + b.data_[5] * u;
a.data_[6] = a.data_[6]*v + b.data_[6] * u;
a.data_[7] = a.data_[7]*v + b.data_[7] * u;
a.data_[8] = a.data_[8]*v + b.data_[8] * u;
return a;
this->data_[1] = this->data_[1]*v + other.data_[1] * u;
this->data_[2] = this->data_[2]*v + other.data_[2] * u;
this->data_[3] = this->data_[3]*v + other.data_[3] * u;
this->data_[4] = this->data_[4]*v + other.data_[4] * u;
this->data_[5] = this->data_[5]*v + other.data_[5] * u;
this->data_[6] = this->data_[6]*v + other.data_[6] * u;
this->data_[7] = this->data_[7]*v + other.data_[7] * u;
this->data_[8] = this->data_[8]*v + other.data_[8] * u;
return *this;
}
template <class RhsType>
static inline Eval& scalarMulEq(Eval& a, const RhsType& other)
// m(c*u)' = c*u'
template <class RhsValueType>
Evaluation& operator*=(const RhsValueType& other)
{
a.data_[0] *= other;
a.data_[1] *= other;
a.data_[2] *= other;
a.data_[3] *= other;
a.data_[4] *= other;
a.data_[5] *= other;
a.data_[6] *= other;
a.data_[7] *= other;
a.data_[8] *= other;
return a;
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;
return *this;
}
static inline Eval& divEq(Eval& a, const Eval& b)
// 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 ValueT v_vv = 1.0 / b.value();
const ValueT u_vv = a.value() * v_vv * v_vv;
const ValueT v_vv = 1.0 / other.value();
const ValueT u_vv = value() * v_vv * v_vv;
// value
a.data_[0] *= v_vv;
data_[valuepos_] *= v_vv;
// derivatives
a.data_[1] = a.data_[1]*v_vv - b.data_[1]*u_vv;
a.data_[2] = a.data_[2]*v_vv - b.data_[2]*u_vv;
a.data_[3] = a.data_[3]*v_vv - b.data_[3]*u_vv;
a.data_[4] = a.data_[4]*v_vv - b.data_[4]*u_vv;
a.data_[5] = a.data_[5]*v_vv - b.data_[5]*u_vv;
a.data_[6] = a.data_[6]*v_vv - b.data_[6]*u_vv;
a.data_[7] = a.data_[7]*v_vv - b.data_[7]*u_vv;
a.data_[8] = a.data_[8]*v_vv - b.data_[8]*u_vv;
return a;
}
template <class RhsType>
static inline Eval& scalarDivEq(Eval& a, const RhsType& other)
{
a.data_[0] /= other;
a.data_[1] /= other;
a.data_[2] /= other;
a.data_[3] /= other;
a.data_[4] /= other;
a.data_[5] /= other;
a.data_[6] /= other;
a.data_[7] /= other;
a.data_[8] /= other;
return a;
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;
return *this;
}
// divide value and derivatives by value of other
template <class RhsValueType>
static inline Eval divide(const RhsValueType& a, const Eval& b)
Evaluation& operator/=(const RhsValueType& other)
{
Eval result;
result.setValue( a/b.value() );
const ValueT df_dg = - result.value()/b.value();
ValueType tmp = 1.0/other;
data_[0] *= tmp;
data_[1] *= tmp;
data_[2] *= tmp;
data_[3] *= tmp;
data_[4] *= tmp;
data_[5] *= tmp;
data_[6] *= tmp;
data_[7] *= tmp;
data_[8] *= tmp;
return *this;
}
// division of a constant by an Evaluation
template <class RhsValueType>
static inline Evaluation divide(const RhsValueType& a, const Evaluation& b)
{
Evaluation result;
ValueType tmp = 1.0/b.value();
result.setValue( a*tmp );
const ValueT df_dg = - result.value()*tmp;
result.data_[1] = df_dg*b.data_[1];
result.data_[2] = df_dg*b.data_[2];
@ -213,8 +350,191 @@ public:
return result;
}
// 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);
return (result -= other);
}
// subtract constant from evaluation object
template <class RhsValueType>
Evaluation operator-(const RhsValueType& other) const
{
Evaluation result(*this);
return (result -= other);
}
// negation (unary minus) operator
Evaluation operator-() const
{
Evaluation result;
// set value and derivatives to negative
result.data_[0] = - data_[0];
result.data_[1] = - data_[1];
result.data_[2] = - data_[2];
result.data_[3] = - data_[3];
result.data_[4] = - data_[4];
result.data_[5] = - data_[5];
result.data_[6] = - data_[6];
result.data_[7] = - data_[7];
result.data_[8] = - data_[8];
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_[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];
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
template <class RhsValueType>
void setValue(const RhsValueType& val)
{ data_[valuepos_] = val; }
// return varIdx'th derivative
const ValueType& derivative(int varIdx) const
{
assert(0 <= varIdx && varIdx < size);
return data_[dstart_ + varIdx];
}
// set derivative at position varIdx
void setDerivative(int varIdx, const ValueType& derVal)
{
assert(0 <= varIdx && varIdx < size);
data_[dstart_ + varIdx] = derVal;
}
private:
std::array<ValueT, length_> data_;
};
} } // namespace DenseAd, Opm
#endif // OPM_DENSEAD_EVALUATION8_HPP
#endif // OPM_DENSEAD_EVALUATION8_HPP

View File

@ -23,193 +23,329 @@
/*!
* \file
*
* \brief This file specializes the dense-AD Evaluation class for 9 derivatives.
*
* \attention THIS FILE GETS AUTOMATICALLY GENERATED BY THE "genEvalSpecializations.py"
* SCRIPT. DO NOT EDIT IT MANUALLY!
*/
#ifndef OPM_DENSEAD_EVALUATION9_HPP
#define OPM_DENSEAD_EVALUATION9_HPP
#include "Math.hpp"
#include <opm/common/Valgrind.hpp>
#include <dune/common/version.hh>
#include <array>
#include <cmath>
#include <cassert>
#include <cstring>
#include <iostream>
#include <algorithm>
namespace Opm {
namespace DenseAd {
template <class ValueT>
struct EvaluationOps<ValueT, 9>
class Evaluation<ValueT, 9>
{
private:
typedef Evaluation<ValueT, 9 > Eval;
public:
//! field type
typedef ValueT ValueType;
//! number of derivatives
static constexpr int size = 9;
protected:
//! length of internal data vector
static constexpr int length_ = size + 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:
typedef std::array<ValueT, 10 > DataVector;
static inline void assign(Eval& a, const Eval& b)
//! default constructor
Evaluation() : data_()
{}
//! copy other function evaluation
Evaluation(const Evaluation& other)
{
a.data_[0] = b.data_[0];
a.data_[1] = b.data_[1];
a.data_[2] = b.data_[2];
a.data_[3] = b.data_[3];
a.data_[4] = b.data_[4];
a.data_[5] = b.data_[5];
a.data_[6] = b.data_[6];
a.data_[7] = b.data_[7];
a.data_[8] = b.data_[8];
a.data_[9] = b.data_[9];
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];
}
static inline void assignNegative(Eval& a, const Eval& b)
// 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 < size);
data_[varPos + dstart_] = 1.0;
Valgrind::CheckDefined(data_);
}
// set all derivatives to zero
void clearDerivatives()
{
a.data_[0] = - b.data_[0];
a.data_[1] = - b.data_[1];
a.data_[2] = - b.data_[2];
a.data_[3] = - b.data_[3];
a.data_[4] = - b.data_[4];
a.data_[5] = - b.data_[5];
a.data_[6] = - b.data_[6];
a.data_[7] = - b.data_[7];
a.data_[8] = - b.data_[8];
a.data_[9] = - b.data_[9];
data_[1] = 0.0;
data_[2] = 0.0;
data_[3] = 0.0;
data_[4] = 0.0;
data_[5] = 0.0;
data_[6] = 0.0;
data_[7] = 0.0;
data_[8] = 0.0;
data_[9] = 0.0;
}
static inline void clearDerivatives(Eval& a)
// 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 < size; ++varIdx)
os << " " << derivative(varIdx);
}
// copy all derivatives from other
void copyDerivatives(const Evaluation& other)
{
for (int i = dstart_; 0 < dend_; ++i)
data_[i] = other.data_[i];
}
// add value and derivatives from other to this values and derivatives
Evaluation& operator+=(const Evaluation& other)
{
a.data_[1] = 0.0;
a.data_[2] = 0.0;
a.data_[3] = 0.0;
a.data_[4] = 0.0;
a.data_[5] = 0.0;
a.data_[6] = 0.0;
a.data_[7] = 0.0;
a.data_[8] = 0.0;
a.data_[9] = 0.0;
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];
return *this;
}
static inline Eval& addEq(Eval& a, const Eval& b)
// 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)
{
a.data_[0] += b.data_[0];
a.data_[1] += b.data_[1];
a.data_[2] += b.data_[2];
a.data_[3] += b.data_[3];
a.data_[4] += b.data_[4];
a.data_[5] += b.data_[5];
a.data_[6] += b.data_[6];
a.data_[7] += b.data_[7];
a.data_[8] += b.data_[8];
a.data_[9] += b.data_[9];
return a;
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];
return *this;
}
static inline Eval& subEq(Eval& a, const Eval& b)
// 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;
a.data_[0] -= b.data_[0];
a.data_[1] -= b.data_[1];
a.data_[2] -= b.data_[2];
a.data_[3] -= b.data_[3];
a.data_[4] -= b.data_[4];
a.data_[5] -= b.data_[5];
a.data_[6] -= b.data_[6];
a.data_[7] -= b.data_[7];
a.data_[8] -= b.data_[8];
a.data_[9] -= b.data_[9];
return a;
return *this;
}
static inline Eval& mulEq(Eval& a, const Eval& b)
// 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 ValueT u = a.value();
const ValueT v = b.value();
const ValueT u = this->value();
const ValueT v = other.value();
// value
a.data_[0] *= v ;
this->data_[valuepos_] *= v ;
// derivatives
a.data_[1] = a.data_[1]*v + b.data_[1] * u;
a.data_[2] = a.data_[2]*v + b.data_[2] * u;
a.data_[3] = a.data_[3]*v + b.data_[3] * u;
a.data_[4] = a.data_[4]*v + b.data_[4] * u;
a.data_[5] = a.data_[5]*v + b.data_[5] * u;
a.data_[6] = a.data_[6]*v + b.data_[6] * u;
a.data_[7] = a.data_[7]*v + b.data_[7] * u;
a.data_[8] = a.data_[8]*v + b.data_[8] * u;
a.data_[9] = a.data_[9]*v + b.data_[9] * u;
return a;
this->data_[1] = this->data_[1]*v + other.data_[1] * u;
this->data_[2] = this->data_[2]*v + other.data_[2] * u;
this->data_[3] = this->data_[3]*v + other.data_[3] * u;
this->data_[4] = this->data_[4]*v + other.data_[4] * u;
this->data_[5] = this->data_[5]*v + other.data_[5] * u;
this->data_[6] = this->data_[6]*v + other.data_[6] * u;
this->data_[7] = this->data_[7]*v + other.data_[7] * u;
this->data_[8] = this->data_[8]*v + other.data_[8] * u;
this->data_[9] = this->data_[9]*v + other.data_[9] * u;
return *this;
}
template <class RhsType>
static inline Eval& scalarMulEq(Eval& a, const RhsType& other)
// m(c*u)' = c*u'
template <class RhsValueType>
Evaluation& operator*=(const RhsValueType& other)
{
a.data_[0] *= other;
a.data_[1] *= other;
a.data_[2] *= other;
a.data_[3] *= other;
a.data_[4] *= other;
a.data_[5] *= other;
a.data_[6] *= other;
a.data_[7] *= other;
a.data_[8] *= other;
a.data_[9] *= other;
return a;
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;
return *this;
}
static inline Eval& divEq(Eval& a, const Eval& b)
// 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 ValueT v_vv = 1.0 / b.value();
const ValueT u_vv = a.value() * v_vv * v_vv;
const ValueT v_vv = 1.0 / other.value();
const ValueT u_vv = value() * v_vv * v_vv;
// value
a.data_[0] *= v_vv;
data_[valuepos_] *= v_vv;
// derivatives
a.data_[1] = a.data_[1]*v_vv - b.data_[1]*u_vv;
a.data_[2] = a.data_[2]*v_vv - b.data_[2]*u_vv;
a.data_[3] = a.data_[3]*v_vv - b.data_[3]*u_vv;
a.data_[4] = a.data_[4]*v_vv - b.data_[4]*u_vv;
a.data_[5] = a.data_[5]*v_vv - b.data_[5]*u_vv;
a.data_[6] = a.data_[6]*v_vv - b.data_[6]*u_vv;
a.data_[7] = a.data_[7]*v_vv - b.data_[7]*u_vv;
a.data_[8] = a.data_[8]*v_vv - b.data_[8]*u_vv;
a.data_[9] = a.data_[9]*v_vv - b.data_[9]*u_vv;
return a;
}
template <class RhsType>
static inline Eval& scalarDivEq(Eval& a, const RhsType& other)
{
a.data_[0] /= other;
a.data_[1] /= other;
a.data_[2] /= other;
a.data_[3] /= other;
a.data_[4] /= other;
a.data_[5] /= other;
a.data_[6] /= other;
a.data_[7] /= other;
a.data_[8] /= other;
a.data_[9] /= other;
return a;
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;
return *this;
}
// divide value and derivatives by value of other
template <class RhsValueType>
static inline Eval divide(const RhsValueType& a, const Eval& b)
Evaluation& operator/=(const RhsValueType& other)
{
Eval result;
result.setValue( a/b.value() );
const ValueT df_dg = - result.value()/b.value();
ValueType tmp = 1.0/other;
data_[0] *= tmp;
data_[1] *= tmp;
data_[2] *= tmp;
data_[3] *= tmp;
data_[4] *= tmp;
data_[5] *= tmp;
data_[6] *= tmp;
data_[7] *= tmp;
data_[8] *= tmp;
data_[9] *= tmp;
return *this;
}
// division of a constant by an Evaluation
template <class RhsValueType>
static inline Evaluation divide(const RhsValueType& a, const Evaluation& b)
{
Evaluation result;
ValueType tmp = 1.0/b.value();
result.setValue( a*tmp );
const ValueT df_dg = - result.value()*tmp;
result.data_[1] = df_dg*b.data_[1];
result.data_[2] = df_dg*b.data_[2];
@ -223,8 +359,193 @@ public:
return result;
}
// 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);
return (result -= other);
}
// subtract constant from evaluation object
template <class RhsValueType>
Evaluation operator-(const RhsValueType& other) const
{
Evaluation result(*this);
return (result -= other);
}
// negation (unary minus) operator
Evaluation operator-() const
{
Evaluation result;
// set value and derivatives to negative
result.data_[0] = - data_[0];
result.data_[1] = - data_[1];
result.data_[2] = - data_[2];
result.data_[3] = - data_[3];
result.data_[4] = - data_[4];
result.data_[5] = - data_[5];
result.data_[6] = - data_[6];
result.data_[7] = - data_[7];
result.data_[8] = - data_[8];
result.data_[9] = - data_[9];
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_[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];
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
template <class RhsValueType>
void setValue(const RhsValueType& val)
{ data_[valuepos_] = val; }
// return varIdx'th derivative
const ValueType& derivative(int varIdx) const
{
assert(0 <= varIdx && varIdx < size);
return data_[dstart_ + varIdx];
}
// set derivative at position varIdx
void setDerivative(int varIdx, const ValueType& derVal)
{
assert(0 <= varIdx && varIdx < size);
data_[dstart_ + varIdx] = derVal;
}
private:
std::array<ValueT, length_> data_;
};
} } // namespace DenseAd, Opm
#endif // OPM_DENSEAD_EVALUATION9_HPP
#endif // OPM_DENSEAD_EVALUATION9_HPP