00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00032 #ifndef OPM_DENSEAD_EVALUATION_HPP
00033 #define OPM_DENSEAD_EVALUATION_HPP
00034
00035 #include "Evaluation.hpp"
00036 #include "Math.hpp"
00037
00038 #include <opm/common/Valgrind.hpp>
00039
00040 #include <dune/common/version.hh>
00041
00042 #include <array>
00043 #include <cmath>
00044 #include <cassert>
00045 #include <cstring>
00046 #include <iostream>
00047 #include <algorithm>
00048
00049 namespace Opm {
00050 namespace DenseAd {
00051
00056 template <class ValueT, int numDerivs>
00057 class Evaluation
00058 {
00059 public:
00061 typedef ValueT ValueType;
00062
00064 static constexpr int size = numDerivs;
00065
00066 protected:
00068 static constexpr int length_ = size + 1;
00069
00071 static constexpr int valuepos_ = 0;
00073 static constexpr int dstart_ = 1;
00075 static constexpr int dend_ = length_;
00076
00077 public:
00079 Evaluation() : data_()
00080 {}
00081
00083 Evaluation(const Evaluation& other) = default;
00084
00085
00086
00087
00088
00089 template <class RhsValueType>
00090 Evaluation(const RhsValueType& c)
00091 {
00092 setValue( c );
00093 clearDerivatives();
00094 Valgrind::CheckDefined( data_ );
00095 }
00096
00097
00098
00099
00100
00101 template <class RhsValueType>
00102 Evaluation(const RhsValueType& c, int varPos)
00103 {
00104
00105 assert(0 <= varPos && varPos < size);
00106
00107 setValue( c );
00108 clearDerivatives();
00109
00110 data_[varPos + dstart_] = 1.0;
00111 Valgrind::CheckDefined(data_);
00112 }
00113
00114
00115 void clearDerivatives()
00116 {
00117 for (int i = dstart_; i < dend_; ++i) {
00118 data_[i] = 0.0;
00119 }
00120 }
00121
00122
00123 template <class RhsValueType>
00124 static Evaluation createVariable(const RhsValueType& value, int varPos)
00125 {
00126
00127
00128 return Evaluation( value, varPos );
00129 }
00130
00131
00132
00133 template <class RhsValueType>
00134 static Evaluation createConstant(const RhsValueType& value)
00135 {
00136 return Evaluation( value );
00137 }
00138
00139
00140 void print(std::ostream& os = std::cout) const
00141 {
00142
00143 os << "v: " << value() << " / d:";
00144
00145
00146 for (int varIdx = 0; varIdx < size; ++varIdx) {
00147 os << " " << derivative(varIdx);
00148 }
00149 }
00150
00151
00152 void copyDerivatives(const Evaluation& other)
00153 {
00154 for (int i = dstart_; i < dend_; ++i) {
00155 data_[i] = other.data_[i];
00156 }
00157 }
00158
00159
00160
00161 Evaluation& operator+=(const Evaluation& other)
00162 {
00163 for (int i = 0; i < length_; ++i) {
00164 data_[i] += other.data_[i];
00165 }
00166
00167 return *this;
00168 }
00169
00170
00171 template <class RhsValueType>
00172 Evaluation& operator+=(const RhsValueType& other)
00173 {
00174
00175 data_[valuepos_] += other;
00176
00177 return *this;
00178 }
00179
00180
00181 Evaluation& operator-=(const Evaluation& other)
00182 {
00183 for (int i = 0; i < length_; ++i) {
00184 data_[i] -= other.data_[i];
00185 }
00186
00187 return *this;
00188 }
00189
00190
00191 template <class RhsValueType>
00192 Evaluation& operator-=(const RhsValueType& other)
00193 {
00194
00195 data_[ valuepos_ ] -= other;
00196
00197 return *this;
00198 }
00199
00200
00201 Evaluation& operator*=(const Evaluation& other)
00202 {
00203
00204
00205 const ValueType u = this->value();
00206 const ValueType v = other.value();
00207
00208
00209 data_[valuepos_] *= v ;
00210
00211
00212 for (int i = dstart_; i < dend_; ++i) {
00213 data_[i] = data_[i] * v + other.data_[i] * u;
00214 }
00215
00216 return *this;
00217 }
00218
00219
00220 template <class RhsValueType>
00221 Evaluation& operator*=(const RhsValueType& other)
00222 {
00223 for (int i = 0; i < length_; ++i) {
00224 data_[i] *= other;
00225 }
00226
00227 return *this;
00228 }
00229
00230
00231 Evaluation& operator/=(const Evaluation& other)
00232 {
00233
00234
00235 ValueType& u = data_[ valuepos_ ];
00236 const ValueType& v = other.value();
00237 for (unsigned idx = dstart_; idx < dend_; ++idx) {
00238 const ValueType& uPrime = data_[idx];
00239 const ValueType& vPrime = other.data_[idx];
00240
00241 data_[idx] = (v*uPrime - u*vPrime)/(v*v);
00242 }
00243 u /= v;
00244
00245 return *this;
00246 }
00247
00248
00249 template <class RhsValueType>
00250 Evaluation& operator/=(const RhsValueType& other)
00251 {
00252 const ValueType tmp = 1.0/other;
00253
00254 for (int i = 0; i < length_; ++i) {
00255 data_[i] *= tmp;
00256 }
00257
00258 return *this;
00259 }
00260
00261
00262 Evaluation operator+(const Evaluation& other) const
00263 {
00264 Evaluation result(*this);
00265
00266 result += other;
00267
00268 return result;
00269 }
00270
00271
00272 template <class RhsValueType>
00273 Evaluation operator+(const RhsValueType& other) const
00274 {
00275 Evaluation result(*this);
00276
00277 result += other;
00278
00279 return result;
00280 }
00281
00282
00283 Evaluation operator-(const Evaluation& other) const
00284 {
00285 Evaluation result(*this);
00286
00287 result -= other;
00288
00289 return result;
00290 }
00291
00292
00293 template <class RhsValueType>
00294 Evaluation operator-(const RhsValueType& other) const
00295 {
00296 Evaluation result(*this);
00297
00298 result -= other;
00299
00300 return result;
00301 }
00302
00303
00304 Evaluation operator-() const
00305 {
00306 Evaluation result;
00307
00308
00309 for (int i = 0; i < length_; ++i) {
00310 result.data_[i] = - data_[i];
00311 }
00312
00313 return result;
00314 }
00315
00316 Evaluation operator*(const Evaluation& other) const
00317 {
00318 Evaluation result(*this);
00319
00320 result *= other;
00321
00322 return result;
00323 }
00324
00325 template <class RhsValueType>
00326 Evaluation operator*(const RhsValueType& other) const
00327 {
00328 Evaluation result(*this);
00329
00330 result *= other;
00331
00332 return result;
00333 }
00334
00335 Evaluation operator/(const Evaluation& other) const
00336 {
00337 Evaluation result(*this);
00338
00339 result /= other;
00340
00341 return result;
00342 }
00343
00344 template <class RhsValueType>
00345 Evaluation operator/(const RhsValueType& other) const
00346 {
00347 Evaluation result(*this);
00348
00349 result /= other;
00350
00351 return result;
00352 }
00353
00354 template <class RhsValueType>
00355 Evaluation& operator=(const RhsValueType& other)
00356 {
00357 setValue( other );
00358 clearDerivatives();
00359
00360 return *this;
00361 }
00362
00363
00364 Evaluation& operator=(const Evaluation& other) = default;
00365
00366 template <class RhsValueType>
00367 bool operator==(const RhsValueType& other) const
00368 { return value() == other; }
00369
00370 bool operator==(const Evaluation& other) const
00371 {
00372 for (int idx = 0; idx < length_; ++idx) {
00373 if (data_[idx] != other.data_[idx]) {
00374 return false;
00375 }
00376 }
00377 return true;
00378 }
00379
00380 bool operator!=(const Evaluation& other) const
00381 { return !operator==(other); }
00382
00383 template <class RhsValueType>
00384 bool operator>(RhsValueType other) const
00385 { return value() > other; }
00386
00387 bool operator>(const Evaluation& other) const
00388 { return value() > other.value(); }
00389
00390 template <class RhsValueType>
00391 bool operator<(RhsValueType other) const
00392 { return value() < other; }
00393
00394 bool operator<(const Evaluation& other) const
00395 { return value() < other.value(); }
00396
00397 template <class RhsValueType>
00398 bool operator>=(RhsValueType other) const
00399 { return value() >= other; }
00400
00401 bool operator>=(const Evaluation& other) const
00402 { return value() >= other.value(); }
00403
00404 template <class RhsValueType>
00405 bool operator<=(RhsValueType other) const
00406 { return value() <= other; }
00407
00408 bool operator<=(const Evaluation& other) const
00409 { return value() <= other.value(); }
00410
00411
00412 const ValueType& value() const
00413 { return data_[valuepos_]; }
00414
00415
00416 template <class RhsValueType>
00417 void setValue(const RhsValueType& val)
00418 { data_[valuepos_] = val; }
00419
00420
00421 const ValueType& derivative(int varIdx) const
00422 {
00423 assert(0 <= varIdx && varIdx < size);
00424
00425 return data_[dstart_ + varIdx];
00426 }
00427
00428
00429 void setDerivative(int varIdx, const ValueType& derVal)
00430 {
00431 assert(0 <= varIdx && varIdx < size);
00432
00433 data_[dstart_ + varIdx] = derVal;
00434 }
00435
00436 private:
00437 std::array<ValueT, length_> data_;
00438 };
00439
00440
00441 template <class RhsValueType, class ValueType, int numVars>
00442 bool operator<(const RhsValueType& a, const Evaluation<ValueType, numVars>& b)
00443 { return b > a; }
00444
00445 template <class RhsValueType, class ValueType, int numVars>
00446 bool operator>(const RhsValueType& a, const Evaluation<ValueType, numVars>& b)
00447 { return b < a; }
00448
00449 template <class RhsValueType, class ValueType, int numVars>
00450 bool operator<=(const RhsValueType& a, const Evaluation<ValueType, numVars>& b)
00451 { return b >= a; }
00452
00453 template <class RhsValueType, class ValueType, int numVars>
00454 bool operator>=(const RhsValueType& a, const Evaluation<ValueType, numVars>& b)
00455 { return b <= a; }
00456
00457 template <class RhsValueType, class ValueType, int numVars>
00458 bool operator!=(const RhsValueType& a, const Evaluation<ValueType, numVars>& b)
00459 { return a != b.value(); }
00460
00461 template <class RhsValueType, class ValueType, int numVars>
00462 Evaluation<ValueType, numVars> operator+(const RhsValueType& a, const Evaluation<ValueType, numVars>& b)
00463 {
00464 Evaluation<ValueType, numVars> result(b);
00465 result += a;
00466 return result;
00467 }
00468
00469 template <class RhsValueType, class ValueType, int numVars>
00470 Evaluation<ValueType, numVars> operator-(const RhsValueType& a, const Evaluation<ValueType, numVars>& b)
00471 {
00472 Evaluation<ValueType, numVars> result(a);
00473 result -= b;
00474 return result;
00475 }
00476
00477 template <class RhsValueType, class ValueType, int numVars>
00478 Evaluation<ValueType, numVars> operator/(const RhsValueType& a, const Evaluation<ValueType, numVars>& b)
00479 {
00480 Evaluation<ValueType, numVars> tmp(a);
00481 tmp /= b;
00482 return tmp;
00483 }
00484
00485 template <class RhsValueType, class ValueType, int numVars>
00486 Evaluation<ValueType, numVars> operator*(const RhsValueType& a, const Evaluation<ValueType, numVars>& b)
00487 {
00488 Evaluation<ValueType, numVars> result(b);
00489 result *= a;
00490 return result;
00491 }
00492
00493 template <class ValueType, int numVars>
00494 std::ostream& operator<<(std::ostream& os, const Evaluation<ValueType, numVars>& eval)
00495 {
00496 os << eval.value();
00497 return os;
00498 }
00499 } }
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523 #if !(DUNE_VERSION_NEWER(DUNE_COMMON, 2,4))
00524
00525 namespace Opm {
00526 namespace DenseAd {
00527 template <class ValueType, int numVars>
00528 Evaluation<ValueType, numVars> abs(const Evaluation<ValueType, numVars>&);
00529 }}
00530
00531 namespace std {
00532 template <class ValueType, int numVars>
00533 const Opm::DenseAd::Evaluation<ValueType, numVars> abs(const Opm::DenseAd::Evaluation<ValueType, numVars>& x)
00534 { return Opm::DenseAd::abs(x); }
00535
00536 }
00537
00538 #if defined DUNE_DENSEMATRIX_HH
00539 #warning \
00540 "Due to some C++ peculiarity regarding function overloads, the 'Evaluation.hpp'" \
00541 "header file must be included before Dune's 'densematrix.hh' for Dune < 2.4. " \
00542 "(If Evaluations are to be used in conjunction with a dense matrix.)"
00543 #endif
00544
00545 #endif
00546
00547
00548 #include <dune/common/ftraits.hh>
00549
00550 namespace Dune {
00551 template <class ValueType, int numVars>
00552 struct FieldTraits<Opm::DenseAd::Evaluation<ValueType, numVars> >
00553 {
00554 public:
00555 typedef Opm::DenseAd::Evaluation<ValueType, numVars> field_type;
00556
00557
00558 typedef field_type real_type;
00559 };
00560
00561 }
00562
00563 #include "EvaluationSpecializations.hpp"
00564
00565 #endif // OPM_DENSEAD_EVALUATION_HPP