fvbaselocalresidual.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

  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_LOCAL_RESIDUAL_HH
#define EWOMS_FV_BASE_LOCAL_RESIDUAL_HH

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

#include <dune/istl/bvector.hh>
#include <dune/grid/common/geometry.hh>

#include <dune/common/fvector.hh>

#include <opm/material/common/ClassName.hpp>

#include "fvbaseproperties.hh"

#include <cmath>

namespace Ewoms {
template<class TypeTag>
class FvBaseLocalResidual
{
private:
    typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;

    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
    typedef typename GridView::template Codim<0>::Entity Element;

    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
    typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
    typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector;
    typedef typename GET_PROP_TYPE(TypeTag, Constraints) Constraints;
    typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
    typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
    enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };

    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
    typedef typename GET_PROP_TYPE(TypeTag, BoundaryContext) BoundaryContext;
    typedef typename GET_PROP_TYPE(TypeTag, ConstraintsContext) ConstraintsContext;

    typedef Opm::MathToolbox<Evaluation> Toolbox;
    typedef Dune::FieldVector<Evaluation, numEq> EvalEqVector;
    typedef Dune::FieldVector<Evaluation, numEq> VectorBlock;
    typedef Dune::BlockVector<VectorBlock> LocalBlockVector;

    // copying the local residual class is not a good idea
    FvBaseLocalResidual(const FvBaseLocalResidual &)
    {}

public:
    FvBaseLocalResidual()
    { }

    ~FvBaseLocalResidual()
    { }

    static void registerParameters()
    { }

    const LocalBlockVector &residual() const
    { return internalResidual_; }

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

    const LocalBlockVector &storageTerm() const
    { return internalStorageTerm_; }

    const VectorBlock &storageTerm(int dofIdx) const
    { return internalStorageTerm_[dofIdx]; }

    void eval(const Problem &problem, const Element &element)
    {
        ElementContext elemCtx(problem);
        elemCtx.updateAll(element);
        eval(elemCtx);
    }

    void eval(const ElementContext &elemCtx)
    {
        int numDof = elemCtx.numDof(/*timeIdx=*/0);
        int numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
        internalResidual_.resize(numDof);
        internalStorageTerm_.resize(numPrimaryDof);
        asImp_().eval(internalResidual_, internalStorageTerm_, elemCtx);
    }

    void eval(LocalBlockVector &residual,
              LocalBlockVector &storage,
              const ElementContext &elemCtx) const
    {
        int numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
        int numDof = elemCtx.numDof(/*timeIdx=*/0);

        typedef Opm::MathToolbox<Evaluation> Toolbox;

        residual = Toolbox::createConstant(0.0);
        storage = Toolbox::createConstant(0.0);

        // evaluate the flux terms
        asImp_().evalFluxes(residual, elemCtx, /*timeIdx=*/0);
        assert(static_cast<int>(residual.size()) == numDof);
        assert(static_cast<int>(storage.size()) == numPrimaryDof);

#if !defined NDEBUG
        for (int i=0; i < numDof; i++) {
            for (int j = 0; j < numEq; ++ j) {
                assert(std::isfinite(Toolbox::value(residual[i][j])));
                Valgrind::CheckDefined(residual[i][j]);
            }
        }
#endif

        // evaluate the storage and the source terms
        asImp_().evalVolumeTerms_(residual, storage, elemCtx);
        for (int dofIdx=0; dofIdx < numPrimaryDof; dofIdx++) {
            for (int eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
                storage[dofIdx][eqIdx] /= elemCtx.dofTotalVolume(dofIdx, /*timeIdx=*/0);

                assert(std::isfinite(Toolbox::value(residual[dofIdx][eqIdx])));
                Valgrind::CheckDefined(residual[dofIdx][eqIdx]);
            }
        }

        // evaluate the boundary conditions
        asImp_().evalBoundary_(residual, elemCtx, /*timeIdx=*/0);
#if !defined NDEBUG
        for (int i=0; i < numDof; i++) {
            for (int j = 0; j < numEq; ++ j) {
                assert(std::isfinite(Toolbox::value(residual[i][j])));
                Valgrind::CheckDefined(residual[i][j]);
            }
        }
#endif

        // evaluate the constraint DOFs
        asImp_().evalConstraints_(residual, storage, elemCtx, /*timeIdx=*/0);

        // make the residual volume specific (i.e., make it incorrect mass per cubic
        // meter instead of total mass)
        for (int dofIdx=0; dofIdx < numDof; ++dofIdx) {
            if (elemCtx.dofTotalVolume(dofIdx, /*timeIdx=*/0) > 0) {
                // interior DOF
                Scalar dofVolume = elemCtx.dofTotalVolume(dofIdx, /*timeIdx=*/0);

                for (int eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
                    residual[dofIdx][eqIdx] /= dofVolume;
                    assert(std::isfinite(Toolbox::value(residual[dofIdx][eqIdx])));
                    Valgrind::CheckDefined(residual[dofIdx][eqIdx]);
                }
            }
        }
    }

    void evalStorage(const ElementContext &elemCtx, int timeIdx)
    {
        int numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
        internalStorageTerm_.resize(numPrimaryDof);
        evalStorage(internalStorageTerm_,
                    elemCtx,
                    timeIdx);
    }

    void evalStorage(LocalBlockVector &storage,
                     const ElementContext &elemCtx,
                     int timeIdx) const
    {
        if (timeIdx == 0) {
            // for the most current solution, the storage term depends on the current
            // primary variables.

            // calculate the amount of conservation each quantity inside
            // all primary sub control volumes
            int numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
            for (int dofIdx=0; dofIdx < numPrimaryDof; dofIdx++) {
                storage[dofIdx] = Toolbox::createConstant(0.0);
                asImp_().computeStorage(storage[dofIdx], elemCtx, dofIdx, timeIdx);

                for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
                    storage[dofIdx][eqIdx] *=
                        elemCtx.stencil(timeIdx).subControlVolume(dofIdx).volume()
                        * elemCtx.intensiveQuantities(dofIdx, timeIdx).extrusionFactor();
                }
            }
        }
        else {
            // for all previous solutions, the storage term does _not_ depend on the
            // current primary variables, so we use scalars to store it.

            // calculate the amount of conservation each quantity inside
            // all primary sub control volumes
            Dune::FieldVector<Scalar, numEq> tmp;
            int numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
            for (int dofIdx=0; dofIdx < numPrimaryDof; dofIdx++) {
                tmp = 0;
                asImp_().computeStorage(tmp,
                                        elemCtx,
                                        dofIdx,
                                        timeIdx);
                tmp *=
                    elemCtx.stencil(timeIdx).subControlVolume(dofIdx).volume()
                    * elemCtx.intensiveQuantities(dofIdx, timeIdx).extrusionFactor();


                for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
                    storage[dofIdx][eqIdx] = Toolbox::createConstant(tmp[eqIdx]);
            }
        }
    }

    void evalFluxes(LocalBlockVector &residual,
                    const ElementContext &elemCtx,
                    int timeIdx) const
    {
        RateVector flux;

        const auto &stencil = elemCtx.stencil(timeIdx);
        // calculate the mass flux over the sub-control volume faces
        int numInteriorFaces = elemCtx.numInteriorFaces(timeIdx);
        for (int scvfIdx = 0; scvfIdx < numInteriorFaces; scvfIdx++) {
            const auto &face = stencil.interiorFace(scvfIdx);
            int i = face.interiorIndex();
            int j = face.exteriorIndex();

            Valgrind::SetUndefined(flux);
            asImp_().computeFlux(flux, /*context=*/elemCtx, scvfIdx, timeIdx);

            Scalar alpha = elemCtx.extensiveQuantities(scvfIdx, timeIdx).extrusionFactor();
            alpha *= face.area();
            for (int eqIdx = 0; eqIdx < numEq; ++ eqIdx)
                flux[eqIdx] *= alpha;
            Valgrind::CheckDefined(flux);

            // The balance equation for a finite volume is given by
            //
            // dStorage/dt + Flux = Source
            //
            // where the 'Flux' and the 'Source' terms represent the
            // mass per second which leaves the finite
            // volume. Re-arranging this, we get
            //
            // dStorage/dt + Flux - Source = 0
            //
            // Since the mass flux as calculated by computeFlux() goes out of sub-control
            // volume i and into sub-control volume j, we need to add the flux to finite
            // volume i and subtract it from finite volume j
            residual[i] += flux;
            residual[j] -= flux;
        }
    }

    // The following methods _must_ be overloaded by the actual flow
    // models!

    void computeStorage(EqVector &storage,
                        const ElementContext &elemCtx,
                        int dofIdx,
                        int timeIdx) const
    {
        OPM_THROW(std::logic_error,
                   "Not implemented: The local residual " << Opm::className<Implementation>()
                   << " does not implement the required method 'computeStorage()'");
    }

    void computeFlux(RateVector &flux,
                     const ElementContext &elemCtx,
                     int scvfIdx,
                     int timeIdx) const
    {
        OPM_THROW(std::logic_error,
                  "Not implemented: The local residual " << Opm::className<Implementation>()
                  << " does not implement the required method 'computeFlux()'");
    }

    void computeSource(RateVector &source,
                       const ElementContext &elemCtx,
                       int dofIdx,
                       int timeIdx) const
    {
        OPM_THROW(std::logic_error,
                  "Not implemented: The local residual " << Opm::className<Implementation>()
                  << " does not implement the required method 'computeSource()'");
    }

protected:
    void evalBoundary_(LocalBlockVector &residual,
                       const ElementContext &elemCtx,
                       int timeIdx) const
    {
        if (!elemCtx.onBoundary())
            return;

        BoundaryContext boundaryCtx(elemCtx);

        // evaluate the boundary for all boundary faces of the current context
        int numBoundaryFaces = boundaryCtx.numBoundaryFaces(/*timeIdx=*/0);
        for (int faceIdx = 0; faceIdx < numBoundaryFaces; ++faceIdx) {
            // add the residual of all vertices of the boundary
            // segment
            evalBoundarySegment_(residual,
                                 boundaryCtx,
                                 faceIdx,
                                 timeIdx);
        }
    }

    void evalBoundarySegment_(LocalBlockVector &residual,
                              const BoundaryContext &boundaryCtx,
                              int boundaryFaceIdx,
                              int timeIdx) const
    {
        BoundaryRateVector values;

        Valgrind::SetUndefined(values);
        boundaryCtx.problem().boundary(values,
                                       boundaryCtx,
                                       boundaryFaceIdx,
                                       timeIdx);
        Valgrind::CheckDefined(values);

        const auto &stencil = boundaryCtx.stencil(timeIdx);
        int dofIdx = stencil.boundaryFace(boundaryFaceIdx).interiorIndex();
        const auto &insideIntQuants = boundaryCtx.elementContext().intensiveQuantities(dofIdx, timeIdx);
        for (unsigned eqIdx = 0; eqIdx < values.size(); ++eqIdx)
            values[eqIdx] *=
                stencil.boundaryFace(boundaryFaceIdx).area()
                * insideIntQuants.extrusionFactor();

        for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
            residual[dofIdx][eqIdx] += values[eqIdx];
        }
    }

    void evalConstraints_(LocalBlockVector &residual,
                          LocalBlockVector &storage,
                          const ElementContext &elemCtx,
                          int timeIdx) const
    {
        if (!GET_PROP_VALUE(TypeTag, EnableConstraints))
            return;

        const auto &problem = elemCtx.problem();
        Constraints constraints;
        ConstraintsContext constraintsCtx(elemCtx);
        int numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
        for (int dofIdx = 0; dofIdx < numPrimaryDof; ++dofIdx) {
            // ask the problem for the constraint values
            constraints.reset();
            problem.constraints(constraints, elemCtx, dofIdx, timeIdx);

            if (!constraints.isConstraint())
                continue;

            // enforce the constraints
            const PrimaryVariables &priVars = elemCtx.primaryVars(dofIdx, timeIdx);
            for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
                if (!constraints.isConstraint(eqIdx))
                    continue;

                int pvIdx = constraints.eqToPvIndex(eqIdx);

                assert(0 <= pvIdx && pvIdx < numEq);
                Valgrind::CheckDefined(constraints[pvIdx]);

                residual[dofIdx][eqIdx] = priVars.makeEvaluation(pvIdx, timeIdx) - constraints[pvIdx];
                storage[dofIdx][eqIdx] = 0.0;
            }
        }
    }

    void evalVolumeTerms_(LocalBlockVector &residual,
                          LocalBlockVector &storage,
                          const ElementContext &elemCtx) const
    {
        EvalEqVector tmp;
        EqVector tmp2;
        RateVector sourceRate;

        tmp = Toolbox::createConstant(0.0);
        tmp2 = 0.0;

        // evaluate the volumetric terms (storage + source terms)
        int numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
        for (int dofIdx=0; dofIdx < numPrimaryDof; dofIdx++) {
            Scalar extrusionFactor =
                elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0).extrusionFactor();
            Scalar scvVolume =
               elemCtx.stencil(/*timeIdx=*/0).subControlVolume(dofIdx).volume() * extrusionFactor;
            Valgrind::CheckDefined(scvVolume);

            // mass balance within the element. this is the
            // \f$\frac{m}{\partial t}\f$ term if using implicit
            // euler as time discretization.
            //
            // TODO (?): we might need a more explicit way for
            // doing the time discretization...
            asImp_().computeStorage(tmp,
                                    elemCtx,
                                    dofIdx,
                                    /*timeIdx=*/0);
            Valgrind::CheckDefined(tmp);

            asImp_().computeStorage(tmp2,
                                    elemCtx,
                                    dofIdx,
                                    /*timeIdx=*/1);
            Valgrind::CheckDefined(tmp2);

            tmp -= tmp2;
            for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
                tmp[eqIdx] *= scvVolume / elemCtx.simulator().timeStepSize();

            storage[dofIdx] += tmp;
            residual[dofIdx] += tmp;

            Valgrind::CheckDefined(storage[dofIdx]);
            Valgrind::CheckDefined(residual[dofIdx]);

            // subtract the source term from the residual
            asImp_().computeSource(sourceRate, elemCtx, dofIdx, /*timeIdx=*/0);
            for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
                sourceRate[eqIdx] *= scvVolume;
            residual[dofIdx] -= sourceRate;

            // make sure that only defined quantities were used
            // to calculate the residual.
            Valgrind::CheckDefined(storage[dofIdx]);
            Valgrind::CheckDefined(residual[dofIdx]);
        }
    }


private:
    Implementation &asImp_()
    {
      assert(static_cast<Implementation*>(this) != 0);
      return *static_cast<Implementation*>(this);
    }

    const Implementation &asImp_() const
    {
      assert(static_cast<const Implementation*>(this) != 0);
      return *static_cast<const Implementation*>(this);
    }

    LocalBlockVector internalResidual_;
    LocalBlockVector internalStorageTerm_;
};

} // namespace Ewoms

#endif