My Project
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 
39 
40 #include <array>
41 #include <cmath>
42 #include <cassert>
43 #include <cstring>
44 #include <iostream>
45 #include <algorithm>
46 
47 namespace Opm {
48 namespace DenseAd {
51 static constexpr int DynamicSize = -1;
52 
57 template <class ValueT, int numDerivs, unsigned staticSize = 0>
59 {
60 public:
63  static const int numVars = numDerivs;
64 
66  typedef ValueT ValueType;
67 
69  constexpr int size() const
70  { return numDerivs; }
71 
72 protected:
74  constexpr int length_() const
75  { return size() + 1; }
76 
77 
79  constexpr int valuepos_() const
80  { return 0; }
82  constexpr int dstart_() const
83  { return 1; }
85  constexpr int dend_() const
86  { return length_(); }
87 
90  void checkDefined_() const
91  {
92 #ifndef NDEBUG
93  for (const auto& v: data_)
94  Valgrind::CheckDefined(v);
95 #endif
96  }
97 
98 public:
100  Evaluation() : data_()
101  {}
102 
104  Evaluation(const Evaluation& other) = default;
105 
106 
107  // create an evaluation which represents a constant function
108  //
109  // i.e., f(x) = c. this implies an evaluation with the given value and all
110  // derivatives being zero.
111  template <class RhsValueType>
112  Evaluation(const RhsValueType& c)
113  {
114  setValue(c);
115  clearDerivatives();
116 
117  checkDefined_();
118  }
119 
120  // create an evaluation which represents a constant function
121  //
122  // i.e., f(x) = c. this implies an evaluation with the given value and all
123  // derivatives being zero.
124  template <class RhsValueType>
125  Evaluation(const RhsValueType& c, int varPos)
126  {
127  // The variable position must be in represented by the given variable descriptor
128  assert(0 <= varPos && varPos < size());
129 
130  setValue( c );
131  clearDerivatives();
132 
133  data_[varPos + dstart_()] = 1.0;
134 
135  checkDefined_();
136  }
137 
138  // set all derivatives to zero
139  void clearDerivatives()
140  {
141  for (int i = dstart_(); i < dend_(); ++i)
142  data_[i] = 0.0;
143  }
144 
145  // create an uninitialized Evaluation object that is compatible with the
146  // argument, but not initialized
147  //
148  // This basically boils down to the copy constructor without copying
149  // anything. If the number of derivatives is known at compile time, this
150  // is equivalent to creating an uninitialized object using the default
151  // constructor, while for dynamic evaluations, it creates an Evaluation
152  // object which exhibits the same number of derivatives as the argument.
153  static Evaluation createBlank(const Evaluation&)
154  { return Evaluation(); }
155 
156  // create an Evaluation with value and all the derivatives to be zero
157  static Evaluation createConstantZero(const Evaluation&)
158  { return Evaluation(0.); }
159 
160  // create an Evaluation with value to be one and all the derivatives to be zero
161  static Evaluation createConstantOne(const Evaluation&)
162  { return Evaluation(1.); }
163 
164  // create a function evaluation for a "naked" depending variable (i.e., f(x) = x)
165  template <class RhsValueType>
166  static Evaluation createVariable(const RhsValueType& value, int varPos)
167  {
168  // copy function value and set all derivatives to 0, except for the variable
169  // which is represented by the value (which is set to 1.0)
170  return Evaluation(value, varPos);
171  }
172 
173  template <class RhsValueType>
174  static Evaluation createVariable(int nVars, const RhsValueType& value, int varPos)
175  {
176  if (nVars != 0)
177  throw std::logic_error("This statically-sized evaluation can only represent objects"
178  " with 0 derivatives");
179 
180  // copy function value and set all derivatives to 0, except for the variable
181  // which is represented by the value (which is set to 1.0)
182  return Evaluation(nVars, value, varPos);
183  }
184 
185  template <class RhsValueType>
186  static Evaluation createVariable(const Evaluation&, const RhsValueType& value, int varPos)
187  {
188  // copy function value and set all derivatives to 0, except for the variable
189  // which is represented by the value (which is set to 1.0)
190  return Evaluation(value, varPos);
191  }
192 
193 
194  // "evaluate" a constant function (i.e. a function that does not depend on the set of
195  // relevant variables, f(x) = c).
196  template <class RhsValueType>
197  static Evaluation createConstant(int nVars, const RhsValueType& value)
198  {
199  if (nVars != 0)
200  throw std::logic_error("This statically-sized evaluation can only represent objects"
201  " with 0 derivatives");
202  return Evaluation(value);
203  }
204 
205  // "evaluate" a constant function (i.e. a function that does not depend on the set of
206  // relevant variables, f(x) = c).
207  template <class RhsValueType>
208  static Evaluation createConstant(const RhsValueType& value)
209  {
210  return Evaluation(value);
211  }
212 
213  // "evaluate" a constant function (i.e. a function that does not depend on the set of
214  // relevant variables, f(x) = c).
215  template <class RhsValueType>
216  static Evaluation createConstant(const Evaluation&, const RhsValueType& value)
217  {
218  return Evaluation(value);
219  }
220 
221  // print the value and the derivatives of the function evaluation
222  void print(std::ostream& os = std::cout) const
223  {
224  // print value
225  os << "v: " << value() << " / d:";
226 
227  // print derivatives
228  for (int varIdx = 0; varIdx < size(); ++varIdx) {
229  os << " " << derivative(varIdx);
230  }
231  }
232 
233  // copy all derivatives from other
234  void copyDerivatives(const Evaluation& other)
235  {
236  assert(size() == other.size());
237 
238  for (int i = dstart_(); i < dend_(); ++i)
239  data_[i] = other.data_[i];
240  }
241 
242 
243  // add value and derivatives from other to this values and derivatives
244  Evaluation& operator+=(const Evaluation& other)
245  {
246  assert(size() == other.size());
247 
248  for (int i = 0; i < length_(); ++i)
249  data_[i] += other.data_[i];
250 
251  return *this;
252  }
253 
254  // add value from other to this values
255  template <class RhsValueType>
256  Evaluation& operator+=(const RhsValueType& other)
257  {
258  // value is added, derivatives stay the same
259  data_[valuepos_()] += other;
260 
261  return *this;
262  }
263 
264  // subtract other's value and derivatives from this values
265  Evaluation& operator-=(const Evaluation& other)
266  {
267  assert(size() == other.size());
268 
269  for (int i = 0; i < length_(); ++i)
270  data_[i] -= other.data_[i];
271 
272  return *this;
273  }
274 
275  // subtract other's value from this values
276  template <class RhsValueType>
277  Evaluation& operator-=(const RhsValueType& other)
278  {
279  // for constants, values are subtracted, derivatives stay the same
280  data_[valuepos_()] -= other;
281 
282  return *this;
283  }
284 
285  // multiply values and apply chain rule to derivatives: (u*v)' = (v'u + u'v)
286  Evaluation& operator*=(const Evaluation& other)
287  {
288  assert(size() == other.size());
289 
290  // while the values are multiplied, the derivatives follow the product rule,
291  // i.e., (u*v)' = (v'u + u'v).
292  const ValueType u = this->value();
293  const ValueType v = other.value();
294 
295  // value
296  data_[valuepos_()] *= v ;
297 
298  // derivatives
299  for (int i = dstart_(); i < dend_(); ++i)
300  data_[i] = data_[i] * v + other.data_[i] * u;
301 
302  return *this;
303  }
304 
305  // m(c*u)' = c*u'
306  template <class RhsValueType>
307  Evaluation& operator*=(const RhsValueType& other)
308  {
309  for (int i = 0; i < length_(); ++i)
310  data_[i] *= other;
311 
312  return *this;
313  }
314 
315  // m(u*v)' = (vu' - uv')/v^2
316  Evaluation& operator/=(const Evaluation& other)
317  {
318  assert(size() == other.size());
319 
320  // values are divided, derivatives follow the rule for division, i.e., (u/v)' = (v'u -
321  // u'v)/v^2.
322  ValueType& u = data_[valuepos_()];
323  const ValueType& v = other.value();
324  for (int idx = dstart_(); idx < dend_(); ++idx) {
325  const ValueType& uPrime = data_[idx];
326  const ValueType& vPrime = other.data_[idx];
327 
328  data_[idx] = (v*uPrime - u*vPrime)/(v*v);
329  }
330  u /= v;
331 
332  return *this;
333  }
334 
335  // divide value and derivatives by value of other
336  template <class RhsValueType>
337  Evaluation& operator/=(const RhsValueType& other)
338  {
339  const ValueType tmp = 1.0/other;
340 
341  for (int i = 0; i < length_(); ++i)
342  data_[i] *= tmp;
343 
344  return *this;
345  }
346 
347  // add two evaluation objects
348  Evaluation operator+(const Evaluation& other) const
349  {
350  assert(size() == other.size());
351 
352  Evaluation result(*this);
353 
354  result += other;
355 
356  return result;
357  }
358 
359  // add constant to this object
360  template <class RhsValueType>
361  Evaluation operator+(const RhsValueType& other) const
362  {
363  Evaluation result(*this);
364 
365  result += other;
366 
367  return result;
368  }
369 
370  // subtract two evaluation objects
371  Evaluation operator-(const Evaluation& other) const
372  {
373  assert(size() == other.size());
374 
375  Evaluation result(*this);
376 
377  result -= other;
378 
379  return result;
380  }
381 
382  // subtract constant from evaluation object
383  template <class RhsValueType>
384  Evaluation operator-(const RhsValueType& other) const
385  {
386  Evaluation result(*this);
387 
388  result -= other;
389 
390  return result;
391  }
392 
393  // negation (unary minus) operator
394  Evaluation operator-() const
395  {
396  Evaluation result;
397 
398  // set value and derivatives to negative
399  for (int i = 0; i < length_(); ++i)
400  result.data_[i] = - data_[i];
401 
402  return result;
403  }
404 
405  Evaluation operator*(const Evaluation& other) const
406  {
407  assert(size() == other.size());
408 
409  Evaluation result(*this);
410 
411  result *= other;
412 
413  return result;
414  }
415 
416  template <class RhsValueType>
417  Evaluation operator*(const RhsValueType& other) const
418  {
419  Evaluation result(*this);
420 
421  result *= other;
422 
423  return result;
424  }
425 
426  Evaluation operator/(const Evaluation& other) const
427  {
428  assert(size() == other.size());
429 
430  Evaluation result(*this);
431 
432  result /= other;
433 
434  return result;
435  }
436 
437  template <class RhsValueType>
438  Evaluation operator/(const RhsValueType& other) const
439  {
440  Evaluation result(*this);
441 
442  result /= other;
443 
444  return result;
445  }
446 
447  template <class RhsValueType>
448  Evaluation& operator=(const RhsValueType& other)
449  {
450  setValue( other );
451  clearDerivatives();
452 
453  return *this;
454  }
455 
456  // copy assignment from evaluation
457  Evaluation& operator=(const Evaluation& other) = default;
458 
459  template <class RhsValueType>
460  bool operator==(const RhsValueType& other) const
461  { return value() == other; }
462 
463  bool operator==(const Evaluation& other) const
464  {
465  assert(size() == other.size());
466 
467  for (int idx = 0; idx < length_(); ++idx) {
468  if (data_[idx] != other.data_[idx]) {
469  return false;
470  }
471  }
472  return true;
473  }
474 
475  bool operator!=(const Evaluation& other) const
476  { return !operator==(other); }
477 
478  template <class RhsValueType>
479  bool operator!=(const RhsValueType& other) const
480  { return !operator==(other); }
481 
482  template <class RhsValueType>
483  bool operator>(RhsValueType other) const
484  { return value() > other; }
485 
486  bool operator>(const Evaluation& other) const
487  {
488  assert(size() == other.size());
489 
490  return value() > other.value();
491  }
492 
493  template <class RhsValueType>
494  bool operator<(RhsValueType other) const
495  { return value() < other; }
496 
497  bool operator<(const Evaluation& other) const
498  {
499  assert(size() == other.size());
500 
501  return value() < other.value();
502  }
503 
504  template <class RhsValueType>
505  bool operator>=(RhsValueType other) const
506  { return value() >= other; }
507 
508  bool operator>=(const Evaluation& other) const
509  {
510  assert(size() == other.size());
511 
512  return value() >= other.value();
513  }
514 
515  template <class RhsValueType>
516  bool operator<=(RhsValueType other) const
517  { return value() <= other; }
518 
519  bool operator<=(const Evaluation& other) const
520  {
521  assert(size() == other.size());
522 
523  return value() <= other.value();
524  }
525 
526  // return value of variable
527  const ValueType& value() const
528  { return data_[valuepos_()]; }
529 
530  // set value of variable
531  template <class RhsValueType>
532  void setValue(const RhsValueType& val)
533  { data_[valuepos_()] = val; }
534 
535  // return varIdx'th derivative
536  const ValueType& derivative(int varIdx) const
537  {
538  assert(0 <= varIdx && varIdx < size());
539 
540  return data_[dstart_() + varIdx];
541  }
542 
543  // set derivative at position varIdx
544  void setDerivative(int varIdx, const ValueType& derVal)
545  {
546  assert(0 <= varIdx && varIdx < size());
547 
548  data_[dstart_() + varIdx] = derVal;
549  }
550 
551 private:
552 
553  std::array<ValueT, numDerivs + 1> data_;
554 };
555 
556 // the generic operators are only required for the unspecialized case
557 template <class RhsValueType, class ValueType, int numVars, unsigned staticSize>
558 bool operator<(const RhsValueType& a, const Evaluation<ValueType, numVars, staticSize>& b)
559 { return b > a; }
560 
561 template <class RhsValueType, class ValueType, int numVars, unsigned staticSize>
562 bool operator>(const RhsValueType& a, const Evaluation<ValueType, numVars, staticSize>& b)
563 { return b < a; }
564 
565 template <class RhsValueType, class ValueType, int numVars, unsigned staticSize>
566 bool operator<=(const RhsValueType& a, const Evaluation<ValueType, numVars, staticSize>& b)
567 { return b >= a; }
568 
569 template <class RhsValueType, class ValueType, int numVars, unsigned staticSize>
570 bool operator>=(const RhsValueType& a, const Evaluation<ValueType, numVars, staticSize>& b)
571 { return b <= a; }
572 
573 template <class RhsValueType, class ValueType, int numVars, unsigned staticSize>
574 bool operator!=(const RhsValueType& a, const Evaluation<ValueType, numVars, staticSize>& b)
575 { return a != b.value(); }
576 
577 template <class RhsValueType, class ValueType, int numVars, unsigned staticSize>
578 Evaluation<ValueType, numVars, staticSize> operator+(const RhsValueType& a, const Evaluation<ValueType, numVars, staticSize>& b)
579 {
580  Evaluation<ValueType, numVars, staticSize> result(b);
581  result += a;
582  return result;
583 }
584 
585 template <class RhsValueType, class ValueType, int numVars, unsigned staticSize>
586 Evaluation<ValueType, numVars, staticSize> operator-(const RhsValueType& a, const Evaluation<ValueType, numVars, staticSize>& b)
587 {
588  return -(b - a);
589 }
590 
591 template <class RhsValueType, class ValueType, int numVars, unsigned staticSize>
592 Evaluation<ValueType, numVars, staticSize> operator/(const RhsValueType& a, const Evaluation<ValueType, numVars, staticSize>& b)
593 {
594  Evaluation<ValueType, numVars, staticSize> tmp(a);
595  tmp /= b;
596  return tmp;
597 }
598 
599 template <class RhsValueType, class ValueType, int numVars, unsigned staticSize>
600 Evaluation<ValueType, numVars, staticSize> operator*(const RhsValueType& a, const Evaluation<ValueType, numVars, staticSize>& b)
601 {
602  Evaluation<ValueType, numVars, staticSize> result(b);
603  result *= a;
604  return result;
605 }
606 
607 template <class ValueType, int numVars, unsigned staticSize>
608 std::ostream& operator<<(std::ostream& os, const Evaluation<ValueType, numVars, staticSize>& eval)
609 {
610  os << eval.value();
611  return os;
612 }
613 
614 } // namespace DenseAd
615 } // namespace Opm
616 
617 // this makes the Dune matrix/vector classes happy...
618 #include <dune/common/ftraits.hh>
619 
620 namespace Dune {
621 template <class ValueType, int numVars, unsigned staticSize>
622 struct FieldTraits<Opm::DenseAd::Evaluation<ValueType, numVars, staticSize> >
623 {
624 public:
626  // setting real_type to field_type here potentially leads to slightly worse
627  // performance, but at least it makes things compile.
628  typedef field_type real_type;
629 };
630 
631 } // namespace Dune
632 
634 
635 #endif // OPM_DENSEAD_EVALUATION_HPP
This file includes all specializations for the dense-AD Evaluation class.
Representation of an evaluation of a function and its derivatives w.r.t.
A number of commonly used algebraic functions for the localized OPM automatic differentiation (AD) fr...
Some templates to wrap the valgrind client request macros.
Represents a function evaluation and its derivatives w.r.t.
Definition: Evaluation.hpp:59
Evaluation()
default constructor
Definition: Evaluation.hpp:100
ValueT ValueType
field type
Definition: Evaluation.hpp:66
void checkDefined_() const
instruct valgrind to check that the value and all derivatives of the Evaluation object are well-defin...
Definition: Evaluation.hpp:90
static const int numVars
the template argument which specifies the number of derivatives (-1 == "DynamicSize" means runtime de...
Definition: Evaluation.hpp:63
constexpr int size() const
number of derivatives
Definition: Evaluation.hpp:69
constexpr int valuepos_() const
position index for value
Definition: Evaluation.hpp:79
Evaluation(const Evaluation &other)=default
copy other function evaluation
constexpr int dend_() const
end+1 index for derivatives
Definition: Evaluation.hpp:85
constexpr int length_() const
length of internal data vector
Definition: Evaluation.hpp:74
constexpr int dstart_() const
start index for derivatives
Definition: Evaluation.hpp:82