All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Evaluation.hpp
Go to the documentation of this file.
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 /*
4  This file is part of the Open Porous Media project (OPM).
5 
6  OPM is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 2 of the License, or
9  (at your option) any later version.
10 
11  OPM is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with OPM. If not, see <http://www.gnu.org/licenses/>.
18 
19  Consult the COPYING file in the top-level source directory of this
20  module for the precise wording of the license and the list of
21  copyright holders.
22 */
32 #ifndef OPM_DENSEAD_EVALUATION_HPP
33 #define OPM_DENSEAD_EVALUATION_HPP
34 
35 #include "Evaluation.hpp"
36 #include "Math.hpp"
37 
38 #include <opm/common/Valgrind.hpp>
39 
40 #include <dune/common/version.hh>
41 
42 #include <array>
43 #include <cmath>
44 #include <cassert>
45 #include <cstring>
46 #include <iostream>
47 #include <algorithm>
48 
49 namespace Opm {
50 namespace DenseAd {
51 
56 template <class ValueT, int numDerivs>
58 {
59 public:
61  typedef ValueT ValueType;
62 
64  static constexpr int size = numDerivs;
65 
66 protected:
68  static constexpr int length_ = size + 1;
69 
71  static constexpr int valuepos_ = 0;
73  static constexpr int dstart_ = 1;
75  static constexpr int dend_ = length_;
76 
77 public:
79  Evaluation() : data_()
80  {}
81 
83  Evaluation(const Evaluation& other) = default;
84 
85  // create an evaluation which represents a constant function
86  //
87  // i.e., f(x) = c. this implies an evaluation with the given value and all
88  // derivatives being zero.
89  template <class RhsValueType>
90  Evaluation(const RhsValueType& c)
91  {
92  setValue( c );
93  clearDerivatives();
94  Valgrind::CheckDefined( data_ );
95  }
96 
97  // create an evaluation which represents a constant function
98  //
99  // i.e., f(x) = c. this implies an evaluation with the given value and all
100  // derivatives being zero.
101  template <class RhsValueType>
102  Evaluation(const RhsValueType& c, int varPos)
103  {
104  // The variable position must be in represented by the given variable descriptor
105  assert(0 <= varPos && varPos < size);
106 
107  setValue( c );
108  clearDerivatives();
109 
110  data_[varPos + dstart_] = 1.0;
111  Valgrind::CheckDefined(data_);
112  }
113 
114  // set all derivatives to zero
115  void clearDerivatives()
116  {
117  for (int i = dstart_; i < dend_; ++i) {
118  data_[i] = 0.0;
119  }
120  }
121 
122  // create a function evaluation for a "naked" depending variable (i.e., f(x) = x)
123  template <class RhsValueType>
124  static Evaluation createVariable(const RhsValueType& value, int varPos)
125  {
126  // copy function value and set all derivatives to 0, except for the variable
127  // which is represented by the value (which is set to 1.0)
128  return Evaluation( value, varPos );
129  }
130 
131  // "evaluate" a constant function (i.e. a function that does not depend on the set of
132  // relevant variables, f(x) = c).
133  template <class RhsValueType>
134  static Evaluation createConstant(const RhsValueType& value)
135  {
136  return Evaluation( value );
137  }
138 
139  // print the value and the derivatives of the function evaluation
140  void print(std::ostream& os = std::cout) const
141  {
142  // print value
143  os << "v: " << value() << " / d:";
144 
145  // print derivatives
146  for (int varIdx = 0; varIdx < size; ++varIdx) {
147  os << " " << derivative(varIdx);
148  }
149  }
150 
151  // copy all derivatives from other
152  void copyDerivatives(const Evaluation& other)
153  {
154  for (int i = dstart_; i < dend_; ++i) {
155  data_[i] = other.data_[i];
156  }
157  }
158 
159 
160  // add value and derivatives from other to this values and derivatives
161  Evaluation& operator+=(const Evaluation& other)
162  {
163  for (int i = 0; i < length_; ++i) {
164  data_[i] += other.data_[i];
165  }
166 
167  return *this;
168  }
169 
170  // add value from other to this values
171  template <class RhsValueType>
172  Evaluation& operator+=(const RhsValueType& other)
173  {
174  // value is added, derivatives stay the same
175  data_[valuepos_] += other;
176 
177  return *this;
178  }
179 
180  // subtract other's value and derivatives from this values
181  Evaluation& operator-=(const Evaluation& other)
182  {
183  for (int i = 0; i < length_; ++i) {
184  data_[i] -= other.data_[i];
185  }
186 
187  return *this;
188  }
189 
190  // subtract other's value from this values
191  template <class RhsValueType>
192  Evaluation& operator-=(const RhsValueType& other)
193  {
194  // for constants, values are subtracted, derivatives stay the same
195  data_[ valuepos_ ] -= other;
196 
197  return *this;
198  }
199 
200  // multiply values and apply chain rule to derivatives: (u*v)' = (v'u + u'v)
201  Evaluation& operator*=(const Evaluation& other)
202  {
203  // while the values are multiplied, the derivatives follow the product rule,
204  // i.e., (u*v)' = (v'u + u'v).
205  const ValueType u = this->value();
206  const ValueType v = other.value();
207 
208  // value
209  data_[valuepos_] *= v ;
210 
211  // derivatives
212  for (int i = dstart_; i < dend_; ++i) {
213  data_[i] = data_[i] * v + other.data_[i] * u;
214  }
215 
216  return *this;
217  }
218 
219  // m(c*u)' = c*u'
220  template <class RhsValueType>
221  Evaluation& operator*=(const RhsValueType& other)
222  {
223  for (int i = 0; i < length_; ++i) {
224  data_[i] *= other;
225  }
226 
227  return *this;
228  }
229 
230  // m(u*v)' = (vu' - uv')/v^2
231  Evaluation& operator/=(const Evaluation& other)
232  {
233  // values are divided, derivatives follow the rule for division, i.e., (u/v)' = (v'u -
234  // u'v)/v^2.
235  ValueType& u = data_[ valuepos_ ];
236  const ValueType& v = other.value();
237  for (unsigned idx = dstart_; idx < dend_; ++idx) {
238  const ValueType& uPrime = data_[idx];
239  const ValueType& vPrime = other.data_[idx];
240 
241  data_[idx] = (v*uPrime - u*vPrime)/(v*v);
242  }
243  u /= v;
244 
245  return *this;
246  }
247 
248  // divide value and derivatives by value of other
249  template <class RhsValueType>
250  Evaluation& operator/=(const RhsValueType& other)
251  {
252  const ValueType tmp = 1.0/other;
253 
254  for (int i = 0; i < length_; ++i) {
255  data_[i] *= tmp;
256  }
257 
258  return *this;
259  }
260 
261  // add two evaluation objects
262  Evaluation operator+(const Evaluation& other) const
263  {
264  Evaluation result(*this);
265 
266  result += other;
267 
268  return result;
269  }
270 
271  // add constant to this object
272  template <class RhsValueType>
273  Evaluation operator+(const RhsValueType& other) const
274  {
275  Evaluation result(*this);
276 
277  result += other;
278 
279  return result;
280  }
281 
282  // subtract two evaluation objects
283  Evaluation operator-(const Evaluation& other) const
284  {
285  Evaluation result(*this);
286 
287  result -= other;
288 
289  return result;
290  }
291 
292  // subtract constant from evaluation object
293  template <class RhsValueType>
294  Evaluation operator-(const RhsValueType& other) const
295  {
296  Evaluation result(*this);
297 
298  result -= other;
299 
300  return result;
301  }
302 
303  // negation (unary minus) operator
304  Evaluation operator-() const
305  {
306  Evaluation result;
307 
308  // set value and derivatives to negative
309  for (int i = 0; i < length_; ++i) {
310  result.data_[i] = - data_[i];
311  }
312 
313  return result;
314  }
315 
316  Evaluation operator*(const Evaluation& other) const
317  {
318  Evaluation result(*this);
319 
320  result *= other;
321 
322  return result;
323  }
324 
325  template <class RhsValueType>
326  Evaluation operator*(const RhsValueType& other) const
327  {
328  Evaluation result(*this);
329 
330  result *= other;
331 
332  return result;
333  }
334 
335  Evaluation operator/(const Evaluation& other) const
336  {
337  Evaluation result(*this);
338 
339  result /= other;
340 
341  return result;
342  }
343 
344  template <class RhsValueType>
345  Evaluation operator/(const RhsValueType& other) const
346  {
347  Evaluation result(*this);
348 
349  result /= other;
350 
351  return result;
352  }
353 
354  template <class RhsValueType>
355  Evaluation& operator=(const RhsValueType& other)
356  {
357  setValue( other );
358  clearDerivatives();
359 
360  return *this;
361  }
362 
363  // copy assignment from evaluation
364  Evaluation& operator=(const Evaluation& other) = default;
365 
366  template <class RhsValueType>
367  bool operator==(const RhsValueType& other) const
368  { return value() == other; }
369 
370  bool operator==(const Evaluation& other) const
371  {
372  for (int idx = 0; idx < length_; ++idx) {
373  if (data_[idx] != other.data_[idx]) {
374  return false;
375  }
376  }
377  return true;
378  }
379 
380  bool operator!=(const Evaluation& other) const
381  { return !operator==(other); }
382 
383  template <class RhsValueType>
384  bool operator>(RhsValueType other) const
385  { return value() > other; }
386 
387  bool operator>(const Evaluation& other) const
388  { return value() > other.value(); }
389 
390  template <class RhsValueType>
391  bool operator<(RhsValueType other) const
392  { return value() < other; }
393 
394  bool operator<(const Evaluation& other) const
395  { return value() < other.value(); }
396 
397  template <class RhsValueType>
398  bool operator>=(RhsValueType other) const
399  { return value() >= other; }
400 
401  bool operator>=(const Evaluation& other) const
402  { return value() >= other.value(); }
403 
404  template <class RhsValueType>
405  bool operator<=(RhsValueType other) const
406  { return value() <= other; }
407 
408  bool operator<=(const Evaluation& other) const
409  { return value() <= other.value(); }
410 
411  // return value of variable
412  const ValueType& value() const
413  { return data_[valuepos_]; }
414 
415  // set value of variable
416  template <class RhsValueType>
417  void setValue(const RhsValueType& val)
418  { data_[valuepos_] = val; }
419 
420  // return varIdx'th derivative
421  const ValueType& derivative(int varIdx) const
422  {
423  assert(0 <= varIdx && varIdx < size);
424 
425  return data_[dstart_ + varIdx];
426  }
427 
428  // set derivative at position varIdx
429  void setDerivative(int varIdx, const ValueType& derVal)
430  {
431  assert(0 <= varIdx && varIdx < size);
432 
433  data_[dstart_ + varIdx] = derVal;
434  }
435 
436 private:
437  std::array<ValueT, length_> data_;
438 };
439 
440 // the generic operators are only required for the unspecialized case
441 template <class RhsValueType, class ValueType, int numVars>
442 bool operator<(const RhsValueType& a, const Evaluation<ValueType, numVars>& b)
443 { return b > a; }
444 
445 template <class RhsValueType, class ValueType, int numVars>
446 bool operator>(const RhsValueType& a, const Evaluation<ValueType, numVars>& b)
447 { return b < a; }
448 
449 template <class RhsValueType, class ValueType, int numVars>
450 bool operator<=(const RhsValueType& a, const Evaluation<ValueType, numVars>& b)
451 { return b >= a; }
452 
453 template <class RhsValueType, class ValueType, int numVars>
454 bool operator>=(const RhsValueType& a, const Evaluation<ValueType, numVars>& b)
455 { return b <= a; }
456 
457 template <class RhsValueType, class ValueType, int numVars>
458 bool operator!=(const RhsValueType& a, const Evaluation<ValueType, numVars>& b)
459 { return a != b.value(); }
460 
461 template <class RhsValueType, class ValueType, int numVars>
462 Evaluation<ValueType, numVars> operator+(const RhsValueType& a, const Evaluation<ValueType, numVars>& b)
463 {
464  Evaluation<ValueType, numVars> result(b);
465  result += a;
466  return result;
467 }
468 
469 template <class RhsValueType, class ValueType, int numVars>
470 Evaluation<ValueType, numVars> operator-(const RhsValueType& a, const Evaluation<ValueType, numVars>& b)
471 {
472  Evaluation<ValueType, numVars> result(a);
473  result -= b;
474  return result;
475 }
476 
477 template <class RhsValueType, class ValueType, int numVars>
478 Evaluation<ValueType, numVars> operator/(const RhsValueType& a, const Evaluation<ValueType, numVars>& b)
479 {
480  Evaluation<ValueType, numVars> tmp(a);
481  tmp /= b;
482  return tmp;
483 }
484 
485 template <class RhsValueType, class ValueType, int numVars>
486 Evaluation<ValueType, numVars> operator*(const RhsValueType& a, const Evaluation<ValueType, numVars>& b)
487 {
488  Evaluation<ValueType, numVars> result(b);
489  result *= a;
490  return result;
491 }
492 
493 template <class ValueType, int numVars>
494 std::ostream& operator<<(std::ostream& os, const Evaluation<ValueType, numVars>& eval)
495 {
496  os << eval.value();
497  return os;
498 }
499 } } // namespace DenseAd, Opm
500 
501 // In Dune 2.3, the Evaluation.hpp header must be included before the fmatrix.hh
502 // header. Dune 2.4+ does not suffer from this because of some c++-foo.
503 //
504 // for those who are wondering: in C++ function templates cannot be partially
505 // specialized, and function argument overloads must be known _before_ they are used. The
506 // latter is what we do for the 'Dune::fvmeta::absreal()' function.
507 //
508 // consider the following test program:
509 //
510 // double foo(double i)
511 // { return i; }
512 //
513 // void bar()
514 // { std::cout << foo(0) << "\n"; }
515 //
516 // int foo(int i)
517 // { return i + 1; }
518 //
519 // void foobar()
520 // { std::cout << foo(0) << "\n"; }
521 //
522 // this will print '0' for bar() and '1' for foobar()...
523 #if !(DUNE_VERSION_NEWER(DUNE_COMMON, 2,4))
524 
525 namespace Opm {
526 namespace DenseAd {
527 template <class ValueType, int numVars>
528 Evaluation<ValueType, numVars> abs(const Evaluation<ValueType, numVars>&);
529 }}
530 
531 namespace std {
532 template <class ValueType, int numVars>
534 { return Opm::DenseAd::abs(x); }
535 
536 } // namespace std
537 
538 #if defined DUNE_DENSEMATRIX_HH
539 #warning \
540  "Due to some C++ peculiarity regarding function overloads, the 'Evaluation.hpp'" \
541  "header file must be included before Dune's 'densematrix.hh' for Dune < 2.4. " \
542  "(If Evaluations are to be used in conjunction with a dense matrix.)"
543 #endif
544 
545 #endif
546 
547 // this makes the Dune matrix/vector classes happy...
548 #include <dune/common/ftraits.hh>
549 
550 namespace Dune {
551 template <class ValueType, int numVars>
552 struct FieldTraits<Opm::DenseAd::Evaluation<ValueType, numVars> >
553 {
554 public:
556  // setting real_type to field_type here potentially leads to slightly worse
557  // performance, but at least it makes things compile.
558  typedef field_type real_type;
559 };
560 
561 } // namespace Dune
562 
564 
565 #endif // OPM_DENSEAD_EVALUATION_HPP
Evaluation()
default constructor
Definition: Evaluation.hpp:79
A number of commonly used algebraic functions for the localized OPM automatic differentiation (AD) fr...
This file includes all specializations for the dense-AD Evaluation class.
ValueT ValueType
field type
Definition: Evaluation.hpp:61
static constexpr int valuepos_
position index for value
Definition: Evaluation.hpp:71
static constexpr int length_
length of internal data vector
Definition: Evaluation.hpp:68
static constexpr int dstart_
start index for derivatives
Definition: Evaluation.hpp:73
Represents a function evaluation and its derivatives w.r.t.
Definition: Evaluation.hpp:57
static constexpr int size
number of derivatives
Definition: Evaluation.hpp:64
static constexpr int dend_
end+1 index for derivatives
Definition: Evaluation.hpp:75
Representation of an evaluation of a function and its derivatives w.r.t.