multibinarycouplingmanager.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_MULTIDOMAIN_MULTIBINARY_COUPLINGMANAGER_HH
#define DUMUX_MULTIDOMAIN_MULTIBINARY_COUPLINGMANAGER_HH
 
#include <utility>
#include <memory>
#include <dune/common/hybridutilities.hh>
#include <dumux/common/properties.hh>
#include <dumux/multidomain/traits.hh>
 
namespace Dumux {
 
namespace Detail {
 
template <std::size_t, typename Tuple>
struct HasIndex;
 
template <std::size_t i, typename... Indices>
struct HasIndex<i, std::tuple<Indices...>>
: std::disjunction<std::is_same<Dune::index_constant<i>, Indices>...>
{};
 
} // end namespace Detail
 
template<class MDTraits, class CouplingMap, class ...CouplingMgrs>
class MultiBinaryCouplingManager
{
    // the sub domain type tags
    template<std::size_t id>
    using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
 
    template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
    template<std::size_t id> using GridView = typename GridGeometry<id>::GridView;
    template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity;
    template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView;
    template<std::size_t id> using SubControlVolumeFace = typename FVElementGeometry<id>::SubControlVolumeFace;
    template<std::size_t id> using SubControlVolume = typename FVElementGeometry<id>::SubControlVolume;
    using CouplingStencil = std::vector<std::size_t>;
 
    template<std::size_t id>
    using SubCouplingManagerT = typename std::tuple_element_t<id, std::tuple<CouplingMgrs...>>;
 
    using CMIndices = std::make_index_sequence<sizeof...(CouplingMgrs)>;
    using CouplingManagers = typename Detail::MultiDomainTupleSharedPtr<SubCouplingManagerT, CMIndices>::type;
 
    template<std::size_t id>
    using SubSolutionVector = std::decay_t<decltype(std::declval<typename MDTraits::SolutionVector>()[Dune::index_constant<id>()])>;
    using SolutionVectors = typename MDTraits::template TupleOfSharedPtr<SubSolutionVector>;
 
    static constexpr auto couplingManagerMap_ = CouplingMap::managerMap();
 
    template<std::size_t i, std::size_t j>
    static constexpr auto globalToLocal_(Dune::index_constant<i> domainI, Dune::index_constant<j> domainJ)
    {
        static_assert(i <= MDTraits::numSubDomains && j <= MDTraits::numSubDomains);
        return CouplingMap::globalToLocal(domainI, domainJ);
    }
 
    template<class Map, std::size_t i, std::size_t j>
    static constexpr bool isCoupled_(Dune::index_constant<i> domainI, Dune::index_constant<j>)
    { return Detail::HasIndex<j, std::decay_t<decltype(Map::coupledDomains(domainI))>>::value; }
 
    template<std::size_t i, std::size_t j>
    static constexpr auto subCouplingManagerIndex_(Dune::index_constant<i> domainI, Dune::index_constant<j> domainJ)
    {
        static_assert(
            isCoupled_<CouplingMap>(domainI, domainJ),
            "Sub-coupling manager only exists for coupled domains."
        );
        return couplingManagerMap_[i][j];
    }
 
public:
    template<std::size_t i, std::size_t j>
    using SubCouplingManager = SubCouplingManagerT<couplingManagerMap_[i][j]>;
 
    MultiBinaryCouplingManager()
    {
        using namespace Dune::Hybrid;
        forEach(couplingManagers_, [&](auto&& couplingManager)
        {
            couplingManager = std::make_shared<typename std::decay_t<decltype(couplingManager)>::element_type>();
        });
 
        forEach(solutionVectors_, [&](auto&& solutionVector)
        {
            solutionVector = std::make_shared<typename std::decay_t<decltype(solutionVector)>::element_type>();
        });
    }
 
    template<std::size_t i, std::size_t j>
    auto& subCouplingManager(Dune::index_constant<i> domainI,
                             Dune::index_constant<j> domainJ)
    {
        constexpr auto idx = subCouplingManagerIndex_(domainI, domainJ);
        return *std::get<idx>(couplingManagers_);
    }
 
    template<std::size_t i, std::size_t j>
    const auto& subCouplingManager(Dune::index_constant<i> domainI,
                                   Dune::index_constant<j> domainJ) const
    {
        constexpr auto idx = subCouplingManagerIndex_(domainI, domainJ);
        return *std::get<idx>(couplingManagers_);
    }
 
    template<std::size_t i, std::size_t j, class Apply>
    decltype(auto) subApply(Dune::index_constant<i> domainI,
                            Dune::index_constant<j> domainJ,
                            Apply&& apply)
    {
        constexpr auto localIndices = globalToLocal_(domainI, domainJ);
        return apply(subCouplingManager(domainI, domainJ), localIndices.first, localIndices.second);
    }
 
    template<std::size_t i, std::size_t j, class Apply>
    decltype(auto) subApply(Dune::index_constant<i> domainI,
                            Dune::index_constant<j> domainJ,
                            const Apply& apply) const
    {
        constexpr auto localIndices = globalToLocal_(domainI, domainJ);
        return apply(subCouplingManager(domainI, domainJ), localIndices.first, localIndices.second);
    }
 
    template<std::size_t i, class Apply>
    decltype(auto) subApply(Dune::index_constant<i> domainI, Apply&& apply)
    {
        constexpr auto dm = CouplingMap::coupledDomains(domainI);
        static_assert(std::tuple_size_v<std::decay_t<decltype(dm)>> != 0, "subApply on uncoupled domain!");
        constexpr auto domainJ = std::get<0>(dm);
        constexpr auto localIndices = globalToLocal_(domainI, domainJ);
        return apply(subCouplingManager(domainI, domainJ), localIndices.first);
    }
 
    template<std::size_t i, class Apply>
    decltype(auto) subApply(Dune::index_constant<i> domainI, const Apply& apply) const
    {
        constexpr auto dm = CouplingMap::coupledDomains(domainI);
        static_assert(std::tuple_size_v<std::decay_t<decltype(dm)>> != 0, "subApply on uncoupled domain!");
        constexpr auto domainJ = std::get<0>(dm);
        constexpr auto localIndices = globalToLocal_(domainI, domainJ);
        return apply(subCouplingManager(domainI, domainJ), localIndices.first);
    }
 
    void updateSolution(const typename MDTraits::SolutionVector& curSol)
    {
        using namespace Dune::Hybrid;
        forEach(integralRange(Dune::Hybrid::size(solutionVectors_)), [&](const auto id)
        {
            *std::get<id>(solutionVectors_) = curSol[id];
        });
    }
 
    template<std::size_t id, class JacobianPattern>
    void extendJacobianPattern(Dune::index_constant<id> domainI, JacobianPattern& pattern) const
    {}
 
    template<std::size_t i, class Entity, std::size_t j>
    const auto& couplingStencil(Dune::index_constant<i> domainI,
                                const Entity& entity,
                                Dune::index_constant<j> domainJ) const
    {
        // if the domains are coupled according to the map, forward to sub-coupling manager
        if constexpr (isCoupled_<CouplingMap>(domainI, domainJ))
            return subApply(domainI, domainJ, [&](const auto& cm, auto&& ii, auto&& jj) -> const auto& {
                return cm.couplingStencil(ii, entity, jj);
            });
        else
            return emptyStencil_;
    }
 
    // ! evaluate coupling residual for the derivative low dim DOF with respect to bulk DOF
    // ! we only need to evaluate the part of the residual that will be influence by the bulk DOF
    template<std::size_t i, class LocalAssemblerI, std::size_t j>
    decltype(auto) evalCouplingResidual(Dune::index_constant<i> domainI,
                                        const SubControlVolumeFace<i>& scvfI,
                                        const LocalAssemblerI& localAssemblerI,
                                        Dune::index_constant<j> domainJ,
                                        std::size_t dofIdxGlobalJ) const
    {
        if constexpr (i == j || isCoupled_<CouplingMap>(domainI, domainJ))
            return subApply(domainI, domainJ, [&](const auto& cm, auto&& ii, auto&& jj) -> decltype(auto) {
                return cm.evalCouplingResidual(ii, scvfI, localAssemblerI, jj, dofIdxGlobalJ);
            });
        else
        {
            DUNE_THROW(Dune::InvalidStateException,
                "Calling evalCouplingResidual for uncoupled domains " << i << " and " << j
            );
            return localAssemblerI.evalLocalResidual();
        }
    }
 
    template<std::size_t i, class LocalAssemblerI, std::size_t j>
    decltype(auto) evalCouplingResidual(Dune::index_constant<i> domainI,
                                        const LocalAssemblerI& localAssemblerI,
                                        Dune::index_constant<j> domainJ,
                                        std::size_t dofIdxGlobalJ) const
    {
        if constexpr (i == j || isCoupled_<CouplingMap>(domainI, domainJ))
            return subApply(domainI, domainJ, [&](const auto& cm, auto&& ii, auto&& jj) -> decltype(auto) {
                return cm.evalCouplingResidual(ii, localAssemblerI, jj, dofIdxGlobalJ);
            });
        else
        {
            DUNE_THROW(Dune::InvalidStateException,
                "Calling evalCouplingResidual for uncoupled domains " << i << " and " << j
            );
            return localAssemblerI.evalLocalResidual();
        }
    }
 
    template<std::size_t i, class LocalAssemblerI, std::size_t j>
    decltype(auto) evalCouplingResidual(Dune::index_constant<i> domainI,
                                        const LocalAssemblerI& localAssemblerI,
                                        const SubControlVolume<i>& scvI,
                                        Dune::index_constant<j> domainJ,
                                        std::size_t dofIdxGlobalJ) const
    {
        if constexpr (i == j || isCoupled_<CouplingMap>(domainI, domainJ))
            return subApply(domainI, domainJ, [&](const auto& cm, auto&& ii, auto&& jj) -> decltype(auto) {
                return cm.evalCouplingResidual(ii, localAssemblerI, scvI, jj, dofIdxGlobalJ);
            });
        else
        {
            DUNE_THROW(Dune::InvalidStateException,
                "Calling evalCouplingResidual for uncoupled domains " << i << " and " << j
            );
            return localAssemblerI.evalLocalResidual();
        }
    }
 
    template<std::size_t i, class LocalAssemblerI, std::size_t j, class PrimaryVariables>
    void updateCouplingContext(Dune::index_constant<i> domainI,
                               const LocalAssemblerI& localAssemblerI,
                               Dune::index_constant<j> domainJ,
                               const std::size_t dofIdxGlobalJ,
                               const PrimaryVariables& priVars,
                               int pvIdxJ)
    {
        // only one other manager needs to take action if i != j
        if constexpr (i != j)
        {
            if constexpr (isCoupled_<CouplingMap>(domainI, domainJ))
                subApply(domainI, domainJ, [&](auto& cm, auto&& ii, auto&& jj){
                    cm.updateCouplingContext(ii, localAssemblerI, jj, dofIdxGlobalJ, priVars, pvIdxJ);
                });
        }
        else
        {
            // for i == j, we need to call all relevant managers where domainI is involved and make sure that the
            // context is updated with respect to its own domain (I-I coupling context)
            using namespace Dune::Hybrid;
            forEach(CouplingMap::coupledDomains(domainI), [&](const auto domainJ){
                static_assert(domainI != domainJ);
                subApply(domainI, domainJ, [&](auto& cm, auto&& ii, auto&& jj){
                    cm.updateCouplingContext(ii, localAssemblerI, ii, dofIdxGlobalJ, priVars, pvIdxJ);
                });
            });
        }
    }
 
    template<std::size_t i, class Assembler = int>
    void bindCouplingContext(Dune::index_constant<i> domainI, const Element<i>& element, const Assembler& assembler = 0)
    {
        // for the coupling blocks
        using namespace Dune::Hybrid;
        forEach(CouplingMap::coupledDomains(domainI), [&](const auto domainJ){
            subApply(domainI, domainJ, [&](auto& cm, auto&& ii, auto&& jj) -> void {
                cm.bindCouplingContext(ii, element, assembler);
            });
        });
    }
 
    template<std::size_t i>
    decltype(auto) numericEpsilon(Dune::index_constant<i> domainI,
                                  const std::string& paramGroup) const
    {
        return subApply(domainI, [&](const auto& cm, auto&& ii) -> decltype(auto) {
            return cm.numericEpsilon(ii, paramGroup);
        });
    }
 
    template<std::size_t i, class LocalAssemblerI, class JacobianMatrixDiagBlock, class GridVariables>
    void evalAdditionalDomainDerivatives(Dune::index_constant<i> domainI,
                                         const LocalAssemblerI& localAssemblerI,
                                         const typename LocalAssemblerI::LocalResidual::ElementResidualVector& origResiduals,
                                         JacobianMatrixDiagBlock& A,
                                         GridVariables& gridVariables)
    {}
 
    template<std::size_t i, class LocalAssemblerI, class UpdatableElementVolVars, class UpdatableFluxVarCache>
    void updateCoupledVariables(Dune::index_constant<i> domainI,
                                const LocalAssemblerI& localAssemblerI,
                                UpdatableElementVolVars& elemVolVars,
                                UpdatableFluxVarCache& elemFluxVarsCache)
    {
        // for the coupling blocks
        using namespace Dune::Hybrid;
        forEach(CouplingMap::coupledDomains(domainI), [&](const auto domainJ){
            subApply(domainI, domainJ, [&](auto& cm, auto&& ii, auto&& jj){
                cm.updateCoupledVariables(ii, localAssemblerI, elemVolVars, elemFluxVarsCache);
            });
        });
    }
 
protected:
    SolutionVectors& curSol()
    { return solutionVectors_; }
 
    const SolutionVectors& curSol() const
    { return solutionVectors_; }
 
private:
    CouplingManagers couplingManagers_;
    SolutionVectors solutionVectors_;
 
    CouplingStencil emptyStencil_;
};
 
} // end namespace Dumux
 
#endif