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