| Line | Branch | Exec | Source |
|---|---|---|---|
| 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-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder | ||
| 5 | // SPDX-License-Identifier: GPL-3.0-or-later | ||
| 6 | // | ||
| 7 | /*! | ||
| 8 | * \file | ||
| 9 | * \ingroup FacetCoupling | ||
| 10 | * \copydoc Dumux::FacetCouplingManager | ||
| 11 | */ | ||
| 12 | #ifndef DUMUX_CCMPFA_FACETCOUPLING_MANAGER_HH | ||
| 13 | #define DUMUX_CCMPFA_FACETCOUPLING_MANAGER_HH | ||
| 14 | |||
| 15 | #include <algorithm> | ||
| 16 | #include <cassert> | ||
| 17 | |||
| 18 | #include <dumux/common/properties.hh> | ||
| 19 | #include <dumux/common/indextraits.hh> | ||
| 20 | #include <dumux/common/numericdifferentiation.hh> | ||
| 21 | #include <dumux/common/numeqvector.hh> | ||
| 22 | #include <dumux/assembly/numericepsilon.hh> | ||
| 23 | |||
| 24 | #include <dumux/discretization/method.hh> | ||
| 25 | #include <dumux/discretization/elementsolution.hh> | ||
| 26 | #include <dumux/multidomain/couplingmanager.hh> | ||
| 27 | #include <dumux/multidomain/facet/cellcentered/tpfa/couplingmanager.hh> | ||
| 28 | |||
| 29 | namespace Dumux { | ||
| 30 | |||
| 31 | /*! | ||
| 32 | * \ingroup FacetCoupling | ||
| 33 | * \brief Manages the coupling between bulk elements and lower dimensional elements | ||
| 34 | * where the coupling occurs across the facets of the bulk grid. This implementation | ||
| 35 | * is to be used in conjunction with models using the cell-centered mpfa scheme. | ||
| 36 | * | ||
| 37 | * \tparam MDTraits The multidomain traits containing the types on all sub-domains | ||
| 38 | * \tparam CouplingMapper Class containing maps on the coupling between dofs of different grids | ||
| 39 | * \tparam bulkDomainId The domain id of the bulk problem | ||
| 40 | * \tparam lowDimDomainId The domain id of the lower-dimensional problem | ||
| 41 | */ | ||
| 42 | template<class MDTraits, class CouplingMapper, std::size_t bulkDomainId, std::size_t lowDimDomainId> | ||
| 43 | class FacetCouplingManager<MDTraits, CouplingMapper, bulkDomainId, lowDimDomainId, DiscretizationMethods::CCMpfa> | ||
| 44 | : public FacetCouplingManager<MDTraits, CouplingMapper, bulkDomainId, lowDimDomainId, DiscretizationMethods::CCTpfa> | ||
| 45 | { | ||
| 46 | using ParentType = FacetCouplingManager<MDTraits, CouplingMapper, bulkDomainId, lowDimDomainId, DiscretizationMethods::CCTpfa>; | ||
| 47 | |||
| 48 | // domain id instances | ||
| 49 | using BulkIdType = typename MDTraits::template SubDomain<bulkDomainId>::Index; | ||
| 50 | using LowDimIdType = typename MDTraits::template SubDomain<lowDimDomainId>::Index; | ||
| 51 | static constexpr auto bulkId = BulkIdType(); | ||
| 52 | static constexpr auto lowDimId = LowDimIdType(); | ||
| 53 | |||
| 54 | // the sub-domain type tags | ||
| 55 | template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag; | ||
| 56 | |||
| 57 | // further types specific to the sub-problems | ||
| 58 | template<std::size_t id> using Scalar = GetPropType<SubDomainTypeTag<id>, Properties::Scalar>; | ||
| 59 | template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>; | ||
| 60 | template<std::size_t id> using NumEqVector = Dumux::NumEqVector<GetPropType<SubDomainTypeTag<id>, Properties::PrimaryVariables>>; | ||
| 61 | template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>; | ||
| 62 | template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView; | ||
| 63 | template<std::size_t id> using SubControlVolume = typename GridGeometry<id>::SubControlVolume; | ||
| 64 | template<std::size_t id> using SubControlVolumeFace = typename GridGeometry<id>::SubControlVolumeFace; | ||
| 65 | template<std::size_t id> using GridView = typename GridGeometry<id>::GridView; | ||
| 66 | template<std::size_t id> using GridIndexType = typename IndexTraits< GridView<id> >::GridIndex; | ||
| 67 | template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity; | ||
| 68 | template<std::size_t id> using LocalResidual = typename MDTraits::template SubDomain<id>::LocalResidual; | ||
| 69 | |||
| 70 | template<std::size_t id> using GridVariables = GetPropType<SubDomainTypeTag<id>, Properties::GridVariables>; | ||
| 71 | template<std::size_t id> using GridVolumeVariables = typename GridVariables<id>::GridVolumeVariables; | ||
| 72 | template<std::size_t id> using ElementVolumeVariables = typename GridVolumeVariables<id>::LocalView; | ||
| 73 | |||
| 74 | // grid ids | ||
| 75 | static constexpr int bulkDim = GridView<bulkDomainId>::dimension; | ||
| 76 | static constexpr int lowDimDim = GridView<lowDimDomainId>::dimension; | ||
| 77 | static constexpr auto bulkGridId = CouplingMapper::template gridId<bulkDim>(); | ||
| 78 | static constexpr auto lowDimGridId = CouplingMapper::template gridId<lowDimDim>(); | ||
| 79 | |||
| 80 | static constexpr bool lowDimUsesBox = GridGeometry<lowDimId>::discMethod == DiscretizationMethods::box; | ||
| 81 | |||
| 82 | public: | ||
| 83 | |||
| 84 | //! the type of the solution vector | ||
| 85 | using SolutionVector = typename MDTraits::SolutionVector; | ||
| 86 | |||
| 87 | /*! | ||
| 88 | * \brief Initialize the coupling manager. | ||
| 89 | * | ||
| 90 | * \param bulkProblem The problem to be solved on the bulk domain | ||
| 91 | * \param lowDimProblem The problem to be solved on the lower-dimensional domain | ||
| 92 | * \param couplingMapper The mapper object containing the connectivity between the domains | ||
| 93 | * \param curSol The current solution | ||
| 94 | */ | ||
| 95 | 15 | void init(std::shared_ptr< Problem<bulkId> > bulkProblem, | |
| 96 | std::shared_ptr< Problem<lowDimId> > lowDimProblem, | ||
| 97 | std::shared_ptr< CouplingMapper > couplingMapper, | ||
| 98 | const SolutionVector& curSol) | ||
| 99 | { | ||
| 100 | // Initialize the parent class | ||
| 101 |
3/10✓ Branch 4 taken 13 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 13 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 13 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
|
45 | ParentType::init(bulkProblem, lowDimProblem, couplingMapper, curSol); |
| 102 | |||
| 103 | // determine all bulk scvfs that coincide with low dim elements | ||
| 104 | 15 | bulkScvfIsOnFacetElement_.assign(bulkProblem->gridGeometry().numScvf(), false); | |
| 105 | 15 | const auto& bulkMap = couplingMapper->couplingMap(bulkGridId, lowDimGridId); | |
| 106 |
2/2✓ Branch 0 taken 1405 times.
✓ Branch 1 taken 13 times.
|
2248 | for (const auto& entry : bulkMap) |
| 107 |
2/2✓ Branch 0 taken 940 times.
✓ Branch 1 taken 1405 times.
|
3573 | for (const auto& couplingEntry : entry.second.elementToScvfMap) |
| 108 |
2/2✓ Branch 0 taken 1880 times.
✓ Branch 1 taken 940 times.
|
4020 | for (const auto& scvfIdx : couplingEntry.second) |
| 109 | 2680 | bulkScvfIsOnFacetElement_[scvfIdx] = true; | |
| 110 | |||
| 111 | // store pointer to mapper | ||
| 112 | 15 | couplingMapperPtr_ = couplingMapper; | |
| 113 | 15 | } | |
| 114 | |||
| 115 | /*! | ||
| 116 | * \brief returns true if a bulk scvf coincides with a facet element. | ||
| 117 | */ | ||
| 118 | 119663643 | bool isOnInteriorBoundary(const Element<bulkId>& element, | |
| 119 | const SubControlVolumeFace<bulkId>& scvf) const | ||
| 120 |
14/14✓ Branch 0 taken 1373368 times.
✓ Branch 1 taken 33979340 times.
✓ Branch 2 taken 1160280 times.
✓ Branch 3 taken 31783905 times.
✓ Branch 4 taken 673302 times.
✓ Branch 5 taken 16278634 times.
✓ Branch 6 taken 595944 times.
✓ Branch 7 taken 15327744 times.
✓ Branch 8 taken 829000 times.
✓ Branch 9 taken 17504000 times.
✓ Branch 10 taken 6632 times.
✓ Branch 11 taken 140032 times.
✓ Branch 12 taken 400 times.
✓ Branch 13 taken 11062 times.
|
119663643 | { return bulkScvfIsOnFacetElement_[scvf.index()]; } |
| 121 | |||
| 122 | using ParentType::evalCouplingResidual; | ||
| 123 | /*! | ||
| 124 | * \brief Evaluates the coupling element residual of a lower-dimensional domain element | ||
| 125 | * with respect to a dof in the bulk domain (dofIdxGlobalJ). This is essentially | ||
| 126 | * the fluxes across the facets of the neighboring bulk elements. | ||
| 127 | * \note The coupling residual is independent of w.r.t. which bulk dof it is computed | ||
| 128 | */ | ||
| 129 | template< class LowDimLocalAssembler > | ||
| 130 | typename LocalResidual<lowDimId>::ElementResidualVector | ||
| 131 | 87984 | evalCouplingResidual(LowDimIdType, | |
| 132 | const LowDimLocalAssembler& lowDimLocalAssembler, | ||
| 133 | BulkIdType, | ||
| 134 | GridIndexType<bulkId> dofIdxGlobalJ) | ||
| 135 |
1/2✓ Branch 1 taken 26814 times.
✗ Branch 2 not taken.
|
87984 | { return evalCouplingResidual(lowDimId, lowDimLocalAssembler, bulkId); } |
| 136 | |||
| 137 | /*! | ||
| 138 | * \brief Evaluates the coupling element residual of a lower-dimensional domain element | ||
| 139 | * with respect to a dof in the bulk domain (dofIdxGlobalJ). This is essentially | ||
| 140 | * the fluxes across the facets of the neighboring bulk elements. | ||
| 141 | */ | ||
| 142 | template< class LowDimLocalAssembler > | ||
| 143 | typename LocalResidual<lowDimId>::ElementResidualVector | ||
| 144 | 110977 | evalCouplingResidual(LowDimIdType, const LowDimLocalAssembler& lowDimLocalAssembler, BulkIdType) | |
| 145 | { | ||
| 146 | // make sure this is called for the element for which the context was set | ||
| 147 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 110977 times.
|
110977 | assert(this->lowDimCouplingContext().isSet); |
| 148 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 110977 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 110977 times.
|
110977 | assert(this->problem(lowDimId).gridGeometry().elementMapper().index(lowDimLocalAssembler.element()) == this->lowDimCouplingContext().elementIdx); |
| 149 | |||
| 150 | // fill element residual vector with the sources | ||
| 151 | 110977 | typename LowDimLocalAssembler::LocalResidual::ElementResidualVector res(lowDimLocalAssembler.fvGeometry().numScv()); | |
| 152 | 110977 | res = 0.0; | |
| 153 |
2/2✓ Branch 0 taken 110977 times.
✓ Branch 1 taken 110977 times.
|
221954 | for (const auto& scv : scvs(lowDimLocalAssembler.fvGeometry())) |
| 154 | 206303 | res[scv.localDofIndex()] -= evalSourcesFromBulk(lowDimLocalAssembler.element(), | |
| 155 | 110977 | lowDimLocalAssembler.fvGeometry(), | |
| 156 | 110977 | lowDimLocalAssembler.curElemVolVars(), | |
| 157 | scv); | ||
| 158 | 110977 | return res; | |
| 159 | } | ||
| 160 | |||
| 161 | /*! | ||
| 162 | * \brief Computes the sources in a lower-dimensional sub-control volume stemming from the bulk domain. | ||
| 163 | */ | ||
| 164 | 233013 | NumEqVector<lowDimId> evalSourcesFromBulk(const Element<lowDimId>& element, | |
| 165 | const FVElementGeometry<lowDimId>& fvGeometry, | ||
| 166 | const ElementVolumeVariables<lowDimId>& elemVolVars, | ||
| 167 | const SubControlVolume<lowDimId>& scv) const | ||
| 168 | { | ||
| 169 | // make sure the this is called for the element of the context | ||
| 170 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 175249 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 175249 times.
|
233013 | assert(this->problem(lowDimId).gridGeometry().elementMapper().index(element) == this->lowDimCouplingContext().elementIdx); |
| 171 | |||
| 172 | 233013 | NumEqVector<lowDimId> sources(0.0); | |
| 173 | |||
| 174 | 233013 | const auto& map = couplingMapperPtr_->couplingMap(lowDimGridId, bulkGridId); | |
| 175 |
2/2✓ Branch 1 taken 80828 times.
✓ Branch 2 taken 94421 times.
|
233013 | auto it = map.find(this->lowDimCouplingContext().elementIdx); |
| 176 |
2/2✓ Branch 0 taken 80828 times.
✓ Branch 1 taken 94421 times.
|
233013 | if (it == map.end()) |
| 177 | ✗ | return sources; | |
| 178 | |||
| 179 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 175249 times.
|
233013 | assert(this->lowDimCouplingContext().isSet); |
| 180 |
2/2✓ Branch 0 taken 349616 times.
✓ Branch 1 taken 175249 times.
|
698157 | for (const auto& embedment : it->second.embedments) |
| 181 | { | ||
| 182 | // list of scvfs in the bulk domain whose fluxes enter this scv | ||
| 183 | // if low dim domain uses a cc scheme, this is all scvfs lying on this element | ||
| 184 | // if it uses box, it is the one scvf coinciding with the given scv | ||
| 185 | 465144 | const auto& coincidingScvfs = embedment.second; | |
| 186 |
1/2✓ Branch 1 taken 349616 times.
✗ Branch 2 not taken.
|
465144 | const auto& scvfList = lowDimUsesBox ? std::vector<GridIndexType<lowDimId>>{ coincidingScvfs[scv.localDofIndex()] } |
| 187 | : coincidingScvfs; | ||
| 188 | |||
| 189 |
2/4✓ Branch 1 taken 349616 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 349616 times.
✗ Branch 5 not taken.
|
1145944 | sources += this->evalBulkFluxes(this->problem(bulkId).gridGeometry().element(embedment.first), |
| 190 | 465144 | *this->lowDimCouplingContext().bulkFvGeometry, | |
| 191 | 465144 | *this->lowDimCouplingContext().bulkElemVolVars, | |
| 192 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 349616 times.
|
465144 | *this->lowDimCouplingContext().bulkElemFluxVarsCache, |
| 193 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 349616 times.
|
465144 | *this->lowDimCouplingContext().bulkLocalResidual, |
| 194 | scvfList); | ||
| 195 | } | ||
| 196 | |||
| 197 | 152185 | return sources; | |
| 198 | } | ||
| 199 | |||
| 200 | /*! | ||
| 201 | * \brief Extend the jacobian pattern of the diagonal block of the lowdim domain | ||
| 202 | * by the elements that are in the coupling stencil of the neighboring bulk elements | ||
| 203 | */ | ||
| 204 | template<class JacobianPattern> | ||
| 205 | 15 | void extendJacobianPattern(LowDimIdType, JacobianPattern& pattern) const | |
| 206 | { | ||
| 207 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 13 times.
|
15 | const auto& lowDimFVGridGeometry = this->problem(lowDimId).gridGeometry(); |
| 208 |
2/2✓ Branch 2 taken 480 times.
✓ Branch 3 taken 13 times.
|
1375 | for (const auto& element : elements(lowDimFVGridGeometry.gridView())) |
| 209 | { | ||
| 210 | |||
| 211 | 680 | const auto eIdx = lowDimFVGridGeometry.elementMapper().index(element); | |
| 212 | 680 | const auto& map = couplingMapperPtr_->couplingMap(lowDimGridId, bulkGridId); | |
| 213 |
1/2✓ Branch 0 taken 480 times.
✗ Branch 1 not taken.
|
680 | auto it = map.find(eIdx); |
| 214 | |||
| 215 | // if element is coupled, take one of the neighbors and add coupling stencil to pattern | ||
| 216 |
1/2✓ Branch 0 taken 480 times.
✗ Branch 1 not taken.
|
680 | if (it != map.end()) |
| 217 | { | ||
| 218 | // coupling stencil of the first neighbor | ||
| 219 |
1/2✓ Branch 1 taken 480 times.
✗ Branch 2 not taken.
|
680 | const auto bulkElemIdx = it->second.embedments[0].first; |
| 220 |
1/2✓ Branch 1 taken 480 times.
✗ Branch 2 not taken.
|
680 | const auto& bulkMapEntry = couplingMapperPtr_->couplingMap(bulkGridId, lowDimGridId).at(bulkElemIdx); |
| 221 | 680 | const auto& couplingStencil = bulkMapEntry.couplingStencil; | |
| 222 | |||
| 223 |
2/2✓ Branch 0 taken 1458 times.
✓ Branch 1 taken 480 times.
|
2778 | for (auto globalJ : couplingStencil) |
| 224 | { | ||
| 225 | if (lowDimUsesBox) | ||
| 226 | { | ||
| 227 | for (int i = 0; i < element.subEntities(lowDimDim); ++i) | ||
| 228 | pattern.add(lowDimFVGridGeometry.vertexMapper().subIndex(element, i, lowDimDim), globalJ); | ||
| 229 | } | ||
| 230 | else | ||
| 231 |
1/2✓ Branch 1 taken 1458 times.
✗ Branch 2 not taken.
|
2098 | pattern.add(eIdx, globalJ); |
| 232 | } | ||
| 233 | } | ||
| 234 | } | ||
| 235 | 15 | } | |
| 236 | |||
| 237 | //! The bulk domain has no extended jacobian pattern | ||
| 238 | template<class JacobianPattern> | ||
| 239 | void extendJacobianPattern(BulkIdType, JacobianPattern& pattern) const | ||
| 240 | {} | ||
| 241 | |||
| 242 | /*! | ||
| 243 | * \brief evaluate additional derivatives of the element residual of the low-dim domain with respect | ||
| 244 | * to dofs in the same domain that are not in the regular stencil (see extendJacobianPattern) | ||
| 245 | * \note Here, this is the change of the source term with respect to changes in the variables of the | ||
| 246 | * other elements in the coupling stencil of the neighboring bulk elements. | ||
| 247 | */ | ||
| 248 | template<class LowDimLocalAssembler, class JacobianMatrixDiagBlock, class GridVariables> | ||
| 249 | 4273 | void evalAdditionalDomainDerivatives(LowDimIdType, | |
| 250 | const LowDimLocalAssembler& lowDimLocalAssembler, | ||
| 251 | const typename LowDimLocalAssembler::LocalResidual::ElementResidualVector&, | ||
| 252 | JacobianMatrixDiagBlock& A, | ||
| 253 | GridVariables& gridVariables) | ||
| 254 | { | ||
| 255 | // Since coupling only occurs via the fluxes, there are no | ||
| 256 | // additional derivatives for explicit time integration schemes | ||
| 257 | if (!LowDimLocalAssembler::isImplicit()) | ||
| 258 | return; | ||
| 259 | |||
| 260 | // lambda to update the coupling context for a given lowDim element/dofIdx | ||
| 261 | 79153 | auto updateContext = [&] (auto elemIdx, auto dofIdx, auto priVars, auto pvIdx) | |
| 262 | { | ||
| 263 | // deflect the solution | ||
| 264 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 37440 times.
|
37440 | auto& ldSol = this->curSol(lowDimId); |
| 265 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 37440 times.
|
37440 | ldSol[dofIdx][pvIdx] = priVars[pvIdx]; |
| 266 | |||
| 267 | // update the corresponding vol vars in the bulk context | ||
| 268 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 37440 times.
|
37440 | assert(this->bulkCouplingContext().isSet); |
| 269 | 37440 | const auto& bulkMap = couplingMapperPtr_->couplingMap(bulkGridId, lowDimGridId); | |
| 270 |
2/3✗ Branch 1 not taken.
✓ Branch 2 taken 36120 times.
✓ Branch 3 taken 1320 times.
|
37440 | const auto& couplingElementStencil = bulkMap.find(this->bulkCouplingContext().elementIdx)->second.couplingElementStencil; |
| 271 | |||
| 272 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 37440 times.
|
37440 | auto it = std::find(couplingElementStencil.begin(), couplingElementStencil.end(), elemIdx); |
| 273 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 37440 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 37440 times.
|
37440 | assert(it != couplingElementStencil.end()); |
| 274 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 37440 times.
|
37440 | const auto idxInContext = std::distance(couplingElementStencil.begin(), it); |
| 275 | |||
| 276 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 37440 times.
|
37440 | auto& volVars = this->bulkCouplingContext().lowDimVolVars[idxInContext]; |
| 277 | 37440 | const auto& fvGeom = this->bulkCouplingContext().lowDimFvGeometries[idxInContext]; | |
| 278 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 37440 times.
✓ Branch 2 taken 37440 times.
✗ Branch 3 not taken.
|
37440 | const auto& element = this->problem(lowDimId).gridGeometry().element(elemIdx); |
| 279 | |||
| 280 | // if low dim domain uses the box scheme, we have to create interpolated vol vars | ||
| 281 | if (lowDimUsesBox) | ||
| 282 | { | ||
| 283 | const auto elemGeom = element.geometry(); | ||
| 284 | FacetCoupling::makeInterpolatedVolVars(volVars, this->problem(lowDimId), ldSol, fvGeom, element, elemGeom, elemGeom.center()); | ||
| 285 | } | ||
| 286 | // if low dim domain uses a cc scheme we can directly update the vol vars | ||
| 287 | else | ||
| 288 |
1/2✓ Branch 0 taken 37440 times.
✗ Branch 1 not taken.
|
74880 | volVars.update( elementSolution(element, ldSol, this->problem(lowDimId).gridGeometry()), |
| 289 | 37440 | this->problem(lowDimId), | |
| 290 | element, | ||
| 291 | fvGeom.scv(elemIdx) ); | ||
| 292 | |||
| 293 | // update the element flux variables cache (tij depend on low dim values in context) | ||
| 294 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 37440 times.
|
37440 | const auto contextElem = this->problem(bulkId).gridGeometry().element(this->bulkCouplingContext().elementIdx); |
| 295 |
1/2✓ Branch 1 taken 37440 times.
✗ Branch 2 not taken.
|
37440 | this->lowDimCouplingContext().bulkElemFluxVarsCache->update(contextElem, |
| 296 |
1/2✓ Branch 1 taken 37440 times.
✗ Branch 2 not taken.
|
37440 | *this->lowDimCouplingContext().bulkFvGeometry, |
| 297 |
1/2✓ Branch 1 taken 37440 times.
✗ Branch 2 not taken.
|
37440 | *this->lowDimCouplingContext().bulkElemVolVars); |
| 298 | 37440 | }; | |
| 299 | |||
| 300 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4273 times.
|
4273 | const auto eIdx = this->problem(lowDimId).gridGeometry().elementMapper().index(lowDimLocalAssembler.element()); |
| 301 | |||
| 302 | // bug tracking | ||
| 303 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4273 times.
|
4273 | assert(this->lowDimCouplingContext().isSet); |
| 304 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4273 times.
|
4273 | assert(this->lowDimCouplingContext().elementIdx == eIdx); |
| 305 | |||
| 306 | // if the element is coupled, evaluate additional source derivatives | ||
| 307 | 4273 | const auto& map = couplingMapperPtr_->couplingMap(lowDimGridId, bulkGridId); | |
| 308 |
1/2✓ Branch 0 taken 4273 times.
✗ Branch 1 not taken.
|
4273 | auto it = map.find(eIdx); |
| 309 |
1/2✓ Branch 0 taken 4273 times.
✗ Branch 1 not taken.
|
4273 | if (it != map.end()) |
| 310 | 4273 | evalLowDimSourceDerivatives_(updateContext, lowDimLocalAssembler, A); | |
| 311 | } | ||
| 312 | |||
| 313 | //! The bulk domain has no additional derivatives | ||
| 314 | template<class LocalAssemblerI, class JacobianMatrixDiagBlock, class GridVariables> | ||
| 315 | void evalAdditionalDomainDerivatives(BulkIdType, | ||
| 316 | const LocalAssemblerI& localAssemblerI, | ||
| 317 | const typename LocalAssemblerI::LocalResidual::ElementResidualVector& origResiduals, | ||
| 318 | JacobianMatrixDiagBlock& A, | ||
| 319 | GridVariables& gridVariables) | ||
| 320 | {} | ||
| 321 | |||
| 322 | private: | ||
| 323 | //! evaluates the additional source derivatives w.r.t. to neighboring elements | ||
| 324 | template<class UpdateContext, class LowDimLocalAssembler, class JacobianMatrixDiagBlock> | ||
| 325 | 4273 | void evalLowDimSourceDerivatives_(const UpdateContext& updateContext, | |
| 326 | const LowDimLocalAssembler& lowDimLocalAssembler, | ||
| 327 | JacobianMatrixDiagBlock& A) | ||
| 328 | { | ||
| 329 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4273 times.
|
4273 | const auto& lowDimFVGridGeometry = this->problem(lowDimId).gridGeometry(); |
| 330 | 4273 | const auto eIdx = lowDimFVGridGeometry.elementMapper().index(lowDimLocalAssembler.element()); | |
| 331 | |||
| 332 | // coupling stencil of the first neighbor | ||
| 333 | 4273 | const auto bulkElemIdx = this->bulkCouplingContext().elementIdx; | |
| 334 | 4273 | const auto& bulkMapEntry = couplingMapperPtr_->couplingMap(bulkGridId, lowDimGridId).at(bulkElemIdx); | |
| 335 | 4273 | const auto& couplingStencil = bulkMapEntry.couplingStencil; | |
| 336 | 4273 | const auto& couplingElementStencil = bulkMapEntry.couplingElementStencil; | |
| 337 | |||
| 338 | // compute the undeflected residual (reuse coupling residual function) | ||
| 339 | 4273 | const auto origResidual = evalCouplingResidual(lowDimId, lowDimLocalAssembler, bulkId); | |
| 340 | |||
| 341 | // container of dofs within this element | ||
| 342 | 4273 | std::vector< std::decay_t<decltype(couplingStencil[0])> > elemDofs; | |
| 343 |
1/2✓ Branch 1 taken 4273 times.
✗ Branch 2 not taken.
|
4273 | elemDofs.reserve(lowDimLocalAssembler.fvGeometry().numScv()); |
| 344 |
3/4✓ Branch 1 taken 4273 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4273 times.
✓ Branch 4 taken 4273 times.
|
8546 | for (const auto& scv : scvs(lowDimLocalAssembler.fvGeometry())) |
| 345 |
1/2✓ Branch 1 taken 4273 times.
✗ Branch 2 not taken.
|
4273 | elemDofs.push_back(scv.dofIndex()); |
| 346 | |||
| 347 | // compute derivate for all additional dofs in the stencil | ||
| 348 |
2/2✓ Branch 1 taken 12359 times.
✓ Branch 2 taken 4273 times.
|
24718 | for (const auto couplingElemIdx : couplingElementStencil) |
| 349 | { | ||
| 350 | // skip the same element | ||
| 351 |
2/2✓ Branch 0 taken 4273 times.
✓ Branch 1 taken 8086 times.
|
12359 | if (couplingElemIdx == eIdx) |
| 352 | 4273 | continue; | |
| 353 | |||
| 354 | // container of dofs within the other element | ||
| 355 | 8086 | std::vector< std::decay_t<decltype(couplingStencil[0])> > elemDofsJ; | |
| 356 | if (lowDimUsesBox) | ||
| 357 | { | ||
| 358 | const auto& elemJ = lowDimFVGridGeometry.element(couplingElemIdx); | ||
| 359 | for (int i = 0; i < elemJ.subEntities(lowDimDim); ++i) | ||
| 360 | elemDofsJ.push_back(lowDimFVGridGeometry.vertexMapper().subIndex(elemJ, i, lowDimDim)); | ||
| 361 | } | ||
| 362 | else | ||
| 363 |
1/2✓ Branch 1 taken 8086 times.
✗ Branch 2 not taken.
|
8086 | elemDofsJ.push_back(couplingElemIdx); |
| 364 | |||
| 365 |
2/2✓ Branch 0 taken 8086 times.
✓ Branch 1 taken 8086 times.
|
16172 | for (auto dofIndex : elemDofsJ) |
| 366 | { | ||
| 367 | 8086 | auto partialDerivs = origResidual; | |
| 368 | 8086 | const auto origPriVars = this->curSol(lowDimId)[dofIndex]; | |
| 369 | |||
| 370 | // calculate derivatives w.r.t to the privars at the dof at hand | ||
| 371 | static constexpr auto numEq = std::decay_t<decltype(origPriVars)>::dimension; | ||
| 372 |
2/2✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 8086 times.
|
26806 | for (int pvIdx = 0; pvIdx < numEq; pvIdx++) |
| 373 | { | ||
| 374 | // reset partial derivatives | ||
| 375 | 18720 | partialDerivs = 0.0; | |
| 376 | |||
| 377 | 37440 | auto evalResiduals = [&](Scalar<lowDimId> priVar) | |
| 378 | { | ||
| 379 | 16926 | auto priVars = origPriVars; | |
| 380 | 18720 | priVars[pvIdx] = priVar; | |
| 381 | |||
| 382 | // Update context to deflected solution and reevaluate residual | ||
| 383 | 18720 | updateContext(couplingElemIdx, dofIndex, priVars, pvIdx); | |
| 384 | 18720 | return this->evalCouplingResidual(lowDimId, lowDimLocalAssembler, bulkId); | |
| 385 | }; | ||
| 386 | |||
| 387 |
5/8✓ Branch 0 taken 12 times.
✓ Branch 1 taken 18708 times.
✓ Branch 3 taken 12 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 12 times.
✓ Branch 8 taken 12 times.
✗ Branch 9 not taken.
|
18720 | static const int numDiffMethod = getParamFromGroup<int>(this->problem(lowDimId).paramGroup(), "Assembly.NumericDifferenceMethod"); |
| 388 |
5/8✓ Branch 0 taken 12 times.
✓ Branch 1 taken 18708 times.
✓ Branch 3 taken 12 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 12 times.
✓ Branch 8 taken 12 times.
✗ Branch 9 not taken.
|
18720 | static const NumericEpsilon< Scalar<lowDimId>, numEq > eps{this->problem(lowDimId).paramGroup()}; |
| 389 |
1/2✓ Branch 1 taken 18720 times.
✗ Branch 2 not taken.
|
18720 | NumericDifferentiation::partialDerivative(evalResiduals, origPriVars[pvIdx], partialDerivs, |
| 390 | 18720 | origResidual, eps(origPriVars[pvIdx], pvIdx), numDiffMethod); | |
| 391 | |||
| 392 | // update the global stiffness matrix with the current partial derivatives | ||
| 393 | // A[i][col][eqIdx][pvIdx] is the rate of change of the residual of equation | ||
| 394 | // 'eqIdx' at dof 'i' depending on the primary variable 'pvIdx' at dof 'col'. | ||
| 395 |
2/2✓ Branch 0 taken 18720 times.
✓ Branch 1 taken 18720 times.
|
37440 | for (const auto& scv : scvs(lowDimLocalAssembler.fvGeometry())) |
| 396 |
2/2✓ Branch 0 taken 48672 times.
✓ Branch 1 taken 18720 times.
|
67392 | for (int eqIdx = 0; eqIdx < numEq; eqIdx++) |
| 397 |
1/2✓ Branch 1 taken 48672 times.
✗ Branch 2 not taken.
|
48672 | A[scv.dofIndex()][dofIndex][eqIdx][pvIdx] += partialDerivs[scv.indexInElement()][eqIdx]; |
| 398 | |||
| 399 | // restore the original coupling context | ||
| 400 |
1/2✓ Branch 1 taken 18720 times.
✗ Branch 2 not taken.
|
18720 | updateContext(couplingElemIdx, dofIndex, origPriVars, pvIdx); |
| 401 | } | ||
| 402 | } | ||
| 403 | } | ||
| 404 | 4273 | } | |
| 405 | |||
| 406 | //! store which scvfs coincide with facet element | ||
| 407 | std::vector<bool> bulkScvfIsOnFacetElement_; | ||
| 408 | |||
| 409 | //! store shared_ptr to coupling mapper | ||
| 410 | std::shared_ptr< CouplingMapper > couplingMapperPtr_; | ||
| 411 | }; | ||
| 412 | |||
| 413 | } // end namespace Dumux | ||
| 414 | |||
| 415 | #endif | ||
| 416 |