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