From 619f50b31315711e3487c649cabbb2f9ff7eeabd Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Thu, 9 Sep 2021 15:52:10 +0200
Subject: [PATCH] [md] Add meta-coupling manager combining multiple binary
coupling managers
---
.../multidomain/multibinarycouplingmanager.hh | 312 ++++++++++++++++++
1 file changed, 312 insertions(+)
create mode 100644 dumux/multidomain/multibinarycouplingmanager.hh
diff --git a/dumux/multidomain/multibinarycouplingmanager.hh b/dumux/multidomain/multibinarycouplingmanager.hh
new file mode 100644
index 0000000000..9979e8255e
--- /dev/null
+++ b/dumux/multidomain/multibinarycouplingmanager.hh
@@ -0,0 +1,312 @@
+// -*- 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/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup MultiDomain
+ * \copydoc Dumux::MultiBinaryCouplingManager
+ */
+
+#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 {
+
+/*!
+ * \ingroup MultiDomain
+ * \brief Coupling manager that combines an arbitrary number of binary coupling manager (coupling two domains each)
+ */
+template<class MDTraits, class CouplingManagerMaps, 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 Indices = std::make_index_sequence<sizeof...(CouplingMgrs)>;
+ using CouplingManagers = typename Detail::MultiDomainTupleSharedPtr<SubCouplingManagerT, Indices>::type;
+
+ template<std::size_t id>
+ using SolutionVectorT = std::decay_t<decltype(std::declval<typename MDTraits::SolutionVector>()[Dune::index_constant<id>()])>;
+ using SolutionVectors = typename Detail::MultiDomainTupleSharedPtr<SolutionVectorT, Indices>::type;
+
+public:
+
+ static constexpr auto couplingManagerMap = CouplingManagerMaps::makeCouplingManagerMap();
+
+ 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>();
+ });
+ }
+
+ //! Returns the local domain indices used within the binary sub coupling managers.
+ template<std::size_t i, std::size_t j>
+ static constexpr auto globalToLocalDomainIndices(Dune::index_constant<i> domainI, Dune::index_constant<j> domainJ)
+ {
+ static_assert(i <= MDTraits::numSubDomains && j <= MDTraits::numSubDomains);
+ return CouplingManagerMaps::globalToLocalDomainIndices(domainI, domainJ);
+ }
+
+ //! Returns the coupling manager index for a given domain combination
+ template<std::size_t i, std::size_t j>
+ static constexpr auto getCouplingManagerIndex(Dune::index_constant<i>, Dune::index_constant<j>)
+ {
+ static_assert(i != j);
+ return couplingManagerMap[i][j];
+ }
+
+ //! return the binary sub-coupling manager
+ template<std::size_t i, std::size_t j>
+ auto& subCouplingManager(Dune::index_constant<i> domainI,
+ Dune::index_constant<j> domainJ)
+ {
+ constexpr auto idx = getCouplingManagerIndex(domainI, domainJ);
+ return *std::get<idx>(couplingManagers_);
+ }
+
+ //! return the binary sub-coupling manager
+ 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 = getCouplingManagerIndex(domainI, domainJ);
+ return *std::get<idx>(couplingManagers_);
+ }
+
+ //! Update the solution vector before assembly
+ 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];
+ });
+ }
+
+ /*!
+ * \brief extend the jacobian pattern of the diagonal block of domain i
+ * by those entries that are not already in the uncoupled pattern
+ * \note per default we do not add such additional dependencies
+ * \note Such additional dependencies can arise from the coupling, e.g. if a coupling source
+ * term depends on a non-local average of a quantity of the same domain
+ * \warning if you overload this also implement evalAdditionalDomainDerivatives
+ */
+ template<std::size_t id, class JacobianPattern>
+ void extendJacobianPattern(Dune::index_constant<id> domainI, JacobianPattern& pattern) const
+ {}
+
+ /*!
+ * \brief Return the coupling element stencil for a given bulk domain element
+ */
+ 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
+ {
+ return subCouplingManager(domainI, domainJ).couplingStencil(globalToLocalDomainIndices(domainI, domainJ).domainI,
+ entity,
+ globalToLocalDomainIndices(domainI, domainJ).domainJ);
+ }
+
+ // ! 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
+ {
+ return subCouplingManager(domainI, domainJ).evalCouplingResidual(globalToLocalDomainIndices(domainI, domainJ).domainI,
+ scvfI, localAssemblerI,
+ globalToLocalDomainIndices(domainI, domainJ).domainJ,
+ dofIdxGlobalJ);
+ }
+
+ //! 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 LocalAssemblerI& localAssemblerI,
+ Dune::index_constant<j> domainJ,
+ std::size_t dofIdxGlobalJ) const
+ {
+ return subCouplingManager(domainI, domainJ).evalCouplingResidual(globalToLocalDomainIndices(domainI, domainJ).domainI,
+ localAssemblerI,
+ globalToLocalDomainIndices(domainI, domainJ).domainJ,
+ dofIdxGlobalJ);
+ }
+
+ //! 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 LocalAssemblerI& localAssemblerI,
+ const SubControlVolume<i>& scvI,
+ Dune::index_constant<j> domainJ,
+ std::size_t dofIdxGlobalJ) const
+ {
+ return subCouplingManager(domainI, domainJ).evalCouplingResidual(globalToLocalDomainIndices(domainI, domainJ).domainI,
+ localAssemblerI,
+ scvI,
+ globalToLocalDomainIndices(domainI, domainJ).domainJ,
+ dofIdxGlobalJ);
+ }
+
+ /*!
+ * \brief Update the coupling context for the bulk face residual w.r.t to the lowDim dofs
+ */
+ 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)
+ {
+ subCouplingManager(domainI, domainJ).updateCouplingContext(globalToLocalDomainIndices(domainI, domainJ).domainI,
+ localAssemblerI,
+ globalToLocalDomainIndices(domainI, domainJ).domainJ,
+ 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
+ using namespace Dune::Hybrid;
+ forEach(CouplingManagerMaps::coupledDomains(domainI), [&](const auto domainJ)
+ {
+ static_assert(domainI != domainJ);
+ subCouplingManager(domainI, domainJ).updateCouplingContext(globalToLocalDomainIndices(domainI, domainJ).domainI,
+ localAssemblerI,
+ globalToLocalDomainIndices(domainI, domainJ).domainI,
+ dofIdxGlobalJ,
+ priVars,
+ pvIdxJ);
+ });
+ }
+ }
+
+ //! Bind the coupling context for a low dim element TODO remove Assembler
+ 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(integralRange(Dune::index_constant<MDTraits::numSubDomains>{}), [&](const auto domainJ)
+ {
+ if constexpr (domainI != domainJ)
+ this->subCouplingManager(domainI, domainJ).bindCouplingContext(globalToLocalDomainIndices(domainI, domainJ).domainI,
+ element,
+ assembler);
+ });
+ }
+
+ /*!
+ * \brief return the numeric epsilon used for deflecting primary variables of coupled domain i.
+ * \note specialization for free-flow schemes
+ */
+ template<std::size_t i>
+ decltype(auto) numericEpsilon(Dune::index_constant<i> domainI,
+ const std::string& paramGroup) const
+ {
+ static constexpr auto someOtherDomain = std::get<0>(CouplingManagerMaps::coupledDomains(domainI));
+ return subCouplingManager(domainI, someOtherDomain).numericEpsilon(globalToLocalDomainIndices(domainI, someOtherDomain).domainI, paramGroup);
+ }
+
+ /*!
+ * \brief evaluate additional derivatives of the element residual of a domain with respect
+ * to dofs in the same domain that are not in the regular stencil (see CouplingManager::extendJacobianPattern)
+ * \note Such additional dependencies can arise from the coupling, e.g. if a coupling source
+ * term depends on a non-local average of a quantity of the same domain
+ */
+ 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(integralRange(Dune::index_constant<MDTraits::numSubDomains>{}), [&](const auto domainJ)
+ {
+ if constexpr (domainI != domainJ) // TODO is this correct?
+ subCouplingManager(domainI, domainJ).updateCoupledVariables(globalToLocalDomainIndices(domainI, domainJ).domainI,
+ localAssemblerI,
+ elemVolVars,
+ elemFluxVarsCache);
+ });
+ }
+
+protected:
+ SolutionVectors& curSol()
+ { return solutionVectors_; }
+
+ const SolutionVectors& curSol() const
+ { return solutionVectors_; }
+
+private:
+ CouplingManagers couplingManagers_;
+ SolutionVectors solutionVectors_;
+};
+
+} // end namespace Dumux
+
+#endif
--
GitLab