math.hh

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
//
// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
// SPDX-License-Identifier: GPL-3.0-or-later
//
#ifndef DUMUX_MATH_HH
#define DUMUX_MATH_HH
 
#include <algorithm>
#include <numeric>
#include <cmath>
#include <utility>
 
#include <dune/common/typetraits.hh>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/dynmatrix.hh>
#include <dune/common/float_cmp.hh>
 
namespace Dumux {
 
template <class Scalar>
constexpr Scalar arithmeticMean(Scalar x, Scalar y, Scalar wx = 1.0, Scalar wy = 1.0) noexcept
{
    static_assert(Dune::IsNumber<Scalar>::value, "The arguments x, y, wx, and wy have to be numbers!");
 
    if (x*y <= 0)
        return 0;
    return (x * wx + y * wy)/(wx + wy);
}
 
template <class Scalar>
constexpr Scalar harmonicMean(Scalar x, Scalar y, Scalar wx = 1.0, Scalar wy = 1.0) noexcept
{
    static_assert(Dune::IsNumber<Scalar>::value, "The arguments x, y, wx, and wy have to be numbers!");
 
    if (x*y <= 0)
        return 0;
    return (wx + wy) * x * y / (wy * x + wx * y);
}
 
template <class Scalar>
Scalar geometricMean(Scalar x, Scalar y) noexcept
{
    static_assert(Dune::IsNumber<Scalar>::value, "The arguments x and y have to be numbers!");
 
    if (x*y <= 0)
        return 0;
    using std::sqrt;
    return sqrt(x*y)*sign(x);
}
 
template <class Scalar, int m, int n>
void harmonicMeanMatrix(Dune::FieldMatrix<Scalar, m, n> &K,
                        const Dune::FieldMatrix<Scalar, m, n> &Ki,
                        const Dune::FieldMatrix<Scalar, m, n> &Kj)
{
    for (int rowIdx=0; rowIdx < m; rowIdx++){
        for (int colIdx=0; colIdx< n; colIdx++){
            if (Dune::FloatCmp::ne<Scalar>(Ki[rowIdx][colIdx], Kj[rowIdx][colIdx])) {
                K[rowIdx][colIdx] =
                    harmonicMean(Ki[rowIdx][colIdx],
                                 Kj[rowIdx][colIdx]);
            }
            else
                K[rowIdx][colIdx] = Ki[rowIdx][colIdx];
        }
    }
}
 
template <class Scalar, class SolContainer>
int invertLinearPolynomial(SolContainer &sol,
                           Scalar a,
                           Scalar b)
{
    if (a == 0.0)
        return 0;
 
    sol[0] = -b/a;
    return 1;
}
 
template <class Scalar, class SolContainer>
int invertQuadraticPolynomial(SolContainer &sol,
                              Scalar a,
                              Scalar b,
                              Scalar c)
{
    // check for a line
    if (a == 0.0)
        return invertLinearPolynomial(sol, b, c);
 
    // discriminant
    Scalar Delta = b*b - 4*a*c;
    if (Delta < 0)
        return 0; // no real roots
 
    using std::sqrt;
    Delta = sqrt(Delta);
    sol[0] = (- b + Delta)/(2*a);
    sol[1] = (- b - Delta)/(2*a);
 
    // sort the result
    if (sol[0] > sol[1])
    {
        using std::swap;
        swap(sol[0], sol[1]);
    }
    return 2; // two real roots
}
 
template <class Scalar, class SolContainer>
void invertCubicPolynomialPostProcess_(SolContainer &sol,
                                       int numSol,
                                       Scalar a, Scalar b, Scalar c, Scalar d,
                                       std::size_t numIterations = 1)
{
    const auto eval = [&](auto x){ return d + x*(c + x*(b + x*a)); };
    const auto evalDeriv = [&](auto x){ return c + x*(2*b + x*3*a); };
 
    // do numIterations Newton iterations on the analytic solution
    // and update result if the precision is increased
    for (int i = 0; i < numSol; ++i)
    {
        Scalar x = sol[i];
        Scalar fCurrent = eval(x);
        const Scalar fOld = fCurrent;
        for (int j = 0; j < numIterations; ++j)
        {
            const Scalar fPrime = evalDeriv(x);
            if (fPrime == 0.0)
                break;
            x -= fCurrent/fPrime;
            fCurrent = eval(x);
        }
 
        using std::abs;
        if (abs(fCurrent) < abs(fOld))
            sol[i] = x;
    }
}
 
template <class Scalar, class SolContainer>
int invertCubicPolynomial(SolContainer *sol,
                          Scalar a, Scalar b, Scalar c, Scalar d,
                          std::size_t numPostProcessIterations = 1)
{
    // reduces to a quadratic polynomial
    if (a == 0)
        return invertQuadraticPolynomial(sol, b, c, d);
 
    // normalize the polynomial
    b /= a;
    c /= a;
    d /= a;
    a = 1;
 
    // get rid of the quadratic term by substituting x = t - b/3
    Scalar p = c - b*b/3;
    Scalar q = d + (2*b*b*b - 9*b*c)/27;
 
    // now we are at the form t^3 + p*t + q = 0. First we handle some
    // special cases to avoid divisions by zero later...
    if (p == 0.0 && q == 0.0) {
        // t^3 = 0, i.e. triple root at t = 0
        sol[0] = sol[1] = sol[2] = 0.0 - b/3;
        return 3;
    }
    else if (p == 0.0 && q != 0.0) {
        // t^3 + q = 0,
        //
        // i. e. single real root at t=curt(q)
        using std::cbrt;
        Scalar t = cbrt(q);
        sol[0] = t - b/3;
 
        return 1;
    }
    else if (p != 0.0 && q == 0.0) {
        // t^3 + p*t = 0 = t*(t^2 + p),
        //
        // i. e. roots at t = 0, t^2 + p = 0
        if (p > 0) {
            sol[0] = 0.0 - b/3;
            return 1; // only a single real root at t=0
        }
 
        // two additional real roots at t = sqrt(-p) and t = -sqrt(-p)
        using std::sqrt;
        sol[0] = -sqrt(-p) - b/3;
        sol[1] = 0.0 - b/3;
        sol[2] = sqrt(-p) - b/3;
 
        return 3;
    }
 
    // At this point
    //
    // t^3 + p*t + q = 0
    //
    // with p != 0 and q != 0 holds. Introducing the variables u and v
    // with the properties
    //
    //   u + v = t       and       3*u*v + p = 0
    //
    // leads to
    //
    // u^3 + v^3 + q = 0 .
    //
    // multiplying both sides with u^3 and taking advantage of the
    // fact that u*v = -p/3 leads to
    //
    // u^6 + q*u^3 - p^3/27 = 0
    //
    // Now, substituting u^3 = w yields
    //
    // w^2 + q*w - p^3/27 = 0
    //
    // This is a quadratic equation with the solutions
    //
    // w = -q/2 + sqrt(q^2/4 + p^3/27) and
    // w = -q/2 - sqrt(q^2/4 + p^3/27)
    //
    // Since w is equivalent to u^3 it is sufficient to only look at
    // one of the two cases. Then, there are still 2 cases: positive
    // and negative discriminant.
    Scalar wDisc = q*q/4 + p*p*p/27;
    if (wDisc >= 0) { // the positive discriminant case:
        // calculate the cube root of - q/2 + sqrt(q^2/4 + p^3/27)
        using std::cbrt;
        using std::sqrt;
 
        // Choose the root that is safe against the loss of precision
        // that can cause wDisc to be equal to q*q/4 despite p != 0.
        // We do not want u to be zero in that case. Mathematically,
        // we are happy with either root.
        const Scalar u = [&]{
            return q > 0 ? cbrt(-0.5*q - sqrt(wDisc)) : cbrt(-0.5*q + sqrt(wDisc));
        }();
        // at this point, u != 0 since p^3 = 0 is necessary in order
        // for u = 0 to hold, so
        sol[0] = u - p/(3*u) - b/3;
        // does not produce a division by zero. the remaining two
        // roots of u are rotated by +- 2/3*pi in the complex plane
        // and thus not considered here
        invertCubicPolynomialPostProcess_(sol, 1, a, b, c, d, numPostProcessIterations);
        return 1;
    }
    else { // the negative discriminant case:
        Scalar uCubedRe = - q/2;
        using std::sqrt;
        Scalar uCubedIm = sqrt(-wDisc);
        // calculate the cube root of - q/2 + sqrt(q^2/4 + p^3/27)
        using std::cbrt;
        Scalar uAbs  = cbrt(sqrt(uCubedRe*uCubedRe + uCubedIm*uCubedIm));
        using std::atan2;
        Scalar phi = atan2(uCubedIm, uCubedRe)/3;
 
        // calculate the length and the angle of the primitive root
 
        // with the definitions from above it follows that
        //
        // x = u - p/(3*u) - b/3
        //
        // where x and u are complex numbers. Rewritten in polar form
        // this is equivalent to
        //
        // x = |u|*e^(i*phi) - p*e^(-i*phi)/(3*|u|) - b/3 .
        //
        // Factoring out the e^ terms and subtracting the additional
        // terms, yields
        //
        // x = (e^(i*phi) + e^(-i*phi))*(|u| - p/(3*|u|)) - y - b/3
        //
        // with
        //
        // y = - |u|*e^(-i*phi) + p*e^(i*phi)/(3*|u|) .
        //
        // The crucial observation is the fact that y is the conjugate
        // of - x + b/3. This means that after taking advantage of the
        // relation
        //
        // e^(i*phi) + e^(-i*phi) = 2*cos(phi)
        //
        // the equation
        //
        // x = 2*cos(phi)*(|u| - p / (3*|u|)) - conj(x) - 2*b/3
        //
        // holds. Since |u|, p, b and cos(phi) are real numbers, it
        // follows that Im(x) = - Im(x) and thus Im(x) = 0. This
        // implies
        //
        // Re(x) = x = cos(phi)*(|u| - p / (3*|u|)) - b/3 .
        //
        // Considering the fact that u is a cubic root, we have three
        // values for phi which differ by 2/3*pi. This allows to
        // calculate the three real roots of the polynomial:
        for (int i = 0; i < 3; ++i) {
            using std::cos;
            sol[i] = cos(phi)*(uAbs - p/(3*uAbs)) - b/3;
            phi += 2*M_PI/3;
        }
 
        // post process the obtained solution to increase numerical
        // precision
        invertCubicPolynomialPostProcess_(sol, 3, a, b, c, d, numPostProcessIterations);
 
        // sort the result
        using std::sort;
        sort(sol, sol + 3);
 
        return 3;
    }
 
    // NOT REACHABLE!
    return 0;
}
 
template <class Scalar, int dim>
bool isLarger(const Dune::FieldVector<Scalar, dim> &pos,
              const Dune::FieldVector<Scalar, dim> &smallerVec)
{
    for (int i=0; i < dim; i++)
    {
        if (pos[i]<= smallerVec[i])
        {
            return false;
        }
    }
    return true;
}
 
template <class Scalar, int dim>
bool isSmaller(const Dune::FieldVector<Scalar, dim> &pos,
               const Dune::FieldVector<Scalar, dim> &largerVec)
{
    for (int i=0; i < dim; i++)
    {
        if (pos[i]>= largerVec[i])
        {
            return false;
        }
    }
    return true;
}
 
template <class Scalar, int dim>
bool isBetween(const Dune::FieldVector<Scalar, dim> &pos,
              const Dune::FieldVector<Scalar, dim> &smallerVec,
              const Dune::FieldVector<Scalar, dim> &largerVec)
{
   if (isLarger(pos, smallerVec) && isSmaller(pos, largerVec))
       {
           return true;
       }
   else
       return false;
}
 
namespace InterpolationPolicy { struct Linear; }
 
template <class Policy = InterpolationPolicy::Linear, class Scalar, class ... Parameter>
Scalar interpolate(Scalar ip, Parameter&& ... params)
{ return Policy::interpolate(ip, std::forward<Parameter>(params) ...); }
 
namespace InterpolationPolicy  {
 
struct Linear
{
    template<class Scalar>
    static constexpr Scalar interpolate(Scalar ip, const std::array<Scalar, 2>& params)
    {
        return params[0]*(1.0 - ip) + params[1]*ip;
    }
};
 
struct LinearTable
{
    template<class Scalar, class RandomAccessContainer0, class RandomAccessContainer1>
    static constexpr Scalar interpolate(Scalar ip, const RandomAccessContainer0& range, const RandomAccessContainer1& values)
    {
        // check bounds
        if (ip > range.back()) return values.back();
        if (ip < range[0]) return values[0];
 
        // if we are within bounds find the index of the lower bound
        const auto lookUpIndex = std::distance(range.begin(), std::lower_bound(range.begin(), range.end(), ip));
        if (lookUpIndex == 0)
            return values[0];
 
        const auto ipLinear = (ip - range[lookUpIndex-1])/(range[lookUpIndex] - range[lookUpIndex-1]);
        return Dumux::interpolate<Linear>(ipLinear, std::array<Scalar, 2>{{values[lookUpIndex-1], values[lookUpIndex]}});
    }
 
    template<class Scalar, class RandomAccessContainer>
    static constexpr Scalar interpolate(Scalar ip, const std::pair<RandomAccessContainer, RandomAccessContainer>& table)
    {
        const auto& [range, values] = table;
        return interpolate(ip, range, values);
    }
};
 
} // end namespace InterpolationPolicy
 
template <class Scalar>
std::vector<Scalar> linspace(const Scalar begin, const Scalar end,
                             std::size_t samples,
                             bool endPoint = true)
{
    using std::max;
    samples = max(std::size_t{2}, samples); // only makes sense for 2 or more samples
    const Scalar divisor = endPoint ? samples-1 : samples;
    const Scalar delta = (end-begin)/divisor;
    std::vector<Scalar> vec(samples);
    for (std::size_t i = 0; i < samples; ++i)
        vec[i] = begin + i*delta;
    return vec;
}
 
 
template <class Scalar>
Scalar antoine(Scalar temperature,
               Scalar A,
               Scalar B,
               Scalar C)
{
    const Scalar ln10 = 2.3025850929940459;
    using std::exp;
    return exp(ln10*(A - B/(C + temperature)));
}
 
template<class ValueType>
constexpr int sign(const ValueType& value) noexcept
{
    return (ValueType(0) < value) - (value < ValueType(0));
}
 
template <class Scalar>
Dune::FieldVector<Scalar, 3> crossProduct(const Dune::FieldVector<Scalar, 3> &vec1,
                                          const Dune::FieldVector<Scalar, 3> &vec2)
{
    return {vec1[1]*vec2[2]-vec1[2]*vec2[1],
            vec1[2]*vec2[0]-vec1[0]*vec2[2],
            vec1[0]*vec2[1]-vec1[1]*vec2[0]};
}
 
template <class Scalar>
Scalar crossProduct(const Dune::FieldVector<Scalar, 2> &vec1,
                    const Dune::FieldVector<Scalar, 2> &vec2)
{   return vec1[0]*vec2[1]-vec1[1]*vec2[0]; }
 
template <class Scalar>
Scalar tripleProduct(const Dune::FieldVector<Scalar, 3> &vec1,
                     const Dune::FieldVector<Scalar, 3> &vec2,
                     const Dune::FieldVector<Scalar, 3> &vec3)
{   return crossProduct<Scalar>(vec1, vec2)*vec3; }
 
template <class Scalar, int m, int n>
Dune::FieldMatrix<Scalar, n, m> getTransposed(const Dune::FieldMatrix<Scalar, m, n>& M)
{
    Dune::FieldMatrix<Scalar, n, m> T;
    for (std::size_t i = 0; i < m; ++i)
        for (std::size_t j = 0; j < n; ++j)
            T[j][i] = M[i][j];
 
    return T;
}
 
template <class Scalar>
Dune::DynamicMatrix<Scalar> getTransposed(const Dune::DynamicMatrix<Scalar>& M)
{
    std::size_t rows_T = M.M();
    std::size_t cols_T = M.N();
 
    Dune::DynamicMatrix<Scalar> M_T(rows_T, cols_T, 0.0);
 
    for (std::size_t i = 0; i < rows_T; ++i)
        for (std::size_t j = 0; j < cols_T; ++j)
            M_T[i][j] = M[j][i];
 
    return M_T;
}
 
template <class Scalar>
Dune::DynamicMatrix<Scalar> multiplyMatrices(const Dune::DynamicMatrix<Scalar> &M1,
                                             const Dune::DynamicMatrix<Scalar> &M2)
{
    using size_type = typename Dune::DynamicMatrix<Scalar>::size_type;
    const size_type rows = M1.N();
    const size_type cols = M2.M();
 
    DUNE_ASSERT_BOUNDS(M1.M() == M2.N());
 
    Dune::DynamicMatrix<Scalar> result(rows, cols, 0.0);
    for (size_type i = 0; i < rows; i++)
        for (size_type j = 0; j < cols; j++)
            for (size_type k = 0; k < M1.M(); k++)
                result[i][j] += M1[i][k]*M2[k][j];
 
    return result;
}
 
template <class Scalar, int rows1, int cols1, int cols2>
Dune::FieldMatrix<Scalar, rows1, cols2> multiplyMatrices(const Dune::FieldMatrix<Scalar, rows1, cols1> &M1,
                                                         const Dune::FieldMatrix<Scalar, cols1, cols2> &M2)
{
    using size_type = typename Dune::FieldMatrix<Scalar, rows1, cols2>::size_type;
 
    Dune::FieldMatrix<Scalar, rows1, cols2> result(0.0);
    for (size_type i = 0; i < rows1; i++)
        for (size_type j = 0; j < cols2; j++)
            for (size_type k = 0; k < cols1; k++)
                result[i][j] += M1[i][k]*M2[k][j];
 
    return result;
}
 
 
template <class MatrixType>
typename Dune::DenseMatrix<MatrixType>::field_type
trace(const Dune::DenseMatrix<MatrixType>& M)
{
    const auto rows = M.N();
    DUNE_ASSERT_BOUNDS(rows == M.M()); // rows == cols
 
    using MatType = Dune::DenseMatrix<MatrixType>;
    typename MatType::field_type trace = 0.0;
 
    for (typename MatType::size_type i = 0; i < rows; ++i)
        trace += M[i][i];
 
    return trace;
}
 
template <class MAT, class V>
typename Dune::DenseVector<V>::derived_type
mv(const Dune::DenseMatrix<MAT>& M,
   const Dune::DenseVector<V>& v)
{
    typename Dune::DenseVector<V>::derived_type res(v);
    M.mv(v, res);
    return res;
}
 
template <class FieldScalar, class V>
typename std::enable_if_t<Dune::IsNumber<FieldScalar>::value,
                          typename Dune::DenseVector<V>::derived_type>
mv(const FieldScalar m, const Dune::DenseVector<V>& v)
{
    typename Dune::DenseVector<V>::derived_type res(v);
    res *= m;
    return res;
}
 
template <class V1, class MAT, class V2>
typename Dune::DenseMatrix<MAT>::value_type
vtmv(const Dune::DenseVector<V1>& v1,
     const Dune::DenseMatrix<MAT>& M,
     const Dune::DenseVector<V2>& v2)
{
    return v1*mv(M, v2);
}
 
template <class V1, class FieldScalar, class V2>
typename std::enable_if_t<Dune::IsNumber<FieldScalar>::value, FieldScalar>
vtmv(const Dune::DenseVector<V1>& v1,
     const FieldScalar m,
     const Dune::DenseVector<V2>& v2)
{
    return m*(v1*v2);
}
 
template <class Scalar>
std::array<Scalar, 2> linearRegression(const std::vector<Scalar>& x,
                                       const std::vector<Scalar>& y)
{
    if (x.size() != y.size())
        DUNE_THROW(Dune::InvalidStateException, "x and y array must have the same length.");
 
    const Scalar averageX = std::accumulate(x.begin(), x.end(), 0.0)/x.size();
    const Scalar averageY = std::accumulate(y.begin(), y.end(), 0.0)/y.size();
 
    // calculate temporary variables necessary for slope computation
    const Scalar numerator = std::inner_product(
        x.begin(), x.end(), y.begin(), 0.0, std::plus<Scalar>(),
        [&](auto xx, auto yy) { return (xx - averageX) * (yy - averageY); }
    );
    const Scalar denominator = std::inner_product(
        x.begin(), x.end(), x.begin(), 0.0, std::plus<Scalar>(),
        [&](auto xx, auto yy) { return (xx - averageX) * (yy - averageX); }
    );
 
    // compute slope and intercept of the regression line
    const Scalar slope = numerator / denominator;
    const Scalar intercept = averageY - slope * averageX;
 
    return {intercept, slope};
}
 
} // end namespace Dumux
 
#endif