numericepsilon.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-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
// SPDX-License-Identifier: GPL-3.0-or-later
//
#ifndef DUMUX_ASSEMBLY_NUMERIC_EPSILON_HH
#define DUMUX_ASSEMBLY_NUMERIC_EPSILON_HH
 
#include <cmath>
 
#include <dune/common/fvector.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/numericdifferentiation.hh>
 
namespace Dumux {

template<class Scalar, int numEq>
class NumericEpsilon
{
    using Magnitude = std::decay_t<decltype(abs(std::declval<Scalar>()))>;
    using NumEqVector = Dune::FieldVector<Magnitude, numEq>;
 
public:
    explicit NumericEpsilon(const std::string& paramGroup = "")
    {
        // get epsilons from input file with invalid default
        baseEps_ = getParamFromGroup<Scalar>(paramGroup, "Assembly.NumericDifference.BaseEpsilon", 1e-10);
        magnitude_ = getParamFromGroup<NumEqVector>(paramGroup, "Assembly.NumericDifference.PriVarMagnitude", NumEqVector(-1));
    }

    Scalar operator() (Scalar priVar, int priVarIdx) const noexcept
    {
        return magnitude_[priVarIdx] > 0.0 ? Scalar(baseEps_*magnitude_[priVarIdx])
                                           : Scalar(NumericDifferentiation::epsilon(priVar, baseEps_));
    }
 
private:
    Scalar baseEps_;
    NumEqVector magnitude_;
};
 
} // end namespace Dumux
 
#endif