Math.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_LOCAL_AD_MATH_HPP
33 #define OPM_LOCAL_AD_MATH_HPP
34 
35 #include "Evaluation.hpp"
36 
38 
39 namespace Opm {
40 namespace DenseAd {
41 // forward declaration of the Evaluation template class
42 template <class ValueT, int numVars>
43 class Evaluation;
44 
45 // provide some algebraic functions
46 template <class ValueType, int numVars>
47 Evaluation<ValueType, numVars> abs(const Evaluation<ValueType, numVars>& x)
48 { return (x > 0.0)?x:-x; }
49 
50 template <class ValueType, int numVars>
51 Evaluation<ValueType, numVars> min(const Evaluation<ValueType, numVars>& x1,
52  const Evaluation<ValueType, numVars>& x2)
53 { return (x1 < x2)?x1:x2; }
54 
55 template <class Arg1ValueType, class ValueType, int numVars>
56 Evaluation<ValueType, numVars> min(const Arg1ValueType& x1,
57  const Evaluation<ValueType, numVars>& x2)
58 { return (x1 < x2)?x1:x2; }
59 
60 template <class ValueType, int numVars, class Arg2ValueType>
61 Evaluation<ValueType, numVars> min(const Evaluation<ValueType, numVars>& x1,
62  const Arg2ValueType& x2)
63 { return min(x2, x1); }
64 
65 template <class ValueType, int numVars>
66 Evaluation<ValueType, numVars> max(const Evaluation<ValueType, numVars>& x1,
67  const Evaluation<ValueType, numVars>& x2)
68 { return (x1 > x2)?x1:x2; }
69 
70 template <class Arg1ValueType, class ValueType, int numVars>
71 Evaluation<ValueType, numVars> max(const Arg1ValueType& x1,
72  const Evaluation<ValueType, numVars>& x2)
73 { return (x1 > x2)?x1:x2; }
74 
75 template <class ValueType, int numVars, class Arg2ValueType>
76 Evaluation<ValueType, numVars> max(const Evaluation<ValueType, numVars>& x1,
77  const Arg2ValueType& x2)
78 { return max(x2, x1); }
79 
80 template <class ValueType, int numVars>
81 Evaluation<ValueType, numVars> tan(const Evaluation<ValueType, numVars>& x)
82 {
83  typedef MathToolbox<ValueType> ValueTypeToolbox;
84 
85  Evaluation<ValueType, numVars> result;
86 
87  const ValueType& tmp = ValueTypeToolbox::tan(x.value());
88  result.setValue(tmp);
89 
90  // derivatives use the chain rule
91  const ValueType& df_dx = 1 + tmp*tmp;
92  for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
93  result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
94 
95  return result;
96 }
97 
98 template <class ValueType, int numVars>
99 Evaluation<ValueType, numVars> atan(const Evaluation<ValueType, numVars>& x)
100 {
101  typedef MathToolbox<ValueType> ValueTypeToolbox;
102 
103  Evaluation<ValueType, numVars> result;
104 
105  result.setValue(ValueTypeToolbox::atan(x.value()));
106 
107  // derivatives use the chain rule
108  const ValueType& df_dx = 1/(1 + x.value()*x.value());
109  for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
110  result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
111 
112  return result;
113 }
114 
115 template <class ValueType, int numVars>
116 Evaluation<ValueType, numVars> atan2(const Evaluation<ValueType, numVars>& x,
117  const Evaluation<ValueType, numVars>& y)
118 {
119  typedef MathToolbox<ValueType> ValueTypeToolbox;
120 
121  Evaluation<ValueType, numVars> result;
122 
123  result.setValue(ValueTypeToolbox::atan2(x.value(), y.value()));
124 
125  // derivatives use the chain rule
126  const ValueType& alpha = 1/(1 + (x.value()*x.value())/(y.value()*y.value()));
127  for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) {
128  result.setDerivative(curVarIdx,
129  alpha/(y.value()*y.value())
130  *(x.derivative(curVarIdx)*y.value() - x.value()*y.derivative(curVarIdx)));
131  }
132 
133  return result;
134 }
135 
136 template <class ValueType, int numVars>
137 Evaluation<ValueType, numVars> sin(const Evaluation<ValueType, numVars>& x)
138 {
139  typedef MathToolbox<ValueType> ValueTypeToolbox;
140 
141  Evaluation<ValueType, numVars> result;
142 
143  result.setValue(ValueTypeToolbox::sin(x.value()));
144 
145  // derivatives use the chain rule
146  const ValueType& df_dx = ValueTypeToolbox::cos(x.value());
147  for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
148  result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
149 
150  return result;
151 }
152 
153 template <class ValueType, int numVars>
154 Evaluation<ValueType, numVars> asin(const Evaluation<ValueType, numVars>& x)
155 {
156  typedef MathToolbox<ValueType> ValueTypeToolbox;
157 
158  Evaluation<ValueType, numVars> result;
159 
160  result.setValue(ValueTypeToolbox::asin(x.value()));
161 
162  // derivatives use the chain rule
163  const ValueType& df_dx = 1.0/ValueTypeToolbox::sqrt(1 - x.value()*x.value());
164  for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
165  result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
166 
167  return result;
168 }
169 
170 template <class ValueType, int numVars>
171 Evaluation<ValueType, numVars> cos(const Evaluation<ValueType, numVars>& x)
172 {
173  typedef MathToolbox<ValueType> ValueTypeToolbox;
174 
175  Evaluation<ValueType, numVars> result;
176 
177  result.setValue(ValueTypeToolbox::cos(x.value()));
178 
179  // derivatives use the chain rule
180  const ValueType& df_dx = -ValueTypeToolbox::sin(x.value());
181  for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
182  result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
183 
184  return result;
185 }
186 
187 template <class ValueType, int numVars>
188 Evaluation<ValueType, numVars> acos(const Evaluation<ValueType, numVars>& x)
189 {
190  typedef MathToolbox<ValueType> ValueTypeToolbox;
191 
192  Evaluation<ValueType, numVars> result;
193 
194  result.setValue(ValueTypeToolbox::acos(x.value()));
195 
196  // derivatives use the chain rule
197  const ValueType& df_dx = - 1.0/ValueTypeToolbox::sqrt(1 - x.value()*x.value());
198  for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
199  result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
200 
201  return result;
202 }
203 
204 template <class ValueType, int numVars>
205 Evaluation<ValueType, numVars> sqrt(const Evaluation<ValueType, numVars>& x)
206 {
207  typedef MathToolbox<ValueType> ValueTypeToolbox;
208 
209  Evaluation<ValueType, numVars> result;
210 
211  const ValueType& sqrt_x = ValueTypeToolbox::sqrt(x.value());
212  result.setValue(sqrt_x);
213 
214  // derivatives use the chain rule
215  ValueType df_dx = 0.5/sqrt_x;
216  for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) {
217  result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
218  }
219 
220  return result;
221 }
222 
223 template <class ValueType, int numVars>
224 Evaluation<ValueType, numVars> exp(const Evaluation<ValueType, numVars>& x)
225 {
226  typedef MathToolbox<ValueType> ValueTypeToolbox;
227  Evaluation<ValueType, numVars> result;
228 
229  const ValueType& exp_x = ValueTypeToolbox::exp(x.value());
230  result.setValue(exp_x);
231 
232  // derivatives use the chain rule
233  const ValueType& df_dx = exp_x;
234  for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
235  result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
236 
237  return result;
238 }
239 
240 // exponentiation of arbitrary base with a fixed constant
241 template <class ValueType, int numVars, class ExpType>
242 Evaluation<ValueType, numVars> pow(const Evaluation<ValueType, numVars>& base,
243  const ExpType& exp)
244 {
245  typedef MathToolbox<ValueType> ValueTypeToolbox;
246  Evaluation<ValueType, numVars> result;
247 
248  const ValueType& pow_x = ValueTypeToolbox::pow(base.value(), exp);
249  result.setValue(pow_x);
250 
251  if (base == 0.0) {
252  // we special case the base 0 case because 0.0 is in the valid range of the
253  // base but the generic code leads to NaNs.
254  result = 0.0;
255  }
256  else {
257  // derivatives use the chain rule
258  const ValueType& df_dx = pow_x/base.value()*exp;
259  for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
260  result.setDerivative(curVarIdx, df_dx*base.derivative(curVarIdx));
261  }
262 
263  return result;
264 }
265 
266 // exponentiation of constant base with an arbitrary exponent
267 template <class BaseType, class ValueType, int numVars>
268 Evaluation<ValueType, numVars> pow(const BaseType& base,
269  const Evaluation<ValueType, numVars>& exp)
270 {
271  typedef MathToolbox<ValueType> ValueTypeToolbox;
272 
273  Evaluation<ValueType, numVars> result;
274 
275  if (base == 0.0) {
276  // we special case the base 0 case because 0.0 is in the valid range of the
277  // base but the generic code leads to NaNs.
278  result = 0.0;
279  }
280  else {
281  const ValueType& lnBase = ValueTypeToolbox::log(base);
282  result.setValue(ValueTypeToolbox::exp(lnBase*exp.value()));
283 
284  // derivatives use the chain rule
285  const ValueType& df_dx = lnBase*result.value();
286  for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
287  result.setDerivative(curVarIdx, df_dx*exp.derivative(curVarIdx));
288  }
289 
290  return result;
291 }
292 
293 // this is the most expensive power function. Computationally it is pretty expensive, so
294 // one of the above two variants above should be preferred if possible.
295 template <class ValueType, int numVars>
296 Evaluation<ValueType, numVars> pow(const Evaluation<ValueType, numVars>& base,
297  const Evaluation<ValueType, numVars>& exp)
298 {
299  typedef MathToolbox<ValueType> ValueTypeToolbox;
300 
301  Evaluation<ValueType, numVars> result;
302 
303  if (base == 0.0) {
304  // we special case the base 0 case because 0.0 is in the valid range of the
305  // base but the generic code leads to NaNs.
306  result = 0.0;
307  }
308  else {
309  ValueType valuePow = ValueTypeToolbox::pow(base.value(), exp.value());
310  result.setValue(valuePow);
311 
312  // use the chain rule for the derivatives. since both, the base and the exponent can
313  // potentially depend on the variable set, calculating these is quite elaborate...
314  const ValueType& f = base.value();
315  const ValueType& g = exp.value();
316  const ValueType& logF = ValueTypeToolbox::log(f);
317  for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) {
318  const ValueType& fPrime = base.derivative(curVarIdx);
319  const ValueType& gPrime = exp.derivative(curVarIdx);
320  result.setDerivative(curVarIdx, (g*fPrime/f + logF*gPrime) * valuePow);
321  }
322  }
323 
324  return result;
325 }
326 
327 template <class ValueType, int numVars>
328 Evaluation<ValueType, numVars> log(const Evaluation<ValueType, numVars>& x)
329 {
330  typedef MathToolbox<ValueType> ValueTypeToolbox;
331 
332  Evaluation<ValueType, numVars> result;
333 
334  result.setValue(ValueTypeToolbox::log(x.value()));
335 
336  // derivatives use the chain rule
337  const ValueType& df_dx = 1/x.value();
338  for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
339  result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
340 
341  return result;
342 }
343 
344 } // namespace DenseAd
345 
346 // a kind of traits class for the automatic differentiation case. (The toolbox for the
347 // scalar case is provided by the MathToolbox.hpp header file.)
348 template <class ValueT, int numVars>
349 struct MathToolbox<Opm::DenseAd::Evaluation<ValueT, numVars> >
350 {
351 private:
352 public:
353  typedef ValueT ValueType;
355  typedef typename InnerToolbox::Scalar Scalar;
357 
358  static ValueType value(const Evaluation& eval)
359  { return eval.value(); }
360 
361  static decltype(InnerToolbox::scalarValue(0.0)) scalarValue(const Evaluation& eval)
362  { return InnerToolbox::scalarValue(eval.value()); }
363 
364  static Evaluation createConstant(ValueType value)
365  { return Evaluation::createConstant(value); }
366 
367  static Evaluation createVariable(ValueType value, int varIdx)
368  { return Evaluation::createVariable(value, varIdx); }
369 
370  template <class LhsEval>
371  static typename std::enable_if<std::is_same<Evaluation, LhsEval>::value,
372  LhsEval>::type
373  decay(const Evaluation& eval)
374  { return eval; }
375 
376  template <class LhsEval>
377  static typename std::enable_if<std::is_same<Evaluation, LhsEval>::value,
378  LhsEval>::type
379  decay(const Evaluation&& eval)
380  { return eval; }
381 
382  template <class LhsEval>
383  static typename std::enable_if<std::is_floating_point<LhsEval>::value,
384  LhsEval>::type
385  decay(const Evaluation& eval)
386  { return eval.value(); }
387 
388  // comparison
389  static bool isSame(const Evaluation& a, const Evaluation& b, Scalar tolerance)
390  {
391  typedef MathToolbox<ValueType> ValueTypeToolbox;
392 
393  // make sure that the value of the evaluation is identical
394  if (!ValueTypeToolbox::isSame(a.value(), b.value(), tolerance))
395  return false;
396 
397  // make sure that the derivatives are identical
398  for (int curVarIdx = 0; curVarIdx < numVars; ++curVarIdx)
399  if (!ValueTypeToolbox::isSame(a.derivative(curVarIdx), b.derivative(curVarIdx), tolerance))
400  return false;
401 
402  return true;
403  }
404 
405  // arithmetic functions
406  template <class Arg1Eval, class Arg2Eval>
407  static Evaluation max(const Arg1Eval& arg1, const Arg2Eval& arg2)
408  { return Opm::DenseAd::max(arg1, arg2); }
409 
410  template <class Arg1Eval, class Arg2Eval>
411  static Evaluation min(const Arg1Eval& arg1, const Arg2Eval& arg2)
412  { return Opm::DenseAd::min(arg1, arg2); }
413 
414  static Evaluation abs(const Evaluation& arg)
415  { return Opm::DenseAd::abs(arg); }
416 
417  static Evaluation tan(const Evaluation& arg)
418  { return Opm::DenseAd::tan(arg); }
419 
420  static Evaluation atan(const Evaluation& arg)
421  { return Opm::DenseAd::atan(arg); }
422 
423  static Evaluation atan2(const Evaluation& arg1, const Evaluation& arg2)
424  { return Opm::DenseAd::atan2(arg1, arg2); }
425 
426  static Evaluation sin(const Evaluation& arg)
427  { return Opm::DenseAd::sin(arg); }
428 
429  static Evaluation asin(const Evaluation& arg)
430  { return Opm::DenseAd::asin(arg); }
431 
432  static Evaluation cos(const Evaluation& arg)
433  { return Opm::DenseAd::cos(arg); }
434 
435  static Evaluation acos(const Evaluation& arg)
436  { return Opm::DenseAd::acos(arg); }
437 
438  static Evaluation sqrt(const Evaluation& arg)
439  { return Opm::DenseAd::sqrt(arg); }
440 
441  static Evaluation exp(const Evaluation& arg)
442  { return Opm::DenseAd::exp(arg); }
443 
444  static Evaluation log(const Evaluation& arg)
445  { return Opm::DenseAd::log(arg); }
446 
447  template <class RhsValueType>
448  static Evaluation pow(const Evaluation& arg1, const RhsValueType& arg2)
449  { return Opm::DenseAd::pow(arg1, arg2); }
450 
451  template <class RhsValueType>
452  static Evaluation pow(const RhsValueType& arg1, const Evaluation& arg2)
453  { return Opm::DenseAd::pow(arg1, arg2); }
454 
455  static Evaluation pow(const Evaluation& arg1, const Evaluation& arg2)
456  { return Opm::DenseAd::pow(arg1, arg2); }
457 
458  static bool isfinite(const Evaluation& arg)
459  {
460  if (!InnerToolbox::isfinite(arg.value()))
461  return false;
462 
463  for (int i = 0; i < numVars; ++i)
464  if (!InnerToolbox::isfinite(arg.derivative(i)))
465  return false;
466 
467  return true;
468  }
469 
470  static bool isnan(const Evaluation& arg)
471  {
472  if (InnerToolbox::isnan(arg.value()))
473  return true;
474 
475  for (int i = 0; i < numVars; ++i)
476  if (InnerToolbox::isnan(arg.derivative(i)))
477  return true;
478 
479  return false;
480  }
481 };
482 
483 }
484 
485 #endif
A traits class which provides basic mathematical functions for arbitrary scalar floating point values...
static Scalar createVariable(Scalar value, unsigned)
Given a scalar value, return an evaluation of a linear function.
Definition: MathToolbox.hpp:111
static Scalar atan2(Scalar arg1, Scalar arg2)
The arcus tangens of a value.
Definition: MathToolbox.hpp:170
static Scalar cos(Scalar arg)
The cosine of a value.
Definition: MathToolbox.hpp:182
static bool isfinite(Scalar arg)
Return true iff the argument&#39;s value and all its derivatives are finite values.
Definition: MathToolbox.hpp:206
static Scalar sin(Scalar arg)
The sine of a value.
Definition: MathToolbox.hpp:174
static Scalar atan(Scalar arg)
The arcus tangens of a value.
Definition: MathToolbox.hpp:166
static Scalar acos(Scalar arg)
The arcus cosine of a value.
Definition: MathToolbox.hpp:186
static Scalar min(Scalar arg1, Scalar arg2)
The minimum of two arguments.
Definition: MathToolbox.hpp:154
Definition: Air_Mesitylene.hpp:33
static Scalar pow(Scalar base, Scalar exp)
Exponentiation to an arbitrary base.
Definition: MathToolbox.hpp:202
ScalarT ValueType
The type used to represent values.
Definition: MathToolbox.hpp:65
static Scalar log(Scalar arg)
The natural logarithm of a value.
Definition: MathToolbox.hpp:198
static bool isnan(Scalar arg)
Return true iff the argument&#39;s value or any of its derivatives are NaN values.
Definition: MathToolbox.hpp:210
static LhsEval decay(Scalar value)
Given a function evaluation, constrain it to its value (if necessary).
Definition: MathToolbox.hpp:126
static Scalar tan(Scalar arg)
The tangens of a value.
Definition: MathToolbox.hpp:162
Definition: MathToolbox.hpp:48
ScalarT Scalar
The type used to represent "primitive" scalar values.
Definition: MathToolbox.hpp:51
Opm::MathToolbox< Scalar > InnerToolbox
The toolbox for the type of value objects.
Definition: MathToolbox.hpp:74
static Scalar asin(Scalar arg)
The arcus sine of a value.
Definition: MathToolbox.hpp:178
static Scalar value(Scalar value)
Return the value of the function at a given evaluation point.
Definition: MathToolbox.hpp:82
static Scalar max(Scalar arg1, Scalar arg2)
The maximum of two arguments.
Definition: MathToolbox.hpp:150
static Scalar abs(Scalar arg)
The absolute value.
Definition: MathToolbox.hpp:158
static bool isSame(Scalar a, Scalar b, Scalar tolerance)
Returns true if two values are identical up to a specified tolerance.
Definition: MathToolbox.hpp:137
static Scalar sqrt(Scalar arg)
The square root of a value.
Definition: MathToolbox.hpp:190
static Scalar createConstant(Scalar value)
Given a scalar value, return an evaluation of a constant function.
Definition: MathToolbox.hpp:101
static Scalar exp(Scalar arg)
The natural exponentiation of a value.
Definition: MathToolbox.hpp:194
Represents a function evaluation and its derivatives w.r.t.
Definition: Evaluation.hpp:57
static Scalar scalarValue(Scalar value)
Return the primitive scalar value of a value object.
Definition: MathToolbox.hpp:91
Representation of an evaluation of a function and its derivatives w.r.t.