diff --git a/bin/genEvalSpecializations.py b/bin/genEvalSpecializations.py new file mode 100755 index 000000000..1e22bc916 --- /dev/null +++ b/bin/genEvalSpecializations.py @@ -0,0 +1,763 @@ +#! /usr/bin/python +# +# This script provides "hand loop-unrolled" specializations of the +# Evaluation class of dense automatic differentiation so that the +# compiler can more easily emit SIMD instructions. In an ideal world, +# C++ compilers should be smart enough to do this themselfs, but +# contemporary compilers don't seem to exhibit enough brains. +# +# Usage: In the opm-material top-level source directory, run +# `./bin/genEvalSpecializations.py [MAX_DERIVATIVES]`. The script then +# generates specializations for Evaluations with up to MAX_DERIVATIVES +# derivatives. The default for MAX_DERIVATIVES is 12. To run this +# script, you need a python 2 installation where the Jinja2 module is +# available. +# +import os +import sys +import jinja2 + +maxDerivs = 12 +if len(sys.argv) == 2: + maxDerivs = int(sys.argv[1]) + +fileNames = [] + +specializationTemplate = \ +"""// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \\file + * +{% 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 "Evaluation.hpp" +#include "Math.hpp" + +#include + +#include + +#include +#include +#include +#include +#include +#include + +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 Evaluation +{% else %}\ +template +class Evaluation +{% endif %}\ +{ +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: + //! 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) %}\ + 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 + 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 + Evaluation(const RhsValueType& c, int varPos) + { + // The variable position must be in represented by the given variable descriptor + assert(0 <= varPos && varPos < size); + + setValue( c ); + clearDerivatives(); + + data_[varPos + dstart_] = 1.0; + Valgrind::CheckDefined(data_); + } + + // 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) %}\ + data_[{{i}}] = 0.0; +{% endfor %}\ +{% endif %}\ + } + + // create a function evaluation for a "naked" depending variable (i.e., f(x) = x) + template + 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 + 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) + { +{% if numDerivs < 0 %}\ + for (int i = dstart_; i < dend_; ++i) { + data_[i] = other.data_[i]; + } +{% else %}\ +{% for i in range(1, numDerivs+1) %}\ + data_[{{i}}] = other.data_[{{i}}]; +{% endfor %}\ +{% endif %}\ + } + + + // 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 + 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 + Evaluation& operator-=(const RhsValueType& other) + { + // for constants, values are subtracted, derivatives stay the same + data_[ valuepos_ ] -= other; + + return *this; + } + + // multiply values and apply chain rule to derivatives: (u*v)' = (v'u + u'v) + Evaluation& operator*=(const Evaluation& other) + { + // while the values are multiplied, the derivatives follow the product rule, + // i.e., (u*v)' = (v'u + u'v). + const ValueType u = this->value(); + const ValueType v = other.value(); + + // value + data_[valuepos_] *= v ; + + // derivatives +{% if numDerivs < 0 %}\ + for (int i = dstart_; i < dend_; ++i) { + data_[i] = data_[i] * v + other.data_[i] * u; + } +{% else %}\ +{% for i in range(1, numDerivs+1) %}\ + data_[{{i}}] = data_[{{i}}] * v + other.data_[{{i}}] * u; +{% endfor %}\ +{% endif %}\ + + return *this; + } + + // m(c*u)' = c*u' + template + 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) %}\ + data_[{{i}}] *= other; +{% endfor %}\ +{% endif %}\ + + return *this; + } + + // m(u*v)' = (vu' - uv')/v^2 + Evaluation& operator/=(const Evaluation& other) + { + // values are divided, derivatives follow the rule for division, i.e., (u/v)' = (v'u - + // u'v)/v^2. + ValueType& u = data_[ valuepos_ ]; + const ValueType& v = other.value(); +{% if numDerivs < 0 %}\ + for (unsigned idx = dstart_; idx < dend_; ++idx) { + const ValueType& uPrime = data_[idx]; + const ValueType& vPrime = other.data_[idx]; + + data_[idx] = (v*uPrime - u*vPrime)/(v*v); + } +{% else %}\ +{% for i in range(1, numDerivs+1) %}\ + data_[{{i}}] = (v*data_[{{i}}] - u*other.data_[{{i}}])/(v*v); +{% endfor %}\ +{% endif %}\ + u /= v; + + return *this; + } + + // divide value and derivatives by value of other + template + Evaluation& operator/=(const RhsValueType& other) + { + const 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; + } + + // add two evaluation objects + Evaluation operator+(const Evaluation& other) const + { + Evaluation result(*this); + + result += other; + + return result; + } + + // add constant to this object + template + Evaluation operator+(const RhsValueType& other) const + { + Evaluation result(*this); + + result += other; + + return result; + } + + // subtract two evaluation objects + Evaluation operator-(const Evaluation& other) const + { + Evaluation result(*this); + + result -= other; + + return result; + } + + // subtract constant from evaluation object + template + Evaluation operator-(const RhsValueType& other) const + { + Evaluation result(*this); + + result -= other; + + return result; + } + + // negation (unary minus) operator + Evaluation operator-() const + { + Evaluation result; + + // set value and derivatives to negative +{% 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 + 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 + Evaluation operator/(const RhsValueType& other) const + { + Evaluation result(*this); + + result /= other; + + return result; + } + + template + 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 + 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 + bool operator>(RhsValueType other) const + { return value() > other; } + + bool operator>(const Evaluation& other) const + { return value() > other.value(); } + + template + bool operator<(RhsValueType other) const + { return value() < other; } + + bool operator<(const Evaluation& other) const + { return value() < other.value(); } + + template + bool operator>=(RhsValueType other) const + { return value() >= other; } + + bool operator>=(const Evaluation& other) const + { return value() >= other.value(); } + + template + 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 + 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 data_; +}; + +{% if numDerivs < 0 %}\ +// the generic operators are only required for the unspecialized case +template +bool operator<(const RhsValueType& a, const Evaluation& b) +{ return b > a; } + +template +bool operator>(const RhsValueType& a, const Evaluation& b) +{ return b < a; } + +template +bool operator<=(const RhsValueType& a, const Evaluation& b) +{ return b >= a; } + +template +bool operator>=(const RhsValueType& a, const Evaluation& b) +{ return b <= a; } + +template +bool operator!=(const RhsValueType& a, const Evaluation& b) +{ return a != b.value(); } + +template +Evaluation operator+(const RhsValueType& a, const Evaluation& b) +{ + Evaluation result(b); + result += a; + return result; +} + +template +Evaluation operator-(const RhsValueType& a, const Evaluation& b) +{ + Evaluation result(a); + result -= b; + return result; +} + +template +Evaluation operator/(const RhsValueType& a, const Evaluation& b) +{ + Evaluation tmp(a); + tmp /= b; + return tmp; +} + +template +Evaluation operator*(const RhsValueType& a, const Evaluation& b) +{ + Evaluation result(b); + result *= a; + return result; +} + +template +std::ostream& operator<<(std::ostream& os, const Evaluation& 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 +Evaluation abs(const Evaluation&); +}} + +namespace std { +template +const Opm::DenseAd::Evaluation abs(const Opm::DenseAd::Evaluation& 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 + +namespace Dune { +template +struct FieldTraits > +{ +public: + typedef Opm::DenseAd::Evaluation 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 = \ +"""// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \\file + * + * \\brief This file includes all specializations for the dense-AD Evaluation class. + * + * \\attention THIS FILE GETS AUTOMATICALLY GENERATED BY THE "{{ scriptName }}" + * SCRIPT. DO NOT EDIT IT MANUALLY! + */ +#ifndef OPM_DENSEAD_EVALUATION_SPECIALIZATIONS_HPP +#define OPM_DENSEAD_EVALUATION_SPECIALIZATIONS_HPP + +{% for fileName in fileNames %}\ +#include <{{ fileName }}> +{% endfor %}\ + +#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) + + fileName = "opm/material/densead/Evaluation%d.hpp"%numDerivs + fileNames.append(fileName) + + template = jinja2.Template(specializationTemplate) + fileContents = template.render(numDerivs=numDerivs, scriptName=os.path.basename(sys.argv[0])) + + f = open(fileName, "w") + f.write(fileContents) + f.close() + +template = jinja2.Template(includeSpecializationsTemplate) +fileContents = template.render(fileNames=fileNames, scriptName=os.path.basename(sys.argv[0])) + +f = open("opm/material/densead/EvaluationSpecializations.hpp", "w") +f.write(fileContents) +f.close() diff --git a/opm/material/densead/Evaluation.hpp b/opm/material/densead/Evaluation.hpp index fe574de40..2b06a72c7 100644 --- a/opm/material/densead/Evaluation.hpp +++ b/opm/material/densead/Evaluation.hpp @@ -25,10 +25,14 @@ * * \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 "Evaluation.hpp" #include "Math.hpp" #include @@ -38,16 +42,18 @@ #include #include #include +#include #include #include namespace Opm { namespace DenseAd { + /*! * \brief Represents a function evaluation and its derivatives w.r.t. a fixed set of * variables. */ -template +template class Evaluation { public: @@ -55,28 +61,28 @@ public: typedef ValueT ValueType; //! number of derivatives - static constexpr unsigned size = numVars; + static constexpr int size = numDerivs; protected: //! length of internal data vector - static constexpr unsigned length_ = numVars + 1 ; + static constexpr int length_ = size + 1; //! position index for value - static constexpr unsigned valuepos_ = 0; + static constexpr int valuepos_ = 0; //! start index for derivatives - static constexpr unsigned dstart_ = 1; + static constexpr int dstart_ = 1; //! end+1 index for derivatives - static constexpr unsigned dend_ = length_ ; -public: + static constexpr int dend_ = length_; +public: //! default constructor Evaluation() : data_() {} //! copy other function evaluation - Evaluation(const Evaluation& other) : data_( other.data_ ) - { - } + Evaluation(const Evaluation& other) + : data_(other.data_) + { } // create an evaluation which represents a constant function // @@ -95,12 +101,13 @@ public: // i.e., f(x) = c. this implies an evaluation with the given value and all // derivatives being zero. template - Evaluation(const RhsValueType& c, unsigned varPos) + Evaluation(const RhsValueType& c, int varPos) { + // The variable position must be in represented by the given variable descriptor + assert(0 <= varPos && varPos < size); + setValue( c ); clearDerivatives(); - // The variable position must be in represented by the given variable descriptor - assert(0 <= varPos && varPos < numVars); data_[varPos + dstart_] = 1.0; Valgrind::CheckDefined(data_); @@ -109,13 +116,14 @@ public: // set all derivatives to zero void clearDerivatives() { - for (unsigned i = dstart_; i < dend_; ++i) - data_[ i ] = 0.0; + for (int i = dstart_; i < dend_; ++i) { + data_[i] = 0.0; + } } // create a function evaluation for a "naked" depending variable (i.e., f(x) = x) template - static Evaluation createVariable(const RhsValueType& value, unsigned varPos) + static Evaluation createVariable(const RhsValueType& value, int varPos) { // copy function value and set all derivatives to 0, except for the variable // which is represented by the value (which is set to 1.0) @@ -135,25 +143,28 @@ public: { // print value os << "v: " << value() << " / d:"; + // print derivatives - for (unsigned 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 (unsigned varIdx = dstart_; varIdx < dend_; ++varIdx) - data_[ varIdx ] = other.data_[ varIdx ]; + for (int i = dstart_; i < dend_; ++i) { + data_[i] = other.data_[i]; + } } // add value and derivatives from other to this values and derivatives Evaluation& operator+=(const Evaluation& other) { - // value and derivatives are added - for (unsigned varIdx = 0; varIdx < length_; ++ varIdx) - data_[ varIdx ] += other.data_[ varIdx ]; + for (int i = 0; i < length_; ++i) { + data_[i] += other.data_[i]; + } return *this; } @@ -164,15 +175,16 @@ public: { // value is added, derivatives stay the same data_[valuepos_] += other; + return *this; } // subtract other's value and derivatives from this values Evaluation& operator-=(const Evaluation& other) { - // value and derivatives are subtracted - for (unsigned idx = 0 ; idx < length_ ; ++ idx) - data_[idx] -= other.data_[idx]; + for (int i = 0; i < length_; ++i) { + data_[i] -= other.data_[i]; + } return *this; } @@ -192,31 +204,32 @@ public: { // while the values are multiplied, the derivatives follow the product rule, // i.e., (u*v)' = (v'u + u'v). - ValueType& u = data_[ valuepos_ ]; - const ValueType& v = other.value(); - for (unsigned idx = dstart_; idx < dend_; ++idx) { - const ValueType& uPrime = data_[idx]; - const ValueType& vPrime = other.data_[idx]; + const ValueType u = this->value(); + const ValueType v = other.value(); - data_[idx] = (v*uPrime + u*vPrime); + // value + data_[valuepos_] *= v ; + + // derivatives + for (int i = dstart_; i < dend_; ++i) { + data_[i] = data_[i] * v + other.data_[i] * u; } - u *= v; return *this; } - // m(u*v)' = (v'u + u'v) + // m(c*u)' = c*u' template - Evaluation& operator*=(RhsValueType other) + Evaluation& operator*=(const RhsValueType& other) { - // values and derivatives are multiplied - for (unsigned idx = 0 ; idx < length_ ; ++ idx) - data_[idx] *= other; + for (int i = 0; i < length_; ++i) { + data_[i] *= other; + } return *this; } - // m(u*v)' = (v'u + u'v) + // m(u*v)' = (vu' - uv')/v^2 Evaluation& operator/=(const Evaluation& other) { // values are divided, derivatives follow the rule for division, i.e., (u/v)' = (v'u - @@ -234,14 +247,15 @@ public: return *this; } - // multiply value and derivatives by value of other + // divide value and derivatives by value of other template Evaluation& operator/=(const RhsValueType& other) { - // values and derivatives are divided - ValueType factor = (1.0/other); - for (unsigned idx = 0; idx < length_; ++idx) - data_[idx] *= factor; + const ValueType tmp = 1.0/other; + + for (int i = 0; i < length_; ++i) { + data_[i] *= tmp; + } return *this; } @@ -250,7 +264,9 @@ public: Evaluation operator+(const Evaluation& other) const { Evaluation result(*this); + result += other; + return result; } @@ -259,7 +275,9 @@ public: Evaluation operator+(const RhsValueType& other) const { Evaluation result(*this); + result += other; + return result; } @@ -267,7 +285,9 @@ public: Evaluation operator-(const Evaluation& other) const { Evaluation result(*this); + result -= other; + return result; } @@ -276,7 +296,9 @@ public: Evaluation operator-(const RhsValueType& other) const { Evaluation result(*this); + result -= other; + return result; } @@ -284,9 +306,11 @@ public: Evaluation operator-() const { Evaluation result; + // set value and derivatives to negative - for (unsigned idx = 0; idx < length_; ++idx) - result.data_[idx] = - data_[idx]; + for (int i = 0; i < length_; ++i) { + result.data_[i] = - data_[i]; + } return result; } @@ -294,7 +318,9 @@ public: Evaluation operator*(const Evaluation& other) const { Evaluation result(*this); + result *= other; + return result; } @@ -302,14 +328,18 @@ public: 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; } @@ -317,7 +347,9 @@ public: Evaluation operator/(const RhsValueType& other) const { Evaluation result(*this); + result /= other; + return result; } @@ -326,13 +358,17 @@ public: { setValue( other ); clearDerivatives(); + return *this; } // copy assignment from evaluation Evaluation& operator=(const Evaluation& other) { - data_ = other.data_; + for (int i = 0; i < length_; ++i) { + data_[i] = other.data_[i]; + } + return *this; } @@ -342,10 +378,11 @@ public: bool operator==(const Evaluation& other) const { - for (unsigned idx = 0; idx < length_; ++idx) - if (data_[idx] != other.data_[idx]) + for (int idx = 0; idx < length_; ++idx) { + if (data_[idx] != other.data_[idx]) { return false; - + } + } return true; } @@ -385,105 +422,90 @@ public: { return data_[valuepos_]; } // set value of variable - void setValue(const ValueType& val) + template + void setValue(const RhsValueType& val) { data_[valuepos_] = val; } // return varIdx'th derivative - const ValueType& derivative(unsigned varIdx) const + const ValueType& derivative(int varIdx) const { - assert(varIdx < numVars); - return data_[varIdx + dstart_]; + assert(0 <= varIdx && varIdx < size); + + return data_[dstart_ + varIdx]; } // set derivative at position varIdx - void setDerivative(unsigned varIdx, const ValueType& derVal) + void setDerivative(int varIdx, const ValueType& derVal) { - assert(varIdx < numVars); - data_[varIdx + dstart_] = derVal; + assert(0 <= varIdx && varIdx < size); + + data_[dstart_ + varIdx] = derVal; } -protected: - std::array data_; +private: + std::array data_; }; -template +// the generic operators are only required for the unspecialized case +template bool operator<(const RhsValueType& a, const Evaluation& b) { return b > a; } -template +template bool operator>(const RhsValueType& a, const Evaluation& b) { return b < a; } -template +template bool operator<=(const RhsValueType& a, const Evaluation& b) { return b >= a; } -template +template bool operator>=(const RhsValueType& a, const Evaluation& b) { return b <= a; } -template +template bool operator!=(const RhsValueType& a, const Evaluation& b) { return a != b.value(); } -template +template Evaluation operator+(const RhsValueType& a, const Evaluation& b) { Evaluation result(b); - result += a; - return result; } -template +template Evaluation operator-(const RhsValueType& a, const Evaluation& b) { - Evaluation result; - - result.setValue(a - b.value()); - for (unsigned varIdx = 0; varIdx < numVars; ++varIdx) - result.setDerivative(varIdx, - b.derivative(varIdx)); - + Evaluation result(a); + result -= b; return result; } -template +template Evaluation operator/(const RhsValueType& a, const Evaluation& b) { - Evaluation result; - - result.setValue(a/b.value()); - - // outer derivative - const ValueType& df_dg = - a/(b.value()*b.value()); - for (unsigned varIdx = 0; varIdx < numVars; ++varIdx) - result.setDerivative(varIdx, df_dg*b.derivative(varIdx)); - - return result; + Evaluation tmp(a); + tmp /= b; + return tmp; } -template +template Evaluation operator*(const RhsValueType& a, const Evaluation& b) { - Evaluation result; - - result.setValue(a*b.value()); - for (unsigned varIdx = 0; varIdx < numVars; ++varIdx) - result.setDerivative(varIdx, a*b.derivative(varIdx)); - + Evaluation result(b); + result *= a; return result; } -template +template std::ostream& operator<<(std::ostream& os, const Evaluation& eval) { os << eval.value(); 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. @@ -511,12 +533,12 @@ std::ostream& operator<<(std::ostream& os, const Evaluation& namespace Opm { namespace DenseAd { -template +template Evaluation abs(const Evaluation&); }} namespace std { -template +template const Opm::DenseAd::Evaluation abs(const Opm::DenseAd::Evaluation& x) { return Opm::DenseAd::abs(x); } @@ -535,7 +557,7 @@ const Opm::DenseAd::Evaluation abs(const Opm::DenseAd::Evalu #include namespace Dune { -template +template struct FieldTraits > { public: @@ -547,4 +569,6 @@ public: } // namespace Dune -#endif +#include "EvaluationSpecializations.hpp" + +#endif // OPM_DENSEAD_EVALUATION_HPP diff --git a/opm/material/densead/Evaluation1.hpp b/opm/material/densead/Evaluation1.hpp new file mode 100644 index 000000000..c17353c1c --- /dev/null +++ b/opm/material/densead/Evaluation1.hpp @@ -0,0 +1,431 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * + * \brief 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 "Evaluation.hpp" +#include "Math.hpp" + +#include + +#include + +#include +#include +#include +#include +#include +#include + +namespace Opm { +namespace DenseAd { + +template +class Evaluation +{ +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: + //! default constructor + Evaluation() : data_() + {} + + //! copy other function evaluation + Evaluation(const Evaluation& other) + { + data_[0] = other.data_[0]; + data_[1] = other.data_[1]; + } + + // 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 + 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 + Evaluation(const RhsValueType& c, int varPos) + { + // The variable position must be in represented by the given variable descriptor + assert(0 <= varPos && varPos < size); + + setValue( c ); + clearDerivatives(); + + data_[varPos + dstart_] = 1.0; + Valgrind::CheckDefined(data_); + } + + // set all derivatives to zero + void clearDerivatives() + { + data_[1] = 0.0; + } + + // create a function evaluation for a "naked" depending variable (i.e., f(x) = x) + template + 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 + 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) + { + data_[1] = other.data_[1]; + } + + + // add value and derivatives from other to this values and derivatives + Evaluation& operator+=(const Evaluation& other) + { + data_[0] += other.data_[0]; + data_[1] += other.data_[1]; + + return *this; + } + + // add value from other to this values + template + 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) + { + data_[0] -= other.data_[0]; + data_[1] -= other.data_[1]; + + return *this; + } + + // subtract other's value from this values + template + Evaluation& operator-=(const RhsValueType& other) + { + // for constants, values are subtracted, derivatives stay the same + data_[ valuepos_ ] -= other; + + return *this; + } + + // multiply values and apply chain rule to derivatives: (u*v)' = (v'u + u'v) + Evaluation& operator*=(const Evaluation& other) + { + // while the values are multiplied, the derivatives follow the product rule, + // i.e., (u*v)' = (v'u + u'v). + const ValueType u = this->value(); + const ValueType v = other.value(); + + // value + data_[valuepos_] *= v ; + + // derivatives + data_[1] = data_[1] * v + other.data_[1] * u; + + return *this; + } + + // m(c*u)' = c*u' + template + Evaluation& operator*=(const RhsValueType& other) + { + data_[0] *= other; + data_[1] *= other; + + return *this; + } + + // m(u*v)' = (vu' - uv')/v^2 + Evaluation& operator/=(const Evaluation& other) + { + // values are divided, derivatives follow the rule for division, i.e., (u/v)' = (v'u - + // u'v)/v^2. + ValueType& u = data_[ valuepos_ ]; + const ValueType& v = other.value(); + data_[1] = (v*data_[1] - u*other.data_[1])/(v*v); + u /= v; + + return *this; + } + + // divide value and derivatives by value of other + template + Evaluation& operator/=(const RhsValueType& other) + { + const ValueType tmp = 1.0/other; + + data_[0] *= tmp; + data_[1] *= tmp; + + return *this; + } + + // add two evaluation objects + Evaluation operator+(const Evaluation& other) const + { + Evaluation result(*this); + + result += other; + + return result; + } + + // add constant to this object + template + Evaluation operator+(const RhsValueType& other) const + { + Evaluation result(*this); + + result += other; + + return result; + } + + // subtract two evaluation objects + Evaluation operator-(const Evaluation& other) const + { + Evaluation result(*this); + + result -= other; + + return result; + } + + // subtract constant from evaluation object + template + Evaluation operator-(const RhsValueType& other) const + { + Evaluation result(*this); + + result -= other; + + return result; + } + + // negation (unary minus) operator + Evaluation operator-() const + { + Evaluation result; + + // set value and derivatives to negative + 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 + 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 + Evaluation operator/(const RhsValueType& other) const + { + Evaluation result(*this); + + result /= other; + + return result; + } + + template + 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 + 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 + bool operator>(RhsValueType other) const + { return value() > other; } + + bool operator>(const Evaluation& other) const + { return value() > other.value(); } + + template + bool operator<(RhsValueType other) const + { return value() < other; } + + bool operator<(const Evaluation& other) const + { return value() < other.value(); } + + template + bool operator>=(RhsValueType other) const + { return value() >= other; } + + bool operator>=(const Evaluation& other) const + { return value() >= other.value(); } + + template + 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 + 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 data_; +}; + +} } // namespace DenseAd, Opm + +#endif // OPM_DENSEAD_EVALUATION1_HPP diff --git a/opm/material/densead/Evaluation10.hpp b/opm/material/densead/Evaluation10.hpp new file mode 100644 index 000000000..088a20e57 --- /dev/null +++ b/opm/material/densead/Evaluation10.hpp @@ -0,0 +1,530 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * + * \brief 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 "Evaluation.hpp" +#include "Math.hpp" + +#include + +#include + +#include +#include +#include +#include +#include +#include + +namespace Opm { +namespace DenseAd { + +template +class Evaluation +{ +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: + //! default constructor + Evaluation() : data_() + {} + + //! copy other function evaluation + Evaluation(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]; + } + + // 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 + 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 + Evaluation(const RhsValueType& c, int varPos) + { + // The variable position must be in represented by the given variable descriptor + assert(0 <= varPos && varPos < size); + + setValue( c ); + clearDerivatives(); + + data_[varPos + dstart_] = 1.0; + Valgrind::CheckDefined(data_); + } + + // set all derivatives to zero + void clearDerivatives() + { + 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; + } + + // create a function evaluation for a "naked" depending variable (i.e., f(x) = x) + template + 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 + 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) + { + 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]; + } + + + // add value and derivatives from other to this values and derivatives + Evaluation& operator+=(const Evaluation& other) + { + data_[0] += other.data_[0]; + data_[1] += other.data_[1]; + data_[2] += other.data_[2]; + data_[3] += other.data_[3]; + data_[4] += other.data_[4]; + data_[5] += other.data_[5]; + data_[6] += other.data_[6]; + data_[7] += other.data_[7]; + data_[8] += other.data_[8]; + data_[9] += other.data_[9]; + data_[10] += other.data_[10]; + + return *this; + } + + // add value from other to this values + template + 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) + { + 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; + } + + // subtract other's value from this values + template + Evaluation& operator-=(const RhsValueType& other) + { + // for constants, values are subtracted, derivatives stay the same + data_[ valuepos_ ] -= other; + + return *this; + } + + // multiply values and apply chain rule to derivatives: (u*v)' = (v'u + u'v) + Evaluation& operator*=(const Evaluation& other) + { + // while the values are multiplied, the derivatives follow the product rule, + // i.e., (u*v)' = (v'u + u'v). + const ValueType u = this->value(); + const ValueType v = other.value(); + + // value + data_[valuepos_] *= v ; + + // derivatives + data_[1] = data_[1] * v + other.data_[1] * u; + data_[2] = data_[2] * v + other.data_[2] * u; + data_[3] = data_[3] * v + other.data_[3] * u; + data_[4] = data_[4] * v + other.data_[4] * u; + data_[5] = data_[5] * v + other.data_[5] * u; + data_[6] = data_[6] * v + other.data_[6] * u; + data_[7] = data_[7] * v + other.data_[7] * u; + data_[8] = data_[8] * v + other.data_[8] * u; + data_[9] = data_[9] * v + other.data_[9] * u; + data_[10] = data_[10] * v + other.data_[10] * u; + + return *this; + } + + // m(c*u)' = c*u' + template + Evaluation& operator*=(const RhsValueType& other) + { + data_[0] *= other; + data_[1] *= other; + data_[2] *= other; + data_[3] *= other; + data_[4] *= other; + data_[5] *= other; + data_[6] *= other; + data_[7] *= other; + data_[8] *= other; + data_[9] *= other; + data_[10] *= other; + + return *this; + } + + // m(u*v)' = (vu' - uv')/v^2 + Evaluation& operator/=(const Evaluation& other) + { + // values are divided, derivatives follow the rule for division, i.e., (u/v)' = (v'u - + // u'v)/v^2. + ValueType& u = data_[ valuepos_ ]; + const ValueType& v = other.value(); + data_[1] = (v*data_[1] - u*other.data_[1])/(v*v); + data_[2] = (v*data_[2] - u*other.data_[2])/(v*v); + data_[3] = (v*data_[3] - u*other.data_[3])/(v*v); + data_[4] = (v*data_[4] - u*other.data_[4])/(v*v); + data_[5] = (v*data_[5] - u*other.data_[5])/(v*v); + data_[6] = (v*data_[6] - u*other.data_[6])/(v*v); + data_[7] = (v*data_[7] - u*other.data_[7])/(v*v); + data_[8] = (v*data_[8] - u*other.data_[8])/(v*v); + data_[9] = (v*data_[9] - u*other.data_[9])/(v*v); + data_[10] = (v*data_[10] - u*other.data_[10])/(v*v); + u /= v; + + return *this; + } + + // divide value and derivatives by value of other + template + Evaluation& operator/=(const RhsValueType& other) + { + const 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; + } + + // add two evaluation objects + Evaluation operator+(const Evaluation& other) const + { + Evaluation result(*this); + + result += other; + + return result; + } + + // add constant to this object + template + Evaluation operator+(const RhsValueType& other) const + { + Evaluation result(*this); + + result += other; + + return result; + } + + // subtract two evaluation objects + Evaluation operator-(const Evaluation& other) const + { + Evaluation result(*this); + + result -= other; + + return result; + } + + // subtract constant from evaluation object + template + Evaluation operator-(const RhsValueType& other) const + { + Evaluation result(*this); + + result -= other; + + return result; + } + + // negation (unary minus) operator + Evaluation operator-() const + { + Evaluation result; + + // set value and derivatives to negative + 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 + 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 + Evaluation operator/(const RhsValueType& other) const + { + Evaluation result(*this); + + result /= other; + + return result; + } + + template + 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 + 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 + bool operator>(RhsValueType other) const + { return value() > other; } + + bool operator>(const Evaluation& other) const + { return value() > other.value(); } + + template + bool operator<(RhsValueType other) const + { return value() < other; } + + bool operator<(const Evaluation& other) const + { return value() < other.value(); } + + template + bool operator>=(RhsValueType other) const + { return value() >= other; } + + bool operator>=(const Evaluation& other) const + { return value() >= other.value(); } + + template + 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 + 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 data_; +}; + +} } // namespace DenseAd, Opm + +#endif // OPM_DENSEAD_EVALUATION10_HPP diff --git a/opm/material/densead/Evaluation11.hpp b/opm/material/densead/Evaluation11.hpp new file mode 100644 index 000000000..ce1743c31 --- /dev/null +++ b/opm/material/densead/Evaluation11.hpp @@ -0,0 +1,541 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * + * \brief 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 "Evaluation.hpp" +#include "Math.hpp" + +#include + +#include + +#include +#include +#include +#include +#include +#include + +namespace Opm { +namespace DenseAd { + +template +class Evaluation +{ +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: + //! default constructor + Evaluation() : data_() + {} + + //! copy other function evaluation + Evaluation(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]; + } + + // 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 + 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 + Evaluation(const RhsValueType& c, int varPos) + { + // The variable position must be in represented by the given variable descriptor + assert(0 <= varPos && varPos < size); + + setValue( c ); + clearDerivatives(); + + data_[varPos + dstart_] = 1.0; + Valgrind::CheckDefined(data_); + } + + // set all derivatives to zero + void clearDerivatives() + { + 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; + } + + // create a function evaluation for a "naked" depending variable (i.e., f(x) = x) + template + 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 + 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) + { + 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]; + } + + + // add value and derivatives from other to this values and derivatives + Evaluation& operator+=(const Evaluation& other) + { + data_[0] += other.data_[0]; + data_[1] += other.data_[1]; + data_[2] += other.data_[2]; + data_[3] += other.data_[3]; + data_[4] += other.data_[4]; + data_[5] += other.data_[5]; + data_[6] += other.data_[6]; + data_[7] += other.data_[7]; + data_[8] += other.data_[8]; + data_[9] += other.data_[9]; + data_[10] += other.data_[10]; + data_[11] += other.data_[11]; + + return *this; + } + + // add value from other to this values + template + 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) + { + 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; + } + + // subtract other's value from this values + template + Evaluation& operator-=(const RhsValueType& other) + { + // for constants, values are subtracted, derivatives stay the same + data_[ valuepos_ ] -= other; + + return *this; + } + + // multiply values and apply chain rule to derivatives: (u*v)' = (v'u + u'v) + Evaluation& operator*=(const Evaluation& other) + { + // while the values are multiplied, the derivatives follow the product rule, + // i.e., (u*v)' = (v'u + u'v). + const ValueType u = this->value(); + const ValueType v = other.value(); + + // value + data_[valuepos_] *= v ; + + // derivatives + data_[1] = data_[1] * v + other.data_[1] * u; + data_[2] = data_[2] * v + other.data_[2] * u; + data_[3] = data_[3] * v + other.data_[3] * u; + data_[4] = data_[4] * v + other.data_[4] * u; + data_[5] = data_[5] * v + other.data_[5] * u; + data_[6] = data_[6] * v + other.data_[6] * u; + data_[7] = data_[7] * v + other.data_[7] * u; + data_[8] = data_[8] * v + other.data_[8] * u; + data_[9] = data_[9] * v + other.data_[9] * u; + data_[10] = data_[10] * v + other.data_[10] * u; + data_[11] = data_[11] * v + other.data_[11] * u; + + return *this; + } + + // m(c*u)' = c*u' + template + Evaluation& operator*=(const RhsValueType& other) + { + data_[0] *= other; + data_[1] *= other; + data_[2] *= other; + data_[3] *= other; + data_[4] *= other; + data_[5] *= other; + data_[6] *= other; + data_[7] *= other; + data_[8] *= other; + data_[9] *= other; + data_[10] *= other; + data_[11] *= other; + + return *this; + } + + // m(u*v)' = (vu' - uv')/v^2 + Evaluation& operator/=(const Evaluation& other) + { + // values are divided, derivatives follow the rule for division, i.e., (u/v)' = (v'u - + // u'v)/v^2. + ValueType& u = data_[ valuepos_ ]; + const ValueType& v = other.value(); + data_[1] = (v*data_[1] - u*other.data_[1])/(v*v); + data_[2] = (v*data_[2] - u*other.data_[2])/(v*v); + data_[3] = (v*data_[3] - u*other.data_[3])/(v*v); + data_[4] = (v*data_[4] - u*other.data_[4])/(v*v); + data_[5] = (v*data_[5] - u*other.data_[5])/(v*v); + data_[6] = (v*data_[6] - u*other.data_[6])/(v*v); + data_[7] = (v*data_[7] - u*other.data_[7])/(v*v); + data_[8] = (v*data_[8] - u*other.data_[8])/(v*v); + data_[9] = (v*data_[9] - u*other.data_[9])/(v*v); + data_[10] = (v*data_[10] - u*other.data_[10])/(v*v); + data_[11] = (v*data_[11] - u*other.data_[11])/(v*v); + u /= v; + + return *this; + } + + // divide value and derivatives by value of other + template + Evaluation& operator/=(const RhsValueType& other) + { + const 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; + } + + // add two evaluation objects + Evaluation operator+(const Evaluation& other) const + { + Evaluation result(*this); + + result += other; + + return result; + } + + // add constant to this object + template + Evaluation operator+(const RhsValueType& other) const + { + Evaluation result(*this); + + result += other; + + return result; + } + + // subtract two evaluation objects + Evaluation operator-(const Evaluation& other) const + { + Evaluation result(*this); + + result -= other; + + return result; + } + + // subtract constant from evaluation object + template + Evaluation operator-(const RhsValueType& other) const + { + Evaluation result(*this); + + result -= other; + + return result; + } + + // negation (unary minus) operator + Evaluation operator-() const + { + Evaluation result; + + // set value and derivatives to negative + 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 + 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 + Evaluation operator/(const RhsValueType& other) const + { + Evaluation result(*this); + + result /= other; + + return result; + } + + template + 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 + 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 + bool operator>(RhsValueType other) const + { return value() > other; } + + bool operator>(const Evaluation& other) const + { return value() > other.value(); } + + template + bool operator<(RhsValueType other) const + { return value() < other; } + + bool operator<(const Evaluation& other) const + { return value() < other.value(); } + + template + bool operator>=(RhsValueType other) const + { return value() >= other; } + + bool operator>=(const Evaluation& other) const + { return value() >= other.value(); } + + template + 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 + 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 data_; +}; + +} } // namespace DenseAd, Opm + +#endif // OPM_DENSEAD_EVALUATION11_HPP diff --git a/opm/material/densead/Evaluation12.hpp b/opm/material/densead/Evaluation12.hpp new file mode 100644 index 000000000..15a6fe2de --- /dev/null +++ b/opm/material/densead/Evaluation12.hpp @@ -0,0 +1,552 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * + * \brief 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 "Evaluation.hpp" +#include "Math.hpp" + +#include + +#include + +#include +#include +#include +#include +#include +#include + +namespace Opm { +namespace DenseAd { + +template +class Evaluation +{ +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: + //! default constructor + Evaluation() : data_() + {} + + //! copy other function evaluation + Evaluation(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]; + } + + // 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 + 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 + Evaluation(const RhsValueType& c, int varPos) + { + // The variable position must be in represented by the given variable descriptor + assert(0 <= varPos && varPos < size); + + setValue( c ); + clearDerivatives(); + + data_[varPos + dstart_] = 1.0; + Valgrind::CheckDefined(data_); + } + + // set all derivatives to zero + void clearDerivatives() + { + 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; + } + + // create a function evaluation for a "naked" depending variable (i.e., f(x) = x) + template + 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 + 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) + { + data_[1] = other.data_[1]; + data_[2] = other.data_[2]; + data_[3] = other.data_[3]; + data_[4] = other.data_[4]; + data_[5] = other.data_[5]; + data_[6] = other.data_[6]; + data_[7] = other.data_[7]; + data_[8] = other.data_[8]; + data_[9] = other.data_[9]; + data_[10] = other.data_[10]; + data_[11] = other.data_[11]; + data_[12] = other.data_[12]; + } + + + // add value and derivatives from other to this values and derivatives + Evaluation& operator+=(const Evaluation& other) + { + data_[0] += other.data_[0]; + data_[1] += other.data_[1]; + data_[2] += other.data_[2]; + data_[3] += other.data_[3]; + data_[4] += other.data_[4]; + data_[5] += other.data_[5]; + data_[6] += other.data_[6]; + data_[7] += other.data_[7]; + data_[8] += other.data_[8]; + data_[9] += other.data_[9]; + data_[10] += other.data_[10]; + data_[11] += other.data_[11]; + data_[12] += other.data_[12]; + + return *this; + } + + // add value from other to this values + template + 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) + { + data_[0] -= other.data_[0]; + data_[1] -= other.data_[1]; + data_[2] -= other.data_[2]; + data_[3] -= other.data_[3]; + data_[4] -= other.data_[4]; + data_[5] -= other.data_[5]; + data_[6] -= other.data_[6]; + data_[7] -= other.data_[7]; + data_[8] -= other.data_[8]; + data_[9] -= other.data_[9]; + data_[10] -= other.data_[10]; + data_[11] -= other.data_[11]; + data_[12] -= other.data_[12]; + + return *this; + } + + // subtract other's value from this values + template + Evaluation& operator-=(const RhsValueType& other) + { + // for constants, values are subtracted, derivatives stay the same + data_[ valuepos_ ] -= other; + + return *this; + } + + // multiply values and apply chain rule to derivatives: (u*v)' = (v'u + u'v) + Evaluation& operator*=(const Evaluation& other) + { + // while the values are multiplied, the derivatives follow the product rule, + // i.e., (u*v)' = (v'u + u'v). + const ValueType u = this->value(); + const ValueType v = other.value(); + + // value + data_[valuepos_] *= v ; + + // derivatives + data_[1] = data_[1] * v + other.data_[1] * u; + data_[2] = data_[2] * v + other.data_[2] * u; + data_[3] = data_[3] * v + other.data_[3] * u; + data_[4] = data_[4] * v + other.data_[4] * u; + data_[5] = data_[5] * v + other.data_[5] * u; + data_[6] = data_[6] * v + other.data_[6] * u; + data_[7] = data_[7] * v + other.data_[7] * u; + data_[8] = data_[8] * v + other.data_[8] * u; + data_[9] = data_[9] * v + other.data_[9] * u; + data_[10] = data_[10] * v + other.data_[10] * u; + data_[11] = data_[11] * v + other.data_[11] * u; + data_[12] = data_[12] * v + other.data_[12] * u; + + return *this; + } + + // m(c*u)' = c*u' + template + Evaluation& operator*=(const RhsValueType& other) + { + data_[0] *= other; + data_[1] *= other; + data_[2] *= other; + data_[3] *= other; + data_[4] *= other; + data_[5] *= other; + data_[6] *= other; + data_[7] *= other; + data_[8] *= other; + data_[9] *= other; + data_[10] *= other; + data_[11] *= other; + data_[12] *= other; + + return *this; + } + + // m(u*v)' = (vu' - uv')/v^2 + Evaluation& operator/=(const Evaluation& other) + { + // values are divided, derivatives follow the rule for division, i.e., (u/v)' = (v'u - + // u'v)/v^2. + ValueType& u = data_[ valuepos_ ]; + const ValueType& v = other.value(); + data_[1] = (v*data_[1] - u*other.data_[1])/(v*v); + data_[2] = (v*data_[2] - u*other.data_[2])/(v*v); + data_[3] = (v*data_[3] - u*other.data_[3])/(v*v); + data_[4] = (v*data_[4] - u*other.data_[4])/(v*v); + data_[5] = (v*data_[5] - u*other.data_[5])/(v*v); + data_[6] = (v*data_[6] - u*other.data_[6])/(v*v); + data_[7] = (v*data_[7] - u*other.data_[7])/(v*v); + data_[8] = (v*data_[8] - u*other.data_[8])/(v*v); + data_[9] = (v*data_[9] - u*other.data_[9])/(v*v); + data_[10] = (v*data_[10] - u*other.data_[10])/(v*v); + data_[11] = (v*data_[11] - u*other.data_[11])/(v*v); + data_[12] = (v*data_[12] - u*other.data_[12])/(v*v); + u /= v; + + return *this; + } + + // divide value and derivatives by value of other + template + Evaluation& operator/=(const RhsValueType& other) + { + const 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; + } + + // add two evaluation objects + Evaluation operator+(const Evaluation& other) const + { + Evaluation result(*this); + + result += other; + + return result; + } + + // add constant to this object + template + Evaluation operator+(const RhsValueType& other) const + { + Evaluation result(*this); + + result += other; + + return result; + } + + // subtract two evaluation objects + Evaluation operator-(const Evaluation& other) const + { + Evaluation result(*this); + + result -= other; + + return result; + } + + // subtract constant from evaluation object + template + Evaluation operator-(const RhsValueType& other) const + { + Evaluation result(*this); + + result -= other; + + return result; + } + + // negation (unary minus) operator + Evaluation operator-() const + { + Evaluation result; + + // set value and derivatives to negative + 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 + 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 + Evaluation operator/(const RhsValueType& other) const + { + Evaluation result(*this); + + result /= other; + + return result; + } + + template + 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 + 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 + bool operator>(RhsValueType other) const + { return value() > other; } + + bool operator>(const Evaluation& other) const + { return value() > other.value(); } + + template + bool operator<(RhsValueType other) const + { return value() < other; } + + bool operator<(const Evaluation& other) const + { return value() < other.value(); } + + template + bool operator>=(RhsValueType other) const + { return value() >= other; } + + bool operator>=(const Evaluation& other) const + { return value() >= other.value(); } + + template + 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 + 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 data_; +}; + +} } // namespace DenseAd, Opm + +#endif // OPM_DENSEAD_EVALUATION12_HPP diff --git a/opm/material/densead/Evaluation2.hpp b/opm/material/densead/Evaluation2.hpp new file mode 100644 index 000000000..1c9b9efca --- /dev/null +++ b/opm/material/densead/Evaluation2.hpp @@ -0,0 +1,442 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * + * \brief 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 "Evaluation.hpp" +#include "Math.hpp" + +#include + +#include + +#include +#include +#include +#include +#include +#include + +namespace Opm { +namespace DenseAd { + +template +class Evaluation +{ +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: + //! default constructor + Evaluation() : data_() + {} + + //! copy other function evaluation + Evaluation(const Evaluation& other) + { + data_[0] = other.data_[0]; + data_[1] = other.data_[1]; + data_[2] = other.data_[2]; + } + + // 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 + 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 + Evaluation(const RhsValueType& c, int varPos) + { + // The variable position must be in represented by the given variable descriptor + assert(0 <= varPos && varPos < size); + + setValue( c ); + clearDerivatives(); + + data_[varPos + dstart_] = 1.0; + Valgrind::CheckDefined(data_); + } + + // set all derivatives to zero + void clearDerivatives() + { + data_[1] = 0.0; + data_[2] = 0.0; + } + + // create a function evaluation for a "naked" depending variable (i.e., f(x) = x) + template + 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 + 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) + { + data_[1] = other.data_[1]; + data_[2] = other.data_[2]; + } + + + // add value and derivatives from other to this values and derivatives + Evaluation& operator+=(const Evaluation& other) + { + data_[0] += other.data_[0]; + data_[1] += other.data_[1]; + data_[2] += other.data_[2]; + + return *this; + } + + // add value from other to this values + template + 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) + { + data_[0] -= other.data_[0]; + data_[1] -= other.data_[1]; + data_[2] -= other.data_[2]; + + return *this; + } + + // subtract other's value from this values + template + Evaluation& operator-=(const RhsValueType& other) + { + // for constants, values are subtracted, derivatives stay the same + data_[ valuepos_ ] -= other; + + return *this; + } + + // multiply values and apply chain rule to derivatives: (u*v)' = (v'u + u'v) + Evaluation& operator*=(const Evaluation& other) + { + // while the values are multiplied, the derivatives follow the product rule, + // i.e., (u*v)' = (v'u + u'v). + const ValueType u = this->value(); + const ValueType v = other.value(); + + // value + data_[valuepos_] *= v ; + + // derivatives + data_[1] = data_[1] * v + other.data_[1] * u; + data_[2] = data_[2] * v + other.data_[2] * u; + + return *this; + } + + // m(c*u)' = c*u' + template + Evaluation& operator*=(const RhsValueType& other) + { + data_[0] *= other; + data_[1] *= other; + data_[2] *= other; + + return *this; + } + + // m(u*v)' = (vu' - uv')/v^2 + Evaluation& operator/=(const Evaluation& other) + { + // values are divided, derivatives follow the rule for division, i.e., (u/v)' = (v'u - + // u'v)/v^2. + ValueType& u = data_[ valuepos_ ]; + const ValueType& v = other.value(); + data_[1] = (v*data_[1] - u*other.data_[1])/(v*v); + data_[2] = (v*data_[2] - u*other.data_[2])/(v*v); + u /= v; + + return *this; + } + + // divide value and derivatives by value of other + template + Evaluation& operator/=(const RhsValueType& other) + { + const ValueType tmp = 1.0/other; + + data_[0] *= tmp; + data_[1] *= tmp; + data_[2] *= tmp; + + return *this; + } + + // add two evaluation objects + Evaluation operator+(const Evaluation& other) const + { + Evaluation result(*this); + + result += other; + + return result; + } + + // add constant to this object + template + Evaluation operator+(const RhsValueType& other) const + { + Evaluation result(*this); + + result += other; + + return result; + } + + // subtract two evaluation objects + Evaluation operator-(const Evaluation& other) const + { + Evaluation result(*this); + + result -= other; + + return result; + } + + // subtract constant from evaluation object + template + Evaluation operator-(const RhsValueType& other) const + { + Evaluation result(*this); + + result -= other; + + return result; + } + + // negation (unary minus) operator + Evaluation operator-() const + { + Evaluation result; + + // set value and derivatives to negative + 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 + 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 + Evaluation operator/(const RhsValueType& other) const + { + Evaluation result(*this); + + result /= other; + + return result; + } + + template + 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 + 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 + bool operator>(RhsValueType other) const + { return value() > other; } + + bool operator>(const Evaluation& other) const + { return value() > other.value(); } + + template + bool operator<(RhsValueType other) const + { return value() < other; } + + bool operator<(const Evaluation& other) const + { return value() < other.value(); } + + template + bool operator>=(RhsValueType other) const + { return value() >= other; } + + bool operator>=(const Evaluation& other) const + { return value() >= other.value(); } + + template + 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 + 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 data_; +}; + +} } // namespace DenseAd, Opm + +#endif // OPM_DENSEAD_EVALUATION2_HPP diff --git a/opm/material/densead/Evaluation3.hpp b/opm/material/densead/Evaluation3.hpp new file mode 100644 index 000000000..c715acd08 --- /dev/null +++ b/opm/material/densead/Evaluation3.hpp @@ -0,0 +1,453 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * + * \brief 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 "Evaluation.hpp" +#include "Math.hpp" + +#include + +#include + +#include +#include +#include +#include +#include +#include + +namespace Opm { +namespace DenseAd { + +template +class Evaluation +{ +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: + //! default constructor + Evaluation() : data_() + {} + + //! copy other function evaluation + Evaluation(const Evaluation& other) + { + data_[0] = other.data_[0]; + data_[1] = other.data_[1]; + data_[2] = other.data_[2]; + data_[3] = other.data_[3]; + } + + // create an evaluation which represents a constant function + // + // i.e., f(x) = c. this implies an evaluation with the given value and all + // derivatives being zero. + template + 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 + Evaluation(const RhsValueType& c, int varPos) + { + // The variable position must be in represented by the given variable descriptor + assert(0 <= varPos && varPos < size); + + setValue( c ); + clearDerivatives(); + + data_[varPos + dstart_] = 1.0; + Valgrind::CheckDefined(data_); + } + + // set all derivatives to zero + void clearDerivatives() + { + data_[1] = 0.0; + data_[2] = 0.0; + data_[3] = 0.0; + } + + // create a function evaluation for a "naked" depending variable (i.e., f(x) = x) + template + 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 + 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) + { + data_[1] = other.data_[1]; + data_[2] = other.data_[2]; + data_[3] = other.data_[3]; + } + + + // add value and derivatives from other to this values and derivatives + Evaluation& operator+=(const Evaluation& other) + { + data_[0] += other.data_[0]; + data_[1] += other.data_[1]; + data_[2] += other.data_[2]; + data_[3] += other.data_[3]; + + return *this; + } + + // add value from other to this values + template + 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) + { + data_[0] -= other.data_[0]; + data_[1] -= other.data_[1]; + data_[2] -= other.data_[2]; + data_[3] -= other.data_[3]; + + return *this; + } + + // subtract other's value from this values + template + Evaluation& operator-=(const RhsValueType& other) + { + // for constants, values are subtracted, derivatives stay the same + data_[ valuepos_ ] -= other; + + return *this; + } + + // multiply values and apply chain rule to derivatives: (u*v)' = (v'u + u'v) + Evaluation& operator*=(const Evaluation& other) + { + // while the values are multiplied, the derivatives follow the product rule, + // i.e., (u*v)' = (v'u + u'v). + const ValueType u = this->value(); + const ValueType v = other.value(); + + // value + data_[valuepos_] *= v ; + + // derivatives + data_[1] = data_[1] * v + other.data_[1] * u; + data_[2] = data_[2] * v + other.data_[2] * u; + data_[3] = data_[3] * v + other.data_[3] * u; + + return *this; + } + + // m(c*u)' = c*u' + template + Evaluation& operator*=(const RhsValueType& other) + { + data_[0] *= other; + data_[1] *= other; + data_[2] *= other; + data_[3] *= other; + + return *this; + } + + // m(u*v)' = (vu' - uv')/v^2 + Evaluation& operator/=(const Evaluation& other) + { + // values are divided, derivatives follow the rule for division, i.e., (u/v)' = (v'u - + // u'v)/v^2. + ValueType& u = data_[ valuepos_ ]; + const ValueType& v = other.value(); + data_[1] = (v*data_[1] - u*other.data_[1])/(v*v); + data_[2] = (v*data_[2] - u*other.data_[2])/(v*v); + data_[3] = (v*data_[3] - u*other.data_[3])/(v*v); + u /= v; + + return *this; + } + + // divide value and derivatives by value of other + template + Evaluation& operator/=(const RhsValueType& other) + { + const ValueType tmp = 1.0/other; + + data_[0] *= tmp; + data_[1] *= tmp; + data_[2] *= tmp; + data_[3] *= tmp; + + return *this; + } + + // add two evaluation objects + Evaluation operator+(const Evaluation& other) const + { + Evaluation result(*this); + + result += other; + + return result; + } + + // add constant to this object + template + Evaluation operator+(const RhsValueType& other) const + { + Evaluation result(*this); + + result += other; + + return result; + } + + // subtract two evaluation objects + Evaluation operator-(const Evaluation& other) const + { + Evaluation result(*this); + + result -= other; + + return result; + } + + // subtract constant from evaluation object + template + Evaluation operator-(const RhsValueType& other) const + { + Evaluation result(*this); + + result -= other; + + return result; + } + + // negation (unary minus) operator + Evaluation operator-() const + { + Evaluation result; + + // set value and derivatives to negative + 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 + 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 + Evaluation operator/(const RhsValueType& other) const + { + Evaluation result(*this); + + result /= other; + + return result; + } + + template + 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 + 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 + bool operator>(RhsValueType other) const + { return value() > other; } + + bool operator>(const Evaluation& other) const + { return value() > other.value(); } + + template + bool operator<(RhsValueType other) const + { return value() < other; } + + bool operator<(const Evaluation& other) const + { return value() < other.value(); } + + template + bool operator>=(RhsValueType other) const + { return value() >= other; } + + bool operator>=(const Evaluation& other) const + { return value() >= other.value(); } + + template + 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 + 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 data_; +}; + +} } // namespace DenseAd, Opm + +#endif // OPM_DENSEAD_EVALUATION3_HPP diff --git a/opm/material/densead/Evaluation4.hpp b/opm/material/densead/Evaluation4.hpp new file mode 100644 index 000000000..978f37f4e --- /dev/null +++ b/opm/material/densead/Evaluation4.hpp @@ -0,0 +1,464 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * + * \brief 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 "Evaluation.hpp" +#include "Math.hpp" + +#include + +#include + +#include +#include +#include +#include +#include +#include + +namespace Opm { +namespace DenseAd { + +template +class Evaluation +{ +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: + //! default constructor + Evaluation() : data_() + {} + + //! copy other function evaluation + Evaluation(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]; + } + + // 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 + 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 + Evaluation(const RhsValueType& c, int varPos) + { + // The variable position must be in represented by the given variable descriptor + assert(0 <= varPos && varPos < size); + + setValue( c ); + clearDerivatives(); + + data_[varPos + dstart_] = 1.0; + Valgrind::CheckDefined(data_); + } + + // set all derivatives to zero + void clearDerivatives() + { + data_[1] = 0.0; + data_[2] = 0.0; + data_[3] = 0.0; + data_[4] = 0.0; + } + + // create a function evaluation for a "naked" depending variable (i.e., f(x) = x) + template + 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 + 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) + { + data_[1] = other.data_[1]; + data_[2] = other.data_[2]; + data_[3] = other.data_[3]; + data_[4] = other.data_[4]; + } + + + // add value and derivatives from other to this values and derivatives + Evaluation& operator+=(const Evaluation& other) + { + data_[0] += other.data_[0]; + data_[1] += other.data_[1]; + data_[2] += other.data_[2]; + data_[3] += other.data_[3]; + data_[4] += other.data_[4]; + + return *this; + } + + // add value from other to this values + template + 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) + { + 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; + } + + // subtract other's value from this values + template + Evaluation& operator-=(const RhsValueType& other) + { + // for constants, values are subtracted, derivatives stay the same + data_[ valuepos_ ] -= other; + + return *this; + } + + // multiply values and apply chain rule to derivatives: (u*v)' = (v'u + u'v) + Evaluation& operator*=(const Evaluation& other) + { + // while the values are multiplied, the derivatives follow the product rule, + // i.e., (u*v)' = (v'u + u'v). + const ValueType u = this->value(); + const ValueType v = other.value(); + + // value + data_[valuepos_] *= v ; + + // derivatives + data_[1] = data_[1] * v + other.data_[1] * u; + data_[2] = data_[2] * v + other.data_[2] * u; + data_[3] = data_[3] * v + other.data_[3] * u; + data_[4] = data_[4] * v + other.data_[4] * u; + + return *this; + } + + // m(c*u)' = c*u' + template + Evaluation& operator*=(const RhsValueType& other) + { + data_[0] *= other; + data_[1] *= other; + data_[2] *= other; + data_[3] *= other; + data_[4] *= other; + + return *this; + } + + // m(u*v)' = (vu' - uv')/v^2 + Evaluation& operator/=(const Evaluation& other) + { + // values are divided, derivatives follow the rule for division, i.e., (u/v)' = (v'u - + // u'v)/v^2. + ValueType& u = data_[ valuepos_ ]; + const ValueType& v = other.value(); + data_[1] = (v*data_[1] - u*other.data_[1])/(v*v); + data_[2] = (v*data_[2] - u*other.data_[2])/(v*v); + data_[3] = (v*data_[3] - u*other.data_[3])/(v*v); + data_[4] = (v*data_[4] - u*other.data_[4])/(v*v); + u /= v; + + return *this; + } + + // divide value and derivatives by value of other + template + Evaluation& operator/=(const RhsValueType& other) + { + const ValueType tmp = 1.0/other; + + data_[0] *= tmp; + data_[1] *= tmp; + data_[2] *= tmp; + data_[3] *= tmp; + data_[4] *= tmp; + + return *this; + } + + // add two evaluation objects + Evaluation operator+(const Evaluation& other) const + { + Evaluation result(*this); + + result += other; + + return result; + } + + // add constant to this object + template + Evaluation operator+(const RhsValueType& other) const + { + Evaluation result(*this); + + result += other; + + return result; + } + + // subtract two evaluation objects + Evaluation operator-(const Evaluation& other) const + { + Evaluation result(*this); + + result -= other; + + return result; + } + + // subtract constant from evaluation object + template + Evaluation operator-(const RhsValueType& other) const + { + Evaluation result(*this); + + result -= other; + + return result; + } + + // negation (unary minus) operator + Evaluation operator-() const + { + Evaluation result; + + // set value and derivatives to negative + 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 + 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 + Evaluation operator/(const RhsValueType& other) const + { + Evaluation result(*this); + + result /= other; + + return result; + } + + template + 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 + 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 + bool operator>(RhsValueType other) const + { return value() > other; } + + bool operator>(const Evaluation& other) const + { return value() > other.value(); } + + template + bool operator<(RhsValueType other) const + { return value() < other; } + + bool operator<(const Evaluation& other) const + { return value() < other.value(); } + + template + bool operator>=(RhsValueType other) const + { return value() >= other; } + + bool operator>=(const Evaluation& other) const + { return value() >= other.value(); } + + template + 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 + 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 data_; +}; + +} } // namespace DenseAd, Opm + +#endif // OPM_DENSEAD_EVALUATION4_HPP diff --git a/opm/material/densead/Evaluation5.hpp b/opm/material/densead/Evaluation5.hpp new file mode 100644 index 000000000..36577c46c --- /dev/null +++ b/opm/material/densead/Evaluation5.hpp @@ -0,0 +1,475 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * + * \brief 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 "Evaluation.hpp" +#include "Math.hpp" + +#include + +#include + +#include +#include +#include +#include +#include +#include + +namespace Opm { +namespace DenseAd { + +template +class Evaluation +{ +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: + //! default constructor + Evaluation() : data_() + {} + + //! copy other function evaluation + Evaluation(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]; + } + + // 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 + 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 + Evaluation(const RhsValueType& c, int varPos) + { + // The variable position must be in represented by the given variable descriptor + assert(0 <= varPos && varPos < size); + + setValue( c ); + clearDerivatives(); + + data_[varPos + dstart_] = 1.0; + Valgrind::CheckDefined(data_); + } + + // set all derivatives to zero + void clearDerivatives() + { + data_[1] = 0.0; + data_[2] = 0.0; + data_[3] = 0.0; + data_[4] = 0.0; + data_[5] = 0.0; + } + + // create a function evaluation for a "naked" depending variable (i.e., f(x) = x) + template + 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 + 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) + { + 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]; + } + + + // add value and derivatives from other to this values and derivatives + Evaluation& operator+=(const Evaluation& other) + { + data_[0] += other.data_[0]; + data_[1] += other.data_[1]; + data_[2] += other.data_[2]; + data_[3] += other.data_[3]; + data_[4] += other.data_[4]; + data_[5] += other.data_[5]; + + return *this; + } + + // add value from other to this values + template + 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) + { + 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; + } + + // subtract other's value from this values + template + Evaluation& operator-=(const RhsValueType& other) + { + // for constants, values are subtracted, derivatives stay the same + data_[ valuepos_ ] -= other; + + return *this; + } + + // multiply values and apply chain rule to derivatives: (u*v)' = (v'u + u'v) + Evaluation& operator*=(const Evaluation& other) + { + // while the values are multiplied, the derivatives follow the product rule, + // i.e., (u*v)' = (v'u + u'v). + const ValueType u = this->value(); + const ValueType v = other.value(); + + // value + data_[valuepos_] *= v ; + + // derivatives + data_[1] = data_[1] * v + other.data_[1] * u; + data_[2] = data_[2] * v + other.data_[2] * u; + data_[3] = data_[3] * v + other.data_[3] * u; + data_[4] = data_[4] * v + other.data_[4] * u; + data_[5] = data_[5] * v + other.data_[5] * u; + + return *this; + } + + // m(c*u)' = c*u' + template + Evaluation& operator*=(const RhsValueType& other) + { + data_[0] *= other; + data_[1] *= other; + data_[2] *= other; + data_[3] *= other; + data_[4] *= other; + data_[5] *= other; + + return *this; + } + + // m(u*v)' = (vu' - uv')/v^2 + Evaluation& operator/=(const Evaluation& other) + { + // values are divided, derivatives follow the rule for division, i.e., (u/v)' = (v'u - + // u'v)/v^2. + ValueType& u = data_[ valuepos_ ]; + const ValueType& v = other.value(); + data_[1] = (v*data_[1] - u*other.data_[1])/(v*v); + data_[2] = (v*data_[2] - u*other.data_[2])/(v*v); + data_[3] = (v*data_[3] - u*other.data_[3])/(v*v); + data_[4] = (v*data_[4] - u*other.data_[4])/(v*v); + data_[5] = (v*data_[5] - u*other.data_[5])/(v*v); + u /= v; + + return *this; + } + + // divide value and derivatives by value of other + template + Evaluation& operator/=(const RhsValueType& other) + { + const 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; + } + + // add two evaluation objects + Evaluation operator+(const Evaluation& other) const + { + Evaluation result(*this); + + result += other; + + return result; + } + + // add constant to this object + template + Evaluation operator+(const RhsValueType& other) const + { + Evaluation result(*this); + + result += other; + + return result; + } + + // subtract two evaluation objects + Evaluation operator-(const Evaluation& other) const + { + Evaluation result(*this); + + result -= other; + + return result; + } + + // subtract constant from evaluation object + template + Evaluation operator-(const RhsValueType& other) const + { + Evaluation result(*this); + + result -= other; + + return result; + } + + // negation (unary minus) operator + Evaluation operator-() const + { + Evaluation result; + + // set value and derivatives to negative + 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 + 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 + Evaluation operator/(const RhsValueType& other) const + { + Evaluation result(*this); + + result /= other; + + return result; + } + + template + 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 + 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 + bool operator>(RhsValueType other) const + { return value() > other; } + + bool operator>(const Evaluation& other) const + { return value() > other.value(); } + + template + bool operator<(RhsValueType other) const + { return value() < other; } + + bool operator<(const Evaluation& other) const + { return value() < other.value(); } + + template + bool operator>=(RhsValueType other) const + { return value() >= other; } + + bool operator>=(const Evaluation& other) const + { return value() >= other.value(); } + + template + 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 + 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 data_; +}; + +} } // namespace DenseAd, Opm + +#endif // OPM_DENSEAD_EVALUATION5_HPP diff --git a/opm/material/densead/Evaluation6.hpp b/opm/material/densead/Evaluation6.hpp new file mode 100644 index 000000000..eca218b8d --- /dev/null +++ b/opm/material/densead/Evaluation6.hpp @@ -0,0 +1,486 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * + * \brief 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 "Evaluation.hpp" +#include "Math.hpp" + +#include + +#include + +#include +#include +#include +#include +#include +#include + +namespace Opm { +namespace DenseAd { + +template +class Evaluation +{ +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: + //! default constructor + Evaluation() : data_() + {} + + //! copy other function evaluation + Evaluation(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]; + } + + // 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 + 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 + Evaluation(const RhsValueType& c, int varPos) + { + // The variable position must be in represented by the given variable descriptor + assert(0 <= varPos && varPos < size); + + setValue( c ); + clearDerivatives(); + + data_[varPos + dstart_] = 1.0; + Valgrind::CheckDefined(data_); + } + + // set all derivatives to zero + void clearDerivatives() + { + data_[1] = 0.0; + data_[2] = 0.0; + data_[3] = 0.0; + data_[4] = 0.0; + data_[5] = 0.0; + data_[6] = 0.0; + } + + // create a function evaluation for a "naked" depending variable (i.e., f(x) = x) + template + 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 + 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) + { + data_[1] = other.data_[1]; + data_[2] = other.data_[2]; + data_[3] = other.data_[3]; + data_[4] = other.data_[4]; + data_[5] = other.data_[5]; + data_[6] = other.data_[6]; + } + + + // add value and derivatives from other to this values and derivatives + Evaluation& operator+=(const Evaluation& other) + { + data_[0] += other.data_[0]; + data_[1] += other.data_[1]; + data_[2] += other.data_[2]; + data_[3] += other.data_[3]; + data_[4] += other.data_[4]; + data_[5] += other.data_[5]; + data_[6] += other.data_[6]; + + return *this; + } + + // add value from other to this values + template + 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) + { + data_[0] -= other.data_[0]; + data_[1] -= other.data_[1]; + data_[2] -= other.data_[2]; + data_[3] -= other.data_[3]; + data_[4] -= other.data_[4]; + data_[5] -= other.data_[5]; + data_[6] -= other.data_[6]; + + return *this; + } + + // subtract other's value from this values + template + Evaluation& operator-=(const RhsValueType& other) + { + // for constants, values are subtracted, derivatives stay the same + data_[ valuepos_ ] -= other; + + return *this; + } + + // multiply values and apply chain rule to derivatives: (u*v)' = (v'u + u'v) + Evaluation& operator*=(const Evaluation& other) + { + // while the values are multiplied, the derivatives follow the product rule, + // i.e., (u*v)' = (v'u + u'v). + const ValueType u = this->value(); + const ValueType v = other.value(); + + // value + data_[valuepos_] *= v ; + + // derivatives + data_[1] = data_[1] * v + other.data_[1] * u; + data_[2] = data_[2] * v + other.data_[2] * u; + data_[3] = data_[3] * v + other.data_[3] * u; + data_[4] = data_[4] * v + other.data_[4] * u; + data_[5] = data_[5] * v + other.data_[5] * u; + data_[6] = data_[6] * v + other.data_[6] * u; + + return *this; + } + + // m(c*u)' = c*u' + template + Evaluation& operator*=(const RhsValueType& other) + { + data_[0] *= other; + data_[1] *= other; + data_[2] *= other; + data_[3] *= other; + data_[4] *= other; + data_[5] *= other; + data_[6] *= other; + + return *this; + } + + // m(u*v)' = (vu' - uv')/v^2 + Evaluation& operator/=(const Evaluation& other) + { + // values are divided, derivatives follow the rule for division, i.e., (u/v)' = (v'u - + // u'v)/v^2. + ValueType& u = data_[ valuepos_ ]; + const ValueType& v = other.value(); + data_[1] = (v*data_[1] - u*other.data_[1])/(v*v); + data_[2] = (v*data_[2] - u*other.data_[2])/(v*v); + data_[3] = (v*data_[3] - u*other.data_[3])/(v*v); + data_[4] = (v*data_[4] - u*other.data_[4])/(v*v); + data_[5] = (v*data_[5] - u*other.data_[5])/(v*v); + data_[6] = (v*data_[6] - u*other.data_[6])/(v*v); + u /= v; + + return *this; + } + + // divide value and derivatives by value of other + template + Evaluation& operator/=(const RhsValueType& other) + { + const 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; + } + + // add two evaluation objects + Evaluation operator+(const Evaluation& other) const + { + Evaluation result(*this); + + result += other; + + return result; + } + + // add constant to this object + template + Evaluation operator+(const RhsValueType& other) const + { + Evaluation result(*this); + + result += other; + + return result; + } + + // subtract two evaluation objects + Evaluation operator-(const Evaluation& other) const + { + Evaluation result(*this); + + result -= other; + + return result; + } + + // subtract constant from evaluation object + template + Evaluation operator-(const RhsValueType& other) const + { + Evaluation result(*this); + + result -= other; + + return result; + } + + // negation (unary minus) operator + Evaluation operator-() const + { + Evaluation result; + + // set value and derivatives to negative + 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 + 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 + Evaluation operator/(const RhsValueType& other) const + { + Evaluation result(*this); + + result /= other; + + return result; + } + + template + 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 + 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 + bool operator>(RhsValueType other) const + { return value() > other; } + + bool operator>(const Evaluation& other) const + { return value() > other.value(); } + + template + bool operator<(RhsValueType other) const + { return value() < other; } + + bool operator<(const Evaluation& other) const + { return value() < other.value(); } + + template + bool operator>=(RhsValueType other) const + { return value() >= other; } + + bool operator>=(const Evaluation& other) const + { return value() >= other.value(); } + + template + 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 + 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 data_; +}; + +} } // namespace DenseAd, Opm + +#endif // OPM_DENSEAD_EVALUATION6_HPP diff --git a/opm/material/densead/Evaluation7.hpp b/opm/material/densead/Evaluation7.hpp new file mode 100644 index 000000000..8703ef515 --- /dev/null +++ b/opm/material/densead/Evaluation7.hpp @@ -0,0 +1,497 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * + * \brief 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 "Evaluation.hpp" +#include "Math.hpp" + +#include + +#include + +#include +#include +#include +#include +#include +#include + +namespace Opm { +namespace DenseAd { + +template +class Evaluation +{ +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: + //! default constructor + Evaluation() : data_() + {} + + //! copy other function evaluation + Evaluation(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]; + } + + // 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 + 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 + Evaluation(const RhsValueType& c, int varPos) + { + // The variable position must be in represented by the given variable descriptor + assert(0 <= varPos && varPos < size); + + setValue( c ); + clearDerivatives(); + + data_[varPos + dstart_] = 1.0; + Valgrind::CheckDefined(data_); + } + + // set all derivatives to zero + void clearDerivatives() + { + 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; + } + + // create a function evaluation for a "naked" depending variable (i.e., f(x) = x) + template + 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 + 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) + { + 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]; + } + + + // add value and derivatives from other to this values and derivatives + Evaluation& operator+=(const Evaluation& other) + { + data_[0] += other.data_[0]; + data_[1] += other.data_[1]; + data_[2] += other.data_[2]; + data_[3] += other.data_[3]; + data_[4] += other.data_[4]; + data_[5] += other.data_[5]; + data_[6] += other.data_[6]; + data_[7] += other.data_[7]; + + return *this; + } + + // add value from other to this values + template + 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) + { + 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; + } + + // subtract other's value from this values + template + Evaluation& operator-=(const RhsValueType& other) + { + // for constants, values are subtracted, derivatives stay the same + data_[ valuepos_ ] -= other; + + return *this; + } + + // multiply values and apply chain rule to derivatives: (u*v)' = (v'u + u'v) + Evaluation& operator*=(const Evaluation& other) + { + // while the values are multiplied, the derivatives follow the product rule, + // i.e., (u*v)' = (v'u + u'v). + const ValueType u = this->value(); + const ValueType v = other.value(); + + // value + data_[valuepos_] *= v ; + + // derivatives + data_[1] = data_[1] * v + other.data_[1] * u; + data_[2] = data_[2] * v + other.data_[2] * u; + data_[3] = data_[3] * v + other.data_[3] * u; + data_[4] = data_[4] * v + other.data_[4] * u; + data_[5] = data_[5] * v + other.data_[5] * u; + data_[6] = data_[6] * v + other.data_[6] * u; + data_[7] = data_[7] * v + other.data_[7] * u; + + return *this; + } + + // m(c*u)' = c*u' + template + Evaluation& operator*=(const RhsValueType& other) + { + data_[0] *= other; + data_[1] *= other; + data_[2] *= other; + data_[3] *= other; + data_[4] *= other; + data_[5] *= other; + data_[6] *= other; + data_[7] *= other; + + return *this; + } + + // m(u*v)' = (vu' - uv')/v^2 + Evaluation& operator/=(const Evaluation& other) + { + // values are divided, derivatives follow the rule for division, i.e., (u/v)' = (v'u - + // u'v)/v^2. + ValueType& u = data_[ valuepos_ ]; + const ValueType& v = other.value(); + data_[1] = (v*data_[1] - u*other.data_[1])/(v*v); + data_[2] = (v*data_[2] - u*other.data_[2])/(v*v); + data_[3] = (v*data_[3] - u*other.data_[3])/(v*v); + data_[4] = (v*data_[4] - u*other.data_[4])/(v*v); + data_[5] = (v*data_[5] - u*other.data_[5])/(v*v); + data_[6] = (v*data_[6] - u*other.data_[6])/(v*v); + data_[7] = (v*data_[7] - u*other.data_[7])/(v*v); + u /= v; + + return *this; + } + + // divide value and derivatives by value of other + template + Evaluation& operator/=(const RhsValueType& other) + { + const 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; + } + + // add two evaluation objects + Evaluation operator+(const Evaluation& other) const + { + Evaluation result(*this); + + result += other; + + return result; + } + + // add constant to this object + template + Evaluation operator+(const RhsValueType& other) const + { + Evaluation result(*this); + + result += other; + + return result; + } + + // subtract two evaluation objects + Evaluation operator-(const Evaluation& other) const + { + Evaluation result(*this); + + result -= other; + + return result; + } + + // subtract constant from evaluation object + template + Evaluation operator-(const RhsValueType& other) const + { + Evaluation result(*this); + + result -= other; + + return result; + } + + // negation (unary minus) operator + Evaluation operator-() const + { + Evaluation result; + + // set value and derivatives to negative + 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 + 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 + Evaluation operator/(const RhsValueType& other) const + { + Evaluation result(*this); + + result /= other; + + return result; + } + + template + 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 + 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 + bool operator>(RhsValueType other) const + { return value() > other; } + + bool operator>(const Evaluation& other) const + { return value() > other.value(); } + + template + bool operator<(RhsValueType other) const + { return value() < other; } + + bool operator<(const Evaluation& other) const + { return value() < other.value(); } + + template + bool operator>=(RhsValueType other) const + { return value() >= other; } + + bool operator>=(const Evaluation& other) const + { return value() >= other.value(); } + + template + 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 + 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 data_; +}; + +} } // namespace DenseAd, Opm + +#endif // OPM_DENSEAD_EVALUATION7_HPP diff --git a/opm/material/densead/Evaluation8.hpp b/opm/material/densead/Evaluation8.hpp new file mode 100644 index 000000000..5165f07b5 --- /dev/null +++ b/opm/material/densead/Evaluation8.hpp @@ -0,0 +1,508 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * + * \brief 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 "Evaluation.hpp" +#include "Math.hpp" + +#include + +#include + +#include +#include +#include +#include +#include +#include + +namespace Opm { +namespace DenseAd { + +template +class Evaluation +{ +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: + //! default constructor + Evaluation() : data_() + {} + + //! copy other function evaluation + Evaluation(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]; + } + + // 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 + 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 + Evaluation(const RhsValueType& c, int varPos) + { + // The variable position must be in represented by the given variable descriptor + assert(0 <= varPos && varPos < size); + + setValue( c ); + clearDerivatives(); + + data_[varPos + dstart_] = 1.0; + Valgrind::CheckDefined(data_); + } + + // set all derivatives to zero + void clearDerivatives() + { + 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; + } + + // create a function evaluation for a "naked" depending variable (i.e., f(x) = x) + template + 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 + 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) + { + 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]; + } + + + // add value and derivatives from other to this values and derivatives + Evaluation& operator+=(const Evaluation& other) + { + data_[0] += other.data_[0]; + data_[1] += other.data_[1]; + data_[2] += other.data_[2]; + data_[3] += other.data_[3]; + data_[4] += other.data_[4]; + data_[5] += other.data_[5]; + data_[6] += other.data_[6]; + data_[7] += other.data_[7]; + data_[8] += other.data_[8]; + + return *this; + } + + // add value from other to this values + template + 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) + { + 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; + } + + // subtract other's value from this values + template + Evaluation& operator-=(const RhsValueType& other) + { + // for constants, values are subtracted, derivatives stay the same + data_[ valuepos_ ] -= other; + + return *this; + } + + // multiply values and apply chain rule to derivatives: (u*v)' = (v'u + u'v) + Evaluation& operator*=(const Evaluation& other) + { + // while the values are multiplied, the derivatives follow the product rule, + // i.e., (u*v)' = (v'u + u'v). + const ValueType u = this->value(); + const ValueType v = other.value(); + + // value + data_[valuepos_] *= v ; + + // derivatives + data_[1] = data_[1] * v + other.data_[1] * u; + data_[2] = data_[2] * v + other.data_[2] * u; + data_[3] = data_[3] * v + other.data_[3] * u; + data_[4] = data_[4] * v + other.data_[4] * u; + data_[5] = data_[5] * v + other.data_[5] * u; + data_[6] = data_[6] * v + other.data_[6] * u; + data_[7] = data_[7] * v + other.data_[7] * u; + data_[8] = data_[8] * v + other.data_[8] * u; + + return *this; + } + + // m(c*u)' = c*u' + template + Evaluation& operator*=(const RhsValueType& other) + { + data_[0] *= other; + data_[1] *= other; + data_[2] *= other; + data_[3] *= other; + data_[4] *= other; + data_[5] *= other; + data_[6] *= other; + data_[7] *= other; + data_[8] *= other; + + return *this; + } + + // m(u*v)' = (vu' - uv')/v^2 + Evaluation& operator/=(const Evaluation& other) + { + // values are divided, derivatives follow the rule for division, i.e., (u/v)' = (v'u - + // u'v)/v^2. + ValueType& u = data_[ valuepos_ ]; + const ValueType& v = other.value(); + data_[1] = (v*data_[1] - u*other.data_[1])/(v*v); + data_[2] = (v*data_[2] - u*other.data_[2])/(v*v); + data_[3] = (v*data_[3] - u*other.data_[3])/(v*v); + data_[4] = (v*data_[4] - u*other.data_[4])/(v*v); + data_[5] = (v*data_[5] - u*other.data_[5])/(v*v); + data_[6] = (v*data_[6] - u*other.data_[6])/(v*v); + data_[7] = (v*data_[7] - u*other.data_[7])/(v*v); + data_[8] = (v*data_[8] - u*other.data_[8])/(v*v); + u /= v; + + return *this; + } + + // divide value and derivatives by value of other + template + Evaluation& operator/=(const RhsValueType& other) + { + const 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; + } + + // add two evaluation objects + Evaluation operator+(const Evaluation& other) const + { + Evaluation result(*this); + + result += other; + + return result; + } + + // add constant to this object + template + Evaluation operator+(const RhsValueType& other) const + { + Evaluation result(*this); + + result += other; + + return result; + } + + // subtract two evaluation objects + Evaluation operator-(const Evaluation& other) const + { + Evaluation result(*this); + + result -= other; + + return result; + } + + // subtract constant from evaluation object + template + Evaluation operator-(const RhsValueType& other) const + { + Evaluation result(*this); + + result -= other; + + return result; + } + + // negation (unary minus) operator + Evaluation operator-() const + { + Evaluation result; + + // set value and derivatives to negative + 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 + 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 + Evaluation operator/(const RhsValueType& other) const + { + Evaluation result(*this); + + result /= other; + + return result; + } + + template + 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 + 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 + bool operator>(RhsValueType other) const + { return value() > other; } + + bool operator>(const Evaluation& other) const + { return value() > other.value(); } + + template + bool operator<(RhsValueType other) const + { return value() < other; } + + bool operator<(const Evaluation& other) const + { return value() < other.value(); } + + template + bool operator>=(RhsValueType other) const + { return value() >= other; } + + bool operator>=(const Evaluation& other) const + { return value() >= other.value(); } + + template + 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 + 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 data_; +}; + +} } // namespace DenseAd, Opm + +#endif // OPM_DENSEAD_EVALUATION8_HPP diff --git a/opm/material/densead/Evaluation9.hpp b/opm/material/densead/Evaluation9.hpp new file mode 100644 index 000000000..fe8f7343a --- /dev/null +++ b/opm/material/densead/Evaluation9.hpp @@ -0,0 +1,519 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * + * \brief 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 "Evaluation.hpp" +#include "Math.hpp" + +#include + +#include + +#include +#include +#include +#include +#include +#include + +namespace Opm { +namespace DenseAd { + +template +class Evaluation +{ +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: + //! default constructor + Evaluation() : data_() + {} + + //! copy other function evaluation + Evaluation(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]; + } + + // 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 + 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 + Evaluation(const RhsValueType& c, int varPos) + { + // The variable position must be in represented by the given variable descriptor + assert(0 <= varPos && varPos < size); + + setValue( c ); + clearDerivatives(); + + data_[varPos + dstart_] = 1.0; + Valgrind::CheckDefined(data_); + } + + // set all derivatives to zero + void clearDerivatives() + { + 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; + } + + // create a function evaluation for a "naked" depending variable (i.e., f(x) = x) + template + 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 + 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) + { + 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]; + } + + + // add value and derivatives from other to this values and derivatives + Evaluation& operator+=(const Evaluation& other) + { + data_[0] += other.data_[0]; + data_[1] += other.data_[1]; + data_[2] += other.data_[2]; + data_[3] += other.data_[3]; + data_[4] += other.data_[4]; + data_[5] += other.data_[5]; + data_[6] += other.data_[6]; + data_[7] += other.data_[7]; + data_[8] += other.data_[8]; + data_[9] += other.data_[9]; + + return *this; + } + + // add value from other to this values + template + 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) + { + 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; + } + + // subtract other's value from this values + template + Evaluation& operator-=(const RhsValueType& other) + { + // for constants, values are subtracted, derivatives stay the same + data_[ valuepos_ ] -= other; + + return *this; + } + + // multiply values and apply chain rule to derivatives: (u*v)' = (v'u + u'v) + Evaluation& operator*=(const Evaluation& other) + { + // while the values are multiplied, the derivatives follow the product rule, + // i.e., (u*v)' = (v'u + u'v). + const ValueType u = this->value(); + const ValueType v = other.value(); + + // value + data_[valuepos_] *= v ; + + // derivatives + data_[1] = data_[1] * v + other.data_[1] * u; + data_[2] = data_[2] * v + other.data_[2] * u; + data_[3] = data_[3] * v + other.data_[3] * u; + data_[4] = data_[4] * v + other.data_[4] * u; + data_[5] = data_[5] * v + other.data_[5] * u; + data_[6] = data_[6] * v + other.data_[6] * u; + data_[7] = data_[7] * v + other.data_[7] * u; + data_[8] = data_[8] * v + other.data_[8] * u; + data_[9] = data_[9] * v + other.data_[9] * u; + + return *this; + } + + // m(c*u)' = c*u' + template + Evaluation& operator*=(const RhsValueType& other) + { + data_[0] *= other; + data_[1] *= other; + data_[2] *= other; + data_[3] *= other; + data_[4] *= other; + data_[5] *= other; + data_[6] *= other; + data_[7] *= other; + data_[8] *= other; + data_[9] *= other; + + return *this; + } + + // m(u*v)' = (vu' - uv')/v^2 + Evaluation& operator/=(const Evaluation& other) + { + // values are divided, derivatives follow the rule for division, i.e., (u/v)' = (v'u - + // u'v)/v^2. + ValueType& u = data_[ valuepos_ ]; + const ValueType& v = other.value(); + data_[1] = (v*data_[1] - u*other.data_[1])/(v*v); + data_[2] = (v*data_[2] - u*other.data_[2])/(v*v); + data_[3] = (v*data_[3] - u*other.data_[3])/(v*v); + data_[4] = (v*data_[4] - u*other.data_[4])/(v*v); + data_[5] = (v*data_[5] - u*other.data_[5])/(v*v); + data_[6] = (v*data_[6] - u*other.data_[6])/(v*v); + data_[7] = (v*data_[7] - u*other.data_[7])/(v*v); + data_[8] = (v*data_[8] - u*other.data_[8])/(v*v); + data_[9] = (v*data_[9] - u*other.data_[9])/(v*v); + u /= v; + + return *this; + } + + // divide value and derivatives by value of other + template + Evaluation& operator/=(const RhsValueType& other) + { + const 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; + } + + // add two evaluation objects + Evaluation operator+(const Evaluation& other) const + { + Evaluation result(*this); + + result += other; + + return result; + } + + // add constant to this object + template + Evaluation operator+(const RhsValueType& other) const + { + Evaluation result(*this); + + result += other; + + return result; + } + + // subtract two evaluation objects + Evaluation operator-(const Evaluation& other) const + { + Evaluation result(*this); + + result -= other; + + return result; + } + + // subtract constant from evaluation object + template + Evaluation operator-(const RhsValueType& other) const + { + Evaluation result(*this); + + result -= other; + + return result; + } + + // negation (unary minus) operator + Evaluation operator-() const + { + Evaluation result; + + // set value and derivatives to negative + 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 + 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 + Evaluation operator/(const RhsValueType& other) const + { + Evaluation result(*this); + + result /= other; + + return result; + } + + template + 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 + 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 + bool operator>(RhsValueType other) const + { return value() > other; } + + bool operator>(const Evaluation& other) const + { return value() > other.value(); } + + template + bool operator<(RhsValueType other) const + { return value() < other; } + + bool operator<(const Evaluation& other) const + { return value() < other.value(); } + + template + bool operator>=(RhsValueType other) const + { return value() >= other; } + + bool operator>=(const Evaluation& other) const + { return value() >= other.value(); } + + template + 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 + 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 data_; +}; + +} } // namespace DenseAd, Opm + +#endif // OPM_DENSEAD_EVALUATION9_HPP diff --git a/opm/material/densead/EvaluationSpecializations.hpp b/opm/material/densead/EvaluationSpecializations.hpp new file mode 100644 index 000000000..4d3666387 --- /dev/null +++ b/opm/material/densead/EvaluationSpecializations.hpp @@ -0,0 +1,47 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * + * \brief This file includes all specializations for the dense-AD Evaluation class. + * + * \attention THIS FILE GETS AUTOMATICALLY GENERATED BY THE "genEvalSpecializations.py" + * SCRIPT. DO NOT EDIT IT MANUALLY! + */ +#ifndef OPM_DENSEAD_EVALUATION_SPECIALIZATIONS_HPP +#define OPM_DENSEAD_EVALUATION_SPECIALIZATIONS_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#endif // OPM_DENSEAD_EVALUATION_SPECIALIZATIONS_HPP \ No newline at end of file diff --git a/opm/material/densead/Math.hpp b/opm/material/densead/Math.hpp index 0076f1989..29e8ef589 100644 --- a/opm/material/densead/Math.hpp +++ b/opm/material/densead/Math.hpp @@ -39,45 +39,45 @@ namespace Opm { namespace DenseAd { // forward declaration of the Evaluation template class -template +template class Evaluation; // provide some algebraic functions -template +template Evaluation abs(const Evaluation& x) { return (x > 0.0)?x:-x; } -template +template Evaluation min(const Evaluation& x1, const Evaluation& x2) { return (x1 < x2)?x1:x2; } -template +template Evaluation min(const Arg1ValueType& x1, const Evaluation& x2) { return (x1 < x2)?x1:x2; } -template +template Evaluation min(const Evaluation& x1, const Arg2ValueType& x2) { return min(x2, x1); } -template +template Evaluation max(const Evaluation& x1, const Evaluation& x2) { return (x1 > x2)?x1:x2; } -template +template Evaluation max(const Arg1ValueType& x1, const Evaluation& x2) { return (x1 > x2)?x1:x2; } -template +template Evaluation max(const Evaluation& x1, const Arg2ValueType& x2) { return max(x2, x1); } -template +template Evaluation tan(const Evaluation& x) { typedef MathToolbox ValueTypeToolbox; @@ -89,13 +89,13 @@ Evaluation tan(const Evaluation& x) // derivatives use the chain rule const ValueType& df_dx = 1 + tmp*tmp; - for (unsigned curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) + for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx)); return result; } -template +template Evaluation atan(const Evaluation& x) { typedef MathToolbox ValueTypeToolbox; @@ -106,13 +106,13 @@ Evaluation atan(const Evaluation& x) // derivatives use the chain rule const ValueType& df_dx = 1/(1 + x.value()*x.value()); - for (unsigned curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) + for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx)); return result; } -template +template Evaluation atan2(const Evaluation& x, const Evaluation& y) { @@ -124,7 +124,7 @@ Evaluation atan2(const Evaluation& x, // derivatives use the chain rule const ValueType& alpha = 1/(1 + (x.value()*x.value())/(y.value()*y.value())); - for (unsigned curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) { + for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) { result.setDerivative(curVarIdx, alpha/(y.value()*y.value()) *(x.derivative(curVarIdx)*y.value() - x.value()*y.derivative(curVarIdx))); @@ -133,7 +133,7 @@ Evaluation atan2(const Evaluation& x, return result; } -template +template Evaluation sin(const Evaluation& x) { typedef MathToolbox ValueTypeToolbox; @@ -144,13 +144,13 @@ Evaluation sin(const Evaluation& x) // derivatives use the chain rule const ValueType& df_dx = ValueTypeToolbox::cos(x.value()); - for (unsigned curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) + for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx)); return result; } -template +template Evaluation asin(const Evaluation& x) { typedef MathToolbox ValueTypeToolbox; @@ -161,13 +161,13 @@ Evaluation asin(const Evaluation& x) // derivatives use the chain rule const ValueType& df_dx = 1.0/ValueTypeToolbox::sqrt(1 - x.value()*x.value()); - for (unsigned curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) + for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx)); return result; } -template +template Evaluation cos(const Evaluation& x) { typedef MathToolbox ValueTypeToolbox; @@ -178,13 +178,13 @@ Evaluation cos(const Evaluation& x) // derivatives use the chain rule const ValueType& df_dx = -ValueTypeToolbox::sin(x.value()); - for (unsigned curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) + for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx)); return result; } -template +template Evaluation acos(const Evaluation& x) { typedef MathToolbox ValueTypeToolbox; @@ -195,13 +195,13 @@ Evaluation acos(const Evaluation& x) // derivatives use the chain rule const ValueType& df_dx = - 1.0/ValueTypeToolbox::sqrt(1 - x.value()*x.value()); - for (unsigned curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) + for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx)); return result; } -template +template Evaluation sqrt(const Evaluation& x) { typedef MathToolbox ValueTypeToolbox; @@ -213,14 +213,14 @@ Evaluation sqrt(const Evaluation& x) // derivatives use the chain rule ValueType df_dx = 0.5/sqrt_x; - for (unsigned curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) { + for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) { result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx)); } return result; } -template +template Evaluation exp(const Evaluation& x) { typedef MathToolbox ValueTypeToolbox; @@ -231,14 +231,14 @@ Evaluation exp(const Evaluation& x) // derivatives use the chain rule const ValueType& df_dx = exp_x; - for (unsigned curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) + for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx)); return result; } // exponentiation of arbitrary base with a fixed constant -template +template Evaluation pow(const Evaluation& base, const ExpType& exp) { @@ -256,7 +256,7 @@ Evaluation pow(const Evaluation& base, else { // derivatives use the chain rule const ValueType& df_dx = pow_x/base.value()*exp; - for (unsigned curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) + for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) result.setDerivative(curVarIdx, df_dx*base.derivative(curVarIdx)); } @@ -264,7 +264,7 @@ Evaluation pow(const Evaluation& base, } // exponentiation of constant base with an arbitrary exponent -template +template Evaluation pow(const BaseType& base, const Evaluation& exp) { @@ -283,7 +283,7 @@ Evaluation pow(const BaseType& base, // derivatives use the chain rule const ValueType& df_dx = lnBase*result.value(); - for (unsigned curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) + for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) result.setDerivative(curVarIdx, df_dx*exp.derivative(curVarIdx)); } @@ -292,7 +292,7 @@ Evaluation pow(const BaseType& base, // this is the most expensive power function. Computationally it is pretty expensive, so // one of the above two variants above should be preferred if possible. -template +template Evaluation pow(const Evaluation& base, const Evaluation& exp) { @@ -314,7 +314,7 @@ Evaluation pow(const Evaluation& base, const ValueType& f = base.value(); const ValueType& g = exp.value(); const ValueType& logF = ValueTypeToolbox::log(f); - for (unsigned curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) { + for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) { const ValueType& fPrime = base.derivative(curVarIdx); const ValueType& gPrime = exp.derivative(curVarIdx); result.setDerivative(curVarIdx, (g*fPrime/f + logF*gPrime) * valuePow); @@ -324,7 +324,7 @@ Evaluation pow(const Evaluation& base, return result; } -template +template Evaluation log(const Evaluation& x) { typedef MathToolbox ValueTypeToolbox; @@ -335,7 +335,7 @@ Evaluation log(const Evaluation& x) // derivatives use the chain rule const ValueType& df_dx = 1/x.value(); - for (unsigned curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) + for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx)); return result; @@ -345,7 +345,7 @@ Evaluation log(const Evaluation& x) // a kind of traits class for the automatic differentiation case. (The toolbox for the // scalar case is provided by the MathToolbox.hpp header file.) -template +template struct MathToolbox > { private: @@ -364,7 +364,7 @@ public: static Evaluation createConstant(ValueType value) { return Evaluation::createConstant(value); } - static Evaluation createVariable(ValueType value, unsigned varIdx) + static Evaluation createVariable(ValueType value, int varIdx) { return Evaluation::createVariable(value, varIdx); } template @@ -395,7 +395,7 @@ public: return false; // make sure that the derivatives are identical - for (unsigned curVarIdx = 0; curVarIdx < numVars; ++curVarIdx) + for (int curVarIdx = 0; curVarIdx < numVars; ++curVarIdx) if (!ValueTypeToolbox::isSame(a.derivative(curVarIdx), b.derivative(curVarIdx), tolerance)) return false; @@ -460,7 +460,7 @@ public: if (!InnerToolbox::isfinite(arg.value())) return false; - for (unsigned i = 0; i < numVars; ++i) + for (int i = 0; i < numVars; ++i) if (!InnerToolbox::isfinite(arg.derivative(i))) return false; @@ -472,7 +472,7 @@ public: if (InnerToolbox::isnan(arg.value())) return true; - for (unsigned i = 0; i < numVars; ++i) + for (int i = 0; i < numVars; ++i) if (InnerToolbox::isnan(arg.derivative(i))) return true; diff --git a/tests/test_densead.cpp b/tests/test_densead.cpp index 0d765b97c..5af2caaa9 100644 --- a/tests/test_densead.cpp +++ b/tests/test_densead.cpp @@ -50,9 +50,12 @@ #include #include -static const int numVars = 3; +//static const int numVars = 3; + +template +struct TestEnv +{ -template void testOperators(const Scalar tolerance) { typedef Opm::DenseAd::Evaluation Eval; @@ -256,7 +259,7 @@ void testOperators(const Scalar tolerance) } } -template +template void test1DFunction(AdFn* adFn, ClassicFn* classicFn, Scalar xMin = 1e-6, Scalar xMax = 1000) { typedef Opm::DenseAd::Evaluation Eval; @@ -289,8 +292,7 @@ void test1DFunction(AdFn* adFn, ClassicFn* classicFn, Scalar xMin = 1e-6, Scalar } } -template void test2DFunction1(AdFn* adFn, ClassicFn* classicFn, Scalar xMin, Scalar xMax, Scalar y) { @@ -326,8 +328,7 @@ void test2DFunction1(AdFn* adFn, ClassicFn* classicFn, Scalar xMin, Scalar xMax, } } -template void test2DFunction2(AdFn* adFn, ClassicFn* classicFn, Scalar x, Scalar yMin, Scalar yMax) { @@ -363,7 +364,6 @@ void test2DFunction2(AdFn* adFn, ClassicFn* classicFn, Scalar x, Scalar yMin, Sc } } -template void testPowBase(Scalar baseMin = 1e-2, Scalar baseMax = 100) { typedef Opm::DenseAd::Evaluation Eval; @@ -405,7 +405,6 @@ void testPowBase(Scalar baseMin = 1e-2, Scalar baseMax = 100) } } -template void testPowExp(Scalar expMin = -100, Scalar expMax = 100) { typedef Opm::DenseAd::Evaluation Eval; @@ -447,7 +446,6 @@ void testPowExp(Scalar expMin = -100, Scalar expMax = 100) } } -template void testAtan2() { typedef Opm::DenseAd::Evaluation Eval; @@ -497,16 +495,12 @@ void testAtan2() } // prototypes -double myScalarMin(double a, double b); -double myScalarMax(double a, double b); - -double myScalarMin(double a, double b) +static double myScalarMin(double a, double b) { return std::min(a, b); } -double myScalarMax(double a, double b) +static double myScalarMax(double a, double b) { return std::max(a, b); } -template inline void testAll() { // the following is commented out because it is supposed to produce a compiler @@ -516,82 +510,82 @@ inline void testAll() std::cout << "testing operators and constructors\n"; const Scalar eps = std::numeric_limits::epsilon()*1e3; - testOperators(eps); + testOperators(eps); std::cout << "testing min()\n"; - test2DFunction1(Opm::DenseAd::min, + test2DFunction1(Opm::DenseAd::min, myScalarMin, -1000, 1000, /*p=*/1.234); - test2DFunction2(Opm::DenseAd::min, + test2DFunction2(Opm::DenseAd::min, myScalarMin, /*T=*/1.234, -1000, 1000); std::cout << "testing max()\n"; - test2DFunction1(Opm::DenseAd::max, + test2DFunction1(Opm::DenseAd::max, myScalarMax, -1000, 1000, /*p=*/1.234); - test2DFunction2(Opm::DenseAd::max, + test2DFunction2(Opm::DenseAd::max, myScalarMax, /*T=*/1.234, -1000, 1000); std::cout << "testing pow()\n"; - testPowBase(); - testPowExp(); + testPowBase(); + testPowExp(); std::cout << "testing abs()\n"; - test1DFunction(Opm::DenseAd::abs, + test1DFunction(Opm::DenseAd::abs, static_cast(std::abs)); std::cout << "testing sqrt()\n"; - test1DFunction(Opm::DenseAd::sqrt, + test1DFunction(Opm::DenseAd::sqrt, static_cast(std::sqrt)); std::cout << "testing sin()\n"; - test1DFunction(Opm::DenseAd::sin, + test1DFunction(Opm::DenseAd::sin, static_cast(std::sin), 0, 2*M_PI); std::cout << "testing asin()\n"; - test1DFunction(Opm::DenseAd::asin, + test1DFunction(Opm::DenseAd::asin, static_cast(std::asin), -1.0, 1.0); std::cout << "testing cos()\n"; - test1DFunction(Opm::DenseAd::cos, + test1DFunction(Opm::DenseAd::cos, static_cast(std::cos), 0, 2*M_PI); std::cout << "testing acos()\n"; - test1DFunction(Opm::DenseAd::acos, + test1DFunction(Opm::DenseAd::acos, static_cast(std::acos), -1.0, 1.0); std::cout << "testing tan()\n"; - test1DFunction(Opm::DenseAd::tan, + test1DFunction(Opm::DenseAd::tan, static_cast(std::tan), -M_PI / 2 * 0.95, M_PI / 2 * 0.95); std::cout << "testing atan()\n"; - test1DFunction(Opm::DenseAd::atan, + test1DFunction(Opm::DenseAd::atan, static_cast(std::atan), -10*1000.0, 10*1000.0); std::cout << "testing atan2()\n"; - testAtan2(); + testAtan2(); std::cout << "testing exp()\n"; - test1DFunction(Opm::DenseAd::exp, + test1DFunction(Opm::DenseAd::exp, static_cast(std::exp), -100, 100); std::cout << "testing log()\n"; - test1DFunction(Opm::DenseAd::log, + test1DFunction(Opm::DenseAd::log, static_cast(std::log), 1e-6, 1e9); @@ -651,12 +645,17 @@ inline void testAll() } } +};//TestEnv + + int main(int argc, char **argv) { Dune::MPIHelper::instance(argc, argv); - testAll(); - testAll(); + TestEnv().testAll(); + TestEnv().testAll(); + TestEnv().testAll(); + TestEnv().testAll(); return 0; }