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