All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Evaluation8.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_EVALUATION8_HPP
32 #define OPM_DENSEAD_EVALUATION8_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, 8>
53 {
54 public:
56  typedef ValueT ValueType;
57 
59  static constexpr int size = 8;
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  }
121 
122  // create a function evaluation for a "naked" depending variable (i.e., f(x) = x)
123  template <class RhsValueType>
124  static Evaluation createVariable(const RhsValueType& value, int varPos)
125  {
126  // copy function value and set all derivatives to 0, except for the variable
127  // which is represented by the value (which is set to 1.0)
128  return Evaluation( value, varPos );
129  }
130 
131  // "evaluate" a constant function (i.e. a function that does not depend on the set of
132  // relevant variables, f(x) = c).
133  template <class RhsValueType>
134  static Evaluation createConstant(const RhsValueType& value)
135  {
136  return Evaluation( value );
137  }
138 
139  // print the value and the derivatives of the function evaluation
140  void print(std::ostream& os = std::cout) const
141  {
142  // print value
143  os << "v: " << value() << " / d:";
144 
145  // print derivatives
146  for (int varIdx = 0; varIdx < size; ++varIdx) {
147  os << " " << derivative(varIdx);
148  }
149  }
150 
151  // copy all derivatives from other
152  void copyDerivatives(const Evaluation& other)
153  {
154  data_[1] = other.data_[1];
155  data_[2] = other.data_[2];
156  data_[3] = other.data_[3];
157  data_[4] = other.data_[4];
158  data_[5] = other.data_[5];
159  data_[6] = other.data_[6];
160  data_[7] = other.data_[7];
161  data_[8] = other.data_[8];
162  }
163 
164 
165  // add value and derivatives from other to this values and derivatives
166  Evaluation& operator+=(const Evaluation& other)
167  {
168  data_[0] += other.data_[0];
169  data_[1] += other.data_[1];
170  data_[2] += other.data_[2];
171  data_[3] += other.data_[3];
172  data_[4] += other.data_[4];
173  data_[5] += other.data_[5];
174  data_[6] += other.data_[6];
175  data_[7] += other.data_[7];
176  data_[8] += other.data_[8];
177 
178  return *this;
179  }
180 
181  // add value from other to this values
182  template <class RhsValueType>
183  Evaluation& operator+=(const RhsValueType& other)
184  {
185  // value is added, derivatives stay the same
186  data_[valuepos_] += other;
187 
188  return *this;
189  }
190 
191  // subtract other's value and derivatives from this values
192  Evaluation& operator-=(const Evaluation& other)
193  {
194  data_[0] -= other.data_[0];
195  data_[1] -= other.data_[1];
196  data_[2] -= other.data_[2];
197  data_[3] -= other.data_[3];
198  data_[4] -= other.data_[4];
199  data_[5] -= other.data_[5];
200  data_[6] -= other.data_[6];
201  data_[7] -= other.data_[7];
202  data_[8] -= other.data_[8];
203 
204  return *this;
205  }
206 
207  // subtract other's value from this values
208  template <class RhsValueType>
209  Evaluation& operator-=(const RhsValueType& other)
210  {
211  // for constants, values are subtracted, derivatives stay the same
212  data_[ valuepos_ ] -= other;
213 
214  return *this;
215  }
216 
217  // multiply values and apply chain rule to derivatives: (u*v)' = (v'u + u'v)
218  Evaluation& operator*=(const Evaluation& other)
219  {
220  // while the values are multiplied, the derivatives follow the product rule,
221  // i.e., (u*v)' = (v'u + u'v).
222  const ValueType u = this->value();
223  const ValueType v = other.value();
224 
225  // value
226  data_[valuepos_] *= v ;
227 
228  // derivatives
229  data_[1] = data_[1] * v + other.data_[1] * u;
230  data_[2] = data_[2] * v + other.data_[2] * u;
231  data_[3] = data_[3] * v + other.data_[3] * u;
232  data_[4] = data_[4] * v + other.data_[4] * u;
233  data_[5] = data_[5] * v + other.data_[5] * u;
234  data_[6] = data_[6] * v + other.data_[6] * u;
235  data_[7] = data_[7] * v + other.data_[7] * u;
236  data_[8] = data_[8] * v + other.data_[8] * u;
237 
238  return *this;
239  }
240 
241  // m(c*u)' = c*u'
242  template <class RhsValueType>
243  Evaluation& operator*=(const RhsValueType& other)
244  {
245  data_[0] *= other;
246  data_[1] *= other;
247  data_[2] *= other;
248  data_[3] *= other;
249  data_[4] *= other;
250  data_[5] *= other;
251  data_[6] *= other;
252  data_[7] *= other;
253  data_[8] *= other;
254 
255  return *this;
256  }
257 
258  // m(u*v)' = (vu' - uv')/v^2
259  Evaluation& operator/=(const Evaluation& other)
260  {
261  // values are divided, derivatives follow the rule for division, i.e., (u/v)' = (v'u -
262  // u'v)/v^2.
263  ValueType& u = data_[ valuepos_ ];
264  const ValueType& v = other.value();
265  data_[1] = (v*data_[1] - u*other.data_[1])/(v*v);
266  data_[2] = (v*data_[2] - u*other.data_[2])/(v*v);
267  data_[3] = (v*data_[3] - u*other.data_[3])/(v*v);
268  data_[4] = (v*data_[4] - u*other.data_[4])/(v*v);
269  data_[5] = (v*data_[5] - u*other.data_[5])/(v*v);
270  data_[6] = (v*data_[6] - u*other.data_[6])/(v*v);
271  data_[7] = (v*data_[7] - u*other.data_[7])/(v*v);
272  data_[8] = (v*data_[8] - u*other.data_[8])/(v*v);
273  u /= v;
274 
275  return *this;
276  }
277 
278  // divide value and derivatives by value of other
279  template <class RhsValueType>
280  Evaluation& operator/=(const RhsValueType& other)
281  {
282  const ValueType tmp = 1.0/other;
283 
284  data_[0] *= tmp;
285  data_[1] *= tmp;
286  data_[2] *= tmp;
287  data_[3] *= tmp;
288  data_[4] *= tmp;
289  data_[5] *= tmp;
290  data_[6] *= tmp;
291  data_[7] *= tmp;
292  data_[8] *= tmp;
293 
294  return *this;
295  }
296 
297  // add two evaluation objects
298  Evaluation operator+(const Evaluation& other) const
299  {
300  Evaluation result(*this);
301 
302  result += other;
303 
304  return result;
305  }
306 
307  // add constant to this object
308  template <class RhsValueType>
309  Evaluation operator+(const RhsValueType& other) const
310  {
311  Evaluation result(*this);
312 
313  result += other;
314 
315  return result;
316  }
317 
318  // subtract two evaluation objects
319  Evaluation operator-(const Evaluation& other) const
320  {
321  Evaluation result(*this);
322 
323  result -= other;
324 
325  return result;
326  }
327 
328  // subtract constant from evaluation object
329  template <class RhsValueType>
330  Evaluation operator-(const RhsValueType& other) const
331  {
332  Evaluation result(*this);
333 
334  result -= other;
335 
336  return result;
337  }
338 
339  // negation (unary minus) operator
340  Evaluation operator-() const
341  {
342  Evaluation result;
343 
344  // set value and derivatives to negative
345  result.data_[0] = - data_[0];
346  result.data_[1] = - data_[1];
347  result.data_[2] = - data_[2];
348  result.data_[3] = - data_[3];
349  result.data_[4] = - data_[4];
350  result.data_[5] = - data_[5];
351  result.data_[6] = - data_[6];
352  result.data_[7] = - data_[7];
353  result.data_[8] = - data_[8];
354 
355  return result;
356  }
357 
358  Evaluation operator*(const Evaluation& other) const
359  {
360  Evaluation result(*this);
361 
362  result *= other;
363 
364  return result;
365  }
366 
367  template <class RhsValueType>
368  Evaluation operator*(const RhsValueType& other) const
369  {
370  Evaluation result(*this);
371 
372  result *= other;
373 
374  return result;
375  }
376 
377  Evaluation operator/(const Evaluation& other) const
378  {
379  Evaluation result(*this);
380 
381  result /= other;
382 
383  return result;
384  }
385 
386  template <class RhsValueType>
387  Evaluation operator/(const RhsValueType& other) const
388  {
389  Evaluation result(*this);
390 
391  result /= other;
392 
393  return result;
394  }
395 
396  template <class RhsValueType>
397  Evaluation& operator=(const RhsValueType& other)
398  {
399  setValue( other );
400  clearDerivatives();
401 
402  return *this;
403  }
404 
405  // copy assignment from evaluation
406  Evaluation& operator=(const Evaluation& other) = default;
407 
408  template <class RhsValueType>
409  bool operator==(const RhsValueType& other) const
410  { return value() == other; }
411 
412  bool operator==(const Evaluation& other) const
413  {
414  for (int idx = 0; idx < length_; ++idx) {
415  if (data_[idx] != other.data_[idx]) {
416  return false;
417  }
418  }
419  return true;
420  }
421 
422  bool operator!=(const Evaluation& other) const
423  { return !operator==(other); }
424 
425  template <class RhsValueType>
426  bool operator>(RhsValueType other) const
427  { return value() > other; }
428 
429  bool operator>(const Evaluation& other) const
430  { return value() > other.value(); }
431 
432  template <class RhsValueType>
433  bool operator<(RhsValueType other) const
434  { return value() < other; }
435 
436  bool operator<(const Evaluation& other) const
437  { return value() < other.value(); }
438 
439  template <class RhsValueType>
440  bool operator>=(RhsValueType other) const
441  { return value() >= other; }
442 
443  bool operator>=(const Evaluation& other) const
444  { return value() >= other.value(); }
445 
446  template <class RhsValueType>
447  bool operator<=(RhsValueType other) const
448  { return value() <= other; }
449 
450  bool operator<=(const Evaluation& other) const
451  { return value() <= other.value(); }
452 
453  // return value of variable
454  const ValueType& value() const
455  { return data_[valuepos_]; }
456 
457  // set value of variable
458  template <class RhsValueType>
459  void setValue(const RhsValueType& val)
460  { data_[valuepos_] = val; }
461 
462  // return varIdx'th derivative
463  const ValueType& derivative(int varIdx) const
464  {
465  assert(0 <= varIdx && varIdx < size);
466 
467  return data_[dstart_ + varIdx];
468  }
469 
470  // set derivative at position varIdx
471  void setDerivative(int varIdx, const ValueType& derVal)
472  {
473  assert(0 <= varIdx && varIdx < size);
474 
475  data_[dstart_ + varIdx] = derVal;
476  }
477 
478 private:
479  std::array<ValueT, length_> data_;
480 };
481 
482 } } // namespace DenseAd, Opm
483 
484 #endif // OPM_DENSEAD_EVALUATION8_HPP
Evaluation()
default constructor
Definition: Evaluation.hpp:79
ValueT ValueType
field type
Definition: Evaluation8.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
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
Evaluation()
default constructor
Definition: Evaluation8.hpp:74
Representation of an evaluation of a function and its derivatives w.r.t.