My Project
Loading...
Searching...
No Matches
Evaluation.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_HPP
33#define OPM_DENSEAD_EVALUATION_HPP
34
35#ifndef NDEBUG
37#endif
38
39#include <array>
40#include <cassert>
41#include <iosfwd>
42#include <stdexcept>
43
44namespace Opm {
45namespace DenseAd {
46
49static constexpr int DynamicSize = -1;
50
55template <class ValueT, int numDerivs, unsigned staticSize = 0>
57{
58public:
61 static const int numVars = numDerivs;
62
64 typedef ValueT ValueType;
65
67 constexpr int size() const
68 { return numDerivs; }
69
70protected:
72 constexpr int length_() const
73 { return size() + 1; }
74
75
77 constexpr int valuepos_() const
78 { return 0; }
80 constexpr int dstart_() const
81 { return 1; }
83 constexpr int dend_() const
84 { return length_(); }
85
88 void checkDefined_() const
89 {
90#ifndef NDEBUG
91 for (const auto& v: data_)
92 Valgrind::CheckDefined(v);
93#endif
94 }
95
96public:
98 Evaluation() : data_()
99 {}
100
102 Evaluation(const Evaluation& other) = default;
103
104
105 // create an evaluation which represents a constant function
106 //
107 // i.e., f(x) = c. this implies an evaluation with the given value and all
108 // derivatives being zero.
109 template <class RhsValueType>
110 Evaluation(const RhsValueType& c)
111 {
112 setValue(c);
113 clearDerivatives();
114
116 }
117
118 // create an evaluation which represents a constant function
119 //
120 // i.e., f(x) = c. this implies an evaluation with the given value and all
121 // derivatives being zero.
122 template <class RhsValueType>
123 Evaluation(const RhsValueType& c, int varPos)
124 {
125 // The variable position must be in represented by the given variable descriptor
126 assert(0 <= varPos && varPos < size());
127
128 setValue( c );
129 clearDerivatives();
130
131 data_[varPos + dstart_()] = 1.0;
132
134 }
135
136 // set all derivatives to zero
137 void clearDerivatives()
138 {
139 for (int i = dstart_(); i < dend_(); ++i)
140 data_[i] = 0.0;
141 }
142
143 // create an uninitialized Evaluation object that is compatible with the
144 // argument, but not initialized
145 //
146 // This basically boils down to the copy constructor without copying
147 // anything. If the number of derivatives is known at compile time, this
148 // is equivalent to creating an uninitialized object using the default
149 // constructor, while for dynamic evaluations, it creates an Evaluation
150 // object which exhibits the same number of derivatives as the argument.
151 static Evaluation createBlank(const Evaluation&)
152 { return Evaluation(); }
153
154 // create an Evaluation with value and all the derivatives to be zero
155 static Evaluation createConstantZero(const Evaluation&)
156 { return Evaluation(0.); }
157
158 // create an Evaluation with value to be one and all the derivatives to be zero
159 static Evaluation createConstantOne(const Evaluation&)
160 { return Evaluation(1.); }
161
162 // create a function evaluation for a "naked" depending variable (i.e., f(x) = x)
163 template <class RhsValueType>
164 static Evaluation createVariable(const RhsValueType& value, int varPos)
165 {
166 // copy function value and set all derivatives to 0, except for the variable
167 // which is represented by the value (which is set to 1.0)
168 return Evaluation(value, varPos);
169 }
170
171 template <class RhsValueType>
172 static Evaluation createVariable(int nVars, const RhsValueType& value, int varPos)
173 {
174 if (nVars != 0)
175 throw std::logic_error("This statically-sized evaluation can only represent objects"
176 " with 0 derivatives");
177
178 // copy function value and set all derivatives to 0, except for the variable
179 // which is represented by the value (which is set to 1.0)
180 return Evaluation(nVars, value, varPos);
181 }
182
183 template <class RhsValueType>
184 static Evaluation createVariable(const Evaluation&, const RhsValueType& value, int varPos)
185 {
186 // copy function value and set all derivatives to 0, except for the variable
187 // which is represented by the value (which is set to 1.0)
188 return Evaluation(value, varPos);
189 }
190
191
192 // "evaluate" a constant function (i.e. a function that does not depend on the set of
193 // relevant variables, f(x) = c).
194 template <class RhsValueType>
195 static Evaluation createConstant(int nVars, const RhsValueType& value)
196 {
197 if (nVars != 0)
198 throw std::logic_error("This statically-sized evaluation can only represent objects"
199 " with 0 derivatives");
200 return Evaluation(value);
201 }
202
203 // "evaluate" a constant function (i.e. a function that does not depend on the set of
204 // relevant variables, f(x) = c).
205 template <class RhsValueType>
206 static Evaluation createConstant(const RhsValueType& value)
207 {
208 return Evaluation(value);
209 }
210
211 // "evaluate" a constant function (i.e. a function that does not depend on the set of
212 // relevant variables, f(x) = c).
213 template <class RhsValueType>
214 static Evaluation createConstant(const Evaluation&, const RhsValueType& value)
215 {
216 return Evaluation(value);
217 }
218
219 // copy all derivatives from other
220 void copyDerivatives(const Evaluation& other)
221 {
222 assert(size() == other.size());
223
224 for (int i = dstart_(); i < dend_(); ++i)
225 data_[i] = other.data_[i];
226 }
227
228
229 // add value and derivatives from other to this values and derivatives
230 Evaluation& operator+=(const Evaluation& other)
231 {
232 assert(size() == other.size());
233
234 for (int i = 0; i < length_(); ++i)
235 data_[i] += other.data_[i];
236
237 return *this;
238 }
239
240 // add value from other to this values
241 template <class RhsValueType>
242 Evaluation& operator+=(const RhsValueType& other)
243 {
244 // value is added, derivatives stay the same
245 data_[valuepos_()] += other;
246
247 return *this;
248 }
249
250 // subtract other's value and derivatives from this values
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 // subtract other's value from this values
262 template <class RhsValueType>
263 Evaluation& operator-=(const RhsValueType& other)
264 {
265 // for constants, values are subtracted, derivatives stay the same
266 data_[valuepos_()] -= other;
267
268 return *this;
269 }
270
271 // multiply values and apply chain rule to derivatives: (u*v)' = (v'u + u'v)
272 Evaluation& operator*=(const Evaluation& other)
273 {
274 assert(size() == other.size());
275
276 // while the values are multiplied, the derivatives follow the product rule,
277 // i.e., (u*v)' = (v'u + u'v).
278 const ValueType u = this->value();
279 const ValueType v = other.value();
280
281 // value
282 data_[valuepos_()] *= v ;
283
284 // derivatives
285 for (int i = dstart_(); i < dend_(); ++i)
286 data_[i] = data_[i] * v + other.data_[i] * u;
287
288 return *this;
289 }
290
291 // m(c*u)' = c*u'
292 template <class RhsValueType>
293 Evaluation& operator*=(const RhsValueType& other)
294 {
295 for (int i = 0; i < length_(); ++i)
296 data_[i] *= other;
297
298 return *this;
299 }
300
301 // m(u*v)' = (vu' - uv')/v^2
302 Evaluation& operator/=(const Evaluation& other)
303 {
304 assert(size() == other.size());
305
306 // values are divided, derivatives follow the rule for division, i.e., (u/v)' = (v'u -
307 // u'v)/v^2.
308 ValueType& u = data_[valuepos_()];
309 const ValueType& v = other.value();
310 for (int idx = dstart_(); idx < dend_(); ++idx) {
311 const ValueType& uPrime = data_[idx];
312 const ValueType& vPrime = other.data_[idx];
313
314 data_[idx] = (v*uPrime - u*vPrime)/(v*v);
315 }
316 u /= v;
317
318 return *this;
319 }
320
321 // divide value and derivatives by value of other
322 template <class RhsValueType>
323 Evaluation& operator/=(const RhsValueType& other)
324 {
325 const ValueType tmp = 1.0/other;
326
327 for (int i = 0; i < length_(); ++i)
328 data_[i] *= tmp;
329
330 return *this;
331 }
332
333 // add two evaluation objects
334 Evaluation operator+(const Evaluation& other) const
335 {
336 assert(size() == other.size());
337
338 Evaluation result(*this);
339
340 result += other;
341
342 return result;
343 }
344
345 // add constant to this object
346 template <class RhsValueType>
347 Evaluation operator+(const RhsValueType& other) const
348 {
349 Evaluation result(*this);
350
351 result += other;
352
353 return result;
354 }
355
356 // subtract two evaluation objects
357 Evaluation operator-(const Evaluation& other) const
358 {
359 assert(size() == other.size());
360
361 Evaluation result(*this);
362
363 result -= other;
364
365 return result;
366 }
367
368 // subtract constant from evaluation object
369 template <class RhsValueType>
370 Evaluation operator-(const RhsValueType& other) const
371 {
372 Evaluation result(*this);
373
374 result -= other;
375
376 return result;
377 }
378
379 // negation (unary minus) operator
380 Evaluation operator-() const
381 {
382 Evaluation result;
383
384 // set value and derivatives to negative
385 for (int i = 0; i < length_(); ++i)
386 result.data_[i] = - data_[i];
387
388 return result;
389 }
390
391 Evaluation operator*(const Evaluation& other) const
392 {
393 assert(size() == other.size());
394
395 Evaluation result(*this);
396
397 result *= other;
398
399 return result;
400 }
401
402 template <class RhsValueType>
403 Evaluation operator*(const RhsValueType& other) const
404 {
405 Evaluation result(*this);
406
407 result *= other;
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 template <class RhsValueType>
434 Evaluation& operator=(const RhsValueType& other)
435 {
436 setValue( other );
437 clearDerivatives();
438
439 return *this;
440 }
441
442 // copy assignment from evaluation
443 Evaluation& operator=(const Evaluation& other) = default;
444
445 template <class RhsValueType>
446 bool operator==(const RhsValueType& other) const
447 { return value() == other; }
448
449 bool operator==(const Evaluation& other) const
450 {
451 assert(size() == other.size());
452
453 for (int idx = 0; idx < length_(); ++idx) {
454 if (data_[idx] != other.data_[idx]) {
455 return false;
456 }
457 }
458 return true;
459 }
460
461 bool operator!=(const Evaluation& other) const
462 { return !operator==(other); }
463
464 template <class RhsValueType>
465 bool operator!=(const RhsValueType& other) const
466 { return !operator==(other); }
467
468 template <class RhsValueType>
469 bool operator>(RhsValueType other) const
470 { return value() > other; }
471
472 bool operator>(const Evaluation& other) const
473 {
474 assert(size() == other.size());
475
476 return value() > other.value();
477 }
478
479 template <class RhsValueType>
480 bool operator<(RhsValueType other) const
481 { return value() < other; }
482
483 bool operator<(const Evaluation& other) const
484 {
485 assert(size() == other.size());
486
487 return value() < other.value();
488 }
489
490 template <class RhsValueType>
491 bool operator>=(RhsValueType other) const
492 { return value() >= other; }
493
494 bool operator>=(const Evaluation& other) const
495 {
496 assert(size() == other.size());
497
498 return value() >= other.value();
499 }
500
501 template <class RhsValueType>
502 bool operator<=(RhsValueType other) const
503 { return value() <= other; }
504
505 bool operator<=(const Evaluation& other) const
506 {
507 assert(size() == other.size());
508
509 return value() <= other.value();
510 }
511
512 // return value of variable
513 const ValueType& value() const
514 { return data_[valuepos_()]; }
515
516 // set value of variable
517 template <class RhsValueType>
518 void setValue(const RhsValueType& val)
519 { data_[valuepos_()] = val; }
520
521 // return varIdx'th derivative
522 const ValueType& derivative(int varIdx) const
523 {
524 assert(0 <= varIdx && varIdx < size());
525
526 return data_[dstart_() + varIdx];
527 }
528
529 // set derivative at position varIdx
530 void setDerivative(int varIdx, const ValueType& derVal)
531 {
532 assert(0 <= varIdx && varIdx < size());
533
534 data_[dstart_() + varIdx] = derVal;
535 }
536
537 template<class Serializer>
538 void serializeOp(Serializer& serializer)
539 {
540 serializer(data_);
541 }
542
543private:
544 std::array<ValueT, numDerivs + 1> data_;
545};
546
547// the generic operators are only required for the unspecialized case
548template <class RhsValueType, class ValueType, int numVars, unsigned staticSize>
549bool operator<(const RhsValueType& a, const Evaluation<ValueType, numVars, staticSize>& b)
550{ return b > a; }
551
552template <class RhsValueType, class ValueType, int numVars, unsigned staticSize>
553bool operator>(const RhsValueType& a, const Evaluation<ValueType, numVars, staticSize>& b)
554{ return b < a; }
555
556template <class RhsValueType, class ValueType, int numVars, unsigned staticSize>
557bool operator<=(const RhsValueType& a, const Evaluation<ValueType, numVars, staticSize>& b)
558{ return b >= a; }
559
560template <class RhsValueType, class ValueType, int numVars, unsigned staticSize>
561bool operator>=(const RhsValueType& a, const Evaluation<ValueType, numVars, staticSize>& b)
562{ return b <= a; }
563
564template <class RhsValueType, class ValueType, int numVars, unsigned staticSize>
565bool operator!=(const RhsValueType& a, const Evaluation<ValueType, numVars, staticSize>& b)
566{ return a != b.value(); }
567
568template <class RhsValueType, class ValueType, int numVars, unsigned staticSize>
569Evaluation<ValueType, numVars, staticSize> operator+(const RhsValueType& a, const Evaluation<ValueType, numVars, staticSize>& b)
570{
571 Evaluation<ValueType, numVars, staticSize> result(b);
572 result += a;
573 return result;
574}
575
576template <class RhsValueType, class ValueType, int numVars, unsigned staticSize>
577Evaluation<ValueType, numVars, staticSize> operator-(const RhsValueType& a, const Evaluation<ValueType, numVars, staticSize>& b)
578{
579 return -(b - a);
580}
581
582template <class RhsValueType, class ValueType, int numVars, unsigned staticSize>
583Evaluation<ValueType, numVars, staticSize> operator/(const RhsValueType& a, const Evaluation<ValueType, numVars, staticSize>& b)
584{
585 Evaluation<ValueType, numVars, staticSize> tmp(a);
586 tmp /= b;
587 return tmp;
588}
589
590template <class RhsValueType, class ValueType, int numVars, unsigned staticSize>
591Evaluation<ValueType, numVars, staticSize> operator*(const RhsValueType& a, const Evaluation<ValueType, numVars, staticSize>& b)
592{
593 Evaluation<ValueType, numVars, staticSize> result(b);
594 result *= a;
595 return result;
596}
597
598template<class T>
600{
601 static constexpr bool value = false;
602};
603
604template <class ValueType, int numVars, unsigned staticSize>
605struct is_evaluation<Evaluation<ValueType,numVars,staticSize>>
606{
607 static constexpr bool value = true;
608};
609
610template <class ValueType, int numVars, unsigned staticSize>
611void printEvaluation(std::ostream& os,
613 bool withDer = false);
614
615template <class ValueType, int numVars, unsigned staticSize>
616std::ostream& operator<<(std::ostream& os, const Evaluation<ValueType, numVars, staticSize>& eval)
617{
619 printEvaluation(os, eval.value(), false);
620 else
621 printEvaluation(os, eval, true);
622
623 return os;
624}
625
626} // namespace DenseAd
627} // namespace Opm
628
630
631#endif // OPM_DENSEAD_EVALUATION_HPP
This file includes all specializations for the dense-AD Evaluation class.
Some templates to wrap the valgrind client request macros.
Represents a function evaluation and its derivatives w.r.t.
Definition Evaluation.hpp:57
Evaluation()
default constructor
Definition Evaluation.hpp:98
ValueT ValueType
field type
Definition Evaluation.hpp:64
void checkDefined_() const
instruct valgrind to check that the value and all derivatives of the Evaluation object are well-defin...
Definition Evaluation.hpp:88
static const int numVars
the template argument which specifies the number of derivatives (-1 == "DynamicSize" means runtime de...
Definition Evaluation.hpp:61
constexpr int size() const
number of derivatives
Definition Evaluation.hpp:67
constexpr int valuepos_() const
position index for value
Definition Evaluation.hpp:77
Evaluation(const Evaluation &other)=default
copy other function evaluation
constexpr int dend_() const
end+1 index for derivatives
Definition Evaluation.hpp:83
constexpr int length_() const
length of internal data vector
Definition Evaluation.hpp:72
constexpr int dstart_() const
start index for derivatives
Definition Evaluation.hpp:80
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30
Definition Evaluation.hpp:600