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:
/*****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   This program 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 3 of the License, or       *
 *   (at your option) any later version.                                     *
 *                                                                           *
 *   This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/
#ifndef DUMUX_FV_LOCAL_ASSEMBLER_BASE_HH
#define DUMUX_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/assembly/diffmethod.hh>
 
namespace Dumux {
 
template<class TypeTag, class Assembler, class Implementation, bool useImplicitAssembly>
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::ResidualType;
    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 = GetPropType<TypeTag, Properties::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().prevGridVolVars()),
                           localView(assembler.gridVariables().gridFluxVarsCache()),
                           assembler.localResidual(),
                           element.partitionType() == Dune::GhostEntity)
    {}
 
    explicit FVLocalAssemblerBase(const Assembler& assembler,
                                  const Element& element,
                                  const SolutionVector& curSol,
                                  const FVElementGeometry& fvGeometry,
                                  const ElementVolumeVariables& curElemVolVars,
                                  const ElementVolumeVariables& prevElemVolVars,
                                  const ElementFluxVariablesCache& elemFluxVarsCache,
                                  const LocalResidual& localResidual,
                                  const bool elementIsGhost)
    : assembler_(assembler)
    , element_(element)
    , curSol_(curSol)
    , fvGeometry_(fvGeometry)
    , curElemVolVars_(curElemVolVars)
    , prevElemVolVars_(prevElemVolVars)
    , elemFluxVarsCache_(elemFluxVarsCache)
    , localResidual_(localResidual)
    , elementIsGhost_(elementIsGhost)
    {}
 
    static constexpr bool isImplicit()
    { return useImplicitAssembly; }
 
    ElementResidualVector evalLocalResidual() const
    {
        if (!isImplicit())
            if (this->assembler().isStationaryProblem())
                DUNE_THROW(Dune::InvalidStateException, "Using explicit jacobian assembler with stationary local residual");
 
        if (elementIsGhost())
            return ElementResidualVector(0.0);
 
        return isImplicit() ? evalLocalResidual(curElemVolVars())
                            : evalLocalResidual(prevElemVolVars());
    }
 
    ElementResidualVector evalLocalResidual(const ElementVolumeVariables& elemVolVars) const
    {
        if (!assembler().isStationaryProblem())
        {
            ElementResidualVector residual = evalLocalFluxAndSourceResidual(elemVolVars);
            residual += evalLocalStorageResidual();
            return residual;
        }
        else
            return evalLocalFluxAndSourceResidual(elemVolVars);
    }
 
    ElementResidualVector evalLocalFluxAndSourceResidual() const
    {
        return isImplicit() ? evalLocalFluxAndSourceResidual(curElemVolVars())
                            : evalLocalFluxAndSourceResidual(prevElemVolVars());
     }
 
    ElementResidualVector evalLocalFluxAndSourceResidual(const ElementVolumeVariables& elemVolVars) const
    {
        return localResidual_.evalFluxAndSource(element_, fvGeometry_, elemVolVars, elemFluxVarsCache_, elemBcTypes_);
    }
 
    ElementResidualVector evalLocalStorageResidual() const
    {
        return localResidual_.evalStorage(element_, fvGeometry_, prevElemVolVars_, curElemVolVars_);
    }
 
    void bindLocalViews()
    {
        // get some references for convenience
        const auto& element = this->element();
        const auto& curSol = this->curSol();
        const auto& prevSol = this->assembler().prevSol();
        auto&& fvGeometry = this->fvGeometry();
        auto&& curElemVolVars = this->curElemVolVars();
        auto&& prevElemVolVars = this->prevElemVolVars();
        auto&& elemFluxVarsCache = this->elemFluxVarsCache();
 
        // bind the caches
        fvGeometry.bind(element);
 
        if (isImplicit())
        {
            curElemVolVars.bind(element, fvGeometry, curSol);
            elemFluxVarsCache.bind(element, fvGeometry, curElemVolVars);
            if (!this->assembler().isStationaryProblem())
                prevElemVolVars.bindElement(element, fvGeometry, this->assembler().prevSol());
        }
        else
        {
            curElemVolVars.bindElement(element, fvGeometry, curSol);
            prevElemVolVars.bind(element, fvGeometry, prevSol);
            elemFluxVarsCache.bind(element, fvGeometry, prevElemVolVars);
        }
    }
 
    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_; }
 
    ElementVolumeVariables& prevElemVolVars()
    { return prevElemVolVars_; }
 
    ElementFluxVariablesCache& elemFluxVarsCache()
    { return elemFluxVarsCache_; }
 
    LocalResidual& localResidual()
    { return localResidual_; }
 
    ElementBoundaryTypes& elemBcTypes()
    { return elemBcTypes_; }
 
    const FVElementGeometry& fvGeometry() const
    { return fvGeometry_; }
 
    const ElementVolumeVariables& curElemVolVars() const
    { return curElemVolVars_; }
 
    const ElementVolumeVariables& prevElemVolVars() const
    { return prevElemVolVars_; }
 
    const ElementFluxVariablesCache& elemFluxVarsCache() const
    { return elemFluxVarsCache_; }
 
    const ElementBoundaryTypes& elemBcTypes() const
    { return elemBcTypes_; }
 
    const LocalResidual& localResidual() const
    { return localResidual_; }
 
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_;
    ElementVolumeVariables prevElemVolVars_;
    ElementFluxVariablesCache elemFluxVarsCache_;
    ElementBoundaryTypes elemBcTypes_;
 
    LocalResidual localResidual_; 
    bool elementIsGhost_; 
};
 
 
} // end namespace Dumux
 
#endif