709 lines
27 KiB
C++
709 lines
27 KiB
C++
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
|
|
// vi: set et ts=4 sw=4 sts=4:
|
|
/*
|
|
This file is part of the Open Porous Media project (OPM).
|
|
|
|
OPM is free software: you can redistribute it and/or modify
|
|
it under the terms of the GNU General Public License as published by
|
|
the Free Software Foundation, either version 2 of the License, or
|
|
(at your option) any later version.
|
|
|
|
OPM is distributed in the hope that it will be useful,
|
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
GNU General Public License for more details.
|
|
|
|
You should have received a copy of the GNU General Public License
|
|
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
|
|
|
Consult the COPYING file in the top-level source directory of this
|
|
module for the precise wording of the license and the list of
|
|
copyright holders.
|
|
*/
|
|
/*!
|
|
* \file
|
|
*
|
|
* \brief Low-level tests for the localized automatic differentiation (AD) framework.
|
|
*/
|
|
#include "config.h"
|
|
|
|
// for testing the "!=" and "==" operators, we need to disable the -Wfloat-equal to
|
|
// prevent clang from producing a warning with -Weverything
|
|
#if defined(__GNUC__) || defined(__clang__)
|
|
#pragma GCC diagnostic push
|
|
#pragma GCC diagnostic ignored "-Wfloat-equal"
|
|
#endif
|
|
|
|
#include <opm/material/densead/Evaluation.hpp>
|
|
#include <opm/material/densead/Math.hpp>
|
|
|
|
#include <iostream>
|
|
#include <array>
|
|
#include <cmath>
|
|
#include <algorithm>
|
|
#include <cassert>
|
|
#include <stdexcept>
|
|
|
|
template <class Eval, int numVars, int staticSize, class Scalar, class Implementation>
|
|
struct TestEnvBase
|
|
{
|
|
const Implementation& asImp_() const
|
|
{ return *static_cast<const Implementation*>(this); }
|
|
|
|
void testOperators(const Scalar tolerance)
|
|
{
|
|
// test the constructors of the Opm::DenseAd::Evaluation class
|
|
const Scalar c = 1.234;
|
|
const Scalar x = 4.567;
|
|
const Scalar y = 8.910;
|
|
[[maybe_unused]] Scalar z = Opm::blank(tolerance);
|
|
const Eval cEval = asImp_().createConstant(c);
|
|
[[maybe_unused]] const Eval c2Eval = asImp_().createConstant(c);
|
|
const Eval xEval = asImp_().createVariable(x, 0);
|
|
const Eval yEval = asImp_().createVariable(y, 1);
|
|
[[maybe_unused]] const Eval zEval = Opm::blank(xEval);
|
|
|
|
Eval xyEval = xEval;
|
|
Eval yxEval = yEval;
|
|
|
|
xyEval.copyDerivatives(yxEval);
|
|
yxEval.clearDerivatives();
|
|
|
|
for (int i = 0; i < xyEval.size(); ++i) {
|
|
if (i == 0 && xEval.derivative(i) != 1.0)
|
|
throw std::logic_error("oops: createVariable");
|
|
else if (i == 1 && yEval.derivative(i) != 1.0)
|
|
throw std::logic_error("oops: createVariable");
|
|
else if (i == 1 && xyEval.derivative(i) != 1.0)
|
|
throw std::logic_error("oops: copyDerivatives");
|
|
else continue;
|
|
|
|
// the remaining derivatives must be zero
|
|
if (xEval.derivative(i) != 0.0)
|
|
throw std::logic_error("oops: createVariable");
|
|
if (yEval.derivative(i) != 0.0)
|
|
throw std::logic_error("oops: createVariable");
|
|
if (cEval.derivative(i) != 0.0)
|
|
throw std::logic_error("oops: createConstant");
|
|
if (xyEval.derivative(i) != 0.0)
|
|
throw std::logic_error("oops: copyDerivatives");
|
|
if (xyEval.derivative(i) != 0.0)
|
|
throw std::logic_error("oops: clearDerivatives");
|
|
}
|
|
|
|
// test the non-inplace operators
|
|
{
|
|
Eval a = xEval + yEval;
|
|
if (std::abs(a.value() - (x + y)) > tolerance)
|
|
throw std::logic_error("oops: operator+");
|
|
|
|
Eval b = xEval + c;
|
|
if (std::abs(b.value() - (x + c)) > tolerance)
|
|
throw std::logic_error("oops: operator+");
|
|
|
|
Eval d = xEval + cEval;
|
|
if (std::abs(d.value() - (x + c)) > tolerance)
|
|
throw std::logic_error("oops: operator+");
|
|
}
|
|
|
|
{
|
|
Eval a = xEval - yEval;
|
|
if (std::abs(a.value() - (x - y)) > tolerance)
|
|
throw std::logic_error("oops: operator-");
|
|
|
|
Eval b = xEval - c;
|
|
if (std::abs(b.value() - (x - c)) > tolerance)
|
|
throw std::logic_error("oops: operator-");
|
|
|
|
Eval d = xEval - cEval;
|
|
if (std::abs(d.value() - (x - c)) > tolerance)
|
|
throw std::logic_error("oops: operator-");
|
|
}
|
|
|
|
{
|
|
Eval a = xEval*yEval;
|
|
if (std::abs(a.value() - (x*y)) > tolerance)
|
|
throw std::logic_error("oops: operator*");
|
|
|
|
Eval b = xEval*c;
|
|
if (std::abs(b.value() - (x*c)) > tolerance)
|
|
throw std::logic_error("oops: operator*");
|
|
|
|
Eval d = xEval*cEval;
|
|
if (std::abs(d.value() - (x*c)) > tolerance)
|
|
throw std::logic_error("oops: operator*");
|
|
}
|
|
|
|
{
|
|
Eval a = xEval/yEval;
|
|
if (std::abs(a.value() - (x/y)) > tolerance)
|
|
throw std::logic_error("oops: operator/");
|
|
|
|
Eval b = xEval/c;
|
|
if (std::abs(b.value() - (x/c)) > tolerance)
|
|
throw std::logic_error("oops: operator/");
|
|
|
|
Eval d = xEval/cEval;
|
|
if (std::abs(d.value() - (x/c)) > tolerance)
|
|
throw std::logic_error("oops: operator/");
|
|
}
|
|
|
|
// test the inplace operators
|
|
{
|
|
Eval a = xEval;
|
|
a += yEval;
|
|
if (std::abs(a.value() - (x + y)) > tolerance)
|
|
throw std::logic_error("oops: operator+");
|
|
|
|
Eval b = xEval;
|
|
b += c;
|
|
if (std::abs(b.value() - (x + c)) > tolerance)
|
|
throw std::logic_error("oops: operator+");
|
|
|
|
Eval d = xEval;
|
|
d += cEval;
|
|
if (std::abs(d.value() - (x + c)) > tolerance)
|
|
throw std::logic_error("oops: operator+");
|
|
}
|
|
|
|
{
|
|
Eval a = xEval;
|
|
a -= yEval;
|
|
if (std::abs(a.value() - (x - y)) > tolerance)
|
|
throw std::logic_error("oops: operator-");
|
|
|
|
Eval b = xEval;
|
|
b -= c;
|
|
if (std::abs(b.value() - (x - c)) > tolerance)
|
|
throw std::logic_error("oops: operator-");
|
|
|
|
Eval d = xEval;
|
|
d -= cEval;
|
|
if (std::abs(d.value() - (x - c)) > tolerance)
|
|
throw std::logic_error("oops: operator-");
|
|
}
|
|
|
|
{
|
|
Eval a = xEval;
|
|
a *= yEval;
|
|
if (std::abs(a.value() - (x*y)) > tolerance)
|
|
throw std::logic_error("oops: operator*");
|
|
|
|
Eval b = xEval;
|
|
b *= c;
|
|
if (std::abs(b.value() - (x*c)) > tolerance)
|
|
throw std::logic_error("oops: operator*");
|
|
|
|
Eval d = xEval;
|
|
d *= cEval;
|
|
if (std::abs(d.value() - (x*c)) > tolerance)
|
|
throw std::logic_error("oops: operator*");
|
|
}
|
|
|
|
{
|
|
Eval a = xEval;
|
|
a /= yEval;
|
|
if (std::abs(a.value() - (x/y)) > tolerance)
|
|
throw std::logic_error("oops: operator/");
|
|
|
|
Eval b = xEval;
|
|
b /= c;
|
|
if (std::abs(b.value() - (x/c)) > tolerance)
|
|
throw std::logic_error("oops: operator/");
|
|
|
|
Eval d = xEval;
|
|
d /= cEval;
|
|
if (std::abs(d.value() - (x/c)) > tolerance)
|
|
throw std::logic_error("oops: operator/");
|
|
}
|
|
|
|
{
|
|
Eval a = asImp_().createConstant(1.0);
|
|
Eval b = asImp_().createConstant(2.0);
|
|
if (a >= b)
|
|
throw std::logic_error("oops: operator>=");
|
|
if (a >= 2.0)
|
|
throw std::logic_error("oops: operator>=");
|
|
|
|
if (!(a >= a))
|
|
throw std::logic_error("oops: operator>=");
|
|
if (!(a >= 1.0))
|
|
throw std::logic_error("oops: operator>=");
|
|
if (!(1.0 <= a))
|
|
throw std::logic_error("oops: operator<=");
|
|
|
|
if (b <= a)
|
|
throw std::logic_error("oops: operator<=");
|
|
if (b <= 1.0)
|
|
throw std::logic_error("oops: operator<=");
|
|
|
|
if (!(b <= b))
|
|
throw std::logic_error("oops: operator<=");
|
|
if (!(b <= 2.0))
|
|
throw std::logic_error("oops: operator<=");
|
|
if (!(2.0 >= b))
|
|
throw std::logic_error("oops: operator>=");
|
|
|
|
if (a != a)
|
|
throw std::logic_error("oops: operator!=");
|
|
if (a != 1.0)
|
|
throw std::logic_error("oops: operator!=");
|
|
if (1.0 != a)
|
|
throw std::logic_error("oops: operator!=");
|
|
}
|
|
}
|
|
|
|
template <class AdFn, class ClassicFn>
|
|
void test1DFunction(AdFn* adFn, ClassicFn* classicFn, Scalar xMin = 1e-6, Scalar xMax = 1000)
|
|
{
|
|
int n = 100*1000;
|
|
for (int i = 0; i < n; ++ i) {
|
|
Scalar x = Scalar(i)/(n - 1)*(xMax - xMin) + xMin;
|
|
|
|
const auto& xEval = asImp_().createVariable(x, 0);
|
|
const Eval& yEval = adFn(xEval);
|
|
|
|
Scalar eps = std::numeric_limits<Scalar>::epsilon()*1e7;
|
|
eps = std::max(eps, eps*std::abs(x));
|
|
|
|
Scalar y = classicFn(x);
|
|
Scalar yStar1 = classicFn(x - eps);
|
|
Scalar yStar2 = classicFn(x + eps);
|
|
Scalar yPrime = (yStar2 - yStar1)/(2*eps);
|
|
|
|
if (std::abs(y-yEval.value()) > std::numeric_limits<Scalar>::epsilon()*1e2)
|
|
throw std::logic_error("oops: value");
|
|
|
|
Scalar deltaAbs = std::abs(yPrime - yEval.derivative(0));
|
|
Scalar deltaRel = std::abs(deltaAbs/yPrime);
|
|
if (deltaAbs > 1000*eps && deltaRel > 1000*eps)
|
|
throw std::logic_error("oops: derivative @"+std::to_string((long double) x)+": "
|
|
+ std::to_string((long double) yPrime) + " vs "
|
|
+ std::to_string((long double) yEval.derivative(0))
|
|
+ " delta: " + std::to_string((long double) std::abs(yPrime - yEval.derivative(0))));
|
|
}
|
|
}
|
|
|
|
template <class AdFn,
|
|
class ClassicFn>
|
|
void test2DFunction1(AdFn* adFn, ClassicFn* classicFn, Scalar xMin, Scalar xMax, Scalar y)
|
|
{
|
|
int n = 100*1000;
|
|
for (int i = 0; i < n; ++ i) {
|
|
Scalar x = Scalar(i)/(n - 1)*(xMax - xMin) + xMin;
|
|
|
|
const auto& xEval = asImp_().createVariable(x, 0);
|
|
const auto& yEval = asImp_().createConstant(y);
|
|
const Eval& zEval = adFn(xEval, yEval);
|
|
|
|
Scalar eps = std::numeric_limits<Scalar>::epsilon()*1e7;
|
|
eps = std::max(eps, eps*std::abs(x));
|
|
|
|
Scalar z = classicFn(x, y);
|
|
Scalar zStar1 = classicFn(x - eps, y);
|
|
Scalar zStar2 = classicFn(x + eps, y);
|
|
Scalar zPrime = (zStar2 - zStar1)/(2.*eps);
|
|
|
|
if (std::abs(z - zEval.value())/std::abs(z + zEval.value())
|
|
> std::numeric_limits<Scalar>::epsilon()*1e2)
|
|
throw std::logic_error("oops: value");
|
|
|
|
Scalar deltaAbs = std::abs(zPrime - zEval.derivative(0));
|
|
Scalar deltaRel = std::abs(deltaAbs/zPrime);
|
|
if (deltaAbs > 1000*eps && deltaRel > 1000*eps)
|
|
throw std::logic_error("oops: derivative @"+std::to_string((long double) x)+": "
|
|
+ std::to_string((long double) zPrime) + " vs "
|
|
+ std::to_string((long double) zEval.derivative(0))
|
|
+ " delta: " + std::to_string((long double) std::abs(zPrime - zEval.derivative(0))));
|
|
}
|
|
}
|
|
|
|
template <class AdFn,
|
|
class ClassicFn>
|
|
void test2DFunction2(AdFn* adFn, ClassicFn* classicFn, Scalar x, Scalar yMin, Scalar yMax)
|
|
{
|
|
int n = 100*1000;
|
|
for (int i = 0; i < n; ++ i) {
|
|
Scalar y = Scalar(i)/(n - 1)*(yMax - yMin) + yMin;
|
|
|
|
const auto& xEval = asImp_().createConstant(x);
|
|
const auto& yEval = asImp_().createVariable(y, 1);
|
|
const Eval& zEval = adFn(xEval, yEval);
|
|
|
|
Scalar eps = std::numeric_limits<Scalar>::epsilon()*1e7;
|
|
eps = std::max(eps, eps*std::abs(y));
|
|
|
|
Scalar z = classicFn(x, y);
|
|
Scalar zStar1 = classicFn(x, y - eps);
|
|
Scalar zStar2 = classicFn(x, y + eps);
|
|
Scalar zPrime = (zStar2 - zStar1)/(2*eps);
|
|
|
|
if (std::abs(z - zEval.value())/std::abs(z + zEval.value())
|
|
> std::numeric_limits<Scalar>::epsilon()*1e2)
|
|
throw std::logic_error("oops: value");
|
|
|
|
Scalar deltaAbs = std::abs(zPrime - zEval.derivative(1));
|
|
Scalar deltaRel = std::abs(deltaAbs/zPrime);
|
|
if (deltaAbs > 1000*eps && deltaRel > 1000*eps)
|
|
throw std::logic_error("oops: derivative @"+std::to_string((long double) x)+": "
|
|
+ std::to_string((long double) zPrime) + " vs "
|
|
+ std::to_string((long double) zEval.derivative(1))
|
|
+ " delta: " + std::to_string((long double) std::abs(zPrime - zEval.derivative(1))));
|
|
}
|
|
}
|
|
|
|
void testPowBase(Scalar baseMin = 1e-2, Scalar baseMax = 100)
|
|
{
|
|
typedef Opm::MathToolbox<Eval> EvalToolbox;
|
|
|
|
Scalar exp = 1.234;
|
|
const auto& expEval = asImp_().createConstant(exp);
|
|
|
|
int n = 100*1000;
|
|
for (int i = 0; i < n; ++ i) {
|
|
Scalar base = Scalar(i)/(n - 1)*(baseMax - baseMin) + baseMin;
|
|
|
|
const auto& baseEval = asImp_().createVariable(base, 0);
|
|
const Eval& zEval1 = pow(baseEval, exp);
|
|
const Eval& zEval2 = pow(baseEval, expEval);
|
|
|
|
Scalar eps = std::numeric_limits<Scalar>::epsilon()*1e7;
|
|
eps = std::max(eps, eps*std::abs(base));
|
|
|
|
Scalar z = pow(base, exp);
|
|
Scalar zStar1 = pow(base - eps, exp);
|
|
Scalar zStar2 = pow(base + eps, exp);
|
|
Scalar zPrime = (zStar2 - zStar1)/(2*eps);
|
|
|
|
if (std::abs(z - zEval2.value())/std::abs(z + zEval2.value())
|
|
> std::numeric_limits<Scalar>::epsilon()*1e2)
|
|
throw std::logic_error("oops: value");
|
|
|
|
Scalar deltaAbs = std::abs(zPrime - zEval1.derivative(0));
|
|
Scalar deltaRel = std::abs(deltaAbs/zPrime);
|
|
if (deltaAbs > 1000*eps && deltaRel > 1000*eps)
|
|
throw std::logic_error("oops: derivative @"+std::to_string((long double) base)+": "
|
|
+ std::to_string((long double) zPrime) + " vs "
|
|
+ std::to_string((long double) zEval1.derivative(0))
|
|
+ " delta: " + std::to_string((long double) std::abs(zPrime - zEval1.derivative(0))));
|
|
|
|
if (!EvalToolbox::isSame(zEval1, zEval2, /*tolerance=*/std::numeric_limits<Scalar>::epsilon()*1e3*zEval1.value()))
|
|
throw std::logic_error("oops: pow(Eval, Scalar) != pow(Eval, Eval)");
|
|
}
|
|
}
|
|
|
|
void testPowExp(Scalar expMin = -100, Scalar expMax = 100)
|
|
{
|
|
typedef Opm::MathToolbox<Eval> EvalToolbox;
|
|
|
|
Scalar base = 1.234;
|
|
const auto& baseEval = asImp_().createConstant(base);
|
|
|
|
int n = 100*1000;
|
|
for (int i = 0; i < n; ++ i) {
|
|
Scalar exp = Scalar(i)/(n - 1)*(expMax - expMin) + expMin;
|
|
const auto& expEval = asImp_().createVariable(exp, 1);
|
|
|
|
const Eval& zEval1 = pow(base, expEval);
|
|
const Eval& zEval2 = pow(baseEval, expEval);
|
|
|
|
Scalar eps = std::numeric_limits<Scalar>::epsilon()*1e7;
|
|
eps = std::max(eps, eps*std::abs(exp));
|
|
|
|
Scalar z = pow(base, exp);
|
|
Scalar zStar1 = pow(base, exp - eps);
|
|
Scalar zStar2 = pow(base, exp + eps);
|
|
Scalar zPrime = (zStar2 - zStar1)/(2*eps);
|
|
|
|
if (std::abs(z - zEval2.value())/std::abs(z + zEval2.value())
|
|
> std::numeric_limits<Scalar>::epsilon()*1e2)
|
|
throw std::logic_error("oops: value");
|
|
|
|
Scalar deltaAbs = std::abs(zPrime - zEval1.derivative(1));
|
|
Scalar deltaRel = std::abs(deltaAbs/zPrime);
|
|
if (deltaAbs > 1000*eps && deltaRel > 1000*eps)
|
|
throw std::logic_error("oops: derivative @"+std::to_string((long double) base)+": "
|
|
+ std::to_string((long double) zPrime) + " vs "
|
|
+ std::to_string((long double) zEval1.derivative(1))
|
|
+ " delta: " + std::to_string((long double) std::abs(zPrime - zEval1.derivative(1))));
|
|
|
|
if (!EvalToolbox::isSame(zEval1, zEval2, /*tolerance=*/std::numeric_limits<Scalar>::epsilon()*1e3*zEval1.value()))
|
|
throw std::logic_error("oops: pow(Eval, Scalar) != pow(Eval, Eval)");
|
|
}
|
|
}
|
|
|
|
void testAtan2()
|
|
{
|
|
int n = 1000;
|
|
Scalar maxVal = 10.0;
|
|
for (int i = 1; i < n; ++ i) {
|
|
Scalar x = 2*maxVal*Scalar(i)/n - maxVal;
|
|
if (- 0.05 < x && x < 0.05)
|
|
// avoid numerical problems
|
|
continue;
|
|
|
|
const Eval& xEval = asImp_().createVariable(x, 0);
|
|
|
|
for (int j = 1; j < n; ++ j) {
|
|
Scalar y = 2*maxVal*Scalar(j)/n - maxVal;
|
|
|
|
if (- 0.05 < y && y < 0.05)
|
|
// avoid numerical problems
|
|
continue;
|
|
|
|
const Eval& yEval = asImp_().createVariable(y, 0);
|
|
const Eval& zEval = atan2(xEval, yEval);
|
|
|
|
Scalar eps = std::numeric_limits<Scalar>::epsilon()*1e7;
|
|
eps = std::max(eps, eps*std::abs(x));
|
|
|
|
Scalar z = atan2(x, y);
|
|
Scalar zStar1 = atan2(x - eps, y - eps);
|
|
Scalar zStar2 = atan2(x + eps, y + eps);
|
|
Scalar zPrime = (zStar2 - zStar1)/(2*eps);
|
|
|
|
if (std::abs(z - zEval.value())/std::abs(z + zEval.value())
|
|
> std::numeric_limits<Scalar>::epsilon()*1e2)
|
|
throw std::logic_error("oops: value");
|
|
|
|
Scalar deltaAbs = std::abs(zPrime - zEval.derivative(0));
|
|
Scalar deltaRel = std::abs(deltaAbs/zPrime);
|
|
if (deltaAbs > 1000*eps && deltaRel > 1000*eps)
|
|
throw std::logic_error("oops: derivative @("+std::to_string((long double) x)+","+std::to_string((long double) y)+"): "
|
|
+ std::to_string((long double) zPrime) + " vs "
|
|
+ std::to_string((long double) zEval.derivative(0))
|
|
+ " delta: " + std::to_string((long double) std::abs(zPrime - zEval.derivative(0))));
|
|
|
|
}
|
|
}
|
|
}
|
|
|
|
// prototypes
|
|
static double myScalarMin(double a, double b)
|
|
{ return std::min(a, b); }
|
|
|
|
static double myScalarMax(double a, double b)
|
|
{ return std::max(a, b); }
|
|
|
|
inline void testAll()
|
|
{
|
|
std::cout << " Testing operators and constructors\n";
|
|
const Scalar eps = std::numeric_limits<Scalar>::epsilon()*1e3;
|
|
testOperators(eps);
|
|
|
|
std::cout << " Testing min()\n";
|
|
test2DFunction1(Opm::DenseAd::min<Scalar, numVars, staticSize>,
|
|
myScalarMin,
|
|
-1000, 1000,
|
|
/*p=*/1.234);
|
|
|
|
test2DFunction2(Opm::DenseAd::min<Scalar, numVars, staticSize>,
|
|
myScalarMin,
|
|
/*T=*/1.234,
|
|
-1000, 1000);
|
|
|
|
std::cout << " Testing max()\n";
|
|
test2DFunction1(Opm::DenseAd::max<Scalar, numVars, staticSize>,
|
|
myScalarMax,
|
|
-1000, 1000,
|
|
/*p=*/1.234);
|
|
|
|
test2DFunction2(Opm::DenseAd::max<Scalar, numVars, staticSize>,
|
|
myScalarMax,
|
|
/*T=*/1.234,
|
|
-1000, 1000);
|
|
|
|
std::cout << " Testing pow()\n";
|
|
testPowBase();
|
|
testPowExp();
|
|
|
|
std::cout << " Testing abs()\n";
|
|
test1DFunction(Opm::DenseAd::abs<Scalar, numVars, staticSize>,
|
|
static_cast<Scalar (*)(Scalar)>(std::abs));
|
|
|
|
std::cout << " Testing sqrt()\n";
|
|
test1DFunction(Opm::DenseAd::sqrt<Scalar, numVars, staticSize>,
|
|
static_cast<Scalar (*)(Scalar)>(std::sqrt));
|
|
|
|
std::cout << " Testing sin()\n";
|
|
test1DFunction(Opm::DenseAd::sin<Scalar, numVars, staticSize>,
|
|
static_cast<Scalar (*)(Scalar)>(std::sin),
|
|
0, 2*M_PI);
|
|
|
|
std::cout << " Testing asin()\n";
|
|
test1DFunction(Opm::DenseAd::asin<Scalar, numVars, staticSize>,
|
|
static_cast<Scalar (*)(Scalar)>(std::asin),
|
|
-1.0, 1.0);
|
|
|
|
std::cout << " Testing cos()\n";
|
|
test1DFunction(Opm::DenseAd::cos<Scalar, numVars, staticSize>,
|
|
static_cast<Scalar (*)(Scalar)>(std::cos),
|
|
0, 2*M_PI);
|
|
|
|
std::cout << " Testing acos()\n";
|
|
test1DFunction(Opm::DenseAd::acos<Scalar, numVars, staticSize>,
|
|
static_cast<Scalar (*)(Scalar)>(std::acos),
|
|
-1.0, 1.0);
|
|
|
|
std::cout << " Testing tan()\n";
|
|
test1DFunction(Opm::DenseAd::tan<Scalar, numVars, staticSize>,
|
|
static_cast<Scalar (*)(Scalar)>(std::tan),
|
|
-M_PI / 2 * 0.95, M_PI / 2 * 0.95);
|
|
|
|
std::cout << " Testing atan()\n";
|
|
test1DFunction(Opm::DenseAd::atan<Scalar, numVars, staticSize>,
|
|
static_cast<Scalar (*)(Scalar)>(std::atan),
|
|
-10*1000.0, 10*1000.0);
|
|
|
|
std::cout << " Testing atan2()\n";
|
|
testAtan2();
|
|
|
|
std::cout << " Testing exp()\n";
|
|
test1DFunction(Opm::DenseAd::exp<Scalar, numVars, staticSize>,
|
|
static_cast<Scalar (*)(Scalar)>(std::exp),
|
|
-100, 100);
|
|
|
|
std::cout << " Testing log()\n";
|
|
test1DFunction(Opm::DenseAd::log<Scalar, numVars, staticSize>,
|
|
static_cast<Scalar (*)(Scalar)>(std::log),
|
|
1e-6, 1e9);
|
|
|
|
std::cout << " Testing log10()\n";
|
|
test1DFunction(Opm::DenseAd::log10<Scalar, numVars, staticSize>,
|
|
static_cast<Scalar (*)(Scalar)>(std::log10),
|
|
1e-6, 1e9);
|
|
|
|
while (false) {
|
|
[[maybe_unused]] Scalar val1 = 0.0;
|
|
[[maybe_unused]] Scalar val2 = 1.0;
|
|
[[maybe_unused]] Scalar resultVal;
|
|
[[maybe_unused]] Eval eval1 = asImp_().createConstant(1.0);
|
|
[[maybe_unused]] Eval eval2 = asImp_().createConstant(2.0);
|
|
[[maybe_unused]] Eval resultEval = asImp_().newBlankEval();;
|
|
|
|
// make sure that the convenince functions work (i.e., that everything can be
|
|
// accessed without the MathToolbox<Scalar> detour.)
|
|
resultVal = Opm::constant<Scalar>(val1);
|
|
resultVal = Opm::variable<Scalar>(val1, /*idx=*/0);
|
|
resultVal = Opm::decay<Scalar>(val1);
|
|
resultVal = Opm::scalarValue(val1);
|
|
resultVal = Opm::getValue(val1);
|
|
resultVal = Opm::min(val1, val2);
|
|
resultVal = Opm::max(val1, val2);
|
|
resultVal = Opm::atan2(val1, val2);
|
|
resultVal = Opm::pow(val1, val2);
|
|
resultVal = Opm::abs(val1);
|
|
resultVal = Opm::atan(val1);
|
|
resultVal = Opm::sin(val1);
|
|
resultVal = Opm::asin(val1);
|
|
resultVal = Opm::cos(val1);
|
|
resultVal = Opm::acos(val1);
|
|
resultVal = Opm::sqrt(val1);
|
|
resultVal = Opm::exp(val1);
|
|
resultVal = Opm::log(val1);
|
|
|
|
resultEval = Opm::decay<Eval>(eval1);
|
|
resultVal = Opm::decay<Scalar>(eval1);
|
|
resultVal = Opm::scalarValue(eval1);
|
|
resultVal = Opm::getValue(eval1);
|
|
resultEval = Opm::min(eval1, eval2);
|
|
resultEval = Opm::min(eval1, val2);
|
|
resultEval = Opm::max(eval1, eval2);
|
|
resultEval = Opm::max(eval1, val2);
|
|
resultEval = Opm::atan2(eval1, eval2);
|
|
resultEval = Opm::atan2(eval1, val2);
|
|
resultEval = Opm::pow(eval1, eval2);
|
|
resultEval = Opm::pow(eval1, val2);
|
|
resultEval = Opm::abs(eval1);
|
|
resultEval = Opm::atan(eval1);
|
|
resultEval = Opm::sin(eval1);
|
|
resultEval = Opm::asin(eval1);
|
|
resultEval = Opm::cos(eval1);
|
|
resultEval = Opm::acos(eval1);
|
|
resultEval = Opm::sqrt(eval1);
|
|
resultEval = Opm::exp(eval1);
|
|
resultEval = Opm::log(eval1);
|
|
resultEval = Opm::log10(eval1);
|
|
|
|
}
|
|
}
|
|
|
|
};//StaticTestEnv
|
|
|
|
template <class Scalar, int staticSize>
|
|
struct DynamicTestEnv : public TestEnvBase<Opm::DenseAd::DynamicEvaluation<Scalar, staticSize>,
|
|
-1,
|
|
staticSize,
|
|
Scalar,
|
|
DynamicTestEnv<Scalar, staticSize> >
|
|
{
|
|
typedef Opm::DenseAd::DynamicEvaluation<Scalar, staticSize> Eval;
|
|
DynamicTestEnv(int numDerivs)
|
|
: numDerivs_(numDerivs)
|
|
{}
|
|
|
|
Eval newBlankEval() const
|
|
{ return Eval(); }
|
|
|
|
Eval createConstant(Scalar c) const
|
|
{ return Opm::constant<Scalar, staticSize>(numDerivs_, c); }
|
|
|
|
Eval createVariable(Scalar v, int varIdx) const
|
|
{ return Opm::variable<Scalar, staticSize>(numDerivs_, v, varIdx); }
|
|
|
|
private:
|
|
int numDerivs_;
|
|
};
|
|
|
|
template <class Scalar, int numDerivs, int staticSize = 0>
|
|
struct StaticTestEnv : public TestEnvBase<Opm::DenseAd::Evaluation<Scalar, numDerivs>,
|
|
numDerivs,
|
|
staticSize,
|
|
Scalar,
|
|
StaticTestEnv<Scalar, numDerivs> >
|
|
{
|
|
typedef Opm::DenseAd::Evaluation<Scalar, numDerivs> Eval;
|
|
|
|
StaticTestEnv()
|
|
{}
|
|
|
|
Eval newBlankEval() const
|
|
{ return Eval(); }
|
|
|
|
Eval createConstant(Scalar c) const
|
|
{ return Opm::constant<Eval, Scalar>(c); }
|
|
|
|
Eval createVariable(Scalar v, int varIdx) const
|
|
{ return Opm::variable<Eval, Scalar>(v, varIdx); }
|
|
|
|
private:
|
|
int numDerivs_;
|
|
};
|
|
|
|
int main()
|
|
{
|
|
std::cout << "Testing statically sized evaluations\n";
|
|
std::cout << " -> Scalar == double, n = 15\n";
|
|
StaticTestEnv<double, 15>().testAll();
|
|
std::cout << " -> Scalar == double, n = s\n";
|
|
StaticTestEnv<double, 2>().testAll();
|
|
std::cout << " -> Scalar == float, n = 15\n";
|
|
StaticTestEnv<float, 15>().testAll();
|
|
std::cout << " -> Scalar == float, n = 2\n";
|
|
StaticTestEnv<float, 2>().testAll();
|
|
|
|
std::cout << "Testing dynamically sized evaluations\n";
|
|
std::cout << " -> Scalar == double\n";
|
|
DynamicTestEnv<double, 6>(5).testAll();
|
|
DynamicTestEnv<double, 0>(5).testAll();
|
|
DynamicTestEnv<double, 4>(8).testAll();
|
|
DynamicTestEnv<double, 4>(2).testAll();
|
|
std::cout << " -> Scalar == float\n";
|
|
DynamicTestEnv<float, 6>(5).testAll();
|
|
DynamicTestEnv<float, 0>(5).testAll();
|
|
DynamicTestEnv<float, 4>(8).testAll();
|
|
DynamicTestEnv<float, 4>(2).testAll();
|
|
|
|
return 0;
|
|
}
|