My Project
Evaluation2.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 */
31 #ifndef OPM_DENSEAD_EVALUATION2_HPP
32 #define OPM_DENSEAD_EVALUATION2_HPP
33 
34 #include "Evaluation.hpp"
35 #include "Math.hpp"
36 
38 
39 #include <array>
40 #include <cmath>
41 #include <cassert>
42 #include <cstring>
43 #include <iostream>
44 #include <algorithm>
45 
46 namespace Opm {
47 namespace DenseAd {
48 
49 template <class ValueT>
50 class Evaluation<ValueT, 2>
51 {
52 public:
55  static const int numVars = 2;
56 
58  typedef ValueT ValueType;
59 
61  constexpr int size() const
62  { return 2; };
63 
64 protected:
66  constexpr int length_() const
67  { return size() + 1; }
68 
69 
71  constexpr int valuepos_() const
72  { return 0; }
74  constexpr int dstart_() const
75  { return 1; }
77  constexpr int dend_() const
78  { return length_(); }
79 
82  void checkDefined_() const
83  {
84 #ifndef NDEBUG
85  for (const auto& v: data_)
86  Valgrind::CheckDefined(v);
87 #endif
88  }
89 
90 public:
92  Evaluation() : data_()
93  {}
94 
96  Evaluation(const Evaluation& other) = default;
97 
98 
99  // create an evaluation which represents a constant function
100  //
101  // i.e., f(x) = c. this implies an evaluation with the given value and all
102  // derivatives being zero.
103  template <class RhsValueType>
104  Evaluation(const RhsValueType& c)
105  {
106  setValue(c);
107  clearDerivatives();
108 
109  checkDefined_();
110  }
111 
112  // create an evaluation which represents a constant function
113  //
114  // i.e., f(x) = c. this implies an evaluation with the given value and all
115  // derivatives being zero.
116  template <class RhsValueType>
117  Evaluation(const RhsValueType& c, int varPos)
118  {
119  // The variable position must be in represented by the given variable descriptor
120  assert(0 <= varPos && varPos < size());
121 
122  setValue( c );
123  clearDerivatives();
124 
125  data_[varPos + dstart_()] = 1.0;
126 
127  checkDefined_();
128  }
129 
130  // set all derivatives to zero
131  void clearDerivatives()
132  {
133  data_[1] = 0.0;
134  data_[2] = 0.0;
135  }
136 
137  // create an uninitialized Evaluation object that is compatible with the
138  // argument, but not initialized
139  //
140  // This basically boils down to the copy constructor without copying
141  // anything. If the number of derivatives is known at compile time, this
142  // is equivalent to creating an uninitialized object using the default
143  // constructor, while for dynamic evaluations, it creates an Evaluation
144  // object which exhibits the same number of derivatives as the argument.
145  static Evaluation createBlank(const Evaluation&)
146  { return Evaluation(); }
147 
148  // create an Evaluation with value and all the derivatives to be zero
149  static Evaluation createConstantZero(const Evaluation&)
150  { return Evaluation(0.); }
151 
152  // create an Evaluation with value to be one and all the derivatives to be zero
153  static Evaluation createConstantOne(const Evaluation&)
154  { return Evaluation(1.); }
155 
156  // create a function evaluation for a "naked" depending variable (i.e., f(x) = x)
157  template <class RhsValueType>
158  static Evaluation createVariable(const RhsValueType& value, int varPos)
159  {
160  // copy function value and set all derivatives to 0, except for the variable
161  // which is represented by the value (which is set to 1.0)
162  return Evaluation(value, varPos);
163  }
164 
165  template <class RhsValueType>
166  static Evaluation createVariable(int nVars, const RhsValueType& value, int varPos)
167  {
168  if (nVars != 2)
169  throw std::logic_error("This statically-sized evaluation can only represent objects"
170  " with 2 derivatives");
171 
172  // copy function value and set all derivatives to 0, except for the variable
173  // which is represented by the value (which is set to 1.0)
174  return Evaluation(nVars, value, varPos);
175  }
176 
177  template <class RhsValueType>
178  static Evaluation createVariable(const Evaluation&, const RhsValueType& value, int varPos)
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(value, varPos);
183  }
184 
185 
186  // "evaluate" a constant function (i.e. a function that does not depend on the set of
187  // relevant variables, f(x) = c).
188  template <class RhsValueType>
189  static Evaluation createConstant(int nVars, const RhsValueType& value)
190  {
191  if (nVars != 2)
192  throw std::logic_error("This statically-sized evaluation can only represent objects"
193  " with 2 derivatives");
194  return Evaluation(value);
195  }
196 
197  // "evaluate" a constant function (i.e. a function that does not depend on the set of
198  // relevant variables, f(x) = c).
199  template <class RhsValueType>
200  static Evaluation createConstant(const RhsValueType& value)
201  {
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 Evaluation&, const RhsValueType& value)
209  {
210  return Evaluation(value);
211  }
212 
213  // print the value and the derivatives of the function evaluation
214  void print(std::ostream& os = std::cout) const
215  {
216  // print value
217  os << "v: " << value() << " / d:";
218 
219  // print derivatives
220  for (int varIdx = 0; varIdx < size(); ++varIdx) {
221  os << " " << derivative(varIdx);
222  }
223  }
224 
225  // copy all derivatives from other
226  void copyDerivatives(const Evaluation& other)
227  {
228  assert(size() == other.size());
229 
230  data_[1] = other.data_[1];
231  data_[2] = other.data_[2];
232  }
233 
234 
235  // add value and derivatives from other to this values and derivatives
236  Evaluation& operator+=(const Evaluation& other)
237  {
238  assert(size() == other.size());
239 
240  data_[0] += other.data_[0];
241  data_[1] += other.data_[1];
242  data_[2] += other.data_[2];
243 
244  return *this;
245  }
246 
247  // add value from other to this values
248  template <class RhsValueType>
249  Evaluation& operator+=(const RhsValueType& other)
250  {
251  // value is added, derivatives stay the same
252  data_[valuepos_()] += other;
253 
254  return *this;
255  }
256 
257  // subtract other's value and derivatives from this values
258  Evaluation& operator-=(const Evaluation& other)
259  {
260  assert(size() == other.size());
261 
262  data_[0] -= other.data_[0];
263  data_[1] -= other.data_[1];
264  data_[2] -= other.data_[2];
265 
266  return *this;
267  }
268 
269  // subtract other's value from this values
270  template <class RhsValueType>
271  Evaluation& operator-=(const RhsValueType& other)
272  {
273  // for constants, values are subtracted, derivatives stay the same
274  data_[valuepos_()] -= other;
275 
276  return *this;
277  }
278 
279  // multiply values and apply chain rule to derivatives: (u*v)' = (v'u + u'v)
280  Evaluation& operator*=(const Evaluation& other)
281  {
282  assert(size() == other.size());
283 
284  // while the values are multiplied, the derivatives follow the product rule,
285  // i.e., (u*v)' = (v'u + u'v).
286  const ValueType u = this->value();
287  const ValueType v = other.value();
288 
289  // value
290  data_[valuepos_()] *= v ;
291 
292  // derivatives
293  data_[1] = data_[1] * v + other.data_[1] * u;
294  data_[2] = data_[2] * v + other.data_[2] * u;
295 
296  return *this;
297  }
298 
299  // m(c*u)' = c*u'
300  template <class RhsValueType>
301  Evaluation& operator*=(const RhsValueType& other)
302  {
303  data_[0] *= other;
304  data_[1] *= other;
305  data_[2] *= other;
306 
307  return *this;
308  }
309 
310  // m(u*v)' = (vu' - uv')/v^2
311  Evaluation& operator/=(const Evaluation& other)
312  {
313  assert(size() == other.size());
314 
315  // values are divided, derivatives follow the rule for division, i.e., (u/v)' = (v'u -
316  // u'v)/v^2.
317  ValueType& u = data_[valuepos_()];
318  const ValueType& v = other.value();
319  data_[1] = (v*data_[1] - u*other.data_[1])/(v*v);
320  data_[2] = (v*data_[2] - u*other.data_[2])/(v*v);
321  u /= v;
322 
323  return *this;
324  }
325 
326  // divide value and derivatives by value of other
327  template <class RhsValueType>
328  Evaluation& operator/=(const RhsValueType& other)
329  {
330  const ValueType tmp = 1.0/other;
331 
332  data_[0] *= tmp;
333  data_[1] *= tmp;
334  data_[2] *= tmp;
335 
336  return *this;
337  }
338 
339  // add two evaluation objects
340  Evaluation operator+(const Evaluation& other) const
341  {
342  assert(size() == other.size());
343 
344  Evaluation result(*this);
345 
346  result += other;
347 
348  return result;
349  }
350 
351  // add constant to this object
352  template <class RhsValueType>
353  Evaluation operator+(const RhsValueType& other) const
354  {
355  Evaluation result(*this);
356 
357  result += other;
358 
359  return result;
360  }
361 
362  // subtract two evaluation objects
363  Evaluation operator-(const Evaluation& other) const
364  {
365  assert(size() == other.size());
366 
367  Evaluation result(*this);
368 
369  result -= other;
370 
371  return result;
372  }
373 
374  // subtract constant from evaluation object
375  template <class RhsValueType>
376  Evaluation operator-(const RhsValueType& other) const
377  {
378  Evaluation result(*this);
379 
380  result -= other;
381 
382  return result;
383  }
384 
385  // negation (unary minus) operator
386  Evaluation operator-() const
387  {
388  Evaluation result;
389 
390  // set value and derivatives to negative
391  result.data_[0] = - data_[0];
392  result.data_[1] = - data_[1];
393  result.data_[2] = - data_[2];
394 
395  return result;
396  }
397 
398  Evaluation operator*(const Evaluation& other) const
399  {
400  assert(size() == other.size());
401 
402  Evaluation result(*this);
403 
404  result *= other;
405 
406  return result;
407  }
408 
409  template <class RhsValueType>
410  Evaluation operator*(const RhsValueType& other) const
411  {
412  Evaluation result(*this);
413 
414  result *= other;
415 
416  return result;
417  }
418 
419  Evaluation operator/(const Evaluation& other) const
420  {
421  assert(size() == other.size());
422 
423  Evaluation result(*this);
424 
425  result /= other;
426 
427  return result;
428  }
429 
430  template <class RhsValueType>
431  Evaluation operator/(const RhsValueType& other) const
432  {
433  Evaluation result(*this);
434 
435  result /= other;
436 
437  return result;
438  }
439 
440  template <class RhsValueType>
441  Evaluation& operator=(const RhsValueType& other)
442  {
443  setValue( other );
444  clearDerivatives();
445 
446  return *this;
447  }
448 
449  // copy assignment from evaluation
450  Evaluation& operator=(const Evaluation& other) = default;
451 
452  template <class RhsValueType>
453  bool operator==(const RhsValueType& other) const
454  { return value() == other; }
455 
456  bool operator==(const Evaluation& other) const
457  {
458  assert(size() == other.size());
459 
460  for (int idx = 0; idx < length_(); ++idx) {
461  if (data_[idx] != other.data_[idx]) {
462  return false;
463  }
464  }
465  return true;
466  }
467 
468  bool operator!=(const Evaluation& other) const
469  { return !operator==(other); }
470 
471  template <class RhsValueType>
472  bool operator!=(const RhsValueType& other) const
473  { return !operator==(other); }
474 
475  template <class RhsValueType>
476  bool operator>(RhsValueType other) const
477  { return value() > other; }
478 
479  bool operator>(const Evaluation& other) const
480  {
481  assert(size() == other.size());
482 
483  return value() > other.value();
484  }
485 
486  template <class RhsValueType>
487  bool operator<(RhsValueType other) const
488  { return value() < other; }
489 
490  bool operator<(const Evaluation& other) const
491  {
492  assert(size() == other.size());
493 
494  return value() < other.value();
495  }
496 
497  template <class RhsValueType>
498  bool operator>=(RhsValueType other) const
499  { return value() >= other; }
500 
501  bool operator>=(const Evaluation& other) const
502  {
503  assert(size() == other.size());
504 
505  return value() >= other.value();
506  }
507 
508  template <class RhsValueType>
509  bool operator<=(RhsValueType other) const
510  { return value() <= other; }
511 
512  bool operator<=(const Evaluation& other) const
513  {
514  assert(size() == other.size());
515 
516  return value() <= other.value();
517  }
518 
519  // return value of variable
520  const ValueType& value() const
521  { return data_[valuepos_()]; }
522 
523  // set value of variable
524  template <class RhsValueType>
525  void setValue(const RhsValueType& val)
526  { data_[valuepos_()] = val; }
527 
528  // return varIdx'th derivative
529  const ValueType& derivative(int varIdx) const
530  {
531  assert(0 <= varIdx && varIdx < size());
532 
533  return data_[dstart_() + varIdx];
534  }
535 
536  // set derivative at position varIdx
537  void setDerivative(int varIdx, const ValueType& derVal)
538  {
539  assert(0 <= varIdx && varIdx < size());
540 
541  data_[dstart_() + varIdx] = derVal;
542  }
543 
544 private:
545 
546  std::array<ValueT, 3> data_;
547 };
548 
549 } // namespace DenseAd
550 } // namespace Opm
551 
552 #endif // OPM_DENSEAD_EVALUATION2_HPP
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.
constexpr int valuepos_() const
position index for value
Definition: Evaluation2.hpp:71
ValueT ValueType
field type
Definition: Evaluation2.hpp:58
constexpr int dstart_() const
start index for derivatives
Definition: Evaluation2.hpp:74
constexpr int size() const
number of derivatives
Definition: Evaluation2.hpp:61
Evaluation()
default constructor
Definition: Evaluation2.hpp:92
void checkDefined_() const
instruct valgrind to check that the value and all derivatives of the Evaluation object are well-defin...
Definition: Evaluation2.hpp:82
constexpr int dend_() const
end+1 index for derivatives
Definition: Evaluation2.hpp:77
constexpr int length_() const
length of internal data vector
Definition: Evaluation2.hpp:66
Evaluation(const Evaluation &other)=default
copy other function evaluation
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
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