diff --git a/opm/material/localad/Evaluation.hpp b/opm/material/localad/Evaluation.hpp
new file mode 100644
index 000000000..664f2f7f7
--- /dev/null
+++ b/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
diff --git a/opm/material/localad/Math.hpp b/opm/material/localad/Math.hpp
new file mode 100644
index 000000000..1246b2f61
--- /dev/null
+++ b/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
diff --git a/tests/test_localad.cpp b/tests/test_localad.cpp
new file mode 100644
index 000000000..033bb9f46
--- /dev/null
+++ b/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;
+}