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:
/*
  Copyright (C) 2008-2013 by Andreas Lauser
  Copyright (C) 2011-2012 by Bernd Flemisch
  Copyright (C) 2011 by Klaus Mosthaf

  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/>.
*/
#ifndef EWOMS_FV_BASE_FD_LOCAL_LINEARIZER_HH
#define EWOMS_FV_BASE_FD_LOCAL_LINEARIZER_HH

#include <ewoms/common/propertysystem.hh>
#include <ewoms/common/parametersystem.hh>

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

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

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

namespace Ewoms {
// forward declaration
template<class TypeTag>
class FvBaseFdLocalLinearizer;

namespace Properties {
// declare the property tags required for the finite differences local linearizer
NEW_TYPE_TAG(FiniteDifferenceLocalLinearizer);

NEW_PROP_TAG(LocalLinearizer);
NEW_PROP_TAG(Evaluation);
NEW_PROP_TAG(NumericDifferenceMethod);
NEW_PROP_TAG(BaseEpsilon);

NEW_PROP_TAG(LocalResidual);
NEW_PROP_TAG(Simulator);
NEW_PROP_TAG(Problem);
NEW_PROP_TAG(Model);
NEW_PROP_TAG(PrimaryVariables);
NEW_PROP_TAG(ElementContext);
NEW_PROP_TAG(Scalar);
NEW_PROP_TAG(Evaluation);
NEW_PROP_TAG(GridView);
NEW_PROP_TAG(NumEq);

// set the properties to be spliced in
SET_TYPE_PROP(FiniteDifferenceLocalLinearizer, LocalLinearizer,
              Ewoms::FvBaseFdLocalLinearizer<TypeTag>);

SET_TYPE_PROP(FiniteDifferenceLocalLinearizer, Evaluation,
              typename GET_PROP_TYPE(TypeTag, Scalar));

SET_INT_PROP(FiniteDifferenceLocalLinearizer, NumericDifferenceMethod, +1);

SET_SCALAR_PROP(FiniteDifferenceLocalLinearizer, BaseEpsilon, 0.9123e-10);
}

template<class TypeTag>
class FvBaseFdLocalLinearizer
{
private:
    typedef typename GET_PROP_TYPE(TypeTag, LocalLinearizer) Implementation;
    typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) LocalResidual;
    typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
    typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
    typedef typename GridView::template Codim<0>::Entity Element;

    enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
    typedef Dune::FieldMatrix<Scalar, numEq, numEq> MatrixBlock;
    typedef Dune::Matrix<MatrixBlock> LocalBlockMatrix;

    typedef Dune::FieldVector<Scalar, numEq> VectorBlock;
    typedef Dune::BlockVector<VectorBlock> LocalBlockVector;
    typedef Dune::BlockVector<MatrixBlock> LocalStorageMatrix;

#if __GNUC__ == 4 && __GNUC_MINOR__ <= 6
public:
    // make older GCCs happy by providing a public copy constructor (this is necessary
    // for their implementation of std::vector, although the method is never called...)
    FvBaseFdLocalLinearizer(const FvBaseFdLocalLinearizer&)
        : internalElemContext_(0)
    {}

#else
    // copying local residual objects around is a very bad idea, so we explicitly prevent
    // it...
    FvBaseFdLocalLinearizer(const FvBaseFdLocalLinearizer&) = delete;
#endif
public:
    FvBaseFdLocalLinearizer()
        : internalElemContext_(0)
    { }

    ~FvBaseFdLocalLinearizer()
    { delete internalElemContext_; }

    static void registerParameters()
    {
        EWOMS_REGISTER_PARAM(TypeTag, int, NumericDifferenceMethod,
                             "The method used for numeric differentiation (-1: backward "
                             "differences, 0: central differences, 1: forward differences)");
    }

    void init(Simulator &simulator)
    {
        simulatorPtr_ = &simulator;
        delete internalElemContext_;
        internalElemContext_ = new ElementContext(simulator);
    }

    void linearize(const Element &element)
    {
        internalElemContext_->updateAll(element);

        linearize(*internalElemContext_);
    }

    void linearize(ElementContext &elemCtx)
    {
        // update the weights of the primary variables for the context
        model_().updatePVWeights(elemCtx);

        resize_(elemCtx);
        reset_(elemCtx);

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

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

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

    static Scalar baseEpsilon()
    { return GET_PROP_VALUE(TypeTag, BaseEpsilon); }

    Scalar numericEpsilon(const ElementContext &elemCtx,
                          int dofIdx,
                          int pvIdx) const
    {
        int globalIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
        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 MatrixBlock &jacobian(int domainScvIdx, int rangeScvIdx) const
    { return jacobian_[domainScvIdx][rangeScvIdx]; }

    const MatrixBlock &jacobianStorage(int dofIdx) const
    { return jacobianStorage_[dofIdx]; }

    const VectorBlock &residual(int dofIdx) const
    { return residual_[dofIdx]; }

    const VectorBlock &residualStorage(int dofIdx) const
    { return residualStorage_[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_()
    { return EWOMS_GET_PARAM(TypeTag, int, NumericDifferenceMethod); }

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

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

        residual_.resize(numDof);
        residualStorage_.resize(numPrimaryDof);

        derivResidual_.resize(numDof);
        derivStorage_.resize(numPrimaryDof);
    }

    void reset_(const ElementContext &elemCtx)
    {
        int numDof = elemCtx.numDof(/*timeIdx=*/0);
        int numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
        for (int primaryDofIdx = 0; primaryDofIdx < numPrimaryDof; ++ primaryDofIdx) {
            residualStorage_[primaryDofIdx] = 0.0;

            jacobianStorage_[primaryDofIdx] = 0.0;
            for (int dof2Idx = 0; dof2Idx < numDof; ++ dof2Idx) {
                jacobian_[dof2Idx][primaryDofIdx] = 0.0;
            }
        }

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

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

        PrimaryVariables priVars(elemCtx.primaryVars(dofIdx, /*timeIdx=*/0));
        Scalar eps = asImp_().numericEpsilon(elemCtx, dofIdx, pvIdx);
        Scalar delta = 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 residual
            elemCtx.updateIntensiveQuantities(priVars, dofIdx, /*timeIdx=*/0);
            elemCtx.updateAllExtensiveQuantities();
            localResidual_.eval(derivResidual_, derivStorage_, 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_;
            derivStorage_ = residualStorage_;
        }

        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 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();
            derivStorage_ -= localResidual_.storageTerm();
        }
        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_;
            derivStorage_ -= residualStorage_;
        }

        assert(delta > 0);

        // divide difference in residuals by the magnitude of the
        // deflections between the two function evaluation
        derivResidual_ /= delta;
        derivStorage_ /= 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 updateLocalLinearizer_(const ElementContext &elemCtx,
                              int primaryDofIdx,
                              int pvIdx)
    {
        // store the derivative of the storage term
        for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
            jacobianStorage_[primaryDofIdx][eqIdx][pvIdx] = derivStorage_[primaryDofIdx][eqIdx];

        int numDof = elemCtx.numDof(/*timeIdx=*/0);
        for (int dofIdx = 0; dofIdx < numDof; dofIdx++)
        {
            for (int eqIdx = 0; eqIdx < numEq; eqIdx++) {
                // A[dofIdx][primaryDofIdx][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
                // 'primaryDofIdx'
                jacobian_[dofIdx][primaryDofIdx][eqIdx][pvIdx] = derivResidual_[dofIdx][eqIdx];
                Valgrind::CheckDefined(jacobian_[dofIdx][primaryDofIdx][eqIdx][pvIdx]);
            }
        }
    }

    Simulator *simulatorPtr_;
    Model *modelPtr_;

    ElementContext *internalElemContext_;

    LocalBlockMatrix jacobian_;
    LocalStorageMatrix jacobianStorage_;

    LocalBlockVector residual_;
    LocalBlockVector residualStorage_;

    LocalBlockVector derivResidual_;
    LocalBlockVector derivStorage_;

    LocalResidual localResidual_;
};

} // namespace Ewoms

#endif