From b9d9f893d9c0c5f3a79421e1eceb53c2d1a89ed0 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Thu, 21 May 2015 15:33:01 +0200 Subject: [PATCH] introduce the concept of "math toolboxes" these are "traits" classes and provide a way to access the value of function evaluations, conditional access to its value (i.e., to forward them if the target object of an assignment includes the derivatives or use the function value if not) and some algebraic functions. the main idea is to be able to abstract the differences between plain scalars and function evaluation... --- opm/material/common/MathToolbox.hpp | 207 ++++++++++++++++++++++++++++ opm/material/localad/Math.hpp | 84 +++++++++++ 2 files changed, 291 insertions(+) create mode 100644 opm/material/common/MathToolbox.hpp diff --git a/opm/material/common/MathToolbox.hpp b/opm/material/common/MathToolbox.hpp new file mode 100644 index 000000000..94e0a872a --- /dev/null +++ b/opm/material/common/MathToolbox.hpp @@ -0,0 +1,207 @@ +/* + 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 . +*/ +/*! + * \file + * + * \brief A traits class which provides basic mathematical functions for arbitrary scalar + * floating point values. + * + * The reason why this is done in such a complicated way is to enable other approaches, + * in particular automatic differentiation based ones. + */ +#ifndef OPM_MATERIAL_MATH_TOOLBOX_HPP +#define OPM_MATERIAL_MATH_TOOLBOX_HPP + +#include +#include +#include + +namespace Opm { +template ::value > +struct MathToolbox; + +template ::value, + bool rhsIsScalar = std::is_floating_point::value> +struct ToLhsEvalHelper; + +template +struct ToLhsEvalHelper { + static LhsEval exec(const RhsEval& eval) + { return eval; } +}; + +template +struct ToLhsEvalHelper { + static LhsEval exec(const RhsEval& eval) + { return eval.value; } +}; + +template +struct ToLhsEvalHelper { + static LhsEval exec(const RhsEval& eval) + { return eval; } +}; + +/* + * \brief A traits class which provides basic mathematical functions for arbitrary scalar + * floating point values. + * + * The reason why this is done in such a complicated way is to enable other approaches, + * in particular automatic differentiation based ones. + */ +template +struct MathToolbox +{ +public: + /*! + * \brief The type used to represent scalar values + */ + typedef ScalarT Scalar; + + /*! + * \brief The type used to represent function evaluations + * + * In general, that is the value of the function plus a number of derivatives at the + * evaluation point. In the case of the scalar toolbox, no derivatives will be + * evaluated. + */ + typedef ScalarT Evaluation; + + /*! + * \brief Return the value of the function at a given evaluation point. + * + * For this toolbox, an evaluation is the value, so this method is the identity + * function. + */ + static Scalar value(const Evaluation& eval) + { return eval; } + + /*! + * \brief Given a scalar value, return an evaluation of a constant function. + * + * For this toolbox, an evaluation is the value, so this method is the identity + * function. In general, this returns an evaluation object for which all derivatives + * are zero. + */ + static Scalar createConstant(Scalar value) + { return value; } + + /*! + * \brief Given a scalar value, return an evaluation of a linear function. + * + * i.e., Create an evaluation which represents f(x) = x and the derivatives with + * regard to x. For scalar evaluations (which do not consider derivatives), this + * method does nothing. + */ + static Scalar createVariable(Scalar value, int varIdx) + { return value; } + + /*! + * \brief Given a function evaluation, constrain it to its value (if necessary). + * + * If the left hand side is a scalar and the right hand side is an evaluation, the + * scalar gets the value of the right hand side assigned. Also if both sides are + * scalars, this method returns the identity. The final case (left hand side being an + * evaluation, right hand side is a scalar) yields a compiler error. + * + * The purpose of this method is to be able to transparantly use evaluation objects + * in scalar computations. + */ + template + static LhsEval toLhs(const Evaluation& eval) + { return ToLhsEvalHelper::exec(eval); } + + /*! + * \brief Pass a value through if it is an evaluation, or create a constant + * evaluation if it is a scalar. + * + * In some sense, this method is the opposite of "toLhs()": If the function argument + * is a Scalar, an Evaluation which represents a constant value is returned, if the + * argument is an Evaluation, it is returned as is. This method makes it possible to + * uniformly handle the cases where some condition is either given by a constant + * value or as a value which depends on other variables. (E.g. boundary conditions.) + */ + static Scalar passThroughOrCreateConstant(Scalar value) + { return value; } + + //////////// + // arithmetic functions + //////////// + + //! The maximum of two arguments + static Scalar max(Scalar arg1, Scalar arg2) + { return std::max(arg1, arg2); } + + //! The minimum of two arguments + static Scalar min(Scalar arg1, Scalar arg2) + { return std::min(arg1, arg2); } + + //! The absolute value + static Scalar abs(Scalar arg) + { return std::abs(arg); } + + //! The tangens of a value + static Scalar tan(Scalar arg) + { return std::tan(arg); } + + //! The arcus tangens of a value + static Scalar atan(Scalar arg) + { return std::atan(arg); } + + //! The arcus tangens of a value + static Scalar atan2(Scalar arg1, Scalar arg2) + { return std::atan2(arg1, arg2); } + + //! The sine of a value + static Scalar sin(Scalar arg) + { return std::sin(arg); } + + //! The arcus sine of a value + static Scalar asin(Scalar arg) + { return std::asin(arg); } + + //! The cosine of a value + static Scalar cos(Scalar arg) + { return std::cos(arg); } + + //! The arcus cosine of a value + static Scalar acos(Scalar arg) + { return std::acos(arg); } + + //! The square root of a value + static Scalar sqrt(Scalar arg) + { return std::sqrt(arg); } + + //! The natural exponentiation of a value + static Scalar exp(Scalar arg) + { return std::exp(arg); } + + //! The natural logarithm of a value + static Scalar log(Scalar arg) + { return std::log(arg); } + + //! Exponentiation to an arbitrary base + static Scalar pow(Scalar base, Scalar exp) + { return std::pow(base, exp); } +}; + +} // namespace Opm + +#endif diff --git a/opm/material/localad/Math.hpp b/opm/material/localad/Math.hpp index 1246b2f61..7e18bacf1 100644 --- a/opm/material/localad/Math.hpp +++ b/opm/material/localad/Math.hpp @@ -30,6 +30,8 @@ #include "Evaluation.hpp" +#include + namespace Opm { namespace LocalAd { // provide some algebraic functions @@ -378,6 +380,88 @@ Evaluation log(const Evaluation +struct MathToolbox, false> +{ +private: +public: + typedef ScalarT Scalar; + typedef Opm::LocalAd::Evaluation Evaluation; + + static Scalar value(const Evaluation& eval) + { return eval.value; } + + static Evaluation createConstant(Scalar value) + { return Evaluation::createConstant(value); } + + static Evaluation createVariable(Scalar value, int varIdx) + { return Evaluation::createVariable(value, varIdx); } + + template + static LhsEval toLhs(const Evaluation& eval) + { return ToLhsEvalHelper::exec(eval); } + + static const Evaluation passThroughOrCreateConstant(Scalar value) + { return createConstant(value); } + + static const Evaluation& passThroughOrCreateConstant(const Evaluation& eval) + { return eval; } + + + // arithmetic functions + template + static Evaluation max(const Arg1Eval& arg1, const Arg2Eval& arg2) + { return Opm::LocalAd::max(arg1, arg2); } + + template + static Evaluation min(const Arg1Eval& arg1, const Arg2Eval& arg2) + { return Opm::LocalAd::min(arg1, arg2); } + + static Evaluation abs(const Evaluation& arg) + { return Opm::LocalAd::abs(arg); } + + static Evaluation tan(const Evaluation& arg) + { return Opm::LocalAd::tan(arg); } + + static Evaluation atan(const Evaluation& arg) + { return Opm::LocalAd::atan(arg); } + + static Evaluation atan2(const Evaluation& arg1, const Evaluation& arg2) + { return Opm::LocalAd::atan2(arg1, arg2); } + + static Evaluation sin(const Evaluation& arg) + { return Opm::LocalAd::sin(arg); } + + static Evaluation asin(const Evaluation& arg) + { return Opm::LocalAd::asin(arg); } + + static Evaluation cos(const Evaluation& arg) + { return Opm::LocalAd::cos(arg); } + + static Evaluation acos(const Evaluation& arg) + { return Opm::LocalAd::acos(arg); } + + static Evaluation sqrt(const Evaluation& arg) + { return Opm::LocalAd::sqrt(arg); } + + static Evaluation exp(const Evaluation& arg) + { return Opm::LocalAd::exp(arg); } + + static Evaluation log(const Evaluation& arg) + { return Opm::LocalAd::log(arg); } + + static Evaluation pow(const Evaluation& arg1, typename Evaluation::Scalar arg2) + { return Opm::LocalAd::pow(arg1, arg2); } + + static Evaluation pow(typename Evaluation::Scalar arg1, const Evaluation& arg2) + { return Opm::LocalAd::pow(arg1, arg2); } + + static Evaluation pow(const Evaluation& arg1, const Evaluation& arg2) + { return Opm::LocalAd::pow(arg1, arg2); } +}; + } #endif