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...
This commit is contained in:
Andreas Lauser 2015-05-21 15:33:01 +02:00
parent d34bbc57a9
commit b9d9f893d9
2 changed files with 291 additions and 0 deletions

View File

@ -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 <http://www.gnu.org/licenses/>.
*/
/*!
* \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 <cmath>
#include <algorithm>
#include <type_traits>
namespace Opm {
template <class Evaluation, bool isScalar = std::is_floating_point<Evaluation>::value >
struct MathToolbox;
template <class LhsEval, class RhsEval,
bool lhsIsScalar = std::is_floating_point<LhsEval>::value,
bool rhsIsScalar = std::is_floating_point<RhsEval>::value>
struct ToLhsEvalHelper;
template <class LhsEval, class RhsEval>
struct ToLhsEvalHelper<LhsEval, RhsEval, true, true> {
static LhsEval exec(const RhsEval& eval)
{ return eval; }
};
template <class LhsEval, class RhsEval>
struct ToLhsEvalHelper<LhsEval, RhsEval, true, false> {
static LhsEval exec(const RhsEval& eval)
{ return eval.value; }
};
template <class LhsEval, class RhsEval>
struct ToLhsEvalHelper<LhsEval, RhsEval, false, false> {
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 <class ScalarT>
struct MathToolbox<ScalarT, true>
{
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 <class LhsEval>
static LhsEval toLhs(const Evaluation& eval)
{ return ToLhsEvalHelper<Evaluation, LhsEval>::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

View File

@ -30,6 +30,8 @@
#include "Evaluation.hpp"
#include <opm/material/common/MathToolbox.hpp>
namespace Opm {
namespace LocalAd {
// provide some algebraic functions
@ -378,6 +380,88 @@ Evaluation<Scalar, VarSetTag, numVars> log(const Evaluation<Scalar, VarSetTag, n
} // namespace LocalAd
// a kind of traits class for the automatic differentiation case. (The toolbox for the
// scalar case is provided by the MathToolbox.hpp header file.)
template <class ScalarT, class VariableSetTag, int numVars>
struct MathToolbox<Opm::LocalAd::Evaluation<ScalarT, VariableSetTag, numVars>, false>
{
private:
public:
typedef ScalarT Scalar;
typedef Opm::LocalAd::Evaluation<ScalarT, VariableSetTag, numVars> 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 <class LhsEval>
static LhsEval toLhs(const Evaluation& eval)
{ return ToLhsEvalHelper<LhsEval, Evaluation>::exec(eval); }
static const Evaluation passThroughOrCreateConstant(Scalar value)
{ return createConstant(value); }
static const Evaluation& passThroughOrCreateConstant(const Evaluation& eval)
{ return eval; }
// arithmetic functions
template <class Arg1Eval, class Arg2Eval>
static Evaluation max(const Arg1Eval& arg1, const Arg2Eval& arg2)
{ return Opm::LocalAd::max(arg1, arg2); }
template <class Arg1Eval, class Arg2Eval>
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