staggeredcouplingmanager.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_STAGGERED_COUPLING_MANAGER_HH
#define DUMUX_STAGGERED_COUPLING_MANAGER_HH
 
#include <dumux/multidomain/couplingmanager.hh>
#include <dumux/assembly/numericepsilon.hh>
#include <dumux/common/indextraits.hh>
#include <dumux/discretization/method.hh>
 
namespace Dumux {
 
template<class MDTraits>
class StaggeredCouplingManager: virtual public CouplingManager<MDTraits>
{
    using ParentType = CouplingManager<MDTraits>;
 
    template<std::size_t id> using GridGeometry = typename MDTraits::template SubDomain<id>::GridGeometry;
    template<std::size_t id> using GridView = typename GridGeometry<id>::GridView;
 
    using FVElementGeometry = typename GridGeometry<0>::LocalView;
    using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
    using Element = typename GridView<0>::template Codim<0>::Entity;
 
    using GridIndexType = typename IndexTraits< GridView<0> >::GridIndex;
    using CouplingStencil = std::vector<GridIndexType>;
 
public:
 
    using ParentType::evalCouplingResidual;
 
    using Traits = MDTraits;
 
    static constexpr auto cellCenterIdx = GridGeometry<0>::cellCenterIdx();
    static constexpr auto faceIdx = GridGeometry<0>::faceIdx();
 
    template<std::size_t i, std::size_t j, class LocalAssemblerI, class PriVarsJ>
    void updateCouplingContext(Dune::index_constant<i> domainI,
                               const LocalAssemblerI& localAssemblerI,
                               Dune::index_constant<j> domainJ,
                               const std::size_t dofIdxGlobalJ,
                               const PriVarsJ& priVarsJ,
                               int pvIdxJ)
    {
        auto& curSol = this->curSol(domainJ);
 
        // only proceed if the solution vector has actually been set in the base class
        // (which is not the case if the staggered model is not coupled to another domain)
        if (curSol.size() > 0)
            curSol[dofIdxGlobalJ][pvIdxJ] = priVarsJ[pvIdxJ];
    }
 
    const CouplingStencil& couplingStencil(Dune::index_constant<cellCenterIdx> domainI,
                                           const Element& elementI,
                                           Dune::index_constant<faceIdx> domainJ) const
    {
        const auto& connectivityMap = this->problem(domainI).gridGeometry().connectivityMap();
        const auto eIdx = this->problem(domainI).gridGeometry().elementMapper().index(elementI);
        return connectivityMap(domainI, domainJ, eIdx);
    }
 
    template<std::size_t i, std::size_t j>
    const CouplingStencil couplingStencil(Dune::index_constant<i> domainI,
                                           const SubControlVolumeFace& scvfI,
                                           Dune::index_constant<j> domainJ) const
    {
        static_assert(i != j, "Domain i cannot be coupled to itself!");
        static_assert(AlwaysFalse<Dune::index_constant<i>>::value,
                      "The coupling manager does not implement the couplingStencil() function" );
 
        return CouplingStencil(); // suppress compiler warning of function having no return statement
    }
 
    const CouplingStencil& couplingStencil(Dune::index_constant<faceIdx> domainI,
                                           const SubControlVolumeFace& scvfI,
                                           Dune::index_constant<cellCenterIdx> domainJ) const
    {
        const auto& connectivityMap = this->problem(domainI).gridGeometry().connectivityMap();
        return connectivityMap(domainI, domainJ, scvfI.index());
    }
 
    template<class LocalAssemblerI, std::size_t j>
    decltype(auto) evalCouplingResidual(Dune::index_constant<cellCenterIdx> domainI,
                                        const LocalAssemblerI& localAssemblerI,
                                        Dune::index_constant<j> domainJ,
                                        std::size_t dofIdxGlobalJ) const
    {
        static_assert(domainI != domainJ, "Domain i cannot be coupled to itself!");
        return localAssemblerI.evalLocalResidualForCellCenter();
    }
 
    template<class LocalAssemblerI, std::size_t j>
    decltype(auto) evalCouplingResidual(Dune::index_constant<faceIdx> domainI,
                                        const SubControlVolumeFace& scvfI,
                                        const LocalAssemblerI& localAssemblerI,
                                        Dune::index_constant<j> domainJ,
                                        std::size_t dofIdxGlobalJ) const
    {
        static_assert(domainI != domainJ, "Domain i cannot be coupled to itself!");
        return localAssemblerI.evalLocalResidualForFace(scvfI);
    }
 
    template<std::size_t i, typename std::enable_if_t<(GridGeometry<i>::discMethod != DiscretizationMethods::staggered), int> = 0>
    decltype(auto) numericEpsilon(Dune::index_constant<i> id,
                                  const std::string& paramGroup) const
    {
        return ParentType::numericEpsilon(id, paramGroup);
    }
 
    template<std::size_t i, typename std::enable_if_t<(GridGeometry<i>::discMethod == DiscretizationMethods::staggered), int> = 0>
    decltype(auto) numericEpsilon(Dune::index_constant<i>,
                                  const std::string& paramGroup) const
    {
        constexpr std::size_t numEqCellCenter = Traits::template SubDomain<cellCenterIdx>::PrimaryVariables::dimension;
        constexpr std::size_t numEqFace = Traits::template SubDomain<faceIdx>::PrimaryVariables::dimension;
        constexpr bool isCellCenter = GridGeometry<i>::isCellCenter();
        constexpr std::size_t numEq = isCellCenter ? numEqCellCenter : numEqFace;
        constexpr auto prefix = isCellCenter ? "CellCenter" : "Face";
 
        try {
            if(paramGroup.empty())
                return NumericEpsilon<typename Traits::Scalar, numEq>(prefix);
            else
                return NumericEpsilon<typename Traits::Scalar, numEq>(paramGroup + "." + prefix);
        }
        catch (Dune::RangeError& e)
        {
            DUNE_THROW(Dumux::ParameterException, "For the staggered model, you can to specify \n\n"
                        "  CellCenter.Assembly.NumericDifference.PriVarMagnitude = mCC\n"
                        "  Face.Assembly.NumericDifference.PriVarMagnitude = mFace\n"
                        "  CellCenter.Assembly.NumericDifference.BaseEpsilon = eCC_0 ... eCC_numEqCellCenter-1\n"
                        "  Face.Assembly.NumericDifference.BaseEpsilon = eFace_0 ... eFace_numEqFace-1\n\n"
                        "Wrong number of values set for " << prefix  << " (has " << numEq << " primary variable(s))\n\n" << e);
        }
 
    }
 
};
 
} //end namespace Dumux
 
#endif