3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
multidomain/facet/cellcentered/mpfa/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 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 2 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
24#ifndef DUMUX_CCMPFA_FACETCOUPLING_MANAGER_HH
25#define DUMUX_CCMPFA_FACETCOUPLING_MANAGER_HH
26
27#include <algorithm>
28#include <cassert>
29
34
39
40namespace Dumux {
41
53template<class MDTraits, class CouplingMapper, std::size_t bulkDomainId, std::size_t lowDimDomainId>
54class FacetCouplingManager<MDTraits, CouplingMapper, bulkDomainId, lowDimDomainId, DiscretizationMethod::ccmpfa>
55: public FacetCouplingManager<MDTraits, CouplingMapper, bulkDomainId, lowDimDomainId, DiscretizationMethod::cctpfa>
56{
58
59 // domain id instances
60 using BulkIdType = typename MDTraits::template SubDomain<bulkDomainId>::Index;
61 using LowDimIdType = typename MDTraits::template SubDomain<lowDimDomainId>::Index;
62 static constexpr auto bulkId = BulkIdType();
63 static constexpr auto lowDimId = LowDimIdType();
64
65 // the sub-domain type tags
66 template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
67
68 // further types specific to the sub-problems
69 template<std::size_t id> using Scalar = GetPropType<SubDomainTypeTag<id>, Properties::Scalar>;
70 template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>;
71 template<std::size_t id> using NumEqVector = GetPropType<SubDomainTypeTag<id>, Properties::NumEqVector>;
72 template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
73 template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView;
74 template<std::size_t id> using SubControlVolume = typename GridGeometry<id>::SubControlVolume;
75 template<std::size_t id> using SubControlVolumeFace = typename GridGeometry<id>::SubControlVolumeFace;
76 template<std::size_t id> using GridView = typename GridGeometry<id>::GridView;
77 template<std::size_t id> using GridIndexType = typename IndexTraits< GridView<id> >::GridIndex;
78 template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity;
79 template<std::size_t id> using LocalResidual = GetPropType<SubDomainTypeTag<id>, Properties::LocalResidual>;
80
81 template<std::size_t id> using GridVariables = GetPropType<SubDomainTypeTag<id>, Properties::GridVariables>;
82 template<std::size_t id> using GridVolumeVariables = typename GridVariables<id>::GridVolumeVariables;
83 template<std::size_t id> using ElementVolumeVariables = typename GridVolumeVariables<id>::LocalView;
84
85 // grid ids
86 static constexpr int bulkDim = GridView<bulkDomainId>::dimension;
87 static constexpr int lowDimDim = GridView<lowDimDomainId>::dimension;
88 static constexpr auto bulkGridId = CouplingMapper::template gridId<bulkDim>();
89 static constexpr auto lowDimGridId = CouplingMapper::template gridId<lowDimDim>();
90
91 static constexpr bool lowDimUsesBox = GridGeometry<lowDimId>::discMethod == DiscretizationMethod::box;
92
93public:
94
96 using SolutionVector = typename MDTraits::SolutionVector;
97
106 void init(std::shared_ptr< Problem<bulkId> > bulkProblem,
107 std::shared_ptr< Problem<lowDimId> > lowDimProblem,
108 std::shared_ptr< CouplingMapper > couplingMapper,
109 const SolutionVector& curSol)
110 {
111 // Initialize the parent class
112 ParentType::init(bulkProblem, lowDimProblem, couplingMapper, curSol);
113
114 // determine all bulk scvfs that coincide with low dim elements
115 bulkScvfIsOnFacetElement_.assign(bulkProblem->gridGeometry().numScvf(), false);
116 const auto& bulkMap = couplingMapper->couplingMap(bulkGridId, lowDimGridId);
117 for (const auto& entry : bulkMap)
118 for (const auto& couplingEntry : entry.second.elementToScvfMap)
119 for (const auto& scvfIdx : couplingEntry.second)
120 bulkScvfIsOnFacetElement_[scvfIdx] = true;
121
122 // store pointer to mapper
123 couplingMapperPtr_ = couplingMapper;
124 }
125
129 bool isOnInteriorBoundary(const Element<bulkId>& element,
130 const SubControlVolumeFace<bulkId>& scvf) const
131 { return bulkScvfIsOnFacetElement_[scvf.index()]; }
132
133 using ParentType::evalCouplingResidual;
140 template< class LowDimLocalAssembler >
141 typename LocalResidual<lowDimId>::ElementResidualVector
143 const LowDimLocalAssembler& lowDimLocalAssembler,
144 BulkIdType,
145 GridIndexType<bulkId> dofIdxGlobalJ)
146 { return evalCouplingResidual(lowDimId, lowDimLocalAssembler, bulkId); }
147
153 template< class LowDimLocalAssembler >
154 typename LocalResidual<lowDimId>::ElementResidualVector
155 evalCouplingResidual(LowDimIdType, const LowDimLocalAssembler& lowDimLocalAssembler, BulkIdType)
156 {
157 // make sure this is called for the element for which the context was set
158 assert(this->lowDimCouplingContext().isSet);
159 assert(this->problem(lowDimId).gridGeometry().elementMapper().index(lowDimLocalAssembler.element()) == this->lowDimCouplingContext().elementIdx);
160
161 // fill element residual vector with the sources
162 typename LowDimLocalAssembler::LocalResidual::ElementResidualVector res(lowDimLocalAssembler.fvGeometry().numScv());
163 res = 0.0;
164 for (const auto& scv : scvs(lowDimLocalAssembler.fvGeometry()))
165 res[scv.localDofIndex()] -= evalSourcesFromBulk(lowDimLocalAssembler.element(),
166 lowDimLocalAssembler.fvGeometry(),
167 lowDimLocalAssembler.curElemVolVars(),
168 scv);
169 return res;
170 }
171
175 NumEqVector<lowDimId> evalSourcesFromBulk(const Element<lowDimId>& element,
176 const FVElementGeometry<lowDimId>& fvGeometry,
177 const ElementVolumeVariables<lowDimId>& elemVolVars,
178 const SubControlVolume<lowDimId>& scv)
179 {
180 // make sure the this is called for the element of the context
181 assert(this->problem(lowDimId).gridGeometry().elementMapper().index(element) == this->lowDimCouplingContext().elementIdx);
182
183 NumEqVector<lowDimId> sources(0.0);
184
185 const auto& map = couplingMapperPtr_->couplingMap(lowDimGridId, bulkGridId);
186 auto it = map.find(this->lowDimCouplingContext().elementIdx);
187 if (it == map.end())
188 return sources;
189
190 assert(this->lowDimCouplingContext().isSet);
191 for (const auto& embedment : it->second.embedments)
192 {
193 // list of scvfs in the bulk domain whose fluxes enter this scv
194 // if low dim domain uses a cc scheme, this is all scvfs lying on this element
195 // if it uses box, it is the one scvf coinciding with the given scv
196 const auto& coincidingScvfs = embedment.second;
197 const auto& scvfList = lowDimUsesBox ? std::vector<GridIndexType<lowDimId>>{ coincidingScvfs[scv.localDofIndex()] }
198 : coincidingScvfs;
199
200 sources += this->evalBulkFluxes(this->problem(bulkId).gridGeometry().element(embedment.first),
201 *this->lowDimCouplingContext().bulkFvGeometry,
202 *this->lowDimCouplingContext().bulkElemVolVars,
203 *this->lowDimCouplingContext().bulkElemFluxVarsCache,
204 *this->lowDimCouplingContext().bulkLocalResidual,
205 scvfList);
206 }
207
208 return sources;
209 }
210
215 template<class JacobianPattern>
216 void extendJacobianPattern(LowDimIdType, JacobianPattern& pattern) const
217 {
218 const auto& lowDimFVGridGeometry = this->problem(lowDimId).gridGeometry();
219 for (const auto& element : elements(lowDimFVGridGeometry.gridView()))
220 {
221
222 const auto eIdx = lowDimFVGridGeometry.elementMapper().index(element);
223 const auto& map = couplingMapperPtr_->couplingMap(lowDimGridId, bulkGridId);
224 auto it = map.find(eIdx);
225
226 // if element is coupled, take one of the neighbors and add coupling stencil to pattern
227 if (it != map.end())
228 {
229 // coupling stencil of the first neighbor
230 const auto bulkElemIdx = it->second.embedments[0].first;
231 const auto& bulkMapEntry = couplingMapperPtr_->couplingMap(bulkGridId, lowDimGridId).at(bulkElemIdx);
232 const auto& couplingStencil = bulkMapEntry.couplingStencil;
233
234 for (auto globalJ : couplingStencil)
235 {
236 if (lowDimUsesBox)
237 {
238 for (int i = 0; i < element.subEntities(lowDimDim); ++i)
239 pattern.add(lowDimFVGridGeometry.vertexMapper().subIndex(element, i, lowDimDim), globalJ);
240 }
241 else
242 pattern.add(eIdx, globalJ);
243 }
244 }
245 }
246 }
247
249 template<class JacobianPattern>
250 void extendJacobianPattern(BulkIdType, JacobianPattern& pattern) const
251 {}
252
259 template<class LowDimLocalAssembler, class JacobianMatrixDiagBlock, class GridVariables>
261 const LowDimLocalAssembler& lowDimLocalAssembler,
262 const typename LowDimLocalAssembler::LocalResidual::ElementResidualVector&,
263 JacobianMatrixDiagBlock& A,
264 GridVariables& gridVariables)
265 {
266 // Since coupling only occurs via the fluxes, there are no
267 // additional derivatives for explicit time integration schemes
268 if (!LowDimLocalAssembler::isImplicit())
269 return;
270
271 // lambda to update the coupling context for a given lowDim element/dofIdx
272 auto updateContext = [&] (auto elemIdx, auto dofIdx, auto priVars, auto pvIdx)
273 {
274 // deflect the solution
275 auto& ldSol = this->curSol()[lowDimId];
276 ldSol[dofIdx][pvIdx] = priVars[pvIdx];
277
278 // update the corresponding vol vars in the bulk context
279 assert(this->bulkCouplingContext().isSet);
280 const auto& bulkMap = couplingMapperPtr_->couplingMap(bulkGridId, lowDimGridId);
281 const auto& couplingElementStencil = bulkMap.find(this->bulkCouplingContext().elementIdx)->second.couplingElementStencil;
282
283 auto it = std::find(couplingElementStencil.begin(), couplingElementStencil.end(), elemIdx);
284 assert(it != couplingElementStencil.end());
285 const auto idxInContext = std::distance(couplingElementStencil.begin(), it);
286
287 auto& volVars = this->bulkCouplingContext().lowDimVolVars[idxInContext];
288 const auto& fvGeom = this->bulkCouplingContext().lowDimFvGeometries[idxInContext];
289 const auto& element = this->problem(lowDimId).gridGeometry().element(elemIdx);
290
291 // if low dim domain uses the box scheme, we have to create interpolated vol vars
292 if (lowDimUsesBox)
293 {
294 const auto elemGeom = element.geometry();
295 FacetCoupling::makeInterpolatedVolVars(volVars, this->problem(lowDimId), ldSol, fvGeom, element, elemGeom, elemGeom.center());
296 }
297 // if low dim domain uses a cc scheme we can directly update the vol vars
298 else
299 volVars.update( elementSolution(element, ldSol, this->problem(lowDimId).gridGeometry()),
300 this->problem(lowDimId),
301 element,
302 fvGeom.scv(elemIdx) );
303
304 // update the element flux variables cache (tij depend on low dim values in context)
305 const auto contextElem = this->problem(bulkId).gridGeometry().element(this->bulkCouplingContext().elementIdx);
306 this->lowDimCouplingContext().bulkElemFluxVarsCache->update(contextElem,
307 *this->lowDimCouplingContext().bulkFvGeometry,
308 *this->lowDimCouplingContext().bulkElemVolVars);
309 };
310
311 const auto eIdx = this->problem(lowDimId).gridGeometry().elementMapper().index(lowDimLocalAssembler.element());
312
313 // bug tracking
314 assert(this->lowDimCouplingContext().isSet);
315 assert(this->lowDimCouplingContext().elementIdx == eIdx);
316
317 // if the element is coupled, evaluate additional source derivatives
318 const auto& map = couplingMapperPtr_->couplingMap(lowDimGridId, bulkGridId);
319 auto it = map.find(eIdx);
320 if (it != map.end())
321 evalLowDimSourceDerivatives_(updateContext, lowDimLocalAssembler, A);
322 }
323
325 template<class LocalAssemblerI, class JacobianMatrixDiagBlock, class GridVariables>
327 const LocalAssemblerI& localAssemblerI,
328 const typename LocalAssemblerI::LocalResidual::ElementResidualVector& origResiduals,
329 JacobianMatrixDiagBlock& A,
330 GridVariables& gridVariables)
331 {}
332
333private:
335 template<class UpdateContext, class LowDimLocalAssembler, class JacobianMatrixDiagBlock>
336 void evalLowDimSourceDerivatives_(const UpdateContext& updateContext,
337 const LowDimLocalAssembler& lowDimLocalAssembler,
338 JacobianMatrixDiagBlock& A)
339 {
340 const auto& lowDimFVGridGeometry = this->problem(lowDimId).gridGeometry();
341 const auto eIdx = lowDimFVGridGeometry.elementMapper().index(lowDimLocalAssembler.element());
342
343 // coupling stencil of the first neighbor
344 const auto bulkElemIdx = this->bulkCouplingContext().elementIdx;
345 const auto& bulkMapEntry = couplingMapperPtr_->couplingMap(bulkGridId, lowDimGridId).at(bulkElemIdx);
346 const auto& couplingStencil = bulkMapEntry.couplingStencil;
347 const auto& couplingElementStencil = bulkMapEntry.couplingElementStencil;
348
349 // compute the undeflected residual (reuse coupling residual function)
350 const auto origResidual = evalCouplingResidual(lowDimId, lowDimLocalAssembler, bulkId);
351
352 // container of dofs within this element
353 std::vector< std::decay_t<decltype(couplingStencil[0])> > elemDofs;
354 elemDofs.reserve(lowDimLocalAssembler.fvGeometry().numScv());
355 for (const auto& scv : scvs(lowDimLocalAssembler.fvGeometry()))
356 elemDofs.push_back(scv.dofIndex());
357
358 // compute derivate for all additional dofs in the stencil
359 for (const auto couplingElemIdx : couplingElementStencil)
360 {
361 // skip the same element
362 if (couplingElemIdx == eIdx)
363 continue;
364
365 // container of dofs within the other element
366 std::vector< std::decay_t<decltype(couplingStencil[0])> > elemDofsJ;
367 if (lowDimUsesBox)
368 {
369 const auto& elemJ = lowDimFVGridGeometry.element(couplingElemIdx);
370 for (int i = 0; i < elemJ.subEntities(lowDimDim); ++i)
371 elemDofsJ.push_back(lowDimFVGridGeometry.vertexMapper().subIndex(elemJ, i, lowDimDim));
372 }
373 else
374 elemDofsJ.push_back(couplingElemIdx);
375
376 for (auto dofIndex : elemDofsJ)
377 {
378 auto partialDerivs = origResidual;
379 const auto origPriVars = this->curSol()[lowDimId][dofIndex];
380
381 // calculate derivatives w.r.t to the privars at the dof at hand
382 static constexpr auto numEq = std::decay_t<decltype(origPriVars)>::dimension;
383 for (int pvIdx = 0; pvIdx < numEq; pvIdx++)
384 {
385 // reset partial derivatives
386 partialDerivs = 0.0;
387
388 auto evalResiduals = [&](Scalar<lowDimId> priVar)
389 {
390 auto priVars = origPriVars;
391 priVars[pvIdx] = priVar;
392
393 // Update context to deflected solution and reevaluate residual
394 updateContext(couplingElemIdx, dofIndex, priVars, pvIdx);
395 return this->evalCouplingResidual(lowDimId, lowDimLocalAssembler, bulkId);
396 };
397
398 static const int numDiffMethod = getParamFromGroup<int>(this->problem(lowDimId).paramGroup(), "Assembly.NumericDifferenceMethod");
399 static const NumericEpsilon< Scalar<lowDimId>, numEq > eps{this->problem(lowDimId).paramGroup()};
400 NumericDifferentiation::partialDerivative(evalResiduals, origPriVars[pvIdx], partialDerivs,
401 origResidual, eps(origPriVars[pvIdx], pvIdx), numDiffMethod);
402
403 // update the global stiffness matrix with the current partial derivatives
404 // A[i][col][eqIdx][pvIdx] is the rate of change of the residual of equation
405 // 'eqIdx' at dof 'i' depending on the primary variable 'pvIdx' at dof 'col'.
406 for (const auto& scv : scvs(lowDimLocalAssembler.fvGeometry()))
407 for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
408 A[scv.dofIndex()][dofIndex][eqIdx][pvIdx] += partialDerivs[scv.indexInElement()][eqIdx];
409
410 // restore the original coupling context
411 updateContext(couplingElemIdx, dofIndex, origPriVars, pvIdx);
412 }
413 }
414 }
415 }
416
418 std::vector<bool> bulkScvfIsOnFacetElement_;
419
421 std::shared_ptr< CouplingMapper > couplingMapperPtr_;
422};
423
424} // end namespace Dumux
425
426#endif
An adapter class for local assemblers using numeric differentiation.
A class for numeric differentiation.
Defines the index types used for grid and local indices.
Element solution classes and factory functions.
The available discretization methods in Dumux.
ctype distance(const Dune::FieldVector< ctype, dimWorld > &a, const Dune::FieldVector< ctype, dimWorld > &b)
Compute the shortest distance between two points.
Definition: geometry/distance.hh:138
DiscretizationMethod
The available discretization methods in Dumux.
Definition: method.hh:37
auto elementSolution(const Element &element, const SolutionVector &sol, const GridGeometry &gg) -> std::enable_if_t< GridGeometry::discMethod==DiscretizationMethod::box, BoxElementSolution< typename GridGeometry::LocalView, std::decay_t< decltype(std::declval< SolutionVector >()[0])> > >
Make an element solution for box schemes.
Definition: box/elementsolution.hh:115
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:52
Definition: adapt.hh:29
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:149
Struture to define the index types used for grid and local indices.
Definition: indextraits.hh:38
static void partialDerivative(const Function &function, Scalar x0, FunctionEvalType &derivative, const FunctionEvalType &fx0, const int numericDifferenceMethod=1)
Computes the derivative of a function with repect to a function parameter.
Definition: numericdifferentiation.hh:61
Property to specify the type of scalar values.
Definition: common/properties.hh:43
A vector of size number equations that can be used for Neumann fluxes, sources, residuals,...
Definition: common/properties.hh:51
Property to specify the type of a problem which has to be solved.
Definition: common/properties.hh:57
Definition: common/properties.hh:77
Definition: common/properties.hh:101
The grid variables object managing variable data on the grid (volvars/fluxvars cache)
Definition: common/properties.hh:122
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/mpfa/couplingmanager.hh:106
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/mpfa/couplingmanager.hh:129
void extendJacobianPattern(LowDimIdType, JacobianPattern &pattern) const
Extend the jacobian pattern of the diagonal block of the lowdim domain by the elements that are in th...
Definition: multidomain/facet/cellcentered/mpfa/couplingmanager.hh:216
typename MDTraits::SolutionVector SolutionVector
the type of the solution vector
Definition: multidomain/facet/cellcentered/mpfa/couplingmanager.hh:96
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/mpfa/couplingmanager.hh:155
void evalAdditionalDomainDerivatives(LowDimIdType, const LowDimLocalAssembler &lowDimLocalAssembler, const typename LowDimLocalAssembler::LocalResidual::ElementResidualVector &, JacobianMatrixDiagBlock &A, GridVariables &gridVariables)
evaluate additional derivatives of the element residual of the low-dim domain with respect to dofs in...
Definition: multidomain/facet/cellcentered/mpfa/couplingmanager.hh:260
void evalAdditionalDomainDerivatives(BulkIdType, const LocalAssemblerI &localAssemblerI, const typename LocalAssemblerI::LocalResidual::ElementResidualVector &origResiduals, JacobianMatrixDiagBlock &A, GridVariables &gridVariables)
The bulk domain has no additional derivatives.
Definition: multidomain/facet/cellcentered/mpfa/couplingmanager.hh:326
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/mpfa/couplingmanager.hh:142
void extendJacobianPattern(BulkIdType, JacobianPattern &pattern) const
The bulk domain has no extended jacobian pattern.
Definition: multidomain/facet/cellcentered/mpfa/couplingmanager.hh:250
NumEqVector< lowDimId > evalSourcesFromBulk(const Element< lowDimId > &element, const FVElementGeometry< lowDimId > &fvGeometry, const ElementVolumeVariables< lowDimId > &elemVolVars, const SubControlVolume< lowDimId > &scv)
Computes the sources in a lower-dimensional sub-control volume stemming from the bulk domain.
Definition: multidomain/facet/cellcentered/mpfa/couplingmanager.hh:175
Manages the coupling between bulk elements and lower dimensional elements where the coupling occurs a...
Definition: multidomain/facet/cellcentered/tpfa/couplingmanager.hh:53
Implementation for the coupling manager between two domains of dimension d and (d-1) for models consi...
Definition: multidomain/facet/couplingmanager.hh:95
Declares all properties used in Dumux.
The interface of the coupling manager for multi domain problems.