assembly/localresidual.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_BASE_LOCAL_RESIDUAL_HH
#define DUMUX_BASE_LOCAL_RESIDUAL_HH
 
#include <cassert>
 
#include <dune/common/exceptions.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/extrusion.hh>
#include <dumux/discretization/cvfe/quadraturerules.hh>
 
namespace Dumux::Experimental {
 
template<class TypeTag>
class LocalResidual
{
    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 Extrusion = Extrusion_t<GridGeometry>;
    using NumEqVector = Dumux::NumEqVector<GetPropType<TypeTag, Properties::PrimaryVariables>>;
    using ElementBoundaryTypes = GetPropType<TypeTag, Properties::ElementBoundaryTypes>;
    using GridVariablesCache = typename GridVariables::GridVariablesCache;
    using ElementVariables = typename GridVariablesCache::LocalView;
    using TimeLoop = TimeLoopBase<Scalar>;
 
public:
    using ElementResidualVector = ReservedBlockVector<NumEqVector, Dumux::Detail::LocalDofs::maxNumLocalDofs<FVElementGeometry>()>;
 
    LocalResidual(const Problem* problem,
                  const TimeLoop* timeLoop = nullptr)
    : problem_(problem)
    , timeLoop_(timeLoop)
    {}
 
    // \{
 
    ElementResidualVector evalStorage(const Element& element,
                                      const FVElementGeometry& fvGeometry,
                                      const ElementVariables& prevElemVars,
                                      const ElementVariables& curElemVars) const
    {
        assert(!this->isStationary() && "no time loop set for storage term evaluation");
 
        // initialize the residual vector for all local dofs in this element
        ElementResidualVector residual(Dumux::Detail::LocalDofs::numLocalDofs(fvGeometry));
 
        // evaluate the volume terms (storage + source terms)
        // forward to the local residual specialized for the discretization methods
        for (const auto& scv : scvs(fvGeometry))
            this->asImp().evalStorage(residual, this->problem(), element, fvGeometry, prevElemVars, curElemVars, scv);
 
        // allow for additional contributions (e.g. hybrid CVFE / FE schemes)
        this->asImp().addToElementStorageResidual(residual, this->problem(), element, fvGeometry, prevElemVars, curElemVars);
 
        return residual;
    }
 
    ElementResidualVector evalFluxAndSource(const Element& element,
                                            const FVElementGeometry& fvGeometry,
                                            const ElementVariables& elemVars,
                                            const ElementBoundaryTypes& bcTypes) const
    {
        // initialize the residual vector for all local dofs in this element
        ElementResidualVector residual(Dumux::Detail::LocalDofs::numLocalDofs(fvGeometry));
 
        // evaluate the volume terms (storage + source terms)
        // forward to the local residual specialized for the discretization methods
        for (const auto& scv : scvs(fvGeometry))
            this->asImp().evalSource(residual, this->problem(), element, fvGeometry, elemVars, scv);
 
        // forward to the local residual specialized for the discretization methods
        for (auto&& scvf : scvfs(fvGeometry))
            this->asImp().evalFlux(residual, this->problem(), element, fvGeometry, elemVars, bcTypes, scvf);
 
        // allow for additional contributions (e.g. hybrid CVFE / FE schemes)
        this->asImp().addToElementFluxAndSourceResidual(residual, this->problem(), element, fvGeometry, elemVars, bcTypes);
 
        return residual;
    }
 
    void addToElementStorageResidual(ElementResidualVector& residual,
                                     const Problem& problem,
                                     const Element& element,
                                     const FVElementGeometry& fvGeometry,
                                     const ElementVariables& prevElemVars,
                                     const ElementVariables& curElemVars) const
    {}
 
    void addToElementFluxAndSourceResidual(ElementResidualVector& residual,
                                           const Problem& problem,
                                           const Element& element,
                                           const FVElementGeometry& fvGeometry,
                                           const ElementVariables& curElemVars,
                                           const ElementBoundaryTypes& bcTypes) const
    {}
 
    // \}
 
    // \{
 
    template<class SubControlVolume>
    NumEqVector storageIntegral(const FVElementGeometry& fvGeometry,
                                const ElementVariables& elemVars,
                                const SubControlVolume& scv,
                                bool isPreviousTimeLevel) const
    {
        DUNE_THROW(Dune::NotImplemented, "This model does not implement a storageIntegral method!");
    }
 
    template<class SubControlVolume>
    NumEqVector sourceIntegral(const FVElementGeometry& fvGeometry,
                               const ElementVariables& elemVars,
                               const SubControlVolume& scv) const
    {
        NumEqVector source(0.0);
 
        const auto& problem = this->asImp().problem();
        for (const auto& qpData : CVFE::quadratureRule(fvGeometry, scv))
            source += qpData.weight() * problem.source(fvGeometry, elemVars, qpData.ipData());
 
        // ToDo: point source data with ipData
        // add contribution from possible point sources
        if (!problem.pointSourceMap().empty())
            source += Extrusion::volume(fvGeometry, scv) * problem.scvPointSources(fvGeometry.element(), fvGeometry, elemVars, scv);
 
        source *= elemVars[scv].extrusionFactor();
 
        return source;
    }
 
    template<class SubControlVolumeFace>
    NumEqVector fluxIntegral(const FVElementGeometry& fvGeometry,
                             const ElementVariables& elemVars,
                             const SubControlVolumeFace& scvf) const
    {
        DUNE_THROW(Dune::NotImplemented, "This model does not implement a fluxIntegral method!");
    }
 
    // \}
 
    // \{
 
    template<class SubControlVolume>
    void evalStorage(ElementResidualVector& residual,
                     const Problem& problem,
                     const Element& element,
                     const FVElementGeometry& fvGeometry,
                     const ElementVariables& prevElemVars,
                     const ElementVariables& curElemVars,
                     const SubControlVolume& scv) const
    {
        // mass balance within the element. this is the
        // \f$\frac{\partial S}{\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 = this->asImp().storageIntegral(fvGeometry, prevElemVars, scv, /*previous time level?*/true);
        NumEqVector storage = this->asImp().storageIntegral(fvGeometry, curElemVars, scv, /*previous time level?*/false);
 
        storage -= prevStorage;
        storage /= timeLoop_->timeStepSize();
 
        residual[scv.localDofIndex()] += storage;
    }
 
    template<class SubControlVolume>
    void evalSource(ElementResidualVector& residual,
                    const Problem& problem,
                    const Element& element,
                    const FVElementGeometry& fvGeometry,
                    const ElementVariables& curElemVars,
                    const SubControlVolume& scv) const
    {
        NumEqVector source = this->asImp().sourceIntegral(fvGeometry, curElemVars, scv);
        residual[scv.localDofIndex()] -= source;
    }
 
    template<class SubControlVolumeFace>
    void evalFlux(ElementResidualVector& residual,
                  const Problem& problem,
                  const Element& element,
                  const FVElementGeometry& fvGeometry,
                  const ElementVariables& elemVars,
                  const ElementBoundaryTypes& elemBcTypes,
                  const SubControlVolumeFace& scvf) const {}
 
    template<class SubControlVolumeFace>
    NumEqVector evalFlux(const Problem& problem,
                         const Element& element,
                         const FVElementGeometry& fvGeometry,
                         const ElementVariables& elemVars,
                         const SubControlVolumeFace& scvf) const
    {
        return asImp().evalFlux(problem, element, fvGeometry, elemVars, scvf);
    }
 
    //\}
 
    // \{
 
    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::Experimental
 
#endif