version 3.8
multidomain/boundary/freeflowporenetwork/ffmassporenetwork/couplingmanager.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3//
4// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
5// SPDX-License-Identifier: GPL-3.0-or-later
6//
13#ifndef DUMUX_MULTIDOMAIN_BOUNDARY_FREEFLOWPORENETWORK_FFMASSPORENETWORK_COUPLINGMANAGER_HH
14#define DUMUX_MULTIDOMAIN_BOUNDARY_FREEFLOWPORENETWORK_FFMASSPORENETWORK_COUPLINGMANAGER_HH
15
16#include <utility>
17#include <memory>
18
19#include <dune/common/float_cmp.hh>
20#include <dune/common/exceptions.hh>
25
26namespace Dumux {
27
32template<class MDTraits>
34: public CouplingManager<MDTraits>
35{
36 using Scalar = typename MDTraits::Scalar;
38
39public:
40 static constexpr auto freeFlowMassIndex = typename MDTraits::template SubDomain<0>::Index();
41 static constexpr auto poreNetworkIndex = typename MDTraits::template SubDomain<1>::Index();
42
44private:
45
46 // obtain the type tags of the sub problems
47 using FreeFlowMassTypeTag = typename MDTraits::template SubDomain<freeFlowMassIndex>::TypeTag;
48 using PoreNetworkTypeTag = typename MDTraits::template SubDomain<poreNetworkIndex>::TypeTag;
49
50 using CouplingStencils = std::unordered_map<std::size_t, std::vector<std::size_t> >;
51 using CouplingStencil = CouplingStencils::mapped_type;
52
53 // the sub domain type tags
54 template<std::size_t id>
55 using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
56
57 template<std::size_t id> using GridView = typename GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>::GridView;
58 template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>;
59 template<std::size_t id> using ElementVolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::LocalView;
60 template<std::size_t id> using GridVolumeVariables = GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>;
61 template<std::size_t id> using VolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::VolumeVariables;
62 template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
63 template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView;
64 template<std::size_t id> using GridVariables = GetPropType<SubDomainTypeTag<id>, Properties::GridVariables>;
65 template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity;
66 template<std::size_t id> using PrimaryVariables = GetPropType<SubDomainTypeTag<id>, Properties::PrimaryVariables>;
67 template<std::size_t id> using SubControlVolumeFace = typename FVElementGeometry<id>::SubControlVolumeFace;
68 template<std::size_t id> using SubControlVolume = typename FVElementGeometry<id>::SubControlVolume;
69
70 using VelocityVector = typename Element<freeFlowMassIndex>::Geometry::GlobalCoordinate;
71
72 struct FreeFlowMassCouplingContext
73 {
74 SubControlVolume<poreNetworkIndex> scv;
75 VolumeVariables<poreNetworkIndex> volVars;
76 mutable VelocityVector velocity; // velocity needs to be set externally, not available in this class
77 };
78
79 struct PoreNetworkCouplingContext
80 {
81 SubControlVolume<freeFlowMassIndex> scv;
82 SubControlVolumeFace<freeFlowMassIndex> scvf;
83 VolumeVariables<freeFlowMassIndex> volVars;
84 mutable VelocityVector velocity; // velocity needs to be set externally, not available in this class
85 };
86
87 using CouplingMapper = StaggeredFreeFlowPoreNetworkCouplingMapper;
88
89public:
90
94 // \{
95
97 void init(std::shared_ptr<Problem<freeFlowMassIndex>> freeFlowMassProblem,
98 std::shared_ptr<Problem<poreNetworkIndex>> pnmProblem,
99 std::shared_ptr<CouplingMapper> couplingMapper,
101 {
102 couplingMapper_ = couplingMapper;
103 this->setSubProblems(std::make_tuple(freeFlowMassProblem, pnmProblem));
104 this->attachSolution(curSol);
105 }
106
107 // \}
108
112 // \{
113
117 template<std::size_t i, class Assembler>
118 void bindCouplingContext(Dune::index_constant<i> domainI, const Element<i>& element, const Assembler& assembler) const
119 {
120 bindCouplingContext_(domainI, element);
121 }
122
126 template<std::size_t i, std::size_t j, class LocalAssemblerI>
127 void updateCouplingContext(Dune::index_constant<i> domainI,
128 const LocalAssemblerI& localAssemblerI,
129 Dune::index_constant<j> domainJ,
130 std::size_t dofIdxGlobalJ,
131 const PrimaryVariables<j>& priVarsJ,
132 int pvIdxJ)
133 {
134 this->curSol(domainJ)[dofIdxGlobalJ][pvIdxJ] = priVarsJ[pvIdxJ];
135
136 // we need to update all solution-dependent components of the coupling context
137 // the dof of domain J has been deflected
138 // if domainJ == freeFlowMassIndex: update volvars in the PoreNetworkCouplingContext
139 // if domainJ == poreNetworkIndex: update volvars in the FreeFlowMassCouplingContext
140 // as the update is symmetric we only need to write this once
141 auto& context = std::get<1-j>(couplingContext_);
142 for (auto& c : context)
143 {
144 if (c.scv.dofIndex() == dofIdxGlobalJ)
145 {
146 const auto& problem = this->problem(domainJ);
147 const auto& element = problem.gridGeometry().element(c.scv.elementIndex());
148 const auto elemSol = elementSolution(element, this->curSol(domainJ), problem.gridGeometry());
149 c.volVars.update(elemSol, problem, element, c.scv);
150 }
151 }
152 }
153
154 // \}
155
159 const auto& couplingContext(Dune::index_constant<freeFlowMassIndex> domainI,
160 const FVElementGeometry<freeFlowMassIndex>& fvGeometry,
161 const SubControlVolumeFace<freeFlowMassIndex> scvf) const
162 {
163 auto& contexts = std::get<freeFlowMassIndex>(couplingContext_);
164 const auto eIdx = scvf.insideScvIdx();
165
166 if (contexts.empty() || couplingContextBoundForElement_[freeFlowMassIndex] != eIdx)
167 bindCouplingContext_(freeFlowMassIndex, fvGeometry);
168
169
170 return contexts[0];
171 }
172
176 const auto& couplingContext(Dune::index_constant<poreNetworkIndex> domainI,
177 const FVElementGeometry<poreNetworkIndex>& fvGeometry,
178 const SubControlVolume<poreNetworkIndex> scv) const
179 {
180 auto& contexts = std::get<poreNetworkIndex>(couplingContext_);
181 const auto eIdx = fvGeometry.gridGeometry().elementMapper().index(fvGeometry.element());
182
183 if (contexts.empty() || couplingContextBoundForElement_[poreNetworkIndex] != eIdx)
184 bindCouplingContext_(poreNetworkIndex, fvGeometry);
185
186 return contexts;
187 }
188
192 // \{
193
197 template<std::size_t i, std::size_t j>
198 const CouplingStencil& couplingStencil(Dune::index_constant<i> domainI,
199 const Element<i>& element,
200 Dune::index_constant<j> domainJ) const
201 {
202 const auto eIdx = this->problem(domainI).gridGeometry().elementMapper().index(element);
203 if constexpr (domainI == freeFlowMassIndex)
204 return couplingMapper_->freeFlowMassToPoreNetworkCouplingStencil(eIdx);
205 else
206 return couplingMapper_->poreNetworkToFreeFlowMassCouplingStencil(eIdx);
207 }
208
209 // \}
210
214 bool isCoupled(Dune::index_constant<freeFlowMassIndex> domainI,
215 const SubControlVolumeFace<freeFlowMassIndex>& scvf) const
216 { return couplingMapper_->isCoupledFreeFlowMassScvf(scvf.index()); }
217
223 bool isCoupled(Dune::index_constant<poreNetworkIndex> domainI,
224 const SubControlVolume<poreNetworkIndex>& scv) const
225 { return couplingMapper_->isCoupledPoreNetworkDof(scv.dofIndex()); }
226
227private:
231 template<std::size_t i>
232 bool isCoupledElement_(Dune::index_constant<i> domainI, std::size_t eIdx) const
233 {
234 if constexpr (domainI == freeFlowMassIndex)
235 return couplingMapper_->isCoupledFreeFlowElement(eIdx);
236 else
237 return couplingMapper_->isCoupledPoreNetworkElement(eIdx);
238 }
239
241 template<std::size_t i>
242 VolumeVariables<i> volVars_(Dune::index_constant<i> domainI,
243 const Element<i>& element,
244 const SubControlVolume<i>& scv) const
245 {
246 VolumeVariables<i> volVars;
247 const auto elemSol = elementSolution(
248 element, this->curSol(domainI), this->problem(domainI).gridGeometry()
249 );
250 volVars.update(elemSol, this->problem(domainI), element, scv);
251 return volVars;
252 }
253
257 template<std::size_t i>
258 void bindCouplingContext_(Dune::index_constant<i> domainI, const Element<i>& element) const
259 {
260 const auto fvGeometry = localView(this->problem(domainI).gridGeometry()).bindElement(element);
261 bindCouplingContext_(domainI, fvGeometry);
262 }
263
267 void bindCouplingContext_(Dune::index_constant<poreNetworkIndex> domainI, const FVElementGeometry<poreNetworkIndex>& fvGeometry) const
268 {
269 auto& context = std::get<domainI>(couplingContext_);
270 const auto eIdx = fvGeometry.gridGeometry().elementMapper().index(fvGeometry.element());
271
272 // do nothing if the element is already bound or not coupled to the other domain
273 if ((!context.empty() && couplingContextBoundForElement_[domainI] == eIdx) || !isCoupledElement_(domainI, eIdx))
274 return;
275
276 context.clear();
277 couplingContextBoundForElement_[domainI] = eIdx;
278
279 const auto& stencil = couplingStencil(poreNetworkIndex, fvGeometry.element(), freeFlowMassIndex);
280 const auto& freeFlowElements = couplingMapper_->pnmElementToFreeFlowElementsMap().at(eIdx);
281 auto ffFVGeometry = localView(this->problem(freeFlowMassIndex).gridGeometry());
282
283 for (const auto ffElementIdx : freeFlowElements)
284 {
285 const auto& ffElement = this->problem(freeFlowMassIndex).gridGeometry().element(ffElementIdx);
286 ffFVGeometry.bindElement(ffElement);
287
288 for (const auto& scvf : scvfs(ffFVGeometry))
289 {
290 if (couplingMapper_->isCoupledFreeFlowMassScvf(scvf.index()))
291 {
292 const auto& scv = ffFVGeometry.scv(scvf.insideScvIdx());
293 const auto dofIdx = scv.dofIndex();
294 if (std::any_of(stencil.begin(), stencil.end(), [&](const auto x){ return dofIdx == x; } ))
295 {
296 context.push_back({scv,
297 scvf,
298 volVars_(freeFlowMassIndex, ffElement, scv),
299 VelocityVector{}}
300 );
301 }
302 }
303 }
304 }
305 }
306
310 void bindCouplingContext_(Dune::index_constant<freeFlowMassIndex> domainI, const FVElementGeometry<freeFlowMassIndex>& fvGeometry) const
311 {
312 auto& context = std::get<domainI>(couplingContext_);
313 const auto eIdx = fvGeometry.gridGeometry().elementMapper().index(fvGeometry.element());
314
315 // do nothing if the element is already bound or not coupled to the other domain
316 if ((!context.empty() && couplingContextBoundForElement_[domainI] == eIdx) || !isCoupledElement_(domainI, eIdx))
317 return;
318
319 context.clear();
320 couplingContextBoundForElement_[domainI] = eIdx;
321
322
323 auto poreNetworkFVGeometry = localView(this->problem(poreNetworkIndex).gridGeometry());
324 const auto poreNetworkElemIdx = couplingMapper_->freeFlowElementToPNMElementMap().at(eIdx);
325 const auto& poreNetworkElement = this->problem(poreNetworkIndex).gridGeometry().element(poreNetworkElemIdx);
326 poreNetworkFVGeometry.bindElement(poreNetworkElement);
327
328 auto poreNetworkScv = [&]
329 {
330 SubControlVolume<poreNetworkIndex> result;
331 std::size_t counter = 0;
332 for (auto&& scv : scvs(poreNetworkFVGeometry))
333 {
334 if (couplingMapper_->isCoupledPoreNetworkDof(scv.dofIndex()))
335 {
336 result = scv;
337 ++counter;
338 }
339 }
340
341 if (counter > 1)
342 DUNE_THROW(Dune::InvalidStateException, "Only one pore per throat may be coupled");
343 else
344 return result;
345 }();
346
347 auto volVars = volVars_(poreNetworkIndex, poreNetworkElement, poreNetworkScv);
348
349 context.push_back({std::move(poreNetworkScv),
350 std::move(volVars),
351 VelocityVector{}}
352 );
353 }
354
355 mutable std::tuple<std::vector<FreeFlowMassCouplingContext>, std::vector<PoreNetworkCouplingContext>> couplingContext_;
356 mutable std::array<std::size_t, 2> couplingContextBoundForElement_;
357
358 std::shared_ptr<CouplingMapper> couplingMapper_;
359};
360
361} // end namespace Dumux
362
363#endif
The interface of the coupling manager for multi domain problems.
Definition: multidomain/couplingmanager.hh:48
void attachSolution(SolutionVectorStorage &curSol)
Attach a solution vector stored outside of this class.
Definition: multidomain/couplingmanager.hh:322
void setSubProblems(const std::tuple< std::shared_ptr< SubProblems >... > &problems)
set the pointers to the sub problems
Definition: multidomain/couplingmanager.hh:287
const Problem< i > & problem(Dune::index_constant< i > domainIdx) const
Return a reference to the sub problem.
Definition: multidomain/couplingmanager.hh:309
SubSolutionVector< i > & curSol(Dune::index_constant< i > domainIdx)
the solution vector of the subproblem
Definition: multidomain/couplingmanager.hh:338
typename Traits::template TupleOfSharedPtr< SubSolutionVector > SolutionVectorStorage
the type in which the solution vector is stored in the manager
Definition: multidomain/couplingmanager.hh:71
Coupling manager for free-flow mass and pore-network models.
Definition: multidomain/boundary/freeflowporenetwork/ffmassporenetwork/couplingmanager.hh:35
bool isCoupled(Dune::index_constant< poreNetworkIndex > domainI, const SubControlVolume< poreNetworkIndex > &scv) const
If the boundary entity is on a coupling boundary.
Definition: multidomain/boundary/freeflowporenetwork/ffmassporenetwork/couplingmanager.hh:223
static constexpr auto poreNetworkIndex
Definition: multidomain/boundary/freeflowporenetwork/ffmassporenetwork/couplingmanager.hh:41
typename ParentType::SolutionVectorStorage SolutionVectorStorage
Definition: multidomain/boundary/freeflowporenetwork/ffmassporenetwork/couplingmanager.hh:43
void bindCouplingContext(Dune::index_constant< i > domainI, const Element< i > &element, const Assembler &assembler) const
Methods to be accessed by the assembly.
Definition: multidomain/boundary/freeflowporenetwork/ffmassporenetwork/couplingmanager.hh:118
void init(std::shared_ptr< Problem< freeFlowMassIndex > > freeFlowMassProblem, std::shared_ptr< Problem< poreNetworkIndex > > pnmProblem, std::shared_ptr< CouplingMapper > couplingMapper, SolutionVectorStorage &curSol)
Methods to be accessed by main.
Definition: multidomain/boundary/freeflowporenetwork/ffmassporenetwork/couplingmanager.hh:97
const auto & couplingContext(Dune::index_constant< poreNetworkIndex > domainI, const FVElementGeometry< poreNetworkIndex > &fvGeometry, const SubControlVolume< poreNetworkIndex > scv) const
Access the coupling context needed for the PNM domain.
Definition: multidomain/boundary/freeflowporenetwork/ffmassporenetwork/couplingmanager.hh:176
const auto & couplingContext(Dune::index_constant< freeFlowMassIndex > domainI, const FVElementGeometry< freeFlowMassIndex > &fvGeometry, const SubControlVolumeFace< freeFlowMassIndex > scvf) const
Access the coupling context needed for the Stokes domain.
Definition: multidomain/boundary/freeflowporenetwork/ffmassporenetwork/couplingmanager.hh:159
static constexpr auto freeFlowMassIndex
Definition: multidomain/boundary/freeflowporenetwork/ffmassporenetwork/couplingmanager.hh:40
const CouplingStencil & couplingStencil(Dune::index_constant< i > domainI, const Element< i > &element, Dune::index_constant< j > domainJ) const
The coupling stencils.
Definition: multidomain/boundary/freeflowporenetwork/ffmassporenetwork/couplingmanager.hh:198
bool isCoupled(Dune::index_constant< freeFlowMassIndex > domainI, const SubControlVolumeFace< freeFlowMassIndex > &scvf) const
Returns whether a given scvf is coupled to the other domain.
Definition: multidomain/boundary/freeflowporenetwork/ffmassporenetwork/couplingmanager.hh:214
void updateCouplingContext(Dune::index_constant< i > domainI, const LocalAssemblerI &localAssemblerI, Dune::index_constant< j > domainJ, std::size_t dofIdxGlobalJ, const PrimaryVariables< j > &priVarsJ, int pvIdxJ)
Update the coupling context.
Definition: multidomain/boundary/freeflowporenetwork/ffmassporenetwork/couplingmanager.hh:127
Defines all properties used in Dumux.
GridCache::LocalView localView(const GridCache &gridCache)
Free function to get the local view of a grid cache object.
Definition: localview.hh:26
auto elementSolution(const Element &element, const SolutionVector &sol, const GridGeometry &gg) -> std::enable_if_t< GridGeometry::discMethod==DiscretizationMethods::cctpfa||GridGeometry::discMethod==DiscretizationMethods::ccmpfa, CCElementSolution< typename GridGeometry::LocalView, std::decay_t< decltype(std::declval< SolutionVector >()[0])> > >
Make an element solution for cell-centered schemes.
Definition: cellcentered/elementsolution.hh:101
typename GetProp< TypeTag, Property >::type GetPropType
get the type alias defined in the property
Definition: propertysystem.hh:296
The interface of the coupling manager for multi domain problems.
Definition: adapt.hh:17
The local element solution class for staggered methods.