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

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

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

#include <dune/istl/bvector.hh>

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

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

#include <cassert>
#include <cmath>
#include <cstddef>
#include <stdexcept>
#include <type_traits>

namespace Opm {

template<class TypeTag>
class FvBaseLocalResidual
{
private:
    using Implementation = GetPropType<TypeTag, Properties::LocalResidual>;

    using GridView = GetPropType<TypeTag, Properties::GridView>;
    using Element = typename GridView::template Codim<0>::Entity;

    using Problem = GetPropType<TypeTag, Properties::Problem>;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
    using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
    using RateVector = GetPropType<TypeTag, Properties::RateVector>;
    using EqVector = GetPropType<TypeTag, Properties::EqVector>;
    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
    using BoundaryContext = GetPropType<TypeTag, Properties::BoundaryContext>;

    static constexpr bool useVolumetricResidual = getPropValue<TypeTag, Properties::UseVolumetricResidual>();

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

    using Toolbox = MathToolbox<Evaluation>;
    using EvalVector = Dune::FieldVector<Evaluation, numEq>;

public:
    using LocalEvalBlockVector = Dune::BlockVector<EvalVector, aligned_allocator<EvalVector, alignof(EvalVector)>>;

    FvBaseLocalResidual() = default;

    // copying the local residual class is not a good idea
    FvBaseLocalResidual(const FvBaseLocalResidual&) = delete;

    static void registerParameters()
    {}

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

    const EvalVector& residual(unsigned dofIdx) const
    { return internalResidual_[dofIdx]; }

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

    void eval(ElementContext& elemCtx)
    {
        std::size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
        internalResidual_.resize(numDof);
        asImp_().eval(internalResidual_, elemCtx);
    }

    void eval(LocalEvalBlockVector& residual,
              ElementContext& elemCtx) const
    {
        assert(residual.size() == elemCtx.numDof(/*timeIdx=*/0));

        residual = 0.0;

        // evaluate the flux terms
        asImp_().evalFluxes(residual, elemCtx, /*timeIdx=*/0);

        // evaluate the storage and the source terms
        asImp_().evalVolumeTerms_(residual, elemCtx);

        // evaluate the boundary conditions
        asImp_().evalBoundary_(residual, elemCtx, /*timeIdx=*/0);

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

                    assert(std::isfinite(dofVolume));
                    Valgrind::CheckDefined(dofVolume);

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

    void evalStorage(LocalEvalBlockVector& storage,
                     const ElementContext& elemCtx,
                     unsigned timeIdx) const
    {
        // the derivative of the storage term depends on the current primary variables;
        // for time indices != 0, the storage term is constant (because these solutions
        // are not changed by the Newton method!)
        if (timeIdx == 0) {
            // calculate the amount of conservation each quantity inside
            // all primary sub control volumes
            const std::size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
            for (unsigned dofIdx = 0; dofIdx < numPrimaryDof; ++dofIdx) {
                storage[dofIdx] = 0.0;

                // the volume of the associated DOF
                const Scalar alpha =
                    elemCtx.stencil(timeIdx).subControlVolume(dofIdx).volume() *
                    elemCtx.intensiveQuantities(dofIdx, timeIdx).extrusionFactor();

                // If the degree of freedom which we currently look at is the one at the
                // center of attention, we need to consider the derivatives for the
                // storage term, else the storage term is constant w.r.t. the primary
                // variables of the focused DOF.
                if (dofIdx == elemCtx.focusDofIndex()) {
                    asImp_().computeStorage(storage[dofIdx],
                                            elemCtx,
                                            dofIdx,
                                            timeIdx);

                    for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
                        storage[dofIdx][eqIdx] *= alpha;
                    }
                }
                else {
                    Dune::FieldVector<Scalar, numEq> tmp;
                    asImp_().computeStorage(tmp,
                                            elemCtx,
                                            dofIdx,
                                            timeIdx);

                    for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
                        storage[dofIdx][eqIdx] = tmp[eqIdx] * alpha;
                    }
                }
            }
        }
        else {
            // for all previous solutions, the storage term does _not_ depend on the
            // current primary variables, so we use scalars to store it.
            if (elemCtx.enableStorageCache()) {
                const std::size_t numPrimaryDof = elemCtx.numPrimaryDof(timeIdx);
                for (unsigned dofIdx = 0; dofIdx < numPrimaryDof; ++dofIdx) {
                    const unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
                    const auto& cachedStorage = elemCtx.model().cachedStorage(globalDofIdx, timeIdx);
                    for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
                        storage[dofIdx][eqIdx] = cachedStorage[eqIdx];
                    }
                }
            }
            else {
                // calculate the amount of conservation each quantity inside
                // all primary sub control volumes
                Dune::FieldVector<Scalar, numEq> tmp;
                const std::size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
                for (unsigned dofIdx = 0; dofIdx < numPrimaryDof; ++dofIdx) {
                    tmp = 0.0;
                    asImp_().computeStorage(tmp,
                                            elemCtx,
                                            dofIdx,
                                            timeIdx);
                    tmp *= elemCtx.stencil(timeIdx).subControlVolume(dofIdx).volume() *
                           elemCtx.intensiveQuantities(dofIdx, timeIdx).extrusionFactor();

                    for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx)
                        storage[dofIdx][eqIdx] = tmp[eqIdx];
                }
            }
        }

#ifndef NDEBUG
        const std::size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
        for (unsigned dofIdx = 0; dofIdx < numPrimaryDof; ++dofIdx) {
            for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
                Valgrind::CheckDefined(storage[dofIdx][eqIdx]);
                assert(isfinite(storage[dofIdx][eqIdx]));
            }
        }
#endif
    }

    void evalFluxes(LocalEvalBlockVector& residual,
                    const ElementContext& elemCtx,
                    unsigned timeIdx) const
    {
        RateVector flux;

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

            Valgrind::SetUndefined(flux);
            asImp_().computeFlux(flux, /*context=*/elemCtx, scvfIdx, timeIdx);
            Valgrind::CheckDefined(flux);
#ifndef NDEBUG
            for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
                assert(isfinite(flux[eqIdx]));
            }
#endif

            const Scalar alpha = elemCtx.extensiveQuantities(scvfIdx, timeIdx).extrusionFactor() * face.area();
            Valgrind::CheckDefined(alpha);
            assert(alpha > 0.0);
            assert(isfinite(alpha));

            for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
                flux[eqIdx] *= alpha;
            }

            // 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
            for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
                assert(isfinite(flux[eqIdx]));
                residual[i][eqIdx] += flux[eqIdx];
                residual[j][eqIdx] -= flux[eqIdx];
            }
        }

#if !defined NDEBUG
        // in debug mode, ensure that the residual is well-defined
        const std::size_t numDof = elemCtx.numDof(timeIdx);
        for (unsigned i = 0; i < numDof; ++i) {
            for (unsigned j = 0; j < numEq; ++j) {
                assert(isfinite(residual[i][j]));
                Valgrind::CheckDefined(residual[i][j]);
            }
        }
#endif
    }

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

    void computeStorage(EqVector&,
                        const ElementContext&,
                        unsigned,
                        unsigned) const
    {
        throw std::logic_error("Not implemented: The local residual " + Dune::className<Implementation>() +
                               " does not implement the required method 'computeStorage()'");
    }

    void computeFlux(RateVector&,
                     const ElementContext&,
                     unsigned,
                     unsigned) const
    {
        throw std::logic_error("Not implemented: The local residual " + Dune::className<Implementation>() +
                               " does not implement the required method 'computeFlux()'");
    }

    void computeSource(RateVector&,
                       const ElementContext&,
                       unsigned,
                       unsigned) const
    {
        throw std::logic_error("Not implemented: The local residual " + Dune::className<Implementation>() +
                               " does not implement the required method 'computeSource()'");
    }

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

        BoundaryContext boundaryCtx(elemCtx);
        // move the iterator to the first boundary
        if (boundaryCtx.intersection(0).neighbor()) {
            boundaryCtx.increment();
        }

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

#if !defined NDEBUG
        // in debug mode, ensure that the residual and the storage terms are well-defined
        const std::size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
        for (unsigned i = 0; i < numDof; ++i) {
            for (unsigned j = 0; j < numEq; ++j) {
                assert(isfinite(residual[i][j]));
                Valgrind::CheckDefined(residual[i][j]);
            }
        }
#endif
    }

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

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

        const auto& stencil = boundaryCtx.stencil(timeIdx);
        const unsigned 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();

            Valgrind::CheckDefined(values[eqIdx]);
            assert(isfinite(values[eqIdx]));
        }

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

    void evalVolumeTerms_(LocalEvalBlockVector& residual,
                          ElementContext& elemCtx) const
    {
        EvalVector tmp;
        EqVector tmp2;
        RateVector sourceRate;

        tmp = 0.0;
        tmp2 = 0.0;

        // evaluate the volumetric terms (storage + source terms)
        const std::size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
        for (unsigned dofIdx = 0; dofIdx < numPrimaryDof; ++dofIdx) {
            const Scalar extrusionFactor =
                elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0).extrusionFactor();
            Valgrind::CheckDefined(extrusionFactor);
            assert(isfinite(extrusionFactor));
            assert(extrusionFactor > 0.0);
            const Scalar scvVolume =
               elemCtx.stencil(/*timeIdx=*/0).subControlVolume(dofIdx).volume() * extrusionFactor;
            Valgrind::CheckDefined(scvVolume);
            assert(isfinite(scvVolume));
            assert(scvVolume > 0.0);

            // if the model uses extensive quantities in its storage term, and we use
            // automatic differention and current DOF is also not the one we currently
            // focus on, the storage term does not need any derivatives!
            if (!extensiveStorageTerm &&
                !std::is_same_v<Scalar, Evaluation> &&
                dofIdx != elemCtx.focusDofIndex())
            {
                asImp_().computeStorage(tmp2, elemCtx, dofIdx, /*timeIdx=*/0);
                for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
                    tmp[eqIdx] = tmp2[eqIdx];
                }
            }
            else {
                asImp_().computeStorage(tmp, elemCtx, dofIdx, /*timeIdx=*/0);
            }

#ifndef NDEBUG
            Valgrind::CheckDefined(tmp);
            for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
                assert(isfinite(tmp[eqIdx]));
            }
#endif

            if (elemCtx.enableStorageCache()) {
                const auto& model = elemCtx.model();
                const unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
                if (elemCtx.problem().iterationContext().isFirstGlobalIteration() &&
                    !elemCtx.haveStashedIntensiveQuantities())
                {
                    if (!elemCtx.problem().recycleFirstIterationStorage()) {
                        // we re-calculate the storage term for the solution of the
                        // previous time step from scratch instead of using the one of
                        // the first iteration of the current time step.
                        tmp2 = 0.0;
                        elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/1);
                        asImp_().computeStorage(tmp2, elemCtx,  dofIdx, /*timeIdx=*/1);
                    }
                    else {
                        // if the storage term is cached and we're in the first iteration
                        // of the time step, use the storage term of the first iteration
                        // as the one as the solution of the last time step (this assumes
                        // that the initial guess for the solution at the end of the time
                        // step is the same as the solution at the beginning of the time
                        // step. This is usually true, but some fancy preprocessing
                        // scheme might invalidate that assumption.)
                        for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
                            tmp2[eqIdx] = Toolbox::value(tmp[eqIdx]);
                        }
                    }

                    Valgrind::CheckDefined(tmp2);

                    model.updateCachedStorage(globalDofIdx, /*timeIdx=*/1, tmp2);
                }
                else {
                    // if the mass storage at the beginning of the time step is not cached,
                    // if the storage term is cached and we're not looking at the first
                    // iteration of the time step, we take the cached data.
                    tmp2 = model.cachedStorage(globalDofIdx, /*timeIdx=*/1);
                    Valgrind::CheckDefined(tmp2);
                }
            }
            else {
                // if the mass storage at the beginning of the time step is not cached,
                // we re-calculate it from scratch.
                tmp2 = 0.0;
                asImp_().computeStorage(tmp2, elemCtx, dofIdx, /*timeIdx=*/1);
                Valgrind::CheckDefined(tmp2);
            }

            // Use the implicit Euler time discretization
            for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
                const double dt = elemCtx.simulator().timeStepSize();
                assert(dt > 0);
                tmp[eqIdx] -= tmp2[eqIdx];
                tmp[eqIdx] *= scvVolume / dt;

                residual[dofIdx][eqIdx] += tmp[eqIdx];
            }

            Valgrind::CheckDefined(residual[dofIdx]);

            // deal with the source term
            asImp_().computeSource(sourceRate, elemCtx, dofIdx, /*timeIdx=*/0);

            // if the model uses extensive quantities in its storage term, and we use
            // automatic differention and current DOF is also not the one we currently
            // focus on, the storage term does not need any derivatives!
            if (!extensiveStorageTerm &&
                !std::is_same_v<Scalar, Evaluation> &&
                dofIdx != elemCtx.focusDofIndex())
            {
                for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
                    residual[dofIdx][eqIdx] -= scalarValue(sourceRate[eqIdx])*scvVolume;
                }
            }
            else {
                for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
                    sourceRate[eqIdx] *= scvVolume;
                    residual[dofIdx][eqIdx] -= sourceRate[eqIdx];
                }
            }

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

#if !defined NDEBUG
        // in debug mode, ensure that the residual is well-defined
        std::size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
        for (unsigned i = 0; i < numDof; ++i) {
            for (unsigned j = 0; j < numEq; ++j) {
                assert(isfinite(residual[i][j]));
                Valgrind::CheckDefined(residual[i][j]);
            }
        }
#endif
    }

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

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

    LocalEvalBlockVector internalResidual_;
};

} // namespace Opm

#endif