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