00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00032 #ifndef OPM_LOCAL_AD_MATH_HPP
00033 #define OPM_LOCAL_AD_MATH_HPP
00034
00035 #include "Evaluation.hpp"
00036
00037 #include <opm/material/common/MathToolbox.hpp>
00038
00039 namespace Opm {
00040 namespace DenseAd {
00041
00042 template <class ValueT, int numVars>
00043 class Evaluation;
00044
00045
00046 template <class ValueType, int numVars>
00047 Evaluation<ValueType, numVars> abs(const Evaluation<ValueType, numVars>& x)
00048 { return (x > 0.0)?x:-x; }
00049
00050 template <class ValueType, int numVars>
00051 Evaluation<ValueType, numVars> min(const Evaluation<ValueType, numVars>& x1,
00052 const Evaluation<ValueType, numVars>& x2)
00053 { return (x1 < x2)?x1:x2; }
00054
00055 template <class Arg1ValueType, class ValueType, int numVars>
00056 Evaluation<ValueType, numVars> min(const Arg1ValueType& x1,
00057 const Evaluation<ValueType, numVars>& x2)
00058 { return (x1 < x2)?x1:x2; }
00059
00060 template <class ValueType, int numVars, class Arg2ValueType>
00061 Evaluation<ValueType, numVars> min(const Evaluation<ValueType, numVars>& x1,
00062 const Arg2ValueType& x2)
00063 { return min(x2, x1); }
00064
00065 template <class ValueType, int numVars>
00066 Evaluation<ValueType, numVars> max(const Evaluation<ValueType, numVars>& x1,
00067 const Evaluation<ValueType, numVars>& x2)
00068 { return (x1 > x2)?x1:x2; }
00069
00070 template <class Arg1ValueType, class ValueType, int numVars>
00071 Evaluation<ValueType, numVars> max(const Arg1ValueType& x1,
00072 const Evaluation<ValueType, numVars>& x2)
00073 { return (x1 > x2)?x1:x2; }
00074
00075 template <class ValueType, int numVars, class Arg2ValueType>
00076 Evaluation<ValueType, numVars> max(const Evaluation<ValueType, numVars>& x1,
00077 const Arg2ValueType& x2)
00078 { return max(x2, x1); }
00079
00080 template <class ValueType, int numVars>
00081 Evaluation<ValueType, numVars> tan(const Evaluation<ValueType, numVars>& x)
00082 {
00083 typedef MathToolbox<ValueType> ValueTypeToolbox;
00084
00085 Evaluation<ValueType, numVars> result;
00086
00087 const ValueType& tmp = ValueTypeToolbox::tan(x.value());
00088 result.setValue(tmp);
00089
00090
00091 const ValueType& df_dx = 1 + tmp*tmp;
00092 for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
00093 result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
00094
00095 return result;
00096 }
00097
00098 template <class ValueType, int numVars>
00099 Evaluation<ValueType, numVars> atan(const Evaluation<ValueType, numVars>& x)
00100 {
00101 typedef MathToolbox<ValueType> ValueTypeToolbox;
00102
00103 Evaluation<ValueType, numVars> result;
00104
00105 result.setValue(ValueTypeToolbox::atan(x.value()));
00106
00107
00108 const ValueType& df_dx = 1/(1 + x.value()*x.value());
00109 for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
00110 result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
00111
00112 return result;
00113 }
00114
00115 template <class ValueType, int numVars>
00116 Evaluation<ValueType, numVars> atan2(const Evaluation<ValueType, numVars>& x,
00117 const Evaluation<ValueType, numVars>& y)
00118 {
00119 typedef MathToolbox<ValueType> ValueTypeToolbox;
00120
00121 Evaluation<ValueType, numVars> result;
00122
00123 result.setValue(ValueTypeToolbox::atan2(x.value(), y.value()));
00124
00125
00126 const ValueType& alpha = 1/(1 + (x.value()*x.value())/(y.value()*y.value()));
00127 for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) {
00128 result.setDerivative(curVarIdx,
00129 alpha/(y.value()*y.value())
00130 *(x.derivative(curVarIdx)*y.value() - x.value()*y.derivative(curVarIdx)));
00131 }
00132
00133 return result;
00134 }
00135
00136 template <class ValueType, int numVars>
00137 Evaluation<ValueType, numVars> sin(const Evaluation<ValueType, numVars>& x)
00138 {
00139 typedef MathToolbox<ValueType> ValueTypeToolbox;
00140
00141 Evaluation<ValueType, numVars> result;
00142
00143 result.setValue(ValueTypeToolbox::sin(x.value()));
00144
00145
00146 const ValueType& df_dx = ValueTypeToolbox::cos(x.value());
00147 for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
00148 result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
00149
00150 return result;
00151 }
00152
00153 template <class ValueType, int numVars>
00154 Evaluation<ValueType, numVars> asin(const Evaluation<ValueType, numVars>& x)
00155 {
00156 typedef MathToolbox<ValueType> ValueTypeToolbox;
00157
00158 Evaluation<ValueType, numVars> result;
00159
00160 result.setValue(ValueTypeToolbox::asin(x.value()));
00161
00162
00163 const ValueType& df_dx = 1.0/ValueTypeToolbox::sqrt(1 - x.value()*x.value());
00164 for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
00165 result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
00166
00167 return result;
00168 }
00169
00170 template <class ValueType, int numVars>
00171 Evaluation<ValueType, numVars> cos(const Evaluation<ValueType, numVars>& x)
00172 {
00173 typedef MathToolbox<ValueType> ValueTypeToolbox;
00174
00175 Evaluation<ValueType, numVars> result;
00176
00177 result.setValue(ValueTypeToolbox::cos(x.value()));
00178
00179
00180 const ValueType& df_dx = -ValueTypeToolbox::sin(x.value());
00181 for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
00182 result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
00183
00184 return result;
00185 }
00186
00187 template <class ValueType, int numVars>
00188 Evaluation<ValueType, numVars> acos(const Evaluation<ValueType, numVars>& x)
00189 {
00190 typedef MathToolbox<ValueType> ValueTypeToolbox;
00191
00192 Evaluation<ValueType, numVars> result;
00193
00194 result.setValue(ValueTypeToolbox::acos(x.value()));
00195
00196
00197 const ValueType& df_dx = - 1.0/ValueTypeToolbox::sqrt(1 - x.value()*x.value());
00198 for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
00199 result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
00200
00201 return result;
00202 }
00203
00204 template <class ValueType, int numVars>
00205 Evaluation<ValueType, numVars> sqrt(const Evaluation<ValueType, numVars>& x)
00206 {
00207 typedef MathToolbox<ValueType> ValueTypeToolbox;
00208
00209 Evaluation<ValueType, numVars> result;
00210
00211 const ValueType& sqrt_x = ValueTypeToolbox::sqrt(x.value());
00212 result.setValue(sqrt_x);
00213
00214
00215 ValueType df_dx = 0.5/sqrt_x;
00216 for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) {
00217 result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
00218 }
00219
00220 return result;
00221 }
00222
00223 template <class ValueType, int numVars>
00224 Evaluation<ValueType, numVars> exp(const Evaluation<ValueType, numVars>& x)
00225 {
00226 typedef MathToolbox<ValueType> ValueTypeToolbox;
00227 Evaluation<ValueType, numVars> result;
00228
00229 const ValueType& exp_x = ValueTypeToolbox::exp(x.value());
00230 result.setValue(exp_x);
00231
00232
00233 const ValueType& df_dx = exp_x;
00234 for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
00235 result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
00236
00237 return result;
00238 }
00239
00240
00241 template <class ValueType, int numVars, class ExpType>
00242 Evaluation<ValueType, numVars> pow(const Evaluation<ValueType, numVars>& base,
00243 const ExpType& exp)
00244 {
00245 typedef MathToolbox<ValueType> ValueTypeToolbox;
00246 Evaluation<ValueType, numVars> result;
00247
00248 const ValueType& pow_x = ValueTypeToolbox::pow(base.value(), exp);
00249 result.setValue(pow_x);
00250
00251 if (base == 0.0) {
00252
00253
00254 result = 0.0;
00255 }
00256 else {
00257
00258 const ValueType& df_dx = pow_x/base.value()*exp;
00259 for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
00260 result.setDerivative(curVarIdx, df_dx*base.derivative(curVarIdx));
00261 }
00262
00263 return result;
00264 }
00265
00266
00267 template <class BaseType, class ValueType, int numVars>
00268 Evaluation<ValueType, numVars> pow(const BaseType& base,
00269 const Evaluation<ValueType, numVars>& exp)
00270 {
00271 typedef MathToolbox<ValueType> ValueTypeToolbox;
00272
00273 Evaluation<ValueType, numVars> result;
00274
00275 if (base == 0.0) {
00276
00277
00278 result = 0.0;
00279 }
00280 else {
00281 const ValueType& lnBase = ValueTypeToolbox::log(base);
00282 result.setValue(ValueTypeToolbox::exp(lnBase*exp.value()));
00283
00284
00285 const ValueType& df_dx = lnBase*result.value();
00286 for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
00287 result.setDerivative(curVarIdx, df_dx*exp.derivative(curVarIdx));
00288 }
00289
00290 return result;
00291 }
00292
00293
00294
00295 template <class ValueType, int numVars>
00296 Evaluation<ValueType, numVars> pow(const Evaluation<ValueType, numVars>& base,
00297 const Evaluation<ValueType, numVars>& exp)
00298 {
00299 typedef MathToolbox<ValueType> ValueTypeToolbox;
00300
00301 Evaluation<ValueType, numVars> result;
00302
00303 if (base == 0.0) {
00304
00305
00306 result = 0.0;
00307 }
00308 else {
00309 ValueType valuePow = ValueTypeToolbox::pow(base.value(), exp.value());
00310 result.setValue(valuePow);
00311
00312
00313
00314 const ValueType& f = base.value();
00315 const ValueType& g = exp.value();
00316 const ValueType& logF = ValueTypeToolbox::log(f);
00317 for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx) {
00318 const ValueType& fPrime = base.derivative(curVarIdx);
00319 const ValueType& gPrime = exp.derivative(curVarIdx);
00320 result.setDerivative(curVarIdx, (g*fPrime/f + logF*gPrime) * valuePow);
00321 }
00322 }
00323
00324 return result;
00325 }
00326
00327 template <class ValueType, int numVars>
00328 Evaluation<ValueType, numVars> log(const Evaluation<ValueType, numVars>& x)
00329 {
00330 typedef MathToolbox<ValueType> ValueTypeToolbox;
00331
00332 Evaluation<ValueType, numVars> result;
00333
00334 result.setValue(ValueTypeToolbox::log(x.value()));
00335
00336
00337 const ValueType& df_dx = 1/x.value();
00338 for (int curVarIdx = 0; curVarIdx < result.size; ++curVarIdx)
00339 result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx));
00340
00341 return result;
00342 }
00343
00344 }
00345
00346
00347
00348 template <class ValueT, int numVars>
00349 struct MathToolbox<Opm::DenseAd::Evaluation<ValueT, numVars> >
00350 {
00351 private:
00352 public:
00353 typedef ValueT ValueType;
00354 typedef Opm::MathToolbox<ValueType> InnerToolbox;
00355 typedef typename InnerToolbox::Scalar Scalar;
00356 typedef Opm::DenseAd::Evaluation<ValueType, numVars> Evaluation;
00357
00358 static ValueType value(const Evaluation& eval)
00359 { return eval.value(); }
00360
00361 static decltype(InnerToolbox::scalarValue(0.0)) scalarValue(const Evaluation& eval)
00362 { return InnerToolbox::scalarValue(eval.value()); }
00363
00364 static Evaluation createConstant(ValueType value)
00365 { return Evaluation::createConstant(value); }
00366
00367 static Evaluation createVariable(ValueType value, int varIdx)
00368 { return Evaluation::createVariable(value, varIdx); }
00369
00370 template <class LhsEval>
00371 static typename std::enable_if<std::is_same<Evaluation, LhsEval>::value,
00372 LhsEval>::type
00373 decay(const Evaluation& eval)
00374 { return eval; }
00375
00376 template <class LhsEval>
00377 static typename std::enable_if<std::is_same<Evaluation, LhsEval>::value,
00378 LhsEval>::type
00379 decay(const Evaluation&& eval)
00380 { return eval; }
00381
00382 template <class LhsEval>
00383 static typename std::enable_if<std::is_floating_point<LhsEval>::value,
00384 LhsEval>::type
00385 decay(const Evaluation& eval)
00386 { return eval.value(); }
00387
00388
00389 static bool isSame(const Evaluation& a, const Evaluation& b, Scalar tolerance)
00390 {
00391 typedef MathToolbox<ValueType> ValueTypeToolbox;
00392
00393
00394 if (!ValueTypeToolbox::isSame(a.value(), b.value(), tolerance))
00395 return false;
00396
00397
00398 for (int curVarIdx = 0; curVarIdx < numVars; ++curVarIdx)
00399 if (!ValueTypeToolbox::isSame(a.derivative(curVarIdx), b.derivative(curVarIdx), tolerance))
00400 return false;
00401
00402 return true;
00403 }
00404
00405
00406 template <class Arg1Eval, class Arg2Eval>
00407 static Evaluation max(const Arg1Eval& arg1, const Arg2Eval& arg2)
00408 { return Opm::DenseAd::max(arg1, arg2); }
00409
00410 template <class Arg1Eval, class Arg2Eval>
00411 static Evaluation min(const Arg1Eval& arg1, const Arg2Eval& arg2)
00412 { return Opm::DenseAd::min(arg1, arg2); }
00413
00414 static Evaluation abs(const Evaluation& arg)
00415 { return Opm::DenseAd::abs(arg); }
00416
00417 static Evaluation tan(const Evaluation& arg)
00418 { return Opm::DenseAd::tan(arg); }
00419
00420 static Evaluation atan(const Evaluation& arg)
00421 { return Opm::DenseAd::atan(arg); }
00422
00423 static Evaluation atan2(const Evaluation& arg1, const Evaluation& arg2)
00424 { return Opm::DenseAd::atan2(arg1, arg2); }
00425
00426 static Evaluation sin(const Evaluation& arg)
00427 { return Opm::DenseAd::sin(arg); }
00428
00429 static Evaluation asin(const Evaluation& arg)
00430 { return Opm::DenseAd::asin(arg); }
00431
00432 static Evaluation cos(const Evaluation& arg)
00433 { return Opm::DenseAd::cos(arg); }
00434
00435 static Evaluation acos(const Evaluation& arg)
00436 { return Opm::DenseAd::acos(arg); }
00437
00438 static Evaluation sqrt(const Evaluation& arg)
00439 { return Opm::DenseAd::sqrt(arg); }
00440
00441 static Evaluation exp(const Evaluation& arg)
00442 { return Opm::DenseAd::exp(arg); }
00443
00444 static Evaluation log(const Evaluation& arg)
00445 { return Opm::DenseAd::log(arg); }
00446
00447 template <class RhsValueType>
00448 static Evaluation pow(const Evaluation& arg1, const RhsValueType& arg2)
00449 { return Opm::DenseAd::pow(arg1, arg2); }
00450
00451 template <class RhsValueType>
00452 static Evaluation pow(const RhsValueType& arg1, const Evaluation& arg2)
00453 { return Opm::DenseAd::pow(arg1, arg2); }
00454
00455 static Evaluation pow(const Evaluation& arg1, const Evaluation& arg2)
00456 { return Opm::DenseAd::pow(arg1, arg2); }
00457
00458 static bool isfinite(const Evaluation& arg)
00459 {
00460 if (!InnerToolbox::isfinite(arg.value()))
00461 return false;
00462
00463 for (int i = 0; i < numVars; ++i)
00464 if (!InnerToolbox::isfinite(arg.derivative(i)))
00465 return false;
00466
00467 return true;
00468 }
00469
00470 static bool isnan(const Evaluation& arg)
00471 {
00472 if (InnerToolbox::isnan(arg.value()))
00473 return true;
00474
00475 for (int i = 0; i < numVars; ++i)
00476 if (InnerToolbox::isnan(arg.derivative(i)))
00477 return true;
00478
00479 return false;
00480 }
00481 };
00482
00483 }
00484
00485 #endif