add a "local" AD framework

The basic idea is to replace "plain" scalar values by "function
evaluations": these store the function's value plus a set of
derivatives for it and the chain rule is used to "drag" the
derivatives along.

So far, the framework only implements the "forward" approach to
automatic differentiation and expression templates are not supported
yet. (The latter point may change in the furture, though.)

"local" means that the framework uses static arrays to represent the
derivatives, i.e. the number of derivatives which are required must be
specified at compile time. Compared to dynamic arrays this improves
efficiency considerably if function evaluation objects must be
allocated and deallocated frequently.

opm/material/localad/Evaluation.hpp

@@ -0,0 +1,467 @@
+/*
+  Copyright (C) 2015 by Andreas Lauser
+
+  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/>.
+*/
+/*!
+ * \file
+ *
+ * \brief Representation of an evaluation of a function and its derivatives w.r.t. a set
+ *        of variables in the localized OPM automatic differentiation (AD) framework.
+ */
+#ifndef OPM_LOCAL_AD_EVALUATION_HPP
+#define OPM_LOCAL_AD_EVALUATION_HPP
+
+#include <iostream>
+#include <array>
+#include <cassert>
+#include <opm/material/common/Valgrind.hpp>
+
+namespace Opm {
+namespace LocalAd {
+/*!
+ * \brief Represents a function evaluation and its derivatives w.r.t. a fixed set of
+ *        variables.
+ */
+template <class ScalarT, class VarSetTag, int numVars>
+class Evaluation
+{
+public:
+    typedef ScalarT Scalar;
+
+    enum { size = numVars };
+
+    Evaluation()
+    {};
+
+    // copy other function evaluation
+    Evaluation(const Evaluation& other)
+    {
+        // copy evaluated function value and derivatives
+        value = other.value;
+        std::copy(other.derivatives.begin(), other.derivatives.end(), derivatives.begin());
+    };
+
+    // create an evaluation which represents a constant function
+    //
+    // i.e., f(x) = c. this implies an evaluation with the given value and all
+    // derivatives being zero.
+    Evaluation(Scalar c)
+    {
+        value = c;
+        std::fill(derivatives.begin(), derivatives.end(), 0.0);
+    };
+
+    // create a function evaluation for a "naked" depending variable (i.e., f(x) = x)
+    static Evaluation createVariable(Scalar value, int varPos)
+    {
+        // The variable position must be in represented by the given variable descriptor
+        assert(0 <= varPos && varPos < size);
+
+        Evaluation result;
+
+        // copy function value and set all derivatives to 0, except for the variable
+        // which is represented by the value (which is set to 1.0)
+        result.value = value;
+        std::fill(result.derivatives.begin(), result.derivatives.end(), 0.0);
+        result.derivatives[varPos] = 1.0;
+
+        return result;
+    }
+
+    // "evaluate" a constant function (i.e. a function that does not depend on the set of
+    // relevant variables, f(x) = c).
+    static Evaluation createConstant(Scalar value)
+    {
+        Evaluation result;
+        result.value = value;
+        std::fill(result.derivatives.begin(), result.derivatives.end(), 0.0);
+        Valgrind::CheckDefined(result.value);
+        Valgrind::CheckDefined(result.derivatives);
+        return result;
+    }
+
+    // print the value and the derivatives of the function evaluation
+    void print(std::ostream& os = std::cout) const
+    {
+        os << "v: " << value << " / d:";
+        for (int varIdx = 0; varIdx < derivatives.size(); ++varIdx)
+            os << " " << derivatives[varIdx];
+    }
+
+
+    Evaluation& operator+=(const Evaluation& other)
+    {
+        // value and derivatives are added
+        this->value += other.value;
+        for (int varIdx= 0; varIdx < size; ++varIdx)
+            this->derivatives[varIdx] += other.derivatives[varIdx];
+
+        return *this;
+    }
+
+    Evaluation& operator+=(Scalar other)
+    {
+        // value is added, derivatives stay the same
+        this->value += other;
+
+        return *this;
+    }
+
+    Evaluation& operator-=(const Evaluation& other)
+    {
+        // value and derivatives are subtracted
+        this->value -= other.value;
+        for (int varIdx= 0; varIdx < size; ++varIdx)
+            this->derivatives[varIdx] -= other.derivatives[varIdx];
+
+        return *this;
+    }
+
+    Evaluation& operator-=(Scalar other)
+    {
+        // for constants, values are subtracted, derivatives stay the same
+        this->value -= other;
+
+        return *this;
+    }
+
+    Evaluation& operator*=(const Evaluation& other)
+    {
+        // while the values are multiplied, the derivatives follow the product rule,
+        // i.e., (u*v)' = (v'u + u'v).
+        Scalar u = this->value;
+        Scalar v = other.value;
+        this->value *= v;
+        for (int varIdx= 0; varIdx < size; ++varIdx) {
+            Scalar uPrime = this->derivatives[varIdx];
+            Scalar vPrime = other.derivatives[varIdx];
+
+            this->derivatives[varIdx] = (v*uPrime + u*vPrime);
+        }
+
+        return *this;
+    }
+
+    Evaluation& operator*=(Scalar other)
+    {
+        // values and derivatives are multiplied
+        this->value *= other;
+        for (int varIdx= 0; varIdx < size; ++varIdx)
+            this->derivatives[varIdx] *= other;
+
+        return *this;
+    }
+
+    Evaluation& operator/=(const Evaluation& other)
+    {
+        // values are divided, derivatives follow the rule for division, i.e., (u/v)' = (v'u -
+        // u'v)/v^2.
+        Scalar u = this->value;
+        Scalar v = other.value;
+        this->value /= v;
+        for (int varIdx= 0; varIdx < size; ++varIdx) {
+            Scalar uPrime = this->derivatives[varIdx];
+            Scalar vPrime = other.derivatives[varIdx];
+
+            this->derivatives[varIdx] = (v*uPrime - u*vPrime)/(v*v);
+        }
+
+        return *this;
+    }
+
+    Evaluation& operator/=(Scalar other)
+    {
+        // values and derivatives are divided
+        other = 1.0/other;
+        this->value *= other;
+        for (int varIdx= 0; varIdx < size; ++varIdx)
+            this->derivatives[varIdx] *= other;
+
+        return *this;
+    }
+
+    Evaluation operator+(const Evaluation& other) const
+    {
+        Evaluation result(*this);
+        result += other;
+        return result;
+    }
+
+    Evaluation operator+(Scalar other) const
+    {
+        Evaluation result(*this);
+        result += other;
+        return result;
+    }
+
+    Evaluation operator-(const Evaluation& other) const
+    {
+        Evaluation result(*this);
+        result -= other;
+        return result;
+    }
+
+    Evaluation operator-(Scalar other) const
+    {
+        Evaluation result(*this);
+        result -= other;
+        return result;
+    }
+
+    // negation (unary minus) operator
+    Evaluation operator-() const
+    {
+        Evaluation result;
+        result.value = -this->value;
+        for (int varIdx= 0; varIdx < size; ++varIdx)
+            result.derivatives[varIdx] = - this->derivatives[varIdx];
+
+        return result;
+    }
+
+    Evaluation operator*(const Evaluation& other) const
+    {
+        Evaluation result(*this);
+        result *= other;
+        return result;
+    }
+
+    Evaluation operator*(Scalar other) const
+    {
+        Evaluation result(*this);
+        result *= other;
+        return result;
+    }
+
+    Evaluation operator/(const Evaluation& other) const
+    {
+        Evaluation result(*this);
+        result /= other;
+        return result;
+    }
+
+    Evaluation operator/(Scalar other) const
+    {
+        Evaluation result(*this);
+        result /= other;
+        return result;
+    }
+
+    Evaluation& operator=(Scalar other)
+    {
+        this->value = other;
+        std::fill(this->derivatives.begin(), this->derivatives.end(), 0.0);
+        return *this;
+    }
+
+    // copy assignment from evaluation
+    Evaluation& operator=(const Evaluation& other)
+    {
+        this->value = other.value;
+        std::copy(other.derivatives.begin(), other.derivatives.end(), this->derivatives.begin());
+        return *this;
+    }
+
+    bool operator==(Scalar other) const
+    { return this->value == other; }
+
+    bool operator==(const Evaluation& other) const
+    {
+        if (this->value != other.value)
+            return false;
+
+        for (int varIdx= 0; varIdx < size; ++varIdx)
+            if (this->derivatives[varIdx] != other.derivatives[varIdx])
+                return false;
+
+        return true;
+    }
+
+    bool isSame(const Evaluation& other, Scalar tolerance) const
+    {
+        Scalar tmp = other.value - other.value;
+        if (std::abs(tmp) > tolerance && std::abs(tmp)/tolerance > 1.0)
+            return false;
+
+        for (int varIdx= 0; varIdx < size; ++varIdx) {
+            Scalar tmp = other.derivatives[varIdx] - this->derivatives[varIdx];
+            if (std::abs(tmp) > tolerance && std::abs(tmp)/tolerance > 1.0)
+                return false;
+        }
+
+        return true;
+    }
+
+    bool operator!=(const Evaluation& other) const
+    { return !operator==(other); }
+
+    bool operator>(Scalar other) const
+    { return this->value > other; }
+
+    bool operator>(const Evaluation& other) const
+    { return this->value > other.value; }
+
+    bool operator<(Scalar other) const
+    { return this->value < other; }
+
+    bool operator<(const Evaluation& other) const
+    { return this->value < other.value; }
+
+    bool operator>=(Scalar other) const
+    { return this->value >= other; }
+
+    bool operator>=(const Evaluation& other) const
+    { return this->value >= other.value; }
+
+    bool operator<=(Scalar other) const
+    { return this->value <= other; }
+
+    bool operator<=(const Evaluation& other) const
+    { return this->value <= other.value; }
+
+    // maybe this should be made 'private'...
+    Scalar value;
+    std::array<Scalar, size> derivatives;
+};
+
+template <class ScalarA, class Scalar, class VarSetTag, int numVars>
+bool operator<(const ScalarA& a, const Evaluation<Scalar, VarSetTag, numVars> &b)
+{ return b > a; }
+
+template <class ScalarA, class Scalar, class VarSetTag, int numVars>
+bool operator>(const ScalarA& a, const Evaluation<Scalar, VarSetTag, numVars> &b)
+{ return b < a; }
+
+template <class ScalarA, class Scalar, class VarSetTag, int numVars>
+bool operator<=(const ScalarA& a, const Evaluation<Scalar, VarSetTag, numVars> &b)
+{ return b >= a; }
+
+template <class ScalarA, class Scalar, class VarSetTag, int numVars>
+bool operator>=(const ScalarA& a, const Evaluation<Scalar, VarSetTag, numVars> &b)
+{ return b <= a; }
+
+template <class ScalarA, class Scalar, class VarSetTag, int numVars>
+bool operator!=(const ScalarA& a, const Evaluation<Scalar, VarSetTag, numVars> &b)
+{ return a != b.value; }
+
+template <class ScalarA, class Scalar, class VarSetTag, int numVars>
+Evaluation<Scalar, VarSetTag, numVars> operator+(const ScalarA& a, const Evaluation<Scalar, VarSetTag, numVars> &b)
+{
+    Evaluation<Scalar, VarSetTag, numVars> result(b);
+
+    result += a;
+
+    return result;
+}
+
+template <class ScalarA, class Scalar, class VarSetTag, int numVars>
+Evaluation<Scalar, VarSetTag, numVars> operator-(const ScalarA& a, const Evaluation<Scalar, VarSetTag, numVars> &b)
+{
+    Evaluation<Scalar, VarSetTag, numVars> result;
+
+    result.value = a - b.value;
+    for (int varIdx= 0; varIdx < numVars; ++varIdx)
+        result.derivatives[varIdx] = - b.derivatives[varIdx];
+
+    return result;
+}
+
+template <class ScalarA, class Scalar, class VarSetTag, int numVars>
+Evaluation<Scalar, VarSetTag, numVars> operator/(const ScalarA& a, const Evaluation<Scalar, VarSetTag, numVars> &b)
+{
+    Evaluation<Scalar, VarSetTag, numVars> result;
+
+    result.value = a/b.value;
+
+    // outer derivative
+    Scalar df_dg = - a/(b.value*b.value);
+    for (int varIdx= 0; varIdx < numVars; ++varIdx)
+        result.derivatives[varIdx] = df_dg*b.derivatives[varIdx];
+
+    return result;
+}
+
+template <class ScalarA, class Scalar, class VarSetTag, int numVars>
+Evaluation<Scalar, VarSetTag, numVars> operator*(const ScalarA& a, const Evaluation<Scalar, VarSetTag, numVars> &b)
+{
+    Evaluation<Scalar, VarSetTag, numVars> result;
+
+    result.value = a*b.value;
+    for (int varIdx= 0; varIdx < numVars; ++varIdx)
+        result.derivatives[varIdx] = a*b.derivatives[varIdx];
+
+    return result;
+}
+
+template <class Scalar, class VarSetTag, int numVars>
+std::ostream& operator<<(std::ostream& os, const Evaluation<Scalar, VarSetTag, numVars>& eval)
+{
+    os << eval.value;
+    return os;
+}
+
+} // namespace LocalAd
+} // namespace Opm
+
+// make the Dune matrix/vector classes happy. Obviously, this is not very elegant...
+
+#ifdef DUNE_DENSEMATRIX_HH
+// for those who are wondering: in C++ function templates cannot be partially
+// specialized, and function argument overloads must be known _before_ they are used. The
+// latter is what we do for the 'Dune::fvmeta::absreal()' function.
+//
+// consider the following test program:
+//
+// double foo(double i)
+// { return i; }
+//
+// void bar()
+// { std::cout << foo(0) << "\n"; }
+//
+// int foo(int i)
+// { return i + 1; }
+//
+// void foobar()
+// { std::cout << foo(0) << "\n"; }
+//
+// this will print '0' for bar() and '1' for foobar()...
+//
+#error "Due to some C++ peculiarity regarding function template specialization,"
+#error "the 'evaluation.hh' header must be included before Dune's 'densematrix.hh'!"
+#endif
+
+#include <dune/common/ftraits.hh>
+
+namespace Dune {
+template <class Scalar, class VarSetTag, int numVars>
+struct FieldTraits<Opm::LocalAd::Evaluation<Scalar, VarSetTag, numVars> >
+{
+public:
+    typedef Opm::LocalAd::Evaluation<Scalar, VarSetTag, numVars> field_type;
+    typedef Scalar real_type;
+};
+
+namespace fvmeta {
+template <class Scalar, class VarSetTag, int numVars>
+inline Scalar absreal(const Opm::LocalAd::Evaluation<Scalar, VarSetTag, numVars>& k)
+{
+    return std::abs(k.value);
+}
+
+}} // namespace fvmeta, Dune
+
+#endif

opm/material/localad/Math.hpp

@@ -0,0 +1,383 @@
+/*
+  Copyright (C) 2015 by Andreas Lauser
+
+  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/>.
+*/
+/*!
+ * \file
+ *
+ * \brief A number of commonly used algebraic functions for the localized OPM automatic
+ *        differentiation (AD) framework.
+ *
+ * This file provides AD variants of the the most commonly used functions of the <cmath>
+ * header file.
+ */
+#ifndef OPM_LOCAL_AD_MATH_HPP
+#define OPM_LOCAL_AD_MATH_HPP
+
+#include "Evaluation.hpp"
+
+namespace Opm {
+namespace LocalAd {
+// provide some algebraic functions
+template <class Scalar, class VarSetTag, int numVars>
+Evaluation<Scalar, VarSetTag, numVars> abs(const Evaluation<Scalar, VarSetTag, numVars>& x)
+{
+    Evaluation<Scalar, VarSetTag, numVars> result;
+
+    result.value = std::abs(x.value);
+
+    // derivatives use the chain rule
+    if (x.value < 0.0) {
+        for (unsigned curVarIdx = 0; curVarIdx < result.derivatives.size(); ++curVarIdx)
+            result.derivatives[curVarIdx] = -x.derivatives[curVarIdx];
+    }
+    else {
+        for (unsigned curVarIdx = 0; curVarIdx < result.derivatives.size(); ++curVarIdx)
+            result.derivatives[curVarIdx] = x.derivatives[curVarIdx];
+    }
+
+    return result;
+}
+
+template <class Scalar, class VarSetTag, int numVars>
+Evaluation<Scalar, VarSetTag, numVars> min(const Evaluation<Scalar, VarSetTag, numVars>& x1,
+                                           const Evaluation<Scalar, VarSetTag, numVars>& x2)
+{
+    Evaluation<Scalar, VarSetTag, numVars> result;
+
+    if (x1.value < x2.value) {
+        result.value = x1.value;
+
+        std::copy(x1.derivatives.begin(),
+                  x1.derivatives.end(),
+                  result.derivatives.begin());
+    }
+    else  {
+        result.value = x2.value;
+
+        std::copy(x2.derivatives.begin(),
+                  x2.derivatives.end(),
+                  result.derivatives.begin());
+    }
+
+    return result;
+}
+
+template <class ScalarA, class Scalar, class VarSetTag, int numVars>
+Evaluation<Scalar, VarSetTag, numVars> min(ScalarA x1,
+                                           const Evaluation<Scalar, VarSetTag, numVars>& x2)
+{
+    Evaluation<Scalar, VarSetTag, numVars> result;
+
+    if (x1 < x2.value) {
+        result.value = x1;
+
+        std::fill(result.derivatives.begin(),
+                  result.derivatives.end(),
+                  0.0);
+    }
+    else  {
+        result.value = x2.value;
+
+        std::copy(x2.derivatives.begin(),
+                  x2.derivatives.end(),
+                  result.derivatives.begin());
+    }
+
+    return result;
+}
+
+template <class ScalarB, class Scalar, class VarSetTag, int numVars>
+Evaluation<Scalar, VarSetTag, numVars> min(const Evaluation<Scalar, VarSetTag, numVars>& x2,
+                                           ScalarB x1)
+{ return min(x1, x2); }
+
+template <class Scalar, class VarSetTag, int numVars>
+Evaluation<Scalar, VarSetTag, numVars> max(const Evaluation<Scalar, VarSetTag, numVars>& x1,
+                                           const Evaluation<Scalar, VarSetTag, numVars>& x2)
+{
+    Evaluation<Scalar, VarSetTag, numVars> result;
+
+    if (x1.value > x2.value) {
+        result.value = x1.value;
+
+        std::copy(x1.derivatives.begin(),
+                  x1.derivatives.end(),
+                  result.derivatives.begin());
+    }
+    else  {
+        result.value = x2.value;
+
+        std::copy(x2.derivatives.begin(),
+                  x2.derivatives.end(),
+                  result.derivatives.begin());
+    }
+
+    return result;
+}
+
+template <class ScalarA, class Scalar, class VarSetTag, int numVars>
+Evaluation<Scalar, VarSetTag, numVars> max(ScalarA x1,
+                                           const Evaluation<Scalar, VarSetTag, numVars>& x2)
+{
+    Evaluation<Scalar, VarSetTag, numVars> result;
+
+    if (x1 > x2.value) {
+        result.value = x1;
+
+        std::fill(result.derivatives.begin(),
+                  result.derivatives.end(),
+                  0.0);
+    }
+    else  {
+        result.value = x2.value;
+
+        std::copy(x2.derivatives.begin(),
+                  x2.derivatives.end(),
+                  result.derivatives.begin());
+    }
+
+    return result;
+}
+
+template <class ScalarB, class Scalar, class VarSetTag, int numVars>
+Evaluation<Scalar, VarSetTag, numVars> max(const Evaluation<Scalar, VarSetTag, numVars>& x2,
+                                           ScalarB x1)
+{ return max(x1, x2); }
+
+template <class Scalar, class VarSetTag, int numVars>
+Evaluation<Scalar, VarSetTag, numVars> tan(const Evaluation<Scalar, VarSetTag, numVars>& x)
+{
+    Evaluation<Scalar, VarSetTag, numVars> result;
+
+    Scalar tmp = std::tan(x.value);
+    result.value = tmp;
+
+    // derivatives use the chain rule
+    Scalar df_dx = 1 + tmp*tmp;
+    for (unsigned curVarIdx = 0; curVarIdx < result.derivatives.size(); ++curVarIdx)
+        result.derivatives[curVarIdx] = df_dx*x.derivatives[curVarIdx];
+
+    return result;
+}
+
+template <class Scalar, class VarSetTag, int numVars>
+Evaluation<Scalar, VarSetTag, numVars> atan(const Evaluation<Scalar, VarSetTag, numVars>& x)
+{
+    Evaluation<Scalar, VarSetTag, numVars> result;
+
+    result.value = std::atan(x.value);
+
+    // derivatives use the chain rule
+    Scalar df_dx = 1/(1 + x.value*x.value);
+    for (unsigned curVarIdx = 0; curVarIdx < result.derivatives.size(); ++curVarIdx)
+        result.derivatives[curVarIdx] = df_dx*x.derivatives[curVarIdx];
+
+    return result;
+}
+
+template <class Scalar, class VarSetTag, int numVars>
+Evaluation<Scalar, VarSetTag, numVars> atan2(const Evaluation<Scalar, VarSetTag, numVars>& x,
+                                             const Evaluation<Scalar, VarSetTag, numVars>& y)
+{
+    Evaluation<Scalar, VarSetTag, numVars> result;
+
+    result.value = std::atan2(x.value, y.value);
+
+    // derivatives use the chain rule
+    Scalar alpha = 1/(1 + (x.value*x.value)/(y.value*y.value));
+    for (unsigned curVarIdx = 0; curVarIdx < result.derivatives.size(); ++curVarIdx) {
+        result.derivatives[curVarIdx] =
+            alpha
+            /(y.value*y.value)
+            *(x.derivatives[curVarIdx]*y.value - x.value*y.derivatives[curVarIdx]);
+    }
+
+    return result;
+}
+
+template <class Scalar, class VarSetTag, int numVars>
+Evaluation<Scalar, VarSetTag, numVars> sin(const Evaluation<Scalar, VarSetTag, numVars>& x)
+{
+    Evaluation<Scalar, VarSetTag, numVars> result;
+
+    result.value = std::sin(x.value);
+
+    // derivatives use the chain rule
+    Scalar df_dx = std::cos(x.value);
+    for (unsigned curVarIdx = 0; curVarIdx < result.derivatives.size(); ++curVarIdx)
+        result.derivatives[curVarIdx] = df_dx*x.derivatives[curVarIdx];
+
+    return result;
+}
+
+template <class Scalar, class VarSetTag, int numVars>
+Evaluation<Scalar, VarSetTag, numVars> asin(const Evaluation<Scalar, VarSetTag, numVars>& x)
+{
+    Evaluation<Scalar, VarSetTag, numVars> result;
+
+    result.value = std::asin(x.value);
+
+    // derivatives use the chain rule
+    Scalar df_dx = 1.0/std::sqrt(1 - x.value*x.value);
+    for (unsigned curVarIdx = 0; curVarIdx < result.derivatives.size(); ++curVarIdx)
+        result.derivatives[curVarIdx] = df_dx*x.derivatives[curVarIdx];
+
+    return result;
+}
+
+template <class Scalar, class VarSetTag, int numVars>
+Evaluation<Scalar, VarSetTag, numVars> cos(const Evaluation<Scalar, VarSetTag, numVars>& x)
+{
+    Evaluation<Scalar, VarSetTag, numVars> result;
+
+    result.value = std::cos(x.value);
+
+    // derivatives use the chain rule
+    Scalar df_dx = -std::sin(x.value);
+    for (unsigned curVarIdx = 0; curVarIdx < result.derivatives.size(); ++curVarIdx)
+        result.derivatives[curVarIdx] = df_dx*x.derivatives[curVarIdx];
+
+    return result;
+}
+
+template <class Scalar, class VarSetTag, int numVars>
+Evaluation<Scalar, VarSetTag, numVars> acos(const Evaluation<Scalar, VarSetTag, numVars>& x)
+{
+    Evaluation<Scalar, VarSetTag, numVars> result;
+
+    result.value = std::acos(x.value);
+
+    // derivatives use the chain rule
+    Scalar df_dx = - 1.0/std::sqrt(1 - x.value*x.value);
+    for (unsigned curVarIdx = 0; curVarIdx < result.derivatives.size(); ++curVarIdx)
+        result.derivatives[curVarIdx] = df_dx*x.derivatives[curVarIdx];
+
+    return result;
+}
+
+template <class Scalar, class VarSetTag, int numVars>
+Evaluation<Scalar, VarSetTag, numVars> sqrt(const Evaluation<Scalar, VarSetTag, numVars>& x)
+{
+    Evaluation<Scalar, VarSetTag, numVars> result;
+
+    Scalar sqrt_x = std::sqrt(x.value);
+    result.value = sqrt_x;
+
+    // derivatives use the chain rule
+    Scalar df_dx = 0.5/sqrt_x;
+    for (unsigned curVarIdx = 0; curVarIdx < result.derivatives.size(); ++curVarIdx) {
+        result.derivatives[curVarIdx] = df_dx*x.derivatives[curVarIdx];
+    }
+
+    return result;
+}
+
+template <class Scalar, class VarSetTag, int numVars>
+Evaluation<Scalar, VarSetTag, numVars> exp(const Evaluation<Scalar, VarSetTag, numVars>& x)
+{
+    Evaluation<Scalar, VarSetTag, numVars> result;
+
+    Scalar exp_x = std::exp(x.value);
+    result.value = exp_x;
+
+    // derivatives use the chain rule
+    Scalar df_dx = exp_x;
+    for (unsigned curVarIdx = 0; curVarIdx < result.derivatives.size(); ++curVarIdx)
+        result.derivatives[curVarIdx] = df_dx*x.derivatives[curVarIdx];
+
+    return result;
+}
+
+// exponentiation of arbitrary base with a fixed constant
+template <class Scalar, class VarSetTag, int numVars>
+Evaluation<Scalar, VarSetTag, numVars> pow(const Evaluation<Scalar, VarSetTag, numVars>& base, Scalar exp)
+{
+    Evaluation<Scalar, VarSetTag, numVars> result;
+
+    Scalar pow_x = std::pow(base.value, exp);
+    result.value = pow_x;
+
+    // derivatives use the chain rule
+    Scalar df_dx = pow_x/base.value*exp;
+    for (unsigned curVarIdx = 0; curVarIdx < result.derivatives.size(); ++curVarIdx)
+        result.derivatives[curVarIdx] = df_dx*base.derivatives[curVarIdx];
+
+    return result;
+}
+
+// exponentiation of constant base with an arbitrary exponent
+template <class Scalar, class VarSetTag, int numVars>
+Evaluation<Scalar, VarSetTag, numVars> pow(Scalar base, const Evaluation<Scalar, VarSetTag, numVars>& exp)
+{
+    Evaluation<Scalar, VarSetTag, numVars> result;
+
+    Scalar lnBase = std::log(base);
+    result.value = std::exp(lnBase*exp.value);
+
+    // derivatives use the chain rule
+    Scalar df_dx = lnBase*result.value;
+    for (unsigned curVarIdx = 0; curVarIdx < result.derivatives.size(); ++curVarIdx)
+        result.derivatives[curVarIdx] = df_dx*exp.derivatives[curVarIdx];
+
+    return result;
+}
+
+// this is the most expensive power function. Computationally it is pretty expensive, so
+// one of the above two variants above should be preferred if possible.
+template <class Scalar, class VarSetTag, int numVars>
+Evaluation<Scalar, VarSetTag, numVars> pow(const Evaluation<Scalar, VarSetTag, numVars>& base, const Evaluation<Scalar, VarSetTag, numVars>& exp)
+{
+    Evaluation<Scalar, VarSetTag, numVars> result;
+
+    Scalar valuePow = std::pow(base.value, exp.value);
+    result.value = valuePow;
+
+    // use the chain rule for the derivatives. since both, the base and the exponent can
+    // potentially depend on the variable set, calculating these is quite elaborate...
+    Scalar f = base.value;
+    Scalar g = exp.value;
+    Scalar logF = std::log(f);
+    for (unsigned curVarIdx = 0; curVarIdx < result.derivatives.size(); ++curVarIdx) {
+        Scalar fPrime = base.derivatives[curVarIdx];
+        Scalar gPrime = exp.derivatives[curVarIdx];
+        result.derivatives[curVarIdx] = (g*fPrime/f + logF*gPrime) * valuePow;
+    }
+
+    return result;
+}
+
+template <class Scalar, class VarSetTag, int numVars>
+Evaluation<Scalar, VarSetTag, numVars> log(const Evaluation<Scalar, VarSetTag, numVars>& x)
+{
+    Evaluation<Scalar, VarSetTag, numVars> result;
+
+    result.value = std::log(x.value);
+
+    // derivatives use the chain rule
+    Scalar df_dx = 1/x.value;
+    for (unsigned curVarIdx = 0; curVarIdx < result.derivatives.size(); ++curVarIdx)
+        result.derivatives[curVarIdx] = df_dx*x.derivatives[curVarIdx];
+
+    return result;
+}
+
+} // namespace LocalAd
+
+}
+
+#endif

tests/test_localad.cpp

@@ -0,0 +1,544 @@
+/*
+  Copyright (C) 2015 by Andreas Lauser
+
+  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/>.
+*/
+/*!
+ * \file
+ *
+ * \brief Low-level tests for the localized automatic differentiation (AD) framework.
+ */
+#include "config.h"
+
+#include <iostream>
+#include <array>
+#include <cmath>
+#include <algorithm>
+#include <cassert>
+#include <stdexcept>
+
+#include <opm/material/common/Unused.hpp>
+
+#include <opm/material/localad/Evaluation.hpp>
+#include <opm/material/localad/Math.hpp>
+
+struct TestVariables
+{
+    static const int size = 3;
+
+    static const int temperatureIdx = 0;
+    static const int pressureIdx = 1;
+    static const int saturationIdx = 2;
+};
+
+template <class Scalar, class VariablesDescriptor>
+void testOperators()
+{
+    typedef Opm::LocalAd::Evaluation<Scalar, VariablesDescriptor, VariablesDescriptor::size> Eval;
+
+    // test the constructors of the Opm::LocalAd::Evaluation class
+    const Scalar c = 1.234;
+    const Scalar x = 4.567;
+    const Scalar y = 8.910;
+    const Eval cEval = Eval::createConstant(c);
+    const Eval OPM_UNUSED c2Eval = c;
+    const Eval xEval = Eval::createVariable(x, 0);
+    const Eval yEval = Eval::createVariable(y, 1);
+
+    // test the non-inplace operators
+    {
+        Eval a = xEval + yEval;
+        if (a.value != x + y)
+            throw std::logic_error("oops: operator+");
+
+        Eval b = xEval + c;
+        if (b.value != x + c)
+            throw std::logic_error("oops: operator+");
+
+        Eval d = xEval + cEval;
+        if (d.value != x + c)
+            throw std::logic_error("oops: operator+");
+    }
+
+    {
+        Eval a = xEval - yEval;
+        if (a.value != x - y)
+            throw std::logic_error("oops: operator-");
+
+        Eval b = xEval - c;
+        if (b.value != x - c)
+            throw std::logic_error("oops: operator-");
+
+        Eval d = xEval - cEval;
+        if (d.value != x - c)
+            throw std::logic_error("oops: operator-");
+    }
+
+    {
+        Eval a = xEval*yEval;
+        if (a.value != x*y)
+            throw std::logic_error("oops: operator*");
+
+        Eval b = xEval*c;
+        if (b.value != x*c)
+            throw std::logic_error("oops: operator*");
+
+        Eval d = xEval*cEval;
+        if (d.value != x*c)
+            throw std::logic_error("oops: operator*");
+    }
+
+    {
+        Eval a = xEval/yEval;
+        if (a.value != x/y)
+            throw std::logic_error("oops: operator/");
+
+        Eval b = xEval/c;
+        if (b.value != x/c)
+            throw std::logic_error("oops: operator/");
+
+        Eval d = xEval/cEval;
+        if (d.value != x/c)
+            throw std::logic_error("oops: operator/");
+    }
+
+    // test the inplace operators
+    {
+        Eval a = xEval;
+        a += yEval;
+        if (a.value != x + y)
+            throw std::logic_error("oops: operator+");
+
+        Eval b = xEval;
+        b += c;
+        if (b.value != x + c)
+            throw std::logic_error("oops: operator+");
+
+        Eval d = xEval;
+        d += cEval;
+        if (d.value != x + c)
+            throw std::logic_error("oops: operator+");
+    }
+
+    {
+        Eval a = xEval;
+        a -= yEval;
+        if (a.value != x - y)
+            throw std::logic_error("oops: operator-");
+
+        Eval b = xEval;
+        b -= c;
+        if (b.value != x - c)
+            throw std::logic_error("oops: operator-");
+
+        Eval d = xEval;
+        d -= cEval;
+        if (d.value != x - c)
+            throw std::logic_error("oops: operator-");
+    }
+
+    {
+        Eval a = xEval;
+        a *= yEval;
+        if (a.value != x*y)
+            throw std::logic_error("oops: operator*");
+
+        Eval b = xEval;
+        b *= c;
+        if (b.value != x*c)
+            throw std::logic_error("oops: operator*");
+
+        Eval d = xEval;
+        d *= cEval;
+        if (d.value != x*c)
+            throw std::logic_error("oops: operator*");
+    }
+
+    {
+        Eval a = xEval;
+        a /= yEval;
+        if (a.value != x/y)
+            throw std::logic_error("oops: operator/");
+
+        Eval b = xEval;
+        b /= c;
+        if (b.value != x/c)
+            throw std::logic_error("oops: operator/");
+
+        Eval d = xEval;
+        d /= cEval;
+        if (d.value != x/c)
+            throw std::logic_error("oops: operator/");
+    }
+
+    {
+        Eval a = 1.0;
+        Eval b = 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 Scalar, class VariablesDescriptor, class AdFn, class ClassicFn>
+void test1DFunction(AdFn* adFn, ClassicFn* classicFn, Scalar xMin = 1e-6, Scalar xMax = 1000)
+{
+    typedef Opm::LocalAd::Evaluation<Scalar, VariablesDescriptor, VariablesDescriptor::size> Eval;
+
+    int n = 100*1000;
+    for (int i = 0; i < n; ++ i) {
+        Scalar x = Scalar(i)/(n - 1)*(xMax - xMin) + xMin;
+
+        const auto& xEval = Eval::createVariable(x, 0);
+        const Eval& yEval = adFn(xEval);
+
+        const Scalar eps = 1e-10;
+        Scalar y = classicFn(x);
+        Scalar yStar1 = classicFn(x - eps);
+        Scalar yStar2 = classicFn(x + eps);
+        Scalar yPrime = (yStar2 - yStar1)/(2*eps);
+
+        if (y != yEval.value)
+            throw std::logic_error("oops: value");
+
+        Scalar deltaAbs = std::abs(yPrime - yEval.derivatives[0]);
+        Scalar deltaRel = std::abs(deltaAbs/yPrime);
+        if (deltaAbs > 1e-3 && deltaRel > 1e-3)
+            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.derivatives[0])
+                                   + " delta: " + std::to_string((long double) std::abs(yPrime - yEval.derivatives[0])));
+    }
+}
+
+template <class Scalar,
+          class VariablesDescriptor,
+          class AdFn,
+          class ClassicFn>
+void test2DFunction1(AdFn* adFn, ClassicFn* classicFn, Scalar xMin, Scalar xMax, Scalar y)
+{
+    typedef Opm::LocalAd::Evaluation<Scalar, VariablesDescriptor, VariablesDescriptor::size> Eval;
+
+    int n = 100*1000;
+    for (int i = 0; i < n; ++ i) {
+        Scalar x = Scalar(i)/(n - 1)*(xMax - xMin) + xMin;
+
+        const auto& xEval = Eval::createVariable(x, 0);
+        const auto& yEval = Eval::createConstant(y);
+        const Eval& zEval = adFn(xEval, yEval);
+
+        const Scalar eps = 1e-10;
+        Scalar z = classicFn(x, y);
+        Scalar zStar1 = classicFn(x - eps, y);
+        Scalar zStar2 = classicFn(x + eps, y);
+        Scalar zPrime = (zStar2 - zStar1)/(2*eps);
+
+        if (z != zEval.value)
+            throw std::logic_error("oops: value");
+
+        Scalar deltaAbs = std::abs(zPrime - zEval.derivatives[0]);
+        Scalar deltaRel = std::abs(deltaAbs/zPrime);
+        if (deltaAbs > 1e-3 && deltaRel > 1e-3)
+            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.derivatives[0])
+                                   + " delta: " + std::to_string((long double) std::abs(zPrime - zEval.derivatives[0])));
+    }
+}
+
+template <class Scalar,
+          class VariablesDescriptor,
+          class AdFn,
+          class ClassicFn>
+void test2DFunction2(AdFn* adFn, ClassicFn* classicFn, Scalar x, Scalar yMin, Scalar yMax)
+{
+    typedef Opm::LocalAd::Evaluation<Scalar, VariablesDescriptor, VariablesDescriptor::size> Eval;
+
+    int n = 100*1000;
+    for (int i = 0; i < n; ++ i) {
+        Scalar y = Scalar(i)/(n - 1)*(yMax - yMin) + yMin;
+
+        const auto& xEval = Eval::createConstant(x);
+        const auto& yEval = Eval::createVariable(y, 1);
+        const Eval& zEval = adFn(xEval, yEval);
+
+        const Scalar eps = 1e-10;
+        Scalar z = classicFn(x, y);
+        Scalar zStar1 = classicFn(x, y - eps);
+        Scalar zStar2 = classicFn(x, y + eps);
+        Scalar zPrime = (zStar2 - zStar1)/(2*eps);
+
+        if (z != zEval.value)
+            throw std::logic_error("oops: value");
+
+        Scalar deltaAbs = std::abs(zPrime - zEval.derivatives[1]);
+        Scalar deltaRel = std::abs(deltaAbs/zPrime);
+        if (deltaAbs > 1e-3 && deltaRel > 1e-3)
+            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.derivatives[1])
+                                   + " delta: " + std::to_string((long double) std::abs(zPrime - zEval.derivatives[1])));
+    }
+}
+
+template <class Scalar, class VariablesDescriptor>
+void testPowBase(Scalar baseMin = 1e-2, Scalar baseMax = 100)
+{
+    typedef Opm::LocalAd::Evaluation<Scalar, VariablesDescriptor, VariablesDescriptor::size> Eval;
+
+    Scalar exp = 1.234;
+    const auto& expEval = Eval::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 = Eval::createVariable(base, 0);
+        const Eval& zEval1 = pow(baseEval, exp);
+        const Eval& zEval2 = pow(baseEval, expEval);
+
+        const Scalar eps = 1e-5;
+        Scalar z = pow(base, exp);
+        Scalar zStar1 = pow(base - eps, exp);
+        Scalar zStar2 = pow(base + eps, exp);
+        Scalar zPrime = (zStar2 - zStar1)/(2*eps);
+
+        if (z != zEval2.value)
+            throw std::logic_error("oops: value");
+
+        Scalar deltaAbs = std::abs(zPrime - zEval1.derivatives[0]);
+        Scalar deltaRel = std::abs(deltaAbs/zPrime);
+        if (deltaAbs > 1e-3 && deltaRel > 1e-3)
+            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.derivatives[0])
+                                   + " delta: " + std::to_string((long double) std::abs(zPrime - zEval1.derivatives[0])));
+
+        if (!zEval1.isSame(zEval2, /*tolerance=*/1e-9))
+            throw std::logic_error("oops: pow(Eval, Scalar) != pow(Eval, Eval)");
+    }
+}
+
+template <class Scalar, class VariablesDescriptor>
+void testPowExp(Scalar expMin = -100, Scalar expMax = 100)
+{
+    typedef Opm::LocalAd::Evaluation<Scalar, VariablesDescriptor, VariablesDescriptor::size> Eval;
+
+    Scalar base = 1.234;
+    const auto& baseEval = Eval::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 = Eval::createVariable(exp, 1);
+
+        const Eval& zEval1 = pow(base, expEval);
+        const Eval& zEval2 = pow(baseEval, expEval);
+
+        const Scalar eps = 1e-8;
+        Scalar z = pow(base, exp);
+        Scalar zStar1 = pow(base, exp - eps);
+        Scalar zStar2 = pow(base, exp + eps);
+        Scalar zPrime = (zStar2 - zStar1)/(2*eps);
+
+        if (z != zEval2.value)
+            throw std::logic_error("oops: value");
+
+        Scalar deltaAbs = std::abs(zPrime - zEval1.derivatives[1]);
+        Scalar deltaRel = std::abs(deltaAbs/zPrime);
+        if (deltaAbs > 1e-3 && deltaRel > 1e-3)
+            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.derivatives[1])
+                                   + " delta: " + std::to_string((long double) std::abs(zPrime - zEval1.derivatives[1])));
+
+        if (!zEval1.isSame(zEval2, /*tolerance=*/1e-5))
+            throw std::logic_error("oops: pow(Eval, Scalar) != pow(Eval, Eval)");
+    }
+}
+
+template <class Scalar, class VariablesDescriptor>
+void testAtan2()
+{
+    typedef Opm::LocalAd::Evaluation<Scalar, VariablesDescriptor, VariablesDescriptor::size> Eval;
+
+    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 = Eval::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 = Eval::createVariable(y, 0);
+            const Eval& zEval = atan2(xEval, yEval);
+
+            const Scalar eps = 1e-8;
+            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 (z != zEval.value)
+                throw std::logic_error("oops: value");
+
+            Scalar deltaAbs = std::abs(zPrime - zEval.derivatives[0]);
+            Scalar deltaRel = std::abs(deltaAbs/zPrime);
+            if (deltaAbs > 1e-3 && deltaRel > 1e-3)
+                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.derivatives[0])
+                                       + " delta: " + std::to_string((long double) std::abs(zPrime - zEval.derivatives[0])));
+
+        }
+    }
+}
+
+double myScalarMin(double a, double b)
+{ return std::min(a, b); }
+
+double myScalarMax(double a, double b)
+{ return std::max(a, b); }
+
+int main()
+{
+    typedef double Scalar;
+
+    typedef TestVariables VarsDescriptor;
+
+    // the following is commented out because it is supposed to produce a compiler
+    // error. This is the case since the function does not calculate the derivatives
+    // w.r.t. Pressure but they have been requested...
+    //const auto& result2 = Opm::LocalAd::sqrt(TemperatureEval::createVariable<Pressure>(4.0));
+
+    std::cout << "testing operators and constructors\n";
+    testOperators<Scalar, VarsDescriptor>();
+
+    std::cout << "testing min()\n";
+    test2DFunction1<Scalar, VarsDescriptor>(Opm::LocalAd::min<Scalar, VarsDescriptor, VarsDescriptor::size>,
+                                            myScalarMin,
+                                            -1000, 1000,
+                                            /*p=*/1.234);
+
+    test2DFunction2<Scalar, VarsDescriptor>(Opm::LocalAd::min<Scalar, VarsDescriptor, VarsDescriptor::size>,
+                                            myScalarMin,
+                                            /*T=*/1.234,
+                                            -1000, 1000);
+
+    std::cout << "testing max()\n";
+    test2DFunction1<Scalar, VarsDescriptor>(Opm::LocalAd::max<Scalar, VarsDescriptor, VarsDescriptor::size>,
+                                            myScalarMax,
+                                            -1000, 1000,
+                                            /*p=*/1.234);
+
+    test2DFunction2<Scalar, VarsDescriptor>(Opm::LocalAd::max<Scalar, VarsDescriptor, VarsDescriptor::size>,
+                                            myScalarMax,
+                                            /*T=*/1.234,
+                                            -1000, 1000);
+
+    std::cout << "testing pow()\n";
+    testPowBase<Scalar, VarsDescriptor>();
+    testPowExp<Scalar, VarsDescriptor>();
+
+    std::cout << "testing abs()\n";
+    test1DFunction<Scalar, VarsDescriptor>(Opm::LocalAd::abs<Scalar, VarsDescriptor, VarsDescriptor::size>,
+                                           static_cast<Scalar (*)(Scalar)>(std::abs));
+
+    std::cout << "testing sqrt()\n";
+    test1DFunction<Scalar, VarsDescriptor>(Opm::LocalAd::sqrt<Scalar, VarsDescriptor, VarsDescriptor::size>,
+                                           static_cast<Scalar (*)(Scalar)>(std::sqrt));
+
+    std::cout << "testing sin()\n";
+    test1DFunction<Scalar, VarsDescriptor>(Opm::LocalAd::sin<Scalar, VarsDescriptor, VarsDescriptor::size>,
+                                           static_cast<Scalar (*)(Scalar)>(std::sin),
+                                           0, 2*M_PI);
+
+    std::cout << "testing asin()\n";
+    test1DFunction<Scalar, VarsDescriptor>(Opm::LocalAd::asin<Scalar, VarsDescriptor, VarsDescriptor::size>,
+                                           static_cast<Scalar (*)(Scalar)>(std::asin),
+                                           -1.0, 1.0);
+
+    std::cout << "testing cos()\n";
+    test1DFunction<Scalar, VarsDescriptor>(Opm::LocalAd::cos<Scalar, VarsDescriptor, VarsDescriptor::size>,
+                                           static_cast<Scalar (*)(Scalar)>(std::cos),
+                                           0, 2*M_PI);
+
+    std::cout << "testing acos()\n";
+    test1DFunction<Scalar, VarsDescriptor>(Opm::LocalAd::acos<Scalar, VarsDescriptor, VarsDescriptor::size>,
+                                           static_cast<Scalar (*)(Scalar)>(std::acos),
+                                           -1.0, 1.0);
+
+    std::cout << "testing tan()\n";
+    test1DFunction<Scalar, VarsDescriptor>(Opm::LocalAd::tan<Scalar, VarsDescriptor, VarsDescriptor::size>,
+                                           static_cast<Scalar (*)(Scalar)>(std::tan),
+                                           -M_PI / 2 * 0.95, M_PI / 2 * 0.95);
+
+    std::cout << "testing atan()\n";
+    test1DFunction<Scalar, VarsDescriptor>(Opm::LocalAd::atan<Scalar, VarsDescriptor, VarsDescriptor::size>,
+                                           static_cast<Scalar (*)(Scalar)>(std::atan),
+                                           -10*1000.0, 10*1000.0);
+
+    std::cout << "testing atan2()\n";
+    testAtan2<Scalar, VarsDescriptor>();
+
+    std::cout << "testing exp()\n";
+    test1DFunction<Scalar, VarsDescriptor>(Opm::LocalAd::exp<Scalar, VarsDescriptor, VarsDescriptor::size>,
+                                           static_cast<Scalar (*)(Scalar)>(std::exp),
+                                           -100, 100);
+
+    std::cout << "testing log()\n";
+    test1DFunction<Scalar, VarsDescriptor>(Opm::LocalAd::log<Scalar, VarsDescriptor, VarsDescriptor::size>,
+                                           static_cast<Scalar (*)(Scalar)>(std::log),
+                                           1e-6, 1e9);
+
+    return 0;
+}