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