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