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