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