version 3.10-dev
multidomain/facet/cellcentered/tpfa/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//
12#ifndef DUMUX_CCTPFA_FACETCOUPLING_MANAGER_HH
13#define DUMUX_CCTPFA_FACETCOUPLING_MANAGER_HH
14
15#include <algorithm>
16#include <cassert>
17
25
26namespace Dumux {
27
39template<class MDTraits, class CouplingMapper, std::size_t bulkDomainId, std::size_t lowDimDomainId>
40class FacetCouplingManager<MDTraits, CouplingMapper, bulkDomainId, lowDimDomainId, DiscretizationMethods::CCTpfa>
41: public virtual CouplingManager< MDTraits >
42{
44
45 // convenience aliases and instances of the two domain ids
46 using BulkIdType = typename MDTraits::template SubDomain<bulkDomainId>::Index;
47 using LowDimIdType = typename MDTraits::template SubDomain<lowDimDomainId>::Index;
48 static constexpr auto bulkId = BulkIdType();
49 static constexpr auto lowDimId = LowDimIdType();
50
51 // the sub-domain type tags
52 template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
53
54 // further types specific to the sub-problems
55 template<std::size_t id> using PrimaryVariables = GetPropType<SubDomainTypeTag<id>, Properties::PrimaryVariables>;
56 template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>;
57 template<std::size_t id> using NumEqVector = Dumux::NumEqVector<PrimaryVariables<id>>;
58 template<std::size_t id> using LocalResidual = typename MDTraits::template SubDomain<id>::LocalResidual;
59
60 template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
61 template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView;
62 template<std::size_t id> using SubControlVolume = typename GridGeometry<id>::SubControlVolume;
63 template<std::size_t id> using SubControlVolumeFace = typename GridGeometry<id>::SubControlVolumeFace;
64 template<std::size_t id> using GridView = typename GridGeometry<id>::GridView;
65 template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity;
66 template<std::size_t id> using GridIndexType = typename IndexTraits< GridView<id> >::GridIndex;
67
68 template<std::size_t id> using GridVariables = GetPropType<SubDomainTypeTag<id>, Properties::GridVariables>;
69 template<std::size_t id> using GridVolumeVariables = typename GridVariables<id>::GridVolumeVariables;
70 template<std::size_t id> using ElementVolumeVariables = typename GridVolumeVariables<id>::LocalView;
71 template<std::size_t id> using VolumeVariables = typename ElementVolumeVariables<id>::VolumeVariables;
72 template<std::size_t id> using GridFluxVariablesCache = typename GridVariables<id>::GridFluxVariablesCache;
73 template<std::size_t id> using ElementFluxVariablesCache = typename GridFluxVariablesCache<id>::LocalView;
74
75 // this currently does not work for some grid-wide caches being active
76 static_assert(!GetPropType<SubDomainTypeTag<bulkId>, Properties::GridVariables>::GridFluxVariablesCache::cachingEnabled,
77 "Grid flux variables caching currently not supported in the bulk domain of cc-facet coupling models");
78 static_assert(!GetPropType<SubDomainTypeTag<lowDimId>, Properties::GridVariables>::GridVolumeVariables::cachingEnabled,
79 "Grid volume variables caching currently not supported in the lower-dimensional domain of cc-facet coupling models");
80 static_assert(!GetPropType<SubDomainTypeTag<bulkId>, Properties::GridVariables>::GridVolumeVariables::cachingEnabled,
81 "Grid volume variables caching currently not supported in the bulk domain of cc-facet coupling models");
82
83 // extract corresponding grid ids from the mapper
84 static constexpr int bulkDim = GridView<bulkDomainId>::dimension;
85 static constexpr int lowDimDim = GridView<lowDimDomainId>::dimension;
86 static constexpr auto bulkGridId = CouplingMapper::template gridId<bulkDim>();
87 static constexpr auto lowDimGridId = CouplingMapper::template gridId<lowDimDim>();
88
89 static constexpr bool lowDimUsesBox = GridGeometry<lowDimId>::discMethod == DiscretizationMethods::box;
90
97 struct BulkCouplingContext
98 {
99 bool isSet;
100 GridIndexType< bulkId > elementIdx;
101 std::vector< FVElementGeometry<lowDimId> > lowDimFvGeometries;
102 std::vector< VolumeVariables<lowDimId> > lowDimVolVars;
103
104 void reset()
105 {
106 lowDimFvGeometries.clear();
107 lowDimVolVars.clear();
108 isSet = false;
109 }
110 };
111
123 struct LowDimCouplingContext
124 {
125 bool isSet;
126 GridIndexType< lowDimId > elementIdx;
127 std::unique_ptr< FVElementGeometry<bulkId> > bulkFvGeometry;
128 std::unique_ptr< ElementVolumeVariables<bulkId> > bulkElemVolVars;
129 std::unique_ptr< ElementFluxVariablesCache<bulkId> > bulkElemFluxVarsCache;
130 std::unique_ptr< LocalResidual<bulkId> > bulkLocalResidual;
131
132 void reset()
133 {
134 bulkFvGeometry.reset(nullptr);
135 bulkElemVolVars.reset(nullptr);
136 bulkElemFluxVarsCache.reset(nullptr);
137 isSet = false;
138 }
139 };
140
141public:
142
144 template<std::size_t i, std::size_t j = (i == bulkId) ? lowDimId : bulkId>
145 using CouplingStencilType = typename CouplingMapper::template Stencil< CouplingMapper::template gridId<GridView<j>::dimension>() >;
146
148 using SolutionVector = typename MDTraits::SolutionVector;
149
158 void init(std::shared_ptr< Problem<bulkId> > bulkProblem,
159 std::shared_ptr< Problem<lowDimId> > lowDimProblem,
160 std::shared_ptr< CouplingMapper > couplingMapper,
161 const SolutionVector& curSol)
162 {
163 couplingMapperPtr_ = couplingMapper;
164
165 // set the sub problems
166 this->setSubProblem(bulkProblem, bulkId);
167 this->setSubProblem(lowDimProblem, lowDimId);
168
169 // copy the solution vector
170 ParentType::updateSolution(curSol);
171
172 // determine all bulk elements/scvfs that couple to low dim elements
173 bulkElemIsCoupled_.assign(bulkProblem->gridGeometry().gridView().size(0), false);
174 bulkScvfIsCoupled_.assign(bulkProblem->gridGeometry().numScvf(), false);
175
176 const auto& bulkMap = couplingMapperPtr_->couplingMap(bulkGridId, lowDimGridId);
177 for (const auto& entry : bulkMap)
178 {
179 bulkElemIsCoupled_[entry.first] = true;
180 for (const auto& couplingEntry : entry.second.dofToCouplingScvfMap)
181 for (const auto& scvfIdx : couplingEntry.second)
182 bulkScvfIsCoupled_[scvfIdx] = true;
183 }
184 }
185
190 const Element<bulkId>& element,
191 LowDimIdType domainJ) const
192 {
193 const auto eIdx = this->problem(bulkId).gridGeometry().elementMapper().index(element);
194
195 if (bulkElemIsCoupled_[eIdx])
196 {
197 const auto& map = couplingMapperPtr_->couplingMap(bulkGridId, lowDimGridId);
198 auto it = map.find(eIdx);
199 assert(it != map.end());
200 return it->second.couplingStencil;
201 }
202
203 return getEmptyStencil(lowDimId);
204 }
205
210 const Element<lowDimId>& element,
211 BulkIdType domainJ) const
212 {
213 const auto eIdx = this->problem(lowDimId).gridGeometry().elementMapper().index(element);
214
215 const auto& map = couplingMapperPtr_->couplingMap(lowDimGridId, bulkGridId);
216 auto it = map.find(eIdx);
217 if (it != map.end()) return it->second.couplingStencil;
218 else return getEmptyStencil(bulkId);
219 }
220
224 bool isCoupled(const Element<bulkId>& element,
225 const SubControlVolumeFace<bulkId>& scvf) const
226 { return bulkScvfIsCoupled_[scvf.index()]; }
227
232 bool isOnInteriorBoundary(const Element<bulkId>& element,
233 const SubControlVolumeFace<bulkId>& scvf) const
234 { return isCoupled(element, scvf); }
235
239 const VolumeVariables<lowDimId>& getLowDimVolVars(const Element<bulkId>& element,
240 const SubControlVolumeFace<bulkId>& scvf) const
241 {
242 assert(bulkContext_.isSet);
243
244 const auto lowDimElemIdx = getLowDimElementIndex(element, scvf);
245 const auto& map = couplingMapperPtr_->couplingMap(bulkGridId, lowDimGridId);
246 const auto& s = map.find(bulkContext_.elementIdx)->second.couplingElementStencil;
247 const auto& idxInContext = std::distance( s.begin(), std::find(s.begin(), s.end(), lowDimElemIdx) );
248
249 assert(std::find(s.begin(), s.end(), lowDimElemIdx) != s.end());
250 return bulkContext_.lowDimVolVars[idxInContext];
251 }
252
256 const Element<lowDimId> getLowDimElement(const Element<bulkId>& element,
257 const SubControlVolumeFace<bulkId>& scvf) const
258 {
259 const auto lowDimElemIdx = getLowDimElementIndex(element, scvf);
260 return this->problem(lowDimId).gridGeometry().element(lowDimElemIdx);
261 }
262
266 const GridIndexType<lowDimId> getLowDimElementIndex(const Element<bulkId>& element,
267 const SubControlVolumeFace<bulkId>& scvf) const
268 {
269 assert(bulkScvfIsCoupled_[scvf.index()]);
270
271 const auto& map = couplingMapperPtr_->couplingMap(bulkGridId, lowDimGridId);
272 const auto& couplingData = map.at(scvf.insideScvIdx());
273
274 // search the low dim element idx this scvf is embedded in
275 auto it = std::find_if( couplingData.elementToScvfMap.begin(),
276 couplingData.elementToScvfMap.end(),
277 [&scvf] (auto& dataPair)
278 {
279 const auto& scvfs = dataPair.second;
280 return std::find(scvfs.begin(), scvfs.end(), scvf.index()) != scvfs.end();
281 } );
282
283 assert(it != couplingData.elementToScvfMap.end());
284 return it->first;
285 }
286
292 template< class BulkLocalAssembler >
293 typename LocalResidual<bulkId>::ElementResidualVector
295 const BulkLocalAssembler& bulkLocalAssembler,
296 LowDimIdType,
297 GridIndexType<lowDimId> dofIdxGlobalJ)
298 {
299 const auto& map = couplingMapperPtr_->couplingMap(bulkGridId, lowDimGridId);
300
301 assert(bulkContext_.isSet);
302 assert(bulkElemIsCoupled_[bulkContext_.elementIdx]);
303 assert(map.find(bulkContext_.elementIdx) != map.end());
304 assert(bulkContext_.elementIdx == this->problem(bulkId).gridGeometry().elementMapper().index(bulkLocalAssembler.element()));
305
306 typename LocalResidual<bulkId>::ElementResidualVector res(1);
307 res = 0.0;
308 res[0] = evalBulkFluxes(bulkLocalAssembler.element(),
309 bulkLocalAssembler.fvGeometry(),
310 bulkLocalAssembler.curElemVolVars(),
311 bulkLocalAssembler.elemFluxVarsCache(),
312 bulkLocalAssembler.localResidual(),
313 map.find(bulkContext_.elementIdx)->second.dofToCouplingScvfMap.at(dofIdxGlobalJ));
314 return res;
315 }
316
327 template< class LowDimLocalAssembler >
328 typename LocalResidual<lowDimId>::ElementResidualVector
330 const LowDimLocalAssembler& lowDimLocalAssembler,
331 BulkIdType,
332 GridIndexType<bulkId> dofIdxGlobalJ)
333 { return evalCouplingResidual(lowDimId, lowDimLocalAssembler, bulkId); }
334
340 template< class LowDimLocalAssembler >
341 typename LocalResidual<lowDimId>::ElementResidualVector
342 evalCouplingResidual(LowDimIdType, const LowDimLocalAssembler& lowDimLocalAssembler, BulkIdType)
343 {
344 // make sure this is called for the element for which the context was set
345 assert(lowDimContext_.isSet);
346 assert(this->problem(lowDimId).gridGeometry().elementMapper().index(lowDimLocalAssembler.element()) == lowDimContext_.elementIdx);
347
348 // evaluate sources for the first scv
349 // the sources are element-wise & scv-independent since we use tpfa in bulk domain
350 const auto source = evalSourcesFromBulk(lowDimLocalAssembler.element(),
351 lowDimLocalAssembler.fvGeometry(),
352 lowDimLocalAssembler.curElemVolVars(),
353 *scvs(lowDimLocalAssembler.fvGeometry()).begin());
354
355 // fill element residual vector with the sources
356 typename LocalResidual<lowDimId>::ElementResidualVector res(lowDimLocalAssembler.fvGeometry().numScv());
357 res = 0.0;
358 for (const auto& scv : scvs(lowDimLocalAssembler.fvGeometry()))
359 res[scv.localDofIndex()] -= source;
360
361 return res;
362 }
363
367 NumEqVector<lowDimId> evalSourcesFromBulk(const Element<lowDimId>& element,
368 const FVElementGeometry<lowDimId>& fvGeometry,
369 const ElementVolumeVariables<lowDimId>& elemVolVars,
370 const SubControlVolume<lowDimId>& scv) const
371 {
372 // make sure the this is called for the element of the context
373 assert(this->problem(lowDimId).gridGeometry().elementMapper().index(element) == lowDimContext_.elementIdx);
374
375 NumEqVector<lowDimId> sources(0.0);
376
377 const auto& map = couplingMapperPtr_->couplingMap(lowDimGridId, bulkGridId);
378 auto it = map.find(lowDimContext_.elementIdx);
379 if (it == map.end())
380 return sources;
381
382 assert(lowDimContext_.isSet);
383 for (const auto& embedment : it->second.embedments)
384 sources += evalBulkFluxes(this->problem(bulkId).gridGeometry().element(embedment.first),
385 *lowDimContext_.bulkFvGeometry,
386 *lowDimContext_.bulkElemVolVars,
387 *lowDimContext_.bulkElemFluxVarsCache,
388 *lowDimContext_.bulkLocalResidual,
389 embedment.second);
390
391 // if lowdim domain uses box, we distribute the sources equally among the scvs
392 if (lowDimUsesBox)
393 sources /= fvGeometry.numScv();
394
395 return sources;
396 }
397
403 template< class Assembler >
404 void bindCouplingContext(BulkIdType, const Element<bulkId>& element, const Assembler& assembler)
405 {
406 // clear context
407 bulkContext_.reset();
408
409 // set index in context in any case
410 const auto bulkElemIdx = this->problem(bulkId).gridGeometry().elementMapper().index(element);
411 bulkContext_.elementIdx = bulkElemIdx;
412
413 // if element is coupled, actually set the context
414 if (bulkElemIsCoupled_[bulkElemIdx])
415 {
416 const auto& map = couplingMapperPtr_->couplingMap(bulkGridId, lowDimGridId);
417
418 auto it = map.find(bulkElemIdx); assert(it != map.end());
419 const auto& elementStencil = it->second.couplingElementStencil;
420 bulkContext_.lowDimFvGeometries.reserve(elementStencil.size());
421 bulkContext_.lowDimVolVars.reserve(elementStencil.size());
422
423 const auto& ldGridGeometry = this->problem(lowDimId).gridGeometry();
424 auto fvGeom = localView(ldGridGeometry);
425 for (const auto lowDimElemIdx : elementStencil)
426 {
427 const auto& ldSol = assembler.isImplicit() ? this->curSol(lowDimId) : assembler.prevSol()[lowDimId];
428 const auto& ldProblem = this->problem(lowDimId);
429
430 const auto elemJ = ldGridGeometry.element(lowDimElemIdx);
431 fvGeom.bindElement(elemJ);
432
433 VolumeVariables<lowDimId> volVars;
434
435 // if low dim domain uses the box scheme, we have to create interpolated vol vars
436 if (lowDimUsesBox)
437 {
438 const auto elemGeom = elemJ.geometry();
439 FacetCoupling::makeInterpolatedVolVars(volVars, ldProblem, ldSol, fvGeom, elemJ, elemGeom, elemGeom.center());
440 }
441 // if low dim domain uses a cc scheme we can directly update the vol vars
442 else
443 volVars.update( elementSolution(elemJ, ldSol, ldGridGeometry),
444 ldProblem,
445 elemJ,
446 fvGeom.scv(lowDimElemIdx) );
447
448 bulkContext_.isSet = true;
449 bulkContext_.lowDimFvGeometries.emplace_back( std::move(fvGeom) );
450 bulkContext_.lowDimVolVars.emplace_back( std::move(volVars) );
451 }
452 }
453 }
454
464 template< class Assembler >
465 void bindCouplingContext(LowDimIdType, const Element<lowDimId>& element, const Assembler& assembler)
466 {
467 // reset contexts
468 bulkContext_.reset();
469 lowDimContext_.reset();
470
471 // set index in context in any case
472 const auto lowDimElemIdx = this->problem(lowDimId).gridGeometry().elementMapper().index(element);
473 lowDimContext_.elementIdx = lowDimElemIdx;
474
475 const auto& map = couplingMapperPtr_->couplingMap(lowDimGridId, bulkGridId);
476 auto it = map.find(lowDimElemIdx);
477
478 // if element is coupled, actually set the context
479 if (it != map.end())
480 {
481 // first bind the low dim context for the first neighboring bulk element
482 const auto& bulkGridGeom = this->problem(bulkId).gridGeometry();
483 const auto bulkElem = bulkGridGeom.element(it->second.embedments[0].first);
484 bindCouplingContext(bulkId, bulkElem, assembler);
485
486 // evaluate variables on old/new time level depending on time disc scheme
487 const auto& bulkSol = assembler.isImplicit() ? this->curSol(bulkId) : assembler.prevSol()[bulkId];
488
489 // then simply bind the local views of that first neighbor
490 auto bulkFvGeom = localView(bulkGridGeom).bind(bulkElem);
491 auto bulkElemVolVars = assembler.isImplicit() ? localView(assembler.gridVariables(bulkId).curGridVolVars()).bind(bulkElem, bulkFvGeom, bulkSol)
492 : localView(assembler.gridVariables(bulkId).prevGridVolVars()).bind(bulkElem, bulkFvGeom, bulkSol);
493 auto bulkElemFluxVarsCache = localView(assembler.gridVariables(bulkId).gridFluxVarsCache()).bind(bulkElem, bulkFvGeom, bulkElemVolVars);
494
495 lowDimContext_.isSet = true;
496 lowDimContext_.bulkFvGeometry = std::make_unique< FVElementGeometry<bulkId> >( std::move(bulkFvGeom) );
497 lowDimContext_.bulkElemVolVars = std::make_unique< ElementVolumeVariables<bulkId> >( std::move(bulkElemVolVars) );
498 lowDimContext_.bulkElemFluxVarsCache = std::make_unique< ElementFluxVariablesCache<bulkId> >( std::move(bulkElemFluxVarsCache) );
499 lowDimContext_.bulkLocalResidual = std::make_unique< LocalResidual<bulkId> >(assembler.localResidual(bulkId));
500 }
501 }
502
507 template< class BulkLocalAssembler >
508 void updateCouplingContext(BulkIdType domainI,
509 const BulkLocalAssembler& bulkLocalAssembler,
510 LowDimIdType domainJ,
511 GridIndexType<lowDimId> dofIdxGlobalJ,
512 const PrimaryVariables<lowDimId>& priVarsJ,
513 unsigned int pvIdxJ)
514 {
515 // communicate deflected solution
516 ParentType::updateCouplingContext(domainI, bulkLocalAssembler, domainJ, dofIdxGlobalJ, priVarsJ, pvIdxJ);
517
518 // Since coupling only occurs via the fluxes, the context does not
519 // have to be updated in explicit time discretization schemes, where
520 // they are strictly evaluated on the old time level
521 if (!bulkLocalAssembler.isImplicit())
522 return;
523
524 // skip the rest if context is empty
525 if (bulkContext_.isSet)
526 {
527 const auto& map = couplingMapperPtr_->couplingMap(bulkGridId, lowDimGridId);
528 const auto& couplingElemStencil = map.find(bulkContext_.elementIdx)->second.couplingElementStencil;
529 const auto& ldSol = this->curSol(lowDimId);
530 const auto& ldProblem = this->problem(lowDimId);
531 const auto& ldGridGeometry = this->problem(lowDimId).gridGeometry();
532
533 // find the low-dim elements in coupling stencil, where this dof is contained in
534 const auto couplingElements = [&] ()
535 {
536 if (lowDimUsesBox)
537 {
538 std::vector< Element<lowDimId> > lowDimElems;
539 std::for_each( couplingElemStencil.begin(), couplingElemStencil.end(),
540 [&] (auto lowDimElemIdx)
541 {
542 auto element = ldGridGeometry.element(lowDimElemIdx);
543 for (int i = 0; i < element.geometry().corners(); ++i)
544 {
545 const auto dofIdx = ldGridGeometry.vertexMapper().subIndex(element, i, lowDimDim);
546 if (dofIdxGlobalJ == dofIdx) { lowDimElems.emplace_back( std::move(element) ); break; }
547 }
548 } );
549 return lowDimElems;
550 }
551 // dof index = element index for cc schemes
552 else
553 return std::vector<Element<lowDimId>>( {ldGridGeometry.element(dofIdxGlobalJ)} );
554 } ();
555
556 // update all necessary vol vars in context
557 for (const auto& element : couplingElements)
558 {
559 // find index in coupling context
560 const auto eIdxGlobal = ldGridGeometry.elementMapper().index(element);
561 auto it = std::find(couplingElemStencil.begin(), couplingElemStencil.end(), eIdxGlobal);
562 const auto idxInContext = std::distance(couplingElemStencil.begin(), it);
563 assert(it != couplingElemStencil.end());
564
565 auto& volVars = bulkContext_.lowDimVolVars[idxInContext];
566 const auto& fvGeom = bulkContext_.lowDimFvGeometries[idxInContext];
567 // if low dim domain uses the box scheme, we have to create interpolated vol vars
568 if (lowDimUsesBox)
569 {
570 const auto elemGeom = element.geometry();
571 FacetCoupling::makeInterpolatedVolVars(volVars, ldProblem, ldSol, fvGeom, element, elemGeom, elemGeom.center());
572 }
573 // if low dim domain uses a cc scheme we can directly update the vol vars
574 else
575 volVars.update( elementSolution(element, ldSol, ldGridGeometry),
576 ldProblem,
577 element,
578 fvGeom.scv(eIdxGlobal) );
579 }
580 }
581 }
582
587 template< class BulkLocalAssembler >
588 void updateCouplingContext(BulkIdType domainI,
589 const BulkLocalAssembler& bulkLocalAssembler,
590 BulkIdType domainJ,
591 GridIndexType<bulkId> dofIdxGlobalJ,
592 const PrimaryVariables<bulkId>& priVarsJ,
593 unsigned int pvIdxJ)
594 {
595 // communicate deflected solution
596 ParentType::updateCouplingContext(domainI, bulkLocalAssembler, domainJ, dofIdxGlobalJ, priVarsJ, pvIdxJ);
597 }
598
604 template< class LowDimLocalAssembler >
605 void updateCouplingContext(LowDimIdType domainI,
606 const LowDimLocalAssembler& lowDimLocalAssembler,
607 BulkIdType domainJ,
608 GridIndexType<bulkId> dofIdxGlobalJ,
609 const PrimaryVariables<bulkId>& priVarsJ,
610 unsigned int pvIdxJ)
611 {
612 // communicate deflected solution
613 ParentType::updateCouplingContext(domainI, lowDimLocalAssembler, domainJ, dofIdxGlobalJ, priVarsJ, pvIdxJ);
614
615 // Since coupling only occurs via the fluxes, the context does not
616 // have to be updated in explicit time discretization schemes, where
617 // they are strictly evaluated on the old time level
618 if (!lowDimLocalAssembler.isImplicit())
619 return;
620
621 // skip the rest if context is empty
622 if (lowDimContext_.isSet)
623 {
624 assert(lowDimContext_.elementIdx == this->problem(lowDimId).gridGeometry().elementMapper().index(lowDimLocalAssembler.element()));
625
626 // since we use cc scheme in bulk domain: dof index = element index
627 const auto& bulkGridGeom = this->problem(bulkId).gridGeometry();
628 const auto elementJ = bulkGridGeom.element(dofIdxGlobalJ);
629
630 // update corresponding vol vars in context
631 const auto& scv = lowDimContext_.bulkFvGeometry->scv(dofIdxGlobalJ);
632 const auto elemSol = elementSolution(elementJ, this->curSol(bulkId), bulkGridGeom);
633 (*lowDimContext_.bulkElemVolVars)[dofIdxGlobalJ].update(elemSol, this->problem(bulkId), elementJ, scv);
634
635 // update the element flux variables cache (tij might be solution-dependent)
636 if (dofIdxGlobalJ == bulkContext_.elementIdx)
637 lowDimContext_.bulkElemFluxVarsCache->update( elementJ, *lowDimContext_.bulkFvGeometry, *lowDimContext_.bulkElemVolVars);
638 else
639 lowDimContext_.bulkElemFluxVarsCache->update( this->problem(bulkId).gridGeometry().element(bulkContext_.elementIdx),
640 *lowDimContext_.bulkFvGeometry,
641 *lowDimContext_.bulkElemVolVars );
642 }
643 }
644
651 template< class LowDimLocalAssembler >
652 void updateCouplingContext(LowDimIdType domainI,
653 const LowDimLocalAssembler& lowDimLocalAssembler,
654 LowDimIdType domainJ,
655 GridIndexType<lowDimId> dofIdxGlobalJ,
656 const PrimaryVariables<lowDimId>& priVarsJ,
657 unsigned int pvIdxJ)
658 {
659 // communicate deflected solution
660 ParentType::updateCouplingContext(domainI, lowDimLocalAssembler, domainJ, dofIdxGlobalJ, priVarsJ, pvIdxJ);
661
662 // Since coupling only occurs via the fluxes, the context does not
663 // have to be updated in explicit time discretization schemes, where
664 // they are strictly evaluated on the old time level
665 if (!lowDimLocalAssembler.isImplicit())
666 return;
667
668 // skip the rest if context is empty
669 if (lowDimContext_.isSet)
670 {
671 const auto& ldSol = this->curSol(lowDimId);
672 const auto& ldProblem = this->problem(lowDimId);
673 const auto& ldGridGeometry = this->problem(lowDimId).gridGeometry();
674
675 assert(bulkContext_.isSet);
676 assert(lowDimContext_.elementIdx == ldGridGeometry.elementMapper().index(lowDimLocalAssembler.element()));
677
678 // update the corresponding vol vars in the bulk context
679 const auto& bulkMap = couplingMapperPtr_->couplingMap(bulkGridId, lowDimGridId);
680 const auto& couplingElementStencil = bulkMap.find(bulkContext_.elementIdx)->second.couplingElementStencil;
681 auto it = std::find(couplingElementStencil.begin(), couplingElementStencil.end(), lowDimContext_.elementIdx);
682 assert(it != couplingElementStencil.end());
683 const auto idxInContext = std::distance(couplingElementStencil.begin(), it);
684
685 auto& volVars = bulkContext_.lowDimVolVars[idxInContext];
686 const auto& fvGeom = bulkContext_.lowDimFvGeometries[idxInContext];
687 const auto& element = lowDimLocalAssembler.element();
688 // if low dim domain uses the box scheme, we have to create interpolated vol vars
689 if (lowDimUsesBox)
690 {
691 const auto elemGeom = element.geometry();
692 FacetCoupling::makeInterpolatedVolVars(volVars, ldProblem, ldSol, fvGeom, element, elemGeom, elemGeom.center());
693 }
694 // if low dim domain uses a cc scheme we can directly update the vol vars
695 else
696 volVars.update( elementSolution(element, ldSol, ldGridGeometry),
697 ldProblem,
698 element,
699 fvGeom.scv(lowDimContext_.elementIdx) );
700
701 // update the element flux variables cache (tij depend on low dim values in context)
702 const auto contextElem = this->problem(bulkId).gridGeometry().element(bulkContext_.elementIdx);
703 lowDimContext_.bulkElemFluxVarsCache->update(contextElem, *lowDimContext_.bulkFvGeometry, *lowDimContext_.bulkElemVolVars);
704 }
705 }
706
708 using ParentType::updateCoupledVariables;
709
714 template< class BulkLocalAssembler, class UpdatableFluxVarCache >
715 void updateCoupledVariables(BulkIdType domainI,
716 const BulkLocalAssembler& bulkLocalAssembler,
717 ElementVolumeVariables<bulkId>& elemVolVars,
718 UpdatableFluxVarCache& fluxVarsCache)
719 {
720 // update transmissibilities after low dim context has changed (implicit only)
721 if (bulkLocalAssembler.isImplicit())
722 fluxVarsCache.update(bulkLocalAssembler.element(),
723 bulkLocalAssembler.fvGeometry(),
724 bulkLocalAssembler.curElemVolVars());
725 }
726
731 template< class BulkLocalAssembler, class UpdatableFluxVarCache >
732 void updateCoupledVariables(BulkIdType domainI,
733 const BulkLocalAssembler& bulkLocalAssembler,
734 GridVolumeVariables<bulkId>& gridVolVars,
735 UpdatableFluxVarCache& fluxVarsCache)
736 {
737 // update transmissibilities after low dim context has changed (implicit only)
738 if (bulkLocalAssembler.isImplicit())
739 {
740 const auto elemVolVars = localView(gridVolVars).bind(bulkLocalAssembler.element(), bulkLocalAssembler.fvGeometry(), this->curSol(bulkId));
741 fluxVarsCache.update(bulkLocalAssembler.element(), bulkLocalAssembler.fvGeometry(), elemVolVars);
742 }
743 }
744
746 template<std::size_t id, std::enable_if_t<(id == bulkId || id == lowDimId), int> = 0>
747 const typename CouplingMapper::template Stencil<id>&
748 getEmptyStencil(Dune::index_constant<id>) const
749 { return std::get<(id == bulkId ? 0 : 1)>(emptyStencilTuple_); }
750
751protected:
753 const BulkCouplingContext& bulkCouplingContext() const { return bulkContext_; }
754 const LowDimCouplingContext& lowDimCouplingContext() const { return lowDimContext_; }
755
757 BulkCouplingContext& bulkCouplingContext() { return bulkContext_; }
758 LowDimCouplingContext& lowDimCouplingContext() { return lowDimContext_; }
759
761 template<class BulkScvfIndices>
762 NumEqVector<bulkId> evalBulkFluxes(const Element<bulkId>& elementI,
763 const FVElementGeometry<bulkId>& fvGeometry,
764 const ElementVolumeVariables<bulkId>& elemVolVars,
765 const ElementFluxVariablesCache<bulkId>& elemFluxVarsCache,
766 const LocalResidual<bulkId>& localResidual,
767 const BulkScvfIndices& scvfIndices) const
768 {
769
770 NumEqVector<bulkId> coupledFluxes(0.0);
771 for (const auto& scvfIdx : scvfIndices)
772 coupledFluxes += localResidual.evalFlux(this->problem(bulkId),
773 elementI,
774 fvGeometry,
775 elemVolVars,
776 elemFluxVarsCache,
777 fvGeometry.scvf(scvfIdx));
778 return coupledFluxes;
779 }
780
781private:
782 std::shared_ptr<CouplingMapper> couplingMapperPtr_;
783
786 std::vector<bool> bulkElemIsCoupled_;
787 std::vector<bool> bulkScvfIsCoupled_;
788
790 using BulkStencil = typename CouplingMapper::template Stencil<bulkId>;
791 using LowDimStencil = typename CouplingMapper::template Stencil<lowDimId>;
792 std::tuple<BulkStencil, LowDimStencil> emptyStencilTuple_;
793
795 BulkCouplingContext bulkContext_;
796 LowDimCouplingContext lowDimContext_;
797};
798
799} // end namespace Dumux
800
801#endif
The interface of the coupling manager for multi domain problems.
Definition: multidomain/couplingmanager.hh:37
void bindCouplingContext(LowDimIdType, const Element< lowDimId > &element, const Assembler &assembler)
For the assembly of the element residual of a bulk domain element we need to prepare the local views ...
Definition: multidomain/facet/cellcentered/tpfa/couplingmanager.hh:465
const CouplingStencilType< bulkId > & couplingStencil(BulkIdType domainI, const Element< bulkId > &element, LowDimIdType domainJ) const
The coupling stencil of a given bulk domain element.
Definition: multidomain/facet/cellcentered/tpfa/couplingmanager.hh:189
typename CouplingMapper::template Stencil< CouplingMapper::template gridId< GridView< j >::dimension >() > CouplingStencilType
types used for coupling stencils
Definition: multidomain/facet/cellcentered/tpfa/couplingmanager.hh:145
const CouplingStencilType< lowDimId > & couplingStencil(LowDimIdType domainI, const Element< lowDimId > &element, BulkIdType domainJ) const
The coupling stencil of the lower-dimensional domain with the bulk domain.
Definition: multidomain/facet/cellcentered/tpfa/couplingmanager.hh:209
BulkCouplingContext & bulkCouplingContext()
Return references to the bulk coupling contexts.
Definition: multidomain/facet/cellcentered/tpfa/couplingmanager.hh:757
bool isCoupled(const Element< bulkId > &element, const SubControlVolumeFace< bulkId > &scvf) const
returns true if a bulk scvf flux depends on data in the facet domain.
Definition: multidomain/facet/cellcentered/tpfa/couplingmanager.hh:224
NumEqVector< bulkId > evalBulkFluxes(const Element< bulkId > &elementI, const FVElementGeometry< bulkId > &fvGeometry, const ElementVolumeVariables< bulkId > &elemVolVars, const ElementFluxVariablesCache< bulkId > &elemFluxVarsCache, const LocalResidual< bulkId > &localResidual, const BulkScvfIndices &scvfIndices) const
evaluates the bulk-facet exchange fluxes for a given facet element
Definition: multidomain/facet/cellcentered/tpfa/couplingmanager.hh:762
void updateCoupledVariables(BulkIdType domainI, const BulkLocalAssembler &bulkLocalAssembler, ElementVolumeVariables< bulkId > &elemVolVars, UpdatableFluxVarCache &fluxVarsCache)
Update the transmissibilities in the bulk domain after the coupling context changed.
Definition: multidomain/facet/cellcentered/tpfa/couplingmanager.hh:715
LocalResidual< bulkId >::ElementResidualVector evalCouplingResidual(BulkIdType, const BulkLocalAssembler &bulkLocalAssembler, LowDimIdType, GridIndexType< lowDimId > dofIdxGlobalJ)
Evaluates the coupling element residual of a bulk domain element with respect to a dof in the lower-d...
Definition: multidomain/facet/cellcentered/tpfa/couplingmanager.hh:294
const VolumeVariables< lowDimId > & getLowDimVolVars(const Element< bulkId > &element, const SubControlVolumeFace< bulkId > &scvf) const
returns the vol vars of a lower-dimensional element coinciding with a bulk scvf.
Definition: multidomain/facet/cellcentered/tpfa/couplingmanager.hh:239
LowDimCouplingContext & lowDimCouplingContext()
Definition: multidomain/facet/cellcentered/tpfa/couplingmanager.hh:758
typename MDTraits::SolutionVector SolutionVector
the type of the solution vector
Definition: multidomain/facet/cellcentered/tpfa/couplingmanager.hh:148
void init(std::shared_ptr< Problem< bulkId > > bulkProblem, std::shared_ptr< Problem< lowDimId > > lowDimProblem, std::shared_ptr< CouplingMapper > couplingMapper, const SolutionVector &curSol)
Initialize the coupling manager.
Definition: multidomain/facet/cellcentered/tpfa/couplingmanager.hh:158
const Element< lowDimId > getLowDimElement(const Element< bulkId > &element, const SubControlVolumeFace< bulkId > &scvf) const
returns the lower-dimensional element coinciding with a bulk scvf.
Definition: multidomain/facet/cellcentered/tpfa/couplingmanager.hh:256
LocalResidual< lowDimId >::ElementResidualVector evalCouplingResidual(LowDimIdType, const LowDimLocalAssembler &lowDimLocalAssembler, BulkIdType)
Evaluates the coupling element residual of a lower-dimensional domain element with respect to a dof i...
Definition: multidomain/facet/cellcentered/tpfa/couplingmanager.hh:342
void updateCoupledVariables(BulkIdType domainI, const BulkLocalAssembler &bulkLocalAssembler, GridVolumeVariables< bulkId > &gridVolVars, UpdatableFluxVarCache &fluxVarsCache)
Update the transmissibilities in the bulk domain after the coupling context changed.
Definition: multidomain/facet/cellcentered/tpfa/couplingmanager.hh:732
const LowDimCouplingContext & lowDimCouplingContext() const
Definition: multidomain/facet/cellcentered/tpfa/couplingmanager.hh:754
void updateCouplingContext(BulkIdType domainI, const BulkLocalAssembler &bulkLocalAssembler, LowDimIdType domainJ, GridIndexType< lowDimId > dofIdxGlobalJ, const PrimaryVariables< lowDimId > &priVarsJ, unsigned int pvIdxJ)
After deflecting the solution of the lower-dimensional domain, we have to update the element volume v...
Definition: multidomain/facet/cellcentered/tpfa/couplingmanager.hh:508
bool isOnInteriorBoundary(const Element< bulkId > &element, const SubControlVolumeFace< bulkId > &scvf) const
returns true if a bulk scvf coincides with a facet element.
Definition: multidomain/facet/cellcentered/tpfa/couplingmanager.hh:232
void updateCouplingContext(LowDimIdType domainI, const LowDimLocalAssembler &lowDimLocalAssembler, LowDimIdType domainJ, GridIndexType< lowDimId > dofIdxGlobalJ, const PrimaryVariables< lowDimId > &priVarsJ, unsigned int pvIdxJ)
After deflecting the solution of the lower-dimensional domain has been deflected during the assembly ...
Definition: multidomain/facet/cellcentered/tpfa/couplingmanager.hh:652
void updateCouplingContext(LowDimIdType domainI, const LowDimLocalAssembler &lowDimLocalAssembler, BulkIdType domainJ, GridIndexType< bulkId > dofIdxGlobalJ, const PrimaryVariables< bulkId > &priVarsJ, unsigned int pvIdxJ)
After deflecting the solution of the bulk domain, we have to update the element volume variables and ...
Definition: multidomain/facet/cellcentered/tpfa/couplingmanager.hh:605
const GridIndexType< lowDimId > getLowDimElementIndex(const Element< bulkId > &element, const SubControlVolumeFace< bulkId > &scvf) const
returns the index of the lower-dimensional element coinciding with a bulk scvf.
Definition: multidomain/facet/cellcentered/tpfa/couplingmanager.hh:266
NumEqVector< lowDimId > evalSourcesFromBulk(const Element< lowDimId > &element, const FVElementGeometry< lowDimId > &fvGeometry, const ElementVolumeVariables< lowDimId > &elemVolVars, const SubControlVolume< lowDimId > &scv) const
Computes the sources in a lower-dimensional element stemming from the bulk domain.
Definition: multidomain/facet/cellcentered/tpfa/couplingmanager.hh:367
const BulkCouplingContext & bulkCouplingContext() const
Return const references to the bulk coupling contexts.
Definition: multidomain/facet/cellcentered/tpfa/couplingmanager.hh:753
const CouplingMapper::template Stencil< id > & getEmptyStencil(Dune::index_constant< id >) const
Empty stencil to be returned for elements that aren't coupled.
Definition: multidomain/facet/cellcentered/tpfa/couplingmanager.hh:748
LocalResidual< lowDimId >::ElementResidualVector evalCouplingResidual(LowDimIdType, const LowDimLocalAssembler &lowDimLocalAssembler, BulkIdType, GridIndexType< bulkId > dofIdxGlobalJ)
Evaluates the coupling element residual of a lower-dimensional domain element with respect to a dof i...
Definition: multidomain/facet/cellcentered/tpfa/couplingmanager.hh:329
void bindCouplingContext(BulkIdType, const Element< bulkId > &element, const Assembler &assembler)
For the assembly of the element residual of a bulk domain element we need to prepare all variables of...
Definition: multidomain/facet/cellcentered/tpfa/couplingmanager.hh:404
void updateCouplingContext(BulkIdType domainI, const BulkLocalAssembler &bulkLocalAssembler, BulkIdType domainJ, GridIndexType< bulkId > dofIdxGlobalJ, const PrimaryVariables< bulkId > &priVarsJ, unsigned int pvIdxJ)
Update the coupling context for a derivative bulk -> bulk. Here, we simply have to update the solutio...
Definition: multidomain/facet/cellcentered/tpfa/couplingmanager.hh:588
Implementation for the coupling manager between two domains of dimension d and (d-1) for models consi...
Definition: multidomain/facet/couplingmanager.hh:83
Defines all properties used in Dumux.
Element solution classes and factory functions.
typename NumEqVectorTraits< PrimaryVariables >::type NumEqVector
A vector with the same size as numbers of equations This is the default implementation and has to be ...
Definition: numeqvector.hh:34
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
void makeInterpolatedVolVars(VolumeVariables &volVars, const Problem &problem, const SolutionVector &sol, const FVGeometry &fvGeometry, const typename FVGeometry::GridGeometry::GridView::template Codim< 0 >::Entity &element, const typename FVGeometry::GridGeometry::GridView::template Codim< 0 >::Entity::Geometry &elemGeom, const typename FVGeometry::GridGeometry::GridView::template Codim< 0 >::Entity::Geometry::GlobalCoordinate &pos)
Free function that allows the creation of a volume variables object interpolated to a given position ...
Definition: multidomain/facet/couplingmanager.hh:40
static ctype distance(const Dune::FieldVector< ctype, dimWorld > &a, const Dune::FieldVector< ctype, dimWorld > &b)
Compute the shortest distance between two points.
Definition: distance.hh:282
typename GetProp< TypeTag, Property >::type GetPropType
get the type alias defined in the property
Definition: propertysystem.hh:296
Defines the index types used for grid and local indices.
The available discretization methods in Dumux.
The interface of the coupling manager for multi domain problems.
constexpr Box box
Definition: method.hh:147
Definition: adapt.hh:17
A helper to deduce a vector with the same size as numbers of equations.
Structure to define the index types used for grid and local indices.
Definition: indextraits.hh:26