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