experimental/assembly/fvlocalassemblerbase.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-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
// SPDX-License-Identifier: GPL-3.0-or-later
//
#ifndef DUMUX_EXPERIMENTAL_FV_LOCAL_ASSEMBLER_BASE_HH
#define DUMUX_EXPERIMENTAL_FV_LOCAL_ASSEMBLER_BASE_HH
 
#include <dune/common/reservedvector.hh>
#include <dune/grid/common/gridenums.hh> // for GhostEntity
#include <dune/istl/matrixindexset.hh>
 
#include <dumux/common/reservedblockvector.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/discretization/extrusion.hh>
#include <dumux/assembly/diffmethod.hh>
 
namespace Dumux::Experimental {
 
template<class TypeTag, class Assembler, class Implementation>
class FVLocalAssemblerBase
{
    using Problem = GetPropType<TypeTag, Properties::Problem>;
    using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using JacobianMatrix = GetPropType<TypeTag, Properties::JacobianMatrix>;
    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
    using SolutionVector = typename Assembler::SolutionVector;
    using ElementBoundaryTypes = GetPropType<TypeTag, Properties::ElementBoundaryTypes>;
    using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
    using SubControlVolume = typename FVElementGeometry::SubControlVolume;
    using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
    using GridVolumeVariables = GetPropType<TypeTag, Properties::GridVolumeVariables>;
    using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
    using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
    using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView;
    using Element = typename GridView::template Codim<0>::Entity;
    static constexpr auto numEq = GetPropType<TypeTag, Properties::ModelTraits>::numEq();
 
public:
    using LocalResidual = std::decay_t<decltype(std::declval<Assembler>().localResidual())>;
    using ElementResidualVector = typename LocalResidual::ElementResidualVector;
 
    explicit FVLocalAssemblerBase(const Assembler& assembler,
                                  const Element& element,
                                  const SolutionVector& curSol)
    : FVLocalAssemblerBase(assembler,
                           element,
                           curSol,
                           localView(assembler.gridGeometry()),
                           localView(assembler.gridVariables().curGridVolVars()),
                           localView(assembler.gridVariables().gridFluxVarsCache()),
                           assembler.localResidual(),
                           element.partitionType() == Dune::GhostEntity,
                           assembler.isImplicit())
    {}
 
    explicit FVLocalAssemblerBase(const Assembler& assembler,
                                  const Element& element,
                                  const SolutionVector& curSol,
                                  const FVElementGeometry& fvGeometry,
                                  const ElementVolumeVariables& curElemVolVars,
                                  const ElementFluxVariablesCache& elemFluxVarsCache,
                                  const LocalResidual& localResidual,
                                  const bool elementIsGhost,
                                  const bool isImplicit)
    : assembler_(assembler)
    , element_(element)
    , curSol_(curSol)
    , fvGeometry_(fvGeometry)
    , curElemVolVars_(curElemVolVars)
    , elemFluxVarsCache_(elemFluxVarsCache)
    , localOperator_(localResidual)
    , elementIsGhost_(elementIsGhost)
    , isImplicit_(isImplicit)
    {}
 
    ElementResidualVector evalLocalResidual() const
    {
        if (elementIsGhost())
            return ElementResidualVector(0.0);
 
        return evalLocalResidual(curElemVolVars());
    }
 
    ElementResidualVector evalLocalResidual(const ElementVolumeVariables& elemVolVars) const
    {
        if (!assembler().isStationaryProblem())
        {
            ElementResidualVector residual = evalLocalFluxAndSourceResidual(elemVolVars);
            residual += evalStorage();
            return residual;
        }
        else
            return evalLocalFluxAndSourceResidual(elemVolVars);
    }
 
    ElementResidualVector evalLocalFluxAndSourceResidual() const
    {
        return evalLocalFluxAndSourceResidual(curElemVolVars());
    }
 
    ElementResidualVector evalLocalFluxAndSourceResidual(const ElementVolumeVariables& elemVolVars) const
    {
        return localOperator_.evalFluxAndSource(element_, fvGeometry_, elemVolVars, elemFluxVarsCache_, elemBcTypes_);
    }
 
    ElementResidualVector evalStorage() const
    {
        return localOperator_.evalStorage(fvGeometry_, curElemVolVars_);
    }
 
    void bindLocalViews()
    {
        // get some references for convenience
        const auto& element = this->element();
        const auto& curSol = this->curSol();
        auto&& fvGeometry = this->fvGeometry();
        auto&& curElemVolVars = this->curElemVolVars();
        auto&& elemFluxVarsCache = this->elemFluxVarsCache();
 
        // bind the caches
        fvGeometry.bind(element);
        if (std::abs(this->localResidual().spatialWeight()) < 1e-6)
            curElemVolVars.bindElement(element, fvGeometry, curSol);
        else
        {
            curElemVolVars.bind(element, fvGeometry, curSol);
            elemFluxVarsCache.bind(element, fvGeometry, curElemVolVars);
        }
    }
 
    template<typename ApplyFunction, class P = Problem, typename std::enable_if_t<P::enableInternalDirichletConstraints(), int> = 0>
    void enforceInternalDirichletConstraints(const ApplyFunction& applyDirichlet)
    {
        // enforce Dirichlet constraints strongly by overwriting partial derivatives with 1 or 0
        // and set the residual to (privar - dirichletvalue)
        for (const auto& scvI : scvs(this->fvGeometry()))
        {
            const auto internalDirichletConstraints = asImp_().problem().hasInternalDirichletConstraint(this->element(), scvI);
            if (internalDirichletConstraints.any())
            {
                const auto dirichletValues = asImp_().problem().internalDirichlet(this->element(), scvI);
                // set the Dirichlet conditions in residual and jacobian
                for (int eqIdx = 0; eqIdx < internalDirichletConstraints.size(); ++eqIdx)
                    if (internalDirichletConstraints[eqIdx])
                        applyDirichlet(scvI, dirichletValues, eqIdx, eqIdx);
            }
        }
    }
 
    template<typename ApplyFunction, class P = Problem, typename std::enable_if_t<!P::enableInternalDirichletConstraints(), int> = 0>
    void enforceInternalDirichletConstraints(const ApplyFunction& applyDirichlet)
    {}
 
    const Problem& problem() const
    { return assembler_.problem(); }
 
    const Assembler& assembler() const
    { return assembler_; }
 
    const Element& element() const
    { return element_; }
 
    bool elementIsGhost() const
    { return elementIsGhost_; }
 
    const SolutionVector& curSol() const
    { return curSol_; }
 
    FVElementGeometry& fvGeometry()
    { return fvGeometry_; }
 
    ElementVolumeVariables& curElemVolVars()
    { return curElemVolVars_; }
 
    ElementFluxVariablesCache& elemFluxVarsCache()
    { return elemFluxVarsCache_; }
 
    LocalResidual& localResidual()
    { return localOperator_; }
 
    ElementBoundaryTypes& elemBcTypes()
    { return elemBcTypes_; }
 
    const FVElementGeometry& fvGeometry() const
    { return fvGeometry_; }
 
    const ElementVolumeVariables& curElemVolVars() const
    { return curElemVolVars_; }
 
    const ElementFluxVariablesCache& elemFluxVarsCache() const
    { return elemFluxVarsCache_; }
 
    const ElementBoundaryTypes& elemBcTypes() const
    { return elemBcTypes_; }
 
    const LocalResidual& localResidual() const
    { return localOperator_; }
 
    bool isImplicit() const
    { return isImplicit_; }
 
protected:
    Implementation &asImp_()
    { return *static_cast<Implementation*>(this); }
 
    const Implementation &asImp_() const
    { return *static_cast<const Implementation*>(this); }
 
    template<class T = TypeTag, typename std::enable_if_t<!GetPropType<T, Properties::GridVariables>::GridVolumeVariables::cachingEnabled, int> = 0>
    VolumeVariables& getVolVarAccess(GridVolumeVariables& gridVolVars, ElementVolumeVariables& elemVolVars, const SubControlVolume& scv)
    { return elemVolVars[scv]; }
 
    template<class T = TypeTag, typename std::enable_if_t<GetPropType<T, Properties::GridVariables>::GridVolumeVariables::cachingEnabled, int> = 0>
    VolumeVariables& getVolVarAccess(GridVolumeVariables& gridVolVars, ElementVolumeVariables& elemVolVars, const SubControlVolume& scv)
    { return gridVolVars.volVars(scv); }
 
private:
 
    const Assembler& assembler_; 
    const Element& element_; 
    const SolutionVector& curSol_; 
 
    FVElementGeometry fvGeometry_;
    ElementVolumeVariables curElemVolVars_;
    ElementFluxVariablesCache elemFluxVarsCache_;
    ElementBoundaryTypes elemBcTypes_;
 
    LocalResidual localOperator_; 
    bool elementIsGhost_; 
    bool isImplicit_; 
};
 
} // end namespace Dumux
 
#endif