All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
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 
37 #include <opm/common/Valgrind.hpp>
38 
39 #include <dune/common/version.hh>
40 
41 #include <array>
42 #include <cmath>
43 #include <cassert>
44 #include <cstring>
45 #include <iostream>
46 #include <algorithm>
47 
48 namespace Opm {
49 namespace DenseAd {
50 
51 template <class ValueT>
52 class Evaluation<ValueT, 12>
53 {
54 public:
56  typedef ValueT ValueType;
57 
59  static constexpr int size = 12;
60 
61 protected:
63  static constexpr int length_ = size + 1;
64 
66  static constexpr int valuepos_ = 0;
68  static constexpr int dstart_ = 1;
70  static constexpr int dend_ = length_;
71 
72 public:
74  Evaluation() : data_()
75  {}
76 
78  Evaluation(const Evaluation& other) = default;
79 
80  // create an evaluation which represents a constant function
81  //
82  // i.e., f(x) = c. this implies an evaluation with the given value and all
83  // derivatives being zero.
84  template <class RhsValueType>
85  Evaluation(const RhsValueType& c)
86  {
87  setValue( c );
88  clearDerivatives();
89  Valgrind::CheckDefined( data_ );
90  }
91 
92  // create an evaluation which represents a constant function
93  //
94  // i.e., f(x) = c. this implies an evaluation with the given value and all
95  // derivatives being zero.
96  template <class RhsValueType>
97  Evaluation(const RhsValueType& c, int varPos)
98  {
99  // The variable position must be in represented by the given variable descriptor
100  assert(0 <= varPos && varPos < size);
101 
102  setValue( c );
103  clearDerivatives();
104 
105  data_[varPos + dstart_] = 1.0;
106  Valgrind::CheckDefined(data_);
107  }
108 
109  // set all derivatives to zero
110  void clearDerivatives()
111  {
112  data_[1] = 0.0;
113  data_[2] = 0.0;
114  data_[3] = 0.0;
115  data_[4] = 0.0;
116  data_[5] = 0.0;
117  data_[6] = 0.0;
118  data_[7] = 0.0;
119  data_[8] = 0.0;
120  data_[9] = 0.0;
121  data_[10] = 0.0;
122  data_[11] = 0.0;
123  data_[12] = 0.0;
124  }
125 
126  // create a function evaluation for a "naked" depending variable (i.e., f(x) = x)
127  template <class RhsValueType>
128  static Evaluation createVariable(const RhsValueType& value, int varPos)
129  {
130  // copy function value and set all derivatives to 0, except for the variable
131  // which is represented by the value (which is set to 1.0)
132  return Evaluation( value, varPos );
133  }
134 
135  // "evaluate" a constant function (i.e. a function that does not depend on the set of
136  // relevant variables, f(x) = c).
137  template <class RhsValueType>
138  static Evaluation createConstant(const RhsValueType& value)
139  {
140  return Evaluation( value );
141  }
142 
143  // print the value and the derivatives of the function evaluation
144  void print(std::ostream& os = std::cout) const
145  {
146  // print value
147  os << "v: " << value() << " / d:";
148 
149  // print derivatives
150  for (int varIdx = 0; varIdx < size; ++varIdx) {
151  os << " " << derivative(varIdx);
152  }
153  }
154 
155  // copy all derivatives from other
156  void copyDerivatives(const Evaluation& other)
157  {
158  data_[1] = other.data_[1];
159  data_[2] = other.data_[2];
160  data_[3] = other.data_[3];
161  data_[4] = other.data_[4];
162  data_[5] = other.data_[5];
163  data_[6] = other.data_[6];
164  data_[7] = other.data_[7];
165  data_[8] = other.data_[8];
166  data_[9] = other.data_[9];
167  data_[10] = other.data_[10];
168  data_[11] = other.data_[11];
169  data_[12] = other.data_[12];
170  }
171 
172 
173  // add value and derivatives from other to this values and derivatives
174  Evaluation& operator+=(const Evaluation& other)
175  {
176  data_[0] += other.data_[0];
177  data_[1] += other.data_[1];
178  data_[2] += other.data_[2];
179  data_[3] += other.data_[3];
180  data_[4] += other.data_[4];
181  data_[5] += other.data_[5];
182  data_[6] += other.data_[6];
183  data_[7] += other.data_[7];
184  data_[8] += other.data_[8];
185  data_[9] += other.data_[9];
186  data_[10] += other.data_[10];
187  data_[11] += other.data_[11];
188  data_[12] += other.data_[12];
189 
190  return *this;
191  }
192 
193  // add value from other to this values
194  template <class RhsValueType>
195  Evaluation& operator+=(const RhsValueType& other)
196  {
197  // value is added, derivatives stay the same
198  data_[valuepos_] += other;
199 
200  return *this;
201  }
202 
203  // subtract other's value and derivatives from this values
204  Evaluation& operator-=(const Evaluation& other)
205  {
206  data_[0] -= other.data_[0];
207  data_[1] -= other.data_[1];
208  data_[2] -= other.data_[2];
209  data_[3] -= other.data_[3];
210  data_[4] -= other.data_[4];
211  data_[5] -= other.data_[5];
212  data_[6] -= other.data_[6];
213  data_[7] -= other.data_[7];
214  data_[8] -= other.data_[8];
215  data_[9] -= other.data_[9];
216  data_[10] -= other.data_[10];
217  data_[11] -= other.data_[11];
218  data_[12] -= other.data_[12];
219 
220  return *this;
221  }
222 
223  // subtract other's value from this values
224  template <class RhsValueType>
225  Evaluation& operator-=(const RhsValueType& other)
226  {
227  // for constants, values are subtracted, derivatives stay the same
228  data_[ valuepos_ ] -= other;
229 
230  return *this;
231  }
232 
233  // multiply values and apply chain rule to derivatives: (u*v)' = (v'u + u'v)
234  Evaluation& operator*=(const Evaluation& other)
235  {
236  // while the values are multiplied, the derivatives follow the product rule,
237  // i.e., (u*v)' = (v'u + u'v).
238  const ValueType u = this->value();
239  const ValueType v = other.value();
240 
241  // value
242  data_[valuepos_] *= v ;
243 
244  // derivatives
245  data_[1] = data_[1] * v + other.data_[1] * u;
246  data_[2] = data_[2] * v + other.data_[2] * u;
247  data_[3] = data_[3] * v + other.data_[3] * u;
248  data_[4] = data_[4] * v + other.data_[4] * u;
249  data_[5] = data_[5] * v + other.data_[5] * u;
250  data_[6] = data_[6] * v + other.data_[6] * u;
251  data_[7] = data_[7] * v + other.data_[7] * u;
252  data_[8] = data_[8] * v + other.data_[8] * u;
253  data_[9] = data_[9] * v + other.data_[9] * u;
254  data_[10] = data_[10] * v + other.data_[10] * u;
255  data_[11] = data_[11] * v + other.data_[11] * u;
256  data_[12] = data_[12] * v + other.data_[12] * u;
257 
258  return *this;
259  }
260 
261  // m(c*u)' = c*u'
262  template <class RhsValueType>
263  Evaluation& operator*=(const RhsValueType& other)
264  {
265  data_[0] *= other;
266  data_[1] *= other;
267  data_[2] *= other;
268  data_[3] *= other;
269  data_[4] *= other;
270  data_[5] *= other;
271  data_[6] *= other;
272  data_[7] *= other;
273  data_[8] *= other;
274  data_[9] *= other;
275  data_[10] *= other;
276  data_[11] *= other;
277  data_[12] *= other;
278 
279  return *this;
280  }
281 
282  // m(u*v)' = (vu' - uv')/v^2
283  Evaluation& operator/=(const Evaluation& other)
284  {
285  // values are divided, derivatives follow the rule for division, i.e., (u/v)' = (v'u -
286  // u'v)/v^2.
287  ValueType& u = data_[ valuepos_ ];
288  const ValueType& v = other.value();
289  data_[1] = (v*data_[1] - u*other.data_[1])/(v*v);
290  data_[2] = (v*data_[2] - u*other.data_[2])/(v*v);
291  data_[3] = (v*data_[3] - u*other.data_[3])/(v*v);
292  data_[4] = (v*data_[4] - u*other.data_[4])/(v*v);
293  data_[5] = (v*data_[5] - u*other.data_[5])/(v*v);
294  data_[6] = (v*data_[6] - u*other.data_[6])/(v*v);
295  data_[7] = (v*data_[7] - u*other.data_[7])/(v*v);
296  data_[8] = (v*data_[8] - u*other.data_[8])/(v*v);
297  data_[9] = (v*data_[9] - u*other.data_[9])/(v*v);
298  data_[10] = (v*data_[10] - u*other.data_[10])/(v*v);
299  data_[11] = (v*data_[11] - u*other.data_[11])/(v*v);
300  data_[12] = (v*data_[12] - u*other.data_[12])/(v*v);
301  u /= v;
302 
303  return *this;
304  }
305 
306  // divide value and derivatives by value of other
307  template <class RhsValueType>
308  Evaluation& operator/=(const RhsValueType& other)
309  {
310  const ValueType tmp = 1.0/other;
311 
312  data_[0] *= tmp;
313  data_[1] *= tmp;
314  data_[2] *= tmp;
315  data_[3] *= tmp;
316  data_[4] *= tmp;
317  data_[5] *= tmp;
318  data_[6] *= tmp;
319  data_[7] *= tmp;
320  data_[8] *= tmp;
321  data_[9] *= tmp;
322  data_[10] *= tmp;
323  data_[11] *= tmp;
324  data_[12] *= tmp;
325 
326  return *this;
327  }
328 
329  // add two evaluation objects
330  Evaluation operator+(const Evaluation& other) const
331  {
332  Evaluation result(*this);
333 
334  result += other;
335 
336  return result;
337  }
338 
339  // add constant to this object
340  template <class RhsValueType>
341  Evaluation operator+(const RhsValueType& other) const
342  {
343  Evaluation result(*this);
344 
345  result += other;
346 
347  return result;
348  }
349 
350  // subtract two evaluation objects
351  Evaluation operator-(const Evaluation& other) const
352  {
353  Evaluation result(*this);
354 
355  result -= other;
356 
357  return result;
358  }
359 
360  // subtract constant from evaluation object
361  template <class RhsValueType>
362  Evaluation operator-(const RhsValueType& other) const
363  {
364  Evaluation result(*this);
365 
366  result -= other;
367 
368  return result;
369  }
370 
371  // negation (unary minus) operator
372  Evaluation operator-() const
373  {
374  Evaluation result;
375 
376  // set value and derivatives to negative
377  result.data_[0] = - data_[0];
378  result.data_[1] = - data_[1];
379  result.data_[2] = - data_[2];
380  result.data_[3] = - data_[3];
381  result.data_[4] = - data_[4];
382  result.data_[5] = - data_[5];
383  result.data_[6] = - data_[6];
384  result.data_[7] = - data_[7];
385  result.data_[8] = - data_[8];
386  result.data_[9] = - data_[9];
387  result.data_[10] = - data_[10];
388  result.data_[11] = - data_[11];
389  result.data_[12] = - data_[12];
390 
391  return result;
392  }
393 
394  Evaluation operator*(const Evaluation& other) const
395  {
396  Evaluation result(*this);
397 
398  result *= other;
399 
400  return result;
401  }
402 
403  template <class RhsValueType>
404  Evaluation operator*(const RhsValueType& other) const
405  {
406  Evaluation result(*this);
407 
408  result *= other;
409 
410  return result;
411  }
412 
413  Evaluation operator/(const Evaluation& other) const
414  {
415  Evaluation result(*this);
416 
417  result /= other;
418 
419  return result;
420  }
421 
422  template <class RhsValueType>
423  Evaluation operator/(const RhsValueType& other) const
424  {
425  Evaluation result(*this);
426 
427  result /= other;
428 
429  return result;
430  }
431 
432  template <class RhsValueType>
433  Evaluation& operator=(const RhsValueType& other)
434  {
435  setValue( other );
436  clearDerivatives();
437 
438  return *this;
439  }
440 
441  // copy assignment from evaluation
442  Evaluation& operator=(const Evaluation& other) = default;
443 
444  template <class RhsValueType>
445  bool operator==(const RhsValueType& other) const
446  { return value() == other; }
447 
448  bool operator==(const Evaluation& other) const
449  {
450  for (int idx = 0; idx < length_; ++idx) {
451  if (data_[idx] != other.data_[idx]) {
452  return false;
453  }
454  }
455  return true;
456  }
457 
458  bool operator!=(const Evaluation& other) const
459  { return !operator==(other); }
460 
461  template <class RhsValueType>
462  bool operator>(RhsValueType other) const
463  { return value() > other; }
464 
465  bool operator>(const Evaluation& other) const
466  { return value() > other.value(); }
467 
468  template <class RhsValueType>
469  bool operator<(RhsValueType other) const
470  { return value() < other; }
471 
472  bool operator<(const Evaluation& other) const
473  { return value() < other.value(); }
474 
475  template <class RhsValueType>
476  bool operator>=(RhsValueType other) const
477  { return value() >= other; }
478 
479  bool operator>=(const Evaluation& other) const
480  { return value() >= other.value(); }
481 
482  template <class RhsValueType>
483  bool operator<=(RhsValueType other) const
484  { return value() <= other; }
485 
486  bool operator<=(const Evaluation& other) const
487  { return value() <= other.value(); }
488 
489  // return value of variable
490  const ValueType& value() const
491  { return data_[valuepos_]; }
492 
493  // set value of variable
494  template <class RhsValueType>
495  void setValue(const RhsValueType& val)
496  { data_[valuepos_] = val; }
497 
498  // return varIdx'th derivative
499  const ValueType& derivative(int varIdx) const
500  {
501  assert(0 <= varIdx && varIdx < size);
502 
503  return data_[dstart_ + varIdx];
504  }
505 
506  // set derivative at position varIdx
507  void setDerivative(int varIdx, const ValueType& derVal)
508  {
509  assert(0 <= varIdx && varIdx < size);
510 
511  data_[dstart_ + varIdx] = derVal;
512  }
513 
514 private:
515  std::array<ValueT, length_> data_;
516 };
517 
518 } } // namespace DenseAd, Opm
519 
520 #endif // OPM_DENSEAD_EVALUATION12_HPP
Evaluation()
default constructor
Definition: Evaluation.hpp:79
ValueT ValueType
field type
Definition: Evaluation12.hpp:56
A number of commonly used algebraic functions for the localized OPM automatic differentiation (AD) fr...
ValueT ValueType
field type
Definition: Evaluation.hpp:61
static constexpr int valuepos_
position index for value
Definition: Evaluation.hpp:71
static constexpr int length_
length of internal data vector
Definition: Evaluation.hpp:68
static constexpr int dstart_
start index for derivatives
Definition: Evaluation.hpp:73
Evaluation()
default constructor
Definition: Evaluation12.hpp:74
Represents a function evaluation and its derivatives w.r.t.
Definition: Evaluation.hpp:57
static constexpr int size
number of derivatives
Definition: Evaluation.hpp:64
static constexpr int dend_
end+1 index for derivatives
Definition: Evaluation.hpp:75
Representation of an evaluation of a function and its derivatives w.r.t.