fvbasefdlocallinearizer.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:
/*
  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/>.

  Consult the COPYING file in the top-level source directory of this
  module for the precise wording of the license and the list of
  copyright holders.
*/
#ifndef EWOMS_FV_BASE_FD_LOCAL_LINEARIZER_HH
#define EWOMS_FV_BASE_FD_LOCAL_LINEARIZER_HH

#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>

#include <dune/istl/bvector.hh>
#include <dune/istl/matrix.hh>

#include <opm/material/common/MathToolbox.hpp>
#include <opm/material/common/Valgrind.hpp>

#include <opm/models/discretization/common/fvbaseproperties.hh>
#include <opm/models/utils/parametersystem.hpp>
#include <opm/models/utils/propertysystem.hh>

#include <opm/simulators/linalg/linalgproperties.hh>

#include <algorithm>
#include <cmath>
#include <cstddef>
#include <limits>
#include <memory>

namespace Opm {

// forward declaration
template<class TypeTag>
class FvBaseFdLocalLinearizer;

} // namespace Opm

namespace Opm::Properties {

// declare the property tags required for the finite differences local linearizer

namespace TTag {
struct FiniteDifferenceLocalLinearizer {};
} // namespace TTag

template<class TypeTag, class MyTypeTag>
struct BaseEpsilon { using type = UndefinedProperty; };

// set the properties to be spliced in
template<class TypeTag>
struct LocalLinearizer<TypeTag, TTag::FiniteDifferenceLocalLinearizer>
{ using type = FvBaseFdLocalLinearizer<TypeTag>; };

template<class TypeTag>
struct Evaluation<TypeTag, TTag::FiniteDifferenceLocalLinearizer>
{ using type = GetPropType<TypeTag, Properties::Scalar>; };

template<class TypeTag>
struct BaseEpsilon<TypeTag, TTag::FiniteDifferenceLocalLinearizer>
{
    using type = GetPropType<TypeTag, Properties::Scalar>;
    static constexpr type value = std::max<type>(0.9123e-10,
                                                 std::numeric_limits<type>::epsilon() * 1.23e3);
};

} // namespace Opm::Properties

namespace Opm::Parameters {

struct NumericDifferenceMethod { static constexpr int value = +1; };

} // namespace Opm::Parameters

namespace Opm {

template<class TypeTag>
class FvBaseFdLocalLinearizer
{
private:
    using Implementation = GetPropType<TypeTag, Properties::LocalLinearizer>;
    using LocalResidual = GetPropType<TypeTag, Properties::LocalResidual>;
    using Simulator = GetPropType<TypeTag, Properties::Simulator>;
    using Problem = GetPropType<TypeTag, Properties::Problem>;
    using Model = GetPropType<TypeTag, Properties::Model>;
    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using GridView = GetPropType<TypeTag, Properties::GridView>;
    using Element = typename GridView::template Codim<0>::Entity;

    enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };

    // extract local matrices from jacobian matrix for consistency
    using ScalarMatrixBlock = typename GetPropType<TypeTag, Properties::SparseMatrixAdapter>::MatrixBlock;
    using ScalarVectorBlock = Dune::FieldVector<Scalar, numEq>;

    using ScalarLocalBlockVector = Dune::BlockVector<ScalarVectorBlock>;
    using ScalarLocalBlockMatrix = Dune::Matrix<ScalarMatrixBlock>;

    using LocalEvalBlockVector = typename LocalResidual::LocalEvalBlockVector;

public:
    FvBaseFdLocalLinearizer() = default;

    // Disable copying for efficiency reasons
    FvBaseFdLocalLinearizer(const FvBaseFdLocalLinearizer&) = delete;

    static void registerParameters()
    {
        Parameters::Register<Parameters::NumericDifferenceMethod>
            ("The method used for numeric differentiation (-1: backward "
             "differences, 0: central differences, 1: forward differences)");
    }

    void init(Simulator& simulator)
    {
        simulatorPtr_ = &simulator;
        internalElemContext_ = std::make_unique<ElementContext>(simulator);
    }

    void linearize(const Element& element)
    {
        linearize(*internalElemContext_, element);
    }

    void linearize(ElementContext& elemCtx, const Element& elem)
    {
        elemCtx.updateAll(elem);

        // update the weights of the primary variables for the context
        model_().updatePVWeights(elemCtx);

        resize_(elemCtx);
        reset_(elemCtx);

        // calculate the local residual
        localResidual_.eval(residual_, elemCtx);

        // calculate the local jacobian matrix
        const std::size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
        for (auto dofIdx = 0 * numPrimaryDof; dofIdx < numPrimaryDof; ++dofIdx) {
            for (auto pvIdx = 0 * numEq; pvIdx < numEq; ++pvIdx) {
                asImp_().evalPartialDerivative_(elemCtx, dofIdx, pvIdx);

                // incorporate the partial derivatives into the local Jacobian matrix
                updateLocalJacobian_(elemCtx, dofIdx, pvIdx);
            }
        }
    }

    static Scalar baseEpsilon()
    { return getPropValue<TypeTag, Properties::BaseEpsilon>(); }

    Scalar numericEpsilon(const ElementContext& elemCtx,
                          unsigned dofIdx,
                          unsigned pvIdx) const
    {
        const unsigned globalIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
        const Scalar pvWeight = elemCtx.model().primaryVarWeight(globalIdx, pvIdx);
        assert(pvWeight > 0 && std::isfinite(pvWeight));
        Valgrind::CheckDefined(pvWeight);

        return baseEpsilon() / pvWeight;
    }

    LocalResidual& localResidual()
    { return localResidual_; }

    const LocalResidual& localResidual() const
    { return localResidual_; }

    const ScalarMatrixBlock& jacobian(unsigned domainScvIdx, unsigned rangeScvIdx) const
    { return jacobian_[domainScvIdx][rangeScvIdx]; }

    const ScalarVectorBlock& residual(unsigned dofIdx) const
    { return residual_[dofIdx]; }

protected:
    Implementation& asImp_()
    { return *static_cast<Implementation*>(this); }

    const Implementation& asImp_() const
    { return *static_cast<const Implementation*>(this); }

    const Simulator& simulator_() const
    { return *simulatorPtr_; }

    const Problem& problem_() const
    { return simulatorPtr_->problem(); }

    const Model& model_() const
    { return simulatorPtr_->model(); }

    static int numericDifferenceMethod_()
    {
        static int diff = Parameters::Get<Parameters::NumericDifferenceMethod>();
        return diff;
    }

    void resize_(const ElementContext& elemCtx)
    {
        const std::size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
        const std::size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);

        residual_.resize(numDof);
        jacobian_.setSize(numDof, numPrimaryDof);

        derivResidual_.resize(numDof);
    }

    void reset_(const ElementContext& elemCtx)
    {
        const std::size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
        const std::size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
        for (unsigned primaryDofIdx = 0; primaryDofIdx < numPrimaryDof; ++primaryDofIdx) {
            for (unsigned dof2Idx = 0; dof2Idx < numDof; ++dof2Idx) {
                jacobian_[dof2Idx][primaryDofIdx] = 0.0;
            }
        }

        for (unsigned primaryDofIdx = 0; primaryDofIdx < numDof; ++primaryDofIdx) {
            residual_[primaryDofIdx] = 0.0;
        }
    }

    void evalPartialDerivative_(ElementContext& elemCtx,
                                unsigned dofIdx,
                                unsigned pvIdx)
    {
        // save all quantities which depend on the specified primary
        // variable at the given sub control volume
        elemCtx.stashIntensiveQuantities(dofIdx);

        PrimaryVariables priVars(elemCtx.primaryVars(dofIdx, /*timeIdx=*/0));
        const Scalar eps = asImp_().numericEpsilon(elemCtx, dofIdx, pvIdx);
        Scalar delta = 0.0;

        if (numericDifferenceMethod_() >= 0) {
            // we are not using backward differences, i.e. we need to
            // calculate f(x + \epsilon)

            // deflect primary variables
            priVars[pvIdx] += eps;
            delta += eps;

            // calculate the deflected residual
            elemCtx.updateIntensiveQuantities(priVars, dofIdx, /*timeIdx=*/0);
            elemCtx.updateAllExtensiveQuantities();
            localResidual_.eval(derivResidual_, elemCtx);
        }
        else {
            // we are using backward differences, i.e. we don't need
            // to calculate f(x + \epsilon) and we can recycle the
            // (already calculated) residual f(x)
            derivResidual_ = residual_;
        }

        if (numericDifferenceMethod_() <= 0) {
            // we are not using forward differences, i.e. we don't
            // need to calculate f(x - \epsilon)

            // deflect the primary variables
            priVars[pvIdx] -= delta + eps;
            delta += eps;

            // calculate the deflected residual again, this time we use the local
            // residual's internal storage.
            elemCtx.updateIntensiveQuantities(priVars, dofIdx, /*timeIdx=*/0);
            elemCtx.updateAllExtensiveQuantities();
            localResidual_.eval(elemCtx);

            derivResidual_ -= localResidual_.residual();
        }
        else {
            // we are using forward differences, i.e. we don't need to
            // calculate f(x - \epsilon) and we can recycle the
            // (already calculated) residual f(x)
            derivResidual_ -= residual_;
        }

        assert(delta > 0);

        // divide difference in residuals by the magnitude of the
        // deflections between the two function evaluation
        derivResidual_ /= delta;

        // restore the original state of the element's volume
        // variables
        elemCtx.restoreIntensiveQuantities(dofIdx);

#ifndef NDEBUG
        for (unsigned i = 0; i < derivResidual_.size(); ++i) {
            Valgrind::CheckDefined(derivResidual_[i]);
        }
#endif
    }

    void updateLocalJacobian_(const ElementContext& elemCtx,
                              unsigned focusDofIdx,
                              unsigned pvIdx)
    {
        const std::size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
        for (unsigned dofIdx = 0; dofIdx < numDof; ++dofIdx) {
            for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
                // A[dofIdx][focusDofIdx][eqIdx][pvIdx] is the partial derivative of the
                // residual function 'eqIdx' for the degree of freedom 'dofIdx' with
                // regard to the primary variable 'pvIdx' of the degree of freedom
                // 'focusDofIdx'
                jacobian_[dofIdx][focusDofIdx][eqIdx][pvIdx] = derivResidual_[dofIdx][eqIdx];
                Valgrind::CheckDefined(jacobian_[dofIdx][focusDofIdx][eqIdx][pvIdx]);
            }
        }
    }

    Simulator* simulatorPtr_{};

    std::unique_ptr<ElementContext> internalElemContext_{};

    LocalEvalBlockVector residual_{};
    LocalEvalBlockVector derivResidual_{};
    ScalarLocalBlockMatrix jacobian_{};

    LocalResidual localResidual_{};
};

} // namespace Opm

#endif