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