My Project
DynamicEvaluation.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 */
32 #ifndef OPM_DENSEAD_EVALUATION_DYNAMIC_HPP
33 #define OPM_DENSEAD_EVALUATION_DYNAMIC_HPP
34 
35 #include "Evaluation.hpp"
36 #include "Math.hpp"
37 
39 
41 #include <cmath>
42 #include <cassert>
43 #include <cstring>
44 #include <iostream>
45 #include <algorithm>
46 
47 namespace Opm {
48 namespace DenseAd {
49 
54 template <class ValueT, unsigned staticSize>
55 class Evaluation<ValueT, DynamicSize, staticSize>
56 {
57 public:
60  static const int numVars = DynamicSize;
61 
63  typedef ValueT ValueType;
64 
66  int size() const
67  { return data_.size() - 1; }
68 
69 protected:
71  int length_() const
72  { return data_.size(); }
73 
74 
76  constexpr int valuepos_() const
77  { return 0; }
79  constexpr int dstart_() const
80  { return 1; }
82  int dend_() const
83  { return length_(); }
84 
87  void checkDefined_() const
88  {
89 #ifndef NDEBUG
90  for (int i = dstart_(); i < dend_(); ++i)
91  Valgrind::CheckDefined(data_[i]);
92 #endif
93  }
94 
95 public:
97  Evaluation() : data_()
98  {}
99 
101  Evaluation(const Evaluation& other) = default;
102 
106  : data_(std::move(other.data_))
107  { }
108 
111  {
112  data_ = std::move(other.data_);
113  return *this;
114  }
115 
116  // create a "blank" dynamic evaluation
117  explicit Evaluation(int numDerivatives)
118  : data_(1 + numDerivatives)
119  {}
120 
121  // create a dynamic evaluation which represents a constant function
122  //
123  // i.e., f(x) = c. this implies an evaluation with the given value and all
124  // derivatives being zero.
125  template <class RhsValueType>
126  Evaluation(int numDerivatives, const RhsValueType& c)
127  : data_(1 + numDerivatives, 0.0)
128  {
129  //clearDerivatives();
130  setValue(c);
131 
132  checkDefined_();
133  }
134 
135  // create an evaluation which represents a constant function
136  //
137  // i.e., f(x) = c. this implies an evaluation with the given value and all
138  // derivatives being zero.
139  template <class RhsValueType>
140  Evaluation(int nVars, const RhsValueType& c, int varPos)
141  : data_(1 + nVars, 0.0)
142  {
143  // The variable position must be in represented by the given variable descriptor
144  assert(0 <= varPos && varPos < size());
145 
146  setValue(c);
147 
148  data_[varPos + dstart_()] = 1.0;
149 
150  checkDefined_();
151  }
152 
153  // set all derivatives to zero
154  void clearDerivatives()
155  {
156  for (int i = dstart_(); i < dend_(); ++i)
157  data_[i] = 0.0;
158  }
159 
160  // create an uninitialized Evaluation object that is compatible with the
161  // argument, but not initialized
162  //
163  // This basically boils down to the copy constructor without copying
164  // anything. If the number of derivatives is known at compile time, this
165  // is equivalent to creating an uninitialized object using the default
166  // constructor, while for dynamic evaluations, it creates an Evaluation
167  // object which exhibits the same number of derivatives as the argument.
168  static Evaluation createBlank(const Evaluation& x)
169  { return Evaluation(x.size()); }
170 
171  // create an Evaluation with value and all the derivatives to be zero
172  static Evaluation createConstantZero(const Evaluation& x)
173  { return Evaluation(x.size(), 0.0); }
174 
175  // create an Evaluation with value to be one and all the derivatives to be zero
176  static Evaluation createConstantOne(const Evaluation& x)
177  { return Evaluation(x.size(), 1.); }
178 
179  // create a function evaluation for a "naked" depending variable (i.e., f(x) = x)
180  template <class RhsValueType>
181  static Evaluation createVariable(const RhsValueType&, int)
182  {
183  throw std::logic_error("Dynamically sized evaluations require that the number of "
184  "derivatives is specified when creating an evaluation");
185  }
186 
187  template <class RhsValueType>
188  static Evaluation createVariable(int nVars, 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(nVars, value, varPos);
193  }
194 
195  template <class RhsValueType>
196  static Evaluation createVariable(const Evaluation& x, const RhsValueType& value, int varPos)
197  {
198  // copy function value and set all derivatives to 0, except for the variable
199  // which is represented by the value (which is set to 1.0)
200  return Evaluation(x.size(), value, varPos);
201  }
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(int nVars, const RhsValueType& value)
208  {
209  return Evaluation(nVars, 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 RhsValueType&)
216  {
217  throw std::logic_error("Dynamically-sized evaluation objects require to specify the number of derivatives.");
218  }
219 
220  // "evaluate" a constant function (i.e. a function that does not depend on the set of
221  // relevant variables, f(x) = c).
222  template <class RhsValueType>
223  static Evaluation createConstant(const Evaluation& x, const RhsValueType& value)
224  {
225  return Evaluation(x.size(), value);
226  }
227 
228  // print the value and the derivatives of the function evaluation
229  void print(std::ostream& os = std::cout) const
230  {
231  // print value
232  os << "v: " << value() << " / d:";
233 
234  // print derivatives
235  for (int varIdx = 0; varIdx < size(); ++varIdx) {
236  os << " " << derivative(varIdx);
237  }
238  }
239 
240  // copy all derivatives from other
241  void copyDerivatives(const Evaluation& other)
242  {
243  assert(size() == other.size());
244 
245  for (int i = dstart_(); i < dend_(); ++i)
246  data_[i] = other.data_[i];
247  }
248 
249 
250  // add value and derivatives from other to this values and derivatives
251  Evaluation& operator+=(const Evaluation& other)
252  {
253  assert(size() == other.size());
254 
255  for (int i = 0; i < length_(); ++i)
256  data_[i] += other.data_[i];
257 
258  return *this;
259  }
260 
261  // add value from other to this values
262  template <class RhsValueType>
263  Evaluation& operator+=(const RhsValueType& other)
264  {
265  // value is added, derivatives stay the same
266  data_[valuepos_()] += other;
267 
268  return *this;
269  }
270 
271  // subtract other's value and derivatives from this values
272  Evaluation& operator-=(const Evaluation& other)
273  {
274  assert(size() == other.size());
275 
276  for (int i = 0; i < length_(); ++i)
277  data_[i] -= other.data_[i];
278 
279  return *this;
280  }
281 
282  // subtract other's value from this values
283  template <class RhsValueType>
284  Evaluation& operator-=(const RhsValueType& other)
285  {
286  // for constants, values are subtracted, derivatives stay the same
287  data_[valuepos_()] -= other;
288 
289  return *this;
290  }
291 
292  // multiply values and apply chain rule to derivatives: (u*v)' = (v'u + u'v)
293  Evaluation& operator*=(const Evaluation& other)
294  {
295  assert(size() == other.size());
296 
297  // while the values are multiplied, the derivatives follow the product rule,
298  // i.e., (u*v)' = (v'u + u'v).
299  const ValueType u = this->value();
300  const ValueType v = other.value();
301 
302  // value
303  data_[valuepos_()] *= v ;
304 
305  // derivatives
306  for (int i = dstart_(); i < dend_(); ++i)
307  data_[i] = data_[i] * v + other.data_[i] * u;
308 
309  return *this;
310  }
311 
312  // m(c*u)' = c*u'
313  template <class RhsValueType>
314  Evaluation& operator*=(const RhsValueType& other)
315  {
316  for (int i = 0; i < length_(); ++i)
317  data_[i] *= other;
318 
319  return *this;
320  }
321 
322  // m(u*v)' = (vu' - uv')/v^2
323  Evaluation& operator/=(const Evaluation& other)
324  {
325  assert(size() == other.size());
326 
327  // values are divided, derivatives follow the rule for division, i.e., (u/v)' = (v'u -
328  // u'v)/v^2.
329  ValueType& u = data_[valuepos_()];
330  const ValueType& v = other.value();
331  for (int idx = dstart_(); idx < dend_(); ++idx) {
332  const ValueType& uPrime = data_[idx];
333  const ValueType& vPrime = other.data_[idx];
334 
335  data_[idx] = (v*uPrime - u*vPrime)/(v*v);
336  }
337  u /= v;
338 
339  return *this;
340  }
341 
342  // divide value and derivatives by value of other
343  template <class RhsValueType>
344  Evaluation& operator/=(const RhsValueType& other)
345  {
346  const ValueType tmp = 1.0/other;
347 
348  for (int i = 0; i < length_(); ++i)
349  data_[i] *= tmp;
350 
351  return *this;
352  }
353 
354  // add two evaluation objects
355  Evaluation operator+(const Evaluation& other) const
356  {
357  assert(size() == other.size());
358 
359  Evaluation result(*this);
360 
361  result += other;
362 
363  return result;
364  }
365 
366  // add constant to this object
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  // subtract two evaluation objects
378  Evaluation operator-(const Evaluation& other) const
379  {
380  assert(size() == other.size());
381 
382  Evaluation result(*this);
383 
384  result -= other;
385 
386  return result;
387  }
388 
389  // subtract constant from evaluation object
390  template <class RhsValueType>
391  Evaluation operator-(const RhsValueType& other) const
392  {
393  Evaluation result(*this);
394 
395  result -= other;
396 
397  return result;
398  }
399 
400  // negation (unary minus) operator
401  Evaluation operator-() const
402  {
403  Evaluation result(*this);
404 
405  // set value and derivatives to negative
406  for (int i = 0; i < length_(); ++i)
407  result.data_[i] = - data_[i];
408 
409  return result;
410  }
411 
412  Evaluation operator*(const Evaluation& other) const
413  {
414  assert(size() == other.size());
415 
416  Evaluation result(*this);
417 
418  result *= other;
419 
420  return result;
421  }
422 
423  template <class RhsValueType>
424  Evaluation operator*(const RhsValueType& other) const
425  {
426  Evaluation result(*this);
427 
428  result *= other;
429 
430  return result;
431  }
432 
433  Evaluation operator/(const Evaluation& other) const
434  {
435  assert(size() == other.size());
436 
437  Evaluation result(*this);
438 
439  result /= other;
440 
441  return result;
442  }
443 
444  template <class RhsValueType>
445  Evaluation operator/(const RhsValueType& other) const
446  {
447  Evaluation result(*this);
448 
449  result /= other;
450 
451  return result;
452  }
453 
454  template <class RhsValueType>
455  Evaluation& operator=(const RhsValueType& other)
456  {
457  setValue( other );
458  clearDerivatives();
459 
460  return *this;
461  }
462 
463  // copy assignment from evaluation
464  Evaluation& operator=(const Evaluation& other) = default;
465 
466  template <class RhsValueType>
467  bool operator==(const RhsValueType& other) const
468  { return value() == other; }
469 
470  bool operator==(const Evaluation& other) const
471  {
472  assert(size() == other.size());
473 
474  for (int idx = 0; idx < length_(); ++idx) {
475  if (data_[idx] != other.data_[idx]) {
476  return false;
477  }
478  }
479  return true;
480  }
481 
482  bool operator!=(const Evaluation& other) const
483  { return !operator==(other); }
484 
485  template <class RhsValueType>
486  bool operator!=(const RhsValueType& other) const
487  { return !operator==(other); }
488 
489  template <class RhsValueType>
490  bool operator>(RhsValueType other) const
491  { return value() > other; }
492 
493  bool operator>(const Evaluation& other) const
494  {
495  assert(size() == other.size());
496 
497  return value() > other.value();
498  }
499 
500  template <class RhsValueType>
501  bool operator<(RhsValueType other) const
502  { return value() < other; }
503 
504  bool operator<(const Evaluation& other) const
505  {
506  assert(size() == other.size());
507 
508  return value() < other.value();
509  }
510 
511  template <class RhsValueType>
512  bool operator>=(RhsValueType other) const
513  { return value() >= other; }
514 
515  bool operator>=(const Evaluation& other) const
516  {
517  assert(size() == other.size());
518 
519  return value() >= other.value();
520  }
521 
522  template <class RhsValueType>
523  bool operator<=(RhsValueType other) const
524  { return value() <= other; }
525 
526  bool operator<=(const Evaluation& other) const
527  {
528  assert(size() == other.size());
529 
530  return value() <= other.value();
531  }
532 
533  // return value of variable
534  const ValueType& value() const
535  { return data_[valuepos_()]; }
536 
537  // set value of variable
538  template <class RhsValueType>
539  void setValue(const RhsValueType& val)
540  { data_[valuepos_()] = val; }
541 
542  // return varIdx'th derivative
543  const ValueType& derivative(int varIdx) const
544  {
545  assert(0 <= varIdx && varIdx < size());
546 
547  return data_[dstart_() + varIdx];
548  }
549 
550  // set derivative at position varIdx
551  void setDerivative(int varIdx, const ValueType& derVal)
552  {
553  assert(0 <= varIdx && varIdx < size());
554 
555  data_[dstart_() + varIdx] = derVal;
556  }
557 
558 private:
559 
560  FastSmallVector<ValueT, staticSize> data_;
561 };
562 
563 template <class Scalar, unsigned staticSize = 0>
564 using DynamicEvaluation = Evaluation<Scalar, DynamicSize, staticSize>;
565 
566 } // namespace DenseAd
567 
568 template <class Scalar, unsigned staticSize>
569 DenseAd::Evaluation<Scalar, -1, staticSize> constant(int numDerivatives, const Scalar& value)
570 { return DenseAd::Evaluation<Scalar, -1, staticSize>::createConstant(numDerivatives, value); }
571 
572 template <class Scalar, unsigned staticSize>
573 DenseAd::Evaluation<Scalar, -1, staticSize> variable(int numDerivatives, const Scalar& value, unsigned idx)
574 { return DenseAd::Evaluation<Scalar, -1, staticSize>::createVariable(numDerivatives, value, idx); }
575 
576 } // namespace Opm
577 
578 #endif // OPM_DENSEAD_EVALUATION_DYNAMIC_HPP
Representation of an evaluation of a function and its derivatives w.r.t.
An implementation of vector/array based on small object optimization.
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 dstart_() const
start index for derivatives
Definition: DynamicEvaluation.hpp:79
Evaluation(const Evaluation &other)=default
copy other function evaluation
int length_() const
length of internal data vector
Definition: DynamicEvaluation.hpp:71
Evaluation(Evaluation &&other)
move other function evaluation (this only makes sense for dynamically allocated Evaluations)
Definition: DynamicEvaluation.hpp:105
void checkDefined_() const
instruct valgrind to check that the value and all derivatives of the Evaluation object are well-defin...
Definition: DynamicEvaluation.hpp:87
int size() const
number of derivatives
Definition: DynamicEvaluation.hpp:66
int dend_() const
end+1 index for derivatives
Definition: DynamicEvaluation.hpp:82
Evaluation & operator=(Evaluation &&other)
move assignment
Definition: DynamicEvaluation.hpp:110
ValueT ValueType
field type
Definition: DynamicEvaluation.hpp:63
Evaluation()
default constructor
Definition: DynamicEvaluation.hpp:97
constexpr int valuepos_() const
position index for value
Definition: DynamicEvaluation.hpp:76
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 dend_() const
end+1 index for derivatives
Definition: Evaluation.hpp:85
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