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/common/pointsources.hh>
#include <dumux/discretization/extrusion.hh>
#include <dumux/discretization/cvfe/quadraturerules.hh>
#include <dumux/discretization/concepts.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 ElementDiscretization = 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 GridVariablesCache = typename GridVariables::GridVariablesCache;
    using ElementVariables = typename GridVariablesCache::LocalView;
    using TimeLoop = TimeLoopBase<Scalar>;
 
public:
    using ElementResidualVector = ReservedBlockVector<NumEqVector, Dumux::Detail::LocalDofs::maxNumLocalDofs<ElementDiscretization>()>;

    LocalResidual(const Problem* problem,
                  const TimeLoop* timeLoop = nullptr)
    : problem_(problem)
    , timeLoop_(timeLoop)
    {}

    // \{

    ElementResidualVector evalStorage(const Element& element,
                                      const ElementDiscretization& elemDisc,
                                      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(elemDisc));
 
        // evaluate the volume terms (storage + source terms)
        // forward to the local residual specialized for the discretization methods
        if constexpr (Concepts::FVElementDiscretization<ElementDiscretization>)
        {
            for (const auto& scv : scvs(elemDisc))
                this->asImp().evalStorage(residual, this->problem(), element, elemDisc, prevElemVars, curElemVars, scv);
        }
 
        // allow for additional contributions (e.g. hybrid CVFE / FE schemes)
        this->asImp().addToElementStorageResidual(residual, this->problem(), element, elemDisc, prevElemVars, curElemVars);
 
        return residual;
    }

    ElementResidualVector evalFluxAndSource(const Element& element,
                                            const ElementDiscretization& elemDisc,
                                            const ElementVariables& elemVars) const
    {
        // initialize the residual vector for all local dofs in this element
        ElementResidualVector residual(Dumux::Detail::LocalDofs::numLocalDofs(elemDisc));
 
        if constexpr (Concepts::FVElementDiscretization<ElementDiscretization>)
        {
            // evaluate the volume terms (storage + source terms)
            // forward to the local residual specialized for the discretization methods
            for (const auto& scv : scvs(elemDisc))
                this->asImp().evalSource(residual, this->problem(), element, elemDisc, elemVars, scv);
 
            // forward to the local residual specialized for the discretization methods
            // TODO: Provide scvfs range only over interior scvfs
            for (auto&& scvf : scvfs(elemDisc))
                if(!scvf.boundary())
                    this->asImp().evalFlux(residual, this->problem(), element, elemDisc, elemVars, scvf);
        }
 
        // allow for additional contributions (e.g. hybrid CVFE / FE schemes)
        this->asImp().addToElementFluxAndSourceResidual(residual, this->problem(), element, elemDisc, elemVars);
 
        // add boundary flux contributions
        this->asImp().addBoundaryFluxIntegral(residual, this->problem(), elemDisc, elemVars);
 
        return residual;
    }

    void addToElementStorageResidual(ElementResidualVector& residual,
                                     const Problem& problem,
                                     const Element& element,
                                     const ElementDiscretization& elemDisc,
                                     const ElementVariables& prevElemVars,
                                     const ElementVariables& curElemVars) const
    {}

    void addToElementFluxAndSourceResidual(ElementResidualVector& residual,
                                           const Problem& problem,
                                           const Element& element,
                                           const ElementDiscretization& elemDisc,
                                           const ElementVariables& curElemVars) const
    {}

    void addBoundaryFluxIntegral(ElementResidualVector& residual,
                                 const Problem& problem,
                                 const ElementDiscretization& elemDisc,
                                 const ElementVariables& elemVars) const
    {
        if(!elemDisc.hasBoundaryFaces())
            return;
 
        for (const auto& boundaryFace : boundaryFaces(elemDisc))
        {
            const auto& bcTypes = problem.boundaryTypes(elemDisc, boundaryFace);
            if(!bcTypes.hasFluxBoundary())
                continue;
 
            if constexpr (Concepts::HybridElementDiscretization<ElementDiscretization>
                       || Concepts::FEElementDiscretization<ElementDiscretization>)
                problem.addFEBoundaryFluxIntegral(residual, elemDisc, elemVars, boundaryFace, bcTypes);
 
            if constexpr (Concepts::FVElementDiscretization<ElementDiscretization>)
                for(const auto& scvf : scvfs(elemDisc, boundaryFace))
                    problem.addFVBoundaryFluxIntegral(residual, elemDisc, elemVars, scvf, bcTypes);
        }
    }
 
    // \}

    // \{

    template<class SubControlVolume>
    NumEqVector storageIntegral(const ElementDiscretization& elemDisc,
                                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 ElementDiscretization& elemDisc,
                               const ElementVariables& elemVars,
                               const SubControlVolume& scv) const
    {
        NumEqVector source(0.0);
 
        const auto& problem = this->asImp().problem();
        for (const auto& qpData : CVFE::quadratureRule(elemDisc, scv))
            source += qpData.weight() * problem.source(elemDisc, elemVars, qpData.ipData());
 
        source *= elemVars[scv].extrusionFactor();
 
        // add contribution from possible point sources
        const auto& pointSources = problem.pointSources();
        if (!pointSources.empty())
            for (const auto& context : pointSources.contexts(elemDisc, scv))
            {
                auto psValues = pointSources.eval(elemDisc, elemVars, context);
                source += psValues;
            }
 
        return source;
    }

    template<class SubControlVolumeFace>
    NumEqVector fluxIntegral(const ElementDiscretization& elemDisc,
                             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 ElementDiscretization& elemDisc,
                     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(elemDisc, prevElemVars, scv, /*previous time level?*/true);
        NumEqVector storage = this->asImp().storageIntegral(elemDisc, 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 ElementDiscretization& elemDisc,
                    const ElementVariables& curElemVars,
                    const SubControlVolume& scv) const
    {
        NumEqVector source = this->asImp().sourceIntegral(elemDisc, curElemVars, scv);
        residual[scv.localDofIndex()] -= source;
    }

    template<class SubControlVolumeFace>
    void evalFlux(ElementResidualVector& residual,
                  const Problem& problem,
                  const Element& element,
                  const ElementDiscretization& elemDisc,
                  const ElementVariables& elemVars,
                  const SubControlVolumeFace& scvf) const {}

    template<class SubControlVolumeFace>
    NumEqVector evalFlux(const Problem& problem,
                         const Element& element,
                         const ElementDiscretization& elemDisc,
                         const ElementVariables& elemVars,
                         const SubControlVolumeFace& scvf) const
    {
        return asImp().evalFlux(problem, element, elemDisc, 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