fvlocalresidual.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_FV_LOCAL_RESIDUAL_HH
#define DUMUX_FV_LOCAL_RESIDUAL_HH
 
#include <dune/common/exceptions.hh>
#include <dune/istl/bvector.hh>
 
#include <dumux/common/typetraits/localdofs_.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/timeloop.hh>
#include <dumux/common/reservedblockvector.hh>
#include <dumux/common/numeqvector.hh>
#include <dumux/discretization/method.hh>
#include <dumux/discretization/extrusion.hh>
 
namespace Dumux {
 
template<class TypeTag>
class FVLocalResidual
{
    using Implementation = GetPropType<TypeTag, Properties::LocalResidual>;
    using Problem = GetPropType<TypeTag, Properties::Problem>;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
    using Element = typename GridView::template Codim<0>::Entity;
    using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
    using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
    using SubControlVolume = typename GridGeometry::SubControlVolume;
    using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
    using Extrusion = Extrusion_t<GridGeometry>;
    using NumEqVector = Dumux::NumEqVector<GetPropType<TypeTag, Properties::PrimaryVariables>>;
    using ElementBoundaryTypes = GetPropType<TypeTag, Properties::ElementBoundaryTypes>;
    using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView;
    using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
    using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
    using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
    using TimeLoop = TimeLoopBase<Scalar>;
 
public:
    using ElementResidualVector = ReservedBlockVector<NumEqVector, Detail::LocalDofs::maxNumLocalDofs<FVElementGeometry>()>;
 
    FVLocalResidual(const Problem* problem,
                    const TimeLoop* timeLoop = nullptr)
    : problem_(problem)
    , timeLoop_(timeLoop)
    {}
 
    // \{
 
    ElementResidualVector evalStorage(const Problem& problem,
                                      const Element &element,
                                      const GridGeometry& gridGeometry,
                                      const GridVariables& gridVariables,
                                      const SolutionVector& sol) const
    {
        // make sure FVElementGeometry and volume variables are bound to the element
        const auto fvGeometry = localView(gridGeometry).bind(element);
        const auto elemVolVars = localView(gridVariables.curGridVolVars()).bind(element, fvGeometry, sol);
 
        ElementResidualVector storage(fvGeometry.numScv());
 
        // calculate the amount of conservation each quantity inside
        // all sub control volumes
        for (auto&& scv : scvs(fvGeometry))
        {
            auto localScvIdx = scv.localDofIndex();
            const auto& volVars = elemVolVars[scv];
            storage[localScvIdx] = asImp().computeStorage(problem, scv, volVars);
            storage[localScvIdx] *= Extrusion::volume(fvGeometry, scv) * volVars.extrusionFactor();
        }
 
        return storage;
    }
 
    // \}
 
    // \{
 
    ElementResidualVector evalStorage(const Element& element,
                                      const FVElementGeometry& fvGeometry,
                                      const ElementVolumeVariables& prevElemVolVars,
                                      const ElementVolumeVariables& curElemVolVars) const
    {
        assert(timeLoop_ && "no time loop set for storage term evaluation");
 
        // initialize the residual vector for all scvs in this element
        ElementResidualVector residual(fvGeometry.numScv());
 
        // evaluate the volume terms (storage + source terms)
        // forward to the local residual specialized for the discretization methods
        for (auto&& scv : scvs(fvGeometry))
            asImp().evalStorage(residual, this->problem(), element, fvGeometry, prevElemVolVars, curElemVolVars, scv);
 
        return residual;
    }
 
    ElementResidualVector evalFluxAndSource(const Element& element,
                                            const FVElementGeometry& fvGeometry,
                                            const ElementVolumeVariables& elemVolVars,
                                            const ElementFluxVariablesCache& elemFluxVarsCache,
                                            const ElementBoundaryTypes &bcTypes) const
    {
        // initialize the residual vector for all scvs in this element
        ElementResidualVector residual(fvGeometry.numScv());
 
        // evaluate the volume terms (storage + source terms)
        // forward to the local residual specialized for the discretization methods
        for (auto&& scv : scvs(fvGeometry))
            asImp().evalSource(residual, this->problem(), element, fvGeometry, elemVolVars, scv);
 
        // forward to the local residual specialized for the discretization methods
        for (auto&& scvf : scvfs(fvGeometry))
            asImp().evalFlux(residual, this->problem(), element, fvGeometry, elemVolVars, bcTypes, elemFluxVarsCache, scvf);
 
        return residual;
    }
 
    // \}
 
 
    // \{
 
    NumEqVector computeStorage(const Problem& problem,
                               const SubControlVolume& scv,
                               const VolumeVariables& volVars) const
    {
        DUNE_THROW(Dune::NotImplemented, "This model does not implement a storage method!");
    }
 
    NumEqVector computeSource(const Problem& problem,
                              const Element& element,
                              const FVElementGeometry& fvGeometry,
                              const ElementVolumeVariables& elemVolVars,
                              const SubControlVolume &scv) const
    {
        NumEqVector source(0.0);
 
        // add contributions from volume flux sources
        source += problem.source(element, fvGeometry, elemVolVars, scv);
 
        // add contribution from possible point sources
        if (!problem.pointSourceMap().empty())
            source += problem.scvPointSources(element, fvGeometry, elemVolVars, scv);
 
        return source;
    }
 
    NumEqVector computeFlux(const Problem& problem,
                            const Element& element,
                            const FVElementGeometry& fvGeometry,
                            const ElementVolumeVariables& elemVolVars,
                            const SubControlVolumeFace& scvf,
                            const ElementFluxVariablesCache& elemFluxVarsCache) const
    {
        DUNE_THROW(Dune::NotImplemented, "This model does not implement a flux method!");
    }
 
    // \}
 
    // \{
 
    void evalStorage(ElementResidualVector& residual,
                     const Problem& problem,
                     const Element& element,
                     const FVElementGeometry& fvGeometry,
                     const ElementVolumeVariables& prevElemVolVars,
                     const ElementVolumeVariables& curElemVolVars,
                     const SubControlVolume& scv) const
    {
        const auto& curVolVars = curElemVolVars[scv];
        const auto& prevVolVars = prevElemVolVars[scv];
 
        // mass balance within the element. this is the
        // \f$\frac{m}{\partial t}\f$ term if using implicit or explicit
        // euler as time discretization.
        //
        // TODO: We might need a more explicit way for
        // doing the time discretization...
 
        NumEqVector prevStorage = asImp().computeStorage(problem, scv, prevVolVars);
        NumEqVector storage = asImp().computeStorage(problem, scv, curVolVars);
 
        prevStorage *= prevVolVars.extrusionFactor();
        storage *= curVolVars.extrusionFactor();
 
        storage -= prevStorage;
        storage *= Extrusion::volume(fvGeometry, scv);
        storage /= timeLoop_->timeStepSize();
 
        residual[scv.localDofIndex()] += storage;
    }
 
    void evalSource(ElementResidualVector& residual,
                    const Problem& problem,
                    const Element& element,
                    const FVElementGeometry& fvGeometry,
                    const ElementVolumeVariables& curElemVolVars,
                    const SubControlVolume& scv) const
    {
        const auto& curVolVars = curElemVolVars[scv];
        NumEqVector source = asImp().computeSource(problem, element, fvGeometry, curElemVolVars, scv);
        source *= Extrusion::volume(fvGeometry, scv)*curVolVars.extrusionFactor();
 
        residual[scv.localDofIndex()] -= source;
    }
 
    void evalFlux(ElementResidualVector& residual,
                  const Problem& problem,
                  const Element& element,
                  const FVElementGeometry& fvGeometry,
                  const ElementVolumeVariables& elemVolVars,
                  const ElementBoundaryTypes& elemBcTypes,
                  const ElementFluxVariablesCache& elemFluxVarsCache,
                  const SubControlVolumeFace& scvf) const {}
 
    NumEqVector evalFlux(const Problem& problem,
                         const Element& element,
                         const FVElementGeometry& fvGeometry,
                         const ElementVolumeVariables& elemVolVars,
                         const ElementFluxVariablesCache& elemFluxVarsCache,
                         const SubControlVolumeFace& scvf) const
    {
        return asImp().evalFlux(problem, element, fvGeometry, elemVolVars, elemFluxVarsCache, scvf);
    }
 
    //\}
 
    // \{
 
    template<class PartialDerivativeMatrix>
    void addStorageDerivatives(PartialDerivativeMatrix& partialDerivatives,
                               const Problem& problem,
                               const Element& element,
                               const FVElementGeometry& fvGeometry,
                               const VolumeVariables& curVolVars,
                               const SubControlVolume& scv) const
    {
        DUNE_THROW(Dune::NotImplemented, "analytic storage derivative");
    }
 
    template<class PartialDerivativeMatrix>
    void addSourceDerivatives(PartialDerivativeMatrix& partialDerivatives,
                              const Problem& problem,
                              const Element& element,
                              const FVElementGeometry& fvGeometry,
                              const VolumeVariables& curVolVars,
                              const SubControlVolume& scv) const
    {
        DUNE_THROW(Dune::NotImplemented, "analytic source derivative");
    }
 
    template<class PartialDerivativeMatrices, class T = TypeTag>
    std::enable_if_t<GetPropType<T, Properties::GridGeometry>::discMethod != DiscretizationMethods::box, void>
    addFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices,
                            const Problem& problem,
                            const Element& element,
                            const FVElementGeometry& fvGeometry,
                            const ElementVolumeVariables& curElemVolVars,
                            const ElementFluxVariablesCache& elemFluxVarsCache,
                            const SubControlVolumeFace& scvf) const
    {
        DUNE_THROW(Dune::NotImplemented, "analytic flux derivative for cell-centered models");
    }
 
    template<class JacobianMatrix, class T = TypeTag>
    std::enable_if_t<GetPropType<T, Properties::GridGeometry>::discMethod == DiscretizationMethods::box, void>
    addFluxDerivatives(JacobianMatrix& A,
                            const Problem& problem,
                            const Element& element,
                            const FVElementGeometry& fvGeometry,
                            const ElementVolumeVariables& curElemVolVars,
                            const ElementFluxVariablesCache& elemFluxVarsCache,
                            const SubControlVolumeFace& scvf) const
    {
        DUNE_THROW(Dune::NotImplemented, "analytic flux derivative for box models");
    }
 
    template<class PartialDerivativeMatrices>
    void addCCDirichletFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices,
                                     const Problem& problem,
                                     const Element& element,
                                     const FVElementGeometry& fvGeometry,
                                     const ElementVolumeVariables& curElemVolVars,
                                     const ElementFluxVariablesCache& elemFluxVarsCache,
                                     const SubControlVolumeFace& scvf) const
    {
        DUNE_THROW(Dune::NotImplemented, "analytic Dirichlet flux derivative");
    }
 
    template<class PartialDerivativeMatrices>
    void addRobinFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices,
                                 const Problem& problem,
                                 const Element& element,
                                 const FVElementGeometry& fvGeometry,
                                 const ElementVolumeVariables& curElemVolVars,
                                 const ElementFluxVariablesCache& elemFluxVarsCache,
                                 const SubControlVolumeFace& scvf) const
    {
        DUNE_THROW(Dune::NotImplemented, "analytic Robin flux derivative");
    }
 
    //\}
 
    // \{
 
    const Problem& problem() const
    { return *problem_; }
 
    const TimeLoop& timeLoop() const
    { return *timeLoop_; }
 
    bool isStationary() const
    { return !timeLoop_; }
 
    // \}
protected:
 
    Implementation &asImp()
    { return *static_cast<Implementation*>(this); }
 
    const Implementation &asImp() const
    { return *static_cast<const Implementation*>(this); }
 
private:
    const Problem* problem_; 
    const TimeLoop* timeLoop_; 
};
 
} // end namespace Dumux
 
#endif