| 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 StokesDarcyCoupling | ||
| 10 | * \copydoc Dumux::StokesDarcyCouplingManager | ||
| 11 | */ | ||
| 12 | |||
| 13 | #ifndef DUMUX_STOKES_DARCY_COUPLINGMANAGER_HH | ||
| 14 | #define DUMUX_STOKES_DARCY_COUPLINGMANAGER_HH | ||
| 15 | |||
| 16 | #include <utility> | ||
| 17 | #include <memory> | ||
| 18 | |||
| 19 | #include <dune/common/float_cmp.hh> | ||
| 20 | #include <dune/common/exceptions.hh> | ||
| 21 | #include <dumux/common/properties.hh> | ||
| 22 | #include <dumux/common/numeqvector.hh> | ||
| 23 | #include <dumux/multidomain/staggeredcouplingmanager.hh> | ||
| 24 | #include <dumux/discretization/staggered/elementsolution.hh> | ||
| 25 | |||
| 26 | #include "couplingdata.hh" | ||
| 27 | #include "couplingmapper.hh" | ||
| 28 | |||
| 29 | namespace Dumux { | ||
| 30 | |||
| 31 | /*! | ||
| 32 | * \ingroup StokesDarcyCoupling | ||
| 33 | * \brief Coupling manager for Stokes and Darcy domains with equal dimension. | ||
| 34 | */ | ||
| 35 | template<class MDTraits> | ||
| 36 | class StokesDarcyCouplingManager | ||
| 37 | : public StaggeredCouplingManager<MDTraits> | ||
| 38 | { | ||
| 39 | using Scalar = typename MDTraits::Scalar; | ||
| 40 | using ParentType = StaggeredCouplingManager<MDTraits>; | ||
| 41 | |||
| 42 | public: | ||
| 43 | static constexpr auto stokesFaceIdx = typename MDTraits::template SubDomain<0>::Index(); | ||
| 44 | static constexpr auto stokesCellCenterIdx = typename MDTraits::template SubDomain<1>::Index(); | ||
| 45 | static constexpr auto stokesIdx = stokesFaceIdx; | ||
| 46 | static constexpr auto darcyIdx = typename MDTraits::template SubDomain<2>::Index(); | ||
| 47 | |||
| 48 | private: | ||
| 49 | |||
| 50 | using SolutionVector = typename MDTraits::SolutionVector; | ||
| 51 | |||
| 52 | // obtain the type tags of the sub problems | ||
| 53 | using StokesTypeTag = typename MDTraits::template SubDomain<0>::TypeTag; | ||
| 54 | using DarcyTypeTag = typename MDTraits::template SubDomain<2>::TypeTag; | ||
| 55 | |||
| 56 | using CouplingStencils = std::unordered_map<std::size_t, std::vector<std::size_t> >; | ||
| 57 | using CouplingStencil = CouplingStencils::mapped_type; | ||
| 58 | |||
| 59 | // the sub domain type tags | ||
| 60 | template<std::size_t id> | ||
| 61 | using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag; | ||
| 62 | |||
| 63 | static constexpr bool isCompositional = GetPropType<SubDomainTypeTag<0>, Properties::ModelTraits>::numFluidComponents()> 1; | ||
| 64 | |||
| 65 | template<std::size_t id> using GridView = typename GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>::GridView; | ||
| 66 | template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>; | ||
| 67 | template<std::size_t id> using PrimaryVariables = typename MDTraits::template SubDomain<id>::PrimaryVariables; | ||
| 68 | template<std::size_t id> using NumEqVector = Dumux::NumEqVector<PrimaryVariables<id>>; | ||
| 69 | template<std::size_t id> using ElementVolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::LocalView; | ||
| 70 | template<std::size_t id> using GridVolumeVariables = GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>; | ||
| 71 | template<std::size_t id> using VolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::VolumeVariables; | ||
| 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 ElementBoundaryTypes = GetPropType<SubDomainTypeTag<id>, Properties::ElementBoundaryTypes>; | ||
| 75 | template<std::size_t id> using ElementFluxVariablesCache = typename GetPropType<SubDomainTypeTag<id>, Properties::GridFluxVariablesCache>::LocalView; | ||
| 76 | template<std::size_t id> using GridVariables = GetPropType<SubDomainTypeTag<id>, Properties::GridVariables>; | ||
| 77 | template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity; | ||
| 78 | template<std::size_t id> using SubControlVolumeFace = typename FVElementGeometry<id>::SubControlVolumeFace; | ||
| 79 | |||
| 80 | using CellCenterSolutionVector = GetPropType<StokesTypeTag, Properties::CellCenterSolutionVector>; | ||
| 81 | |||
| 82 | using VelocityVector = typename Element<stokesIdx>::Geometry::GlobalCoordinate; | ||
| 83 | |||
| 84 | using CouplingMapper = StokesDarcyCouplingMapper; | ||
| 85 | |||
| 86 | 22680 | struct StationaryStokesCouplingContext | |
| 87 | { | ||
| 88 | Element<darcyIdx> element; | ||
| 89 | FVElementGeometry<darcyIdx> fvGeometry; | ||
| 90 | std::size_t darcyScvfIdx; | ||
| 91 | std::size_t stokesScvfIdx; | ||
| 92 | VolumeVariables<darcyIdx> volVars; | ||
| 93 | }; | ||
| 94 | |||
| 95 | struct StationaryDarcyCouplingContext | ||
| 96 | { | ||
| 97 | Element<stokesIdx> element; | ||
| 98 | FVElementGeometry<stokesIdx> fvGeometry; | ||
| 99 | std::size_t stokesScvfIdx; | ||
| 100 | std::size_t darcyScvfIdx; | ||
| 101 | VelocityVector velocity; | ||
| 102 | VolumeVariables<stokesIdx> volVars; | ||
| 103 | }; | ||
| 104 | public: | ||
| 105 | |||
| 106 | using CouplingData = StokesDarcyCouplingData<MDTraits, StokesDarcyCouplingManager<MDTraits>>; | ||
| 107 | |||
| 108 | //! Constructor | ||
| 109 | 17 | StokesDarcyCouplingManager(std::shared_ptr<const GridGeometry<stokesIdx>> stokesFvGridGeometry, | |
| 110 | std::shared_ptr<const GridGeometry<darcyIdx>> darcyFvGridGeometry) | ||
| 111 | 17 | { } | |
| 112 | |||
| 113 | /*! | ||
| 114 | * \brief Methods to be accessed by main | ||
| 115 | */ | ||
| 116 | // \{ | ||
| 117 | |||
| 118 | //! Initialize the coupling manager | ||
| 119 | 17 | void init(std::shared_ptr<const Problem<stokesIdx>> stokesProblem, | |
| 120 | std::shared_ptr<const Problem<darcyIdx>> darcyProblem, | ||
| 121 | const SolutionVector& curSol) | ||
| 122 | { | ||
| 123 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 17 times.
|
34 | if (Dune::FloatCmp::ne(stokesProblem->gravity(), darcyProblem->spatialParams().gravity({}))) |
| 124 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Both models must use the same gravity vector"); | |
| 125 | |||
| 126 | 17 | this->setSubProblems(std::make_tuple(stokesProblem, stokesProblem, darcyProblem)); | |
| 127 | 17 | this->updateSolution(curSol); | |
| 128 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 17 times.
|
17 | couplingData_ = std::make_shared<CouplingData>(*this); |
| 129 | 17 | computeStencils(); | |
| 130 | 17 | } | |
| 131 | |||
| 132 | // \} | ||
| 133 | |||
| 134 | //! Prepare the coupling stencils | ||
| 135 | 17 | void computeStencils() | |
| 136 | { | ||
| 137 | 17 | couplingMapper_.computeCouplingMapsAndStencils(*this, | |
| 138 | 17 | darcyToStokesCellCenterCouplingStencils_, | |
| 139 | 17 | darcyToStokesFaceCouplingStencils_, | |
| 140 | 17 | stokesCellCenterCouplingStencils_, | |
| 141 | 17 | stokesFaceCouplingStencils_); | |
| 142 | |||
| 143 |
2/2✓ Branch 0 taken 756 times.
✓ Branch 1 taken 17 times.
|
773 | for(auto&& stencil : darcyToStokesCellCenterCouplingStencils_) |
| 144 | 756 | removeDuplicates_(stencil.second); | |
| 145 |
2/2✓ Branch 0 taken 756 times.
✓ Branch 1 taken 17 times.
|
773 | for(auto&& stencil : darcyToStokesFaceCouplingStencils_) |
| 146 | 756 | removeDuplicates_(stencil.second); | |
| 147 |
2/2✓ Branch 0 taken 756 times.
✓ Branch 1 taken 17 times.
|
773 | for(auto&& stencil : stokesCellCenterCouplingStencils_) |
| 148 | 756 | removeDuplicates_(stencil.second); | |
| 149 |
2/2✓ Branch 0 taken 756 times.
✓ Branch 1 taken 17 times.
|
773 | for(auto&& stencil : stokesFaceCouplingStencils_) |
| 150 | 756 | removeDuplicates_(stencil.second); | |
| 151 | 17 | } | |
| 152 | |||
| 153 | /*! | ||
| 154 | * \brief Methods to be accessed by the assembly | ||
| 155 | */ | ||
| 156 | // \{ | ||
| 157 | |||
| 158 | using ParentType::evalCouplingResidual; | ||
| 159 | |||
| 160 | /*! | ||
| 161 | * \brief prepares all data and variables that are necessary to evaluate the residual of an Darcy element (i.e. Darcy information) | ||
| 162 | */ | ||
| 163 | template<std::size_t i, class Assembler, std::enable_if_t<(i == stokesCellCenterIdx || i == stokesFaceIdx), int> = 0> | ||
| 164 | 826080 | void bindCouplingContext(Dune::index_constant<i> domainI, const Element<stokesCellCenterIdx>& element, const Assembler& assembler) const | |
| 165 | 826080 | { bindCouplingContext(domainI, element); } | |
| 166 | |||
| 167 | /*! | ||
| 168 | * \brief prepares all data and variables that are necessary to evaluate the residual of an Darcy element (i.e. Darcy information) | ||
| 169 | */ | ||
| 170 | template<std::size_t i, std::enable_if_t<(i == stokesCellCenterIdx || i == stokesFaceIdx), int> = 0> | ||
| 171 | 1652160 | void bindCouplingContext(Dune::index_constant<i> domainI, const Element<stokesCellCenterIdx>& element) const | |
| 172 | { | ||
| 173 |
2/2✓ Branch 0 taken 22680 times.
✓ Branch 1 taken 803400 times.
|
1652160 | stokesCouplingContext_.clear(); |
| 174 | |||
| 175 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 826080 times.
|
1652160 | const auto stokesElementIdx = this->problem(stokesIdx).gridGeometry().elementMapper().index(element); |
| 176 | 1652160 | boundStokesElemIdx_ = stokesElementIdx; | |
| 177 | |||
| 178 | // do nothing if the element is not coupled to the other domain | ||
| 179 | 1652160 | if(!couplingMapper_.stokesElementToDarcyElementMap().count(stokesElementIdx)) | |
| 180 | 1606800 | return; | |
| 181 | |||
| 182 | // prepare the coupling context | ||
| 183 | 45360 | const auto& darcyIndices = couplingMapper_.stokesElementToDarcyElementMap().at(stokesElementIdx); | |
| 184 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 22680 times.
|
45360 | auto darcyFvGeometry = localView(this->problem(darcyIdx).gridGeometry()); |
| 185 |
2/2✓ Branch 0 taken 22680 times.
✓ Branch 1 taken 22680 times.
|
90720 | for(auto&& indices : darcyIndices) |
| 186 | { | ||
| 187 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 22680 times.
✓ Branch 3 taken 22680 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 22680 times.
✗ Branch 7 not taken.
|
45360 | const auto& darcyElement = this->problem(darcyIdx).gridGeometry().boundingBoxTree().entitySet().entity(indices.eIdx); |
| 188 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 22680 times.
|
45360 | darcyFvGeometry.bindElement(darcyElement); |
| 189 | 45360 | const auto& scv = (*scvs(darcyFvGeometry).begin()); | |
| 190 | |||
| 191 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 22680 times.
✓ Branch 4 taken 22680 times.
✗ Branch 5 not taken.
|
45360 | const auto darcyElemSol = elementSolution(darcyElement, this->curSol(darcyIdx), this->problem(darcyIdx).gridGeometry()); |
| 192 | 45360 | VolumeVariables<darcyIdx> darcyVolVars; | |
| 193 |
1/2✓ Branch 1 taken 22680 times.
✗ Branch 2 not taken.
|
45360 | darcyVolVars.update(darcyElemSol, this->problem(darcyIdx), darcyElement, scv); |
| 194 | |||
| 195 | // add the context | ||
| 196 |
1/2✓ Branch 1 taken 22680 times.
✗ Branch 2 not taken.
|
90720 | stokesCouplingContext_.push_back({darcyElement, darcyFvGeometry, indices.scvfIdx, indices.flipScvfIdx, darcyVolVars}); |
| 197 | } | ||
| 198 | 90720 | } | |
| 199 | |||
| 200 | /*! | ||
| 201 | * \brief prepares all data and variables that are necessary to evaluate the residual of an Darcy element (i.e. Stokes information) | ||
| 202 | */ | ||
| 203 | template<class Assembler> | ||
| 204 | 393520 | void bindCouplingContext(Dune::index_constant<darcyIdx> domainI, const Element<darcyIdx>& element, const Assembler& assembler) const | |
| 205 | 393520 | { bindCouplingContext(domainI, element); } | |
| 206 | |||
| 207 | /*! | ||
| 208 | * \brief prepares all data and variables that are necessary to evaluate the residual of an Darcy element (i.e. Stokes information) | ||
| 209 | */ | ||
| 210 | 396340 | void bindCouplingContext(Dune::index_constant<darcyIdx> domainI, const Element<darcyIdx>& element) const | |
| 211 | { | ||
| 212 |
2/2✓ Branch 0 taken 14143 times.
✓ Branch 1 taken 382197 times.
|
396340 | darcyCouplingContext_.clear(); |
| 213 | |||
| 214 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 396340 times.
|
396340 | const auto darcyElementIdx = this->problem(darcyIdx).gridGeometry().elementMapper().index(element); |
| 215 | 396340 | boundDarcyElemIdx_ = darcyElementIdx; | |
| 216 | |||
| 217 | // do nothing if the element is not coupled to the other domain | ||
| 218 | 396340 | if(!couplingMapper_.darcyElementToStokesElementMap().count(darcyElementIdx)) | |
| 219 | 382180 | return; | |
| 220 | |||
| 221 | // prepare the coupling context | ||
| 222 | 14160 | const auto& stokesElementIndices = couplingMapper_.darcyElementToStokesElementMap().at(darcyElementIdx); | |
| 223 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 14160 times.
|
14160 | auto stokesFvGeometry = localView(this->problem(stokesIdx).gridGeometry()); |
| 224 | |||
| 225 |
2/2✓ Branch 0 taken 14160 times.
✓ Branch 1 taken 14160 times.
|
28320 | for(auto&& indices : stokesElementIndices) |
| 226 | { | ||
| 227 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 14160 times.
|
14160 | const auto& stokesElement = this->problem(stokesIdx).gridGeometry().boundingBoxTree().entitySet().entity(indices.eIdx); |
| 228 | 14160 | stokesFvGeometry.bindElement(stokesElement); | |
| 229 | |||
| 230 | 14160 | VelocityVector faceVelocity(0.0); | |
| 231 | |||
| 232 |
4/4✓ Branch 0 taken 14160 times.
✓ Branch 1 taken 42480 times.
✓ Branch 2 taken 56640 times.
✓ Branch 3 taken 14160 times.
|
70800 | for(const auto& scvf : scvfs(stokesFvGeometry)) |
| 233 | { | ||
| 234 |
2/2✓ Branch 0 taken 14160 times.
✓ Branch 1 taken 42480 times.
|
56640 | if(scvf.index() == indices.scvfIdx) |
| 235 | 14160 | faceVelocity[scvf.directionIndex()] = this->curSol(stokesFaceIdx)[scvf.dofIndex()]; | |
| 236 | } | ||
| 237 | |||
| 238 | using PriVarsType = typename VolumeVariables<stokesCellCenterIdx>::PrimaryVariables; | ||
| 239 | 14160 | const auto& cellCenterPriVars = this->curSol(stokesCellCenterIdx)[indices.eIdx]; | |
| 240 | 14160 | const auto elemSol = makeElementSolutionFromCellCenterPrivars<PriVarsType>(cellCenterPriVars); | |
| 241 | |||
| 242 | 14160 | VolumeVariables<stokesIdx> stokesVolVars; | |
| 243 |
3/3✓ Branch 0 taken 13120 times.
✓ Branch 1 taken 14160 times.
✓ Branch 2 taken 1040 times.
|
28320 | for(const auto& scv : scvs(stokesFvGeometry)) |
| 244 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 14160 times.
✓ Branch 3 taken 14160 times.
✗ Branch 4 not taken.
|
19320 | stokesVolVars.update(elemSol, this->problem(stokesIdx), stokesElement, scv); |
| 245 | |||
| 246 | // add the context | ||
| 247 |
3/8✓ Branch 0 taken 14160 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 14160 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 14160 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
28320 | darcyCouplingContext_.push_back({stokesElement, stokesFvGeometry, indices.scvfIdx, indices.flipScvfIdx, faceVelocity, stokesVolVars}); |
| 248 | } | ||
| 249 | } | ||
| 250 | |||
| 251 | using StaggeredCouplingManager<MDTraits>::updateCouplingContext; | ||
| 252 | |||
| 253 | //! \copydoc StaggeredCouplingManager::updateCouplingContext | ||
| 254 | template<class LocalAssemblerI> | ||
| 255 | 1161792 | void updateCouplingContext(Dune::index_constant<darcyIdx> domainI, | |
| 256 | const LocalAssemblerI& localAssemblerI, | ||
| 257 | Dune::index_constant<darcyIdx> domainJ, | ||
| 258 | std::size_t dofIdxGlobalJ, | ||
| 259 | const PrimaryVariables<darcyIdx>& priVarsJ, | ||
| 260 | int pvIdxJ) | ||
| 261 | { | ||
| 262 | 1161792 | this->curSol(domainJ)[dofIdxGlobalJ][pvIdxJ] = priVarsJ[pvIdxJ]; | |
| 263 | } | ||
| 264 | |||
| 265 | //! \copydoc StaggeredCouplingManager::updateCouplingContext | ||
| 266 | template<class LocalAssemblerI> | ||
| 267 | 42664 | void updateCouplingContext(Dune::index_constant<darcyIdx> domainI, | |
| 268 | const LocalAssemblerI& localAssemblerI, | ||
| 269 | Dune::index_constant<stokesCellCenterIdx> domainJ, | ||
| 270 | const std::size_t dofIdxGlobalJ, | ||
| 271 | const PrimaryVariables<stokesCellCenterIdx>& priVarsJ, | ||
| 272 | int pvIdxJ) | ||
| 273 | { | ||
| 274 | 42664 | this->curSol(domainJ)[dofIdxGlobalJ] = priVarsJ; | |
| 275 | |||
| 276 |
2/2✓ Branch 0 taken 42664 times.
✓ Branch 1 taken 42664 times.
|
85328 | for (auto& data : darcyCouplingContext_) |
| 277 | { | ||
| 278 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 42664 times.
|
42664 | const auto stokesElemIdx = this->problem(stokesIdx).gridGeometry().elementMapper().index(data.element); |
| 279 | |||
| 280 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 42664 times.
|
42664 | if(stokesElemIdx != dofIdxGlobalJ) |
| 281 | ✗ | continue; | |
| 282 | |||
| 283 | using PriVarsType = typename VolumeVariables<stokesCellCenterIdx>::PrimaryVariables; | ||
| 284 | 42664 | const auto elemSol = makeElementSolutionFromCellCenterPrivars<PriVarsType>(priVarsJ); | |
| 285 | |||
| 286 |
2/2✓ Branch 0 taken 42664 times.
✓ Branch 1 taken 42664 times.
|
85328 | for(const auto& scv : scvs(data.fvGeometry)) |
| 287 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 42664 times.
✓ Branch 3 taken 42664 times.
✗ Branch 4 not taken.
|
49744 | data.volVars.update(elemSol, this->problem(stokesIdx), data.element, scv); |
| 288 | } | ||
| 289 | 42664 | } | |
| 290 | |||
| 291 | //! \copydoc StaggeredCouplingManager::updateCouplingContext | ||
| 292 | template<class LocalAssemblerI> | ||
| 293 | 22680 | void updateCouplingContext(Dune::index_constant<darcyIdx> domainI, | |
| 294 | const LocalAssemblerI& localAssemblerI, | ||
| 295 | Dune::index_constant<stokesFaceIdx> domainJ, | ||
| 296 | const std::size_t dofIdxGlobalJ, | ||
| 297 | const PrimaryVariables<stokesFaceIdx>& priVarsJ, | ||
| 298 | int pvIdxJ) | ||
| 299 | { | ||
| 300 | 22680 | this->curSol(domainJ)[dofIdxGlobalJ] = priVarsJ; | |
| 301 | |||
| 302 |
2/2✓ Branch 0 taken 22680 times.
✓ Branch 1 taken 22680 times.
|
45360 | for (auto& data : darcyCouplingContext_) |
| 303 | { | ||
| 304 |
4/4✓ Branch 0 taken 22680 times.
✓ Branch 1 taken 68040 times.
✓ Branch 2 taken 90720 times.
✓ Branch 3 taken 22680 times.
|
113400 | for(const auto& scvf : scvfs(data.fvGeometry)) |
| 305 | { | ||
| 306 |
2/2✓ Branch 0 taken 22680 times.
✓ Branch 1 taken 68040 times.
|
90720 | if(scvf.dofIndex() == dofIdxGlobalJ) |
| 307 | 22680 | data.velocity[scvf.directionIndex()] = priVarsJ; | |
| 308 | } | ||
| 309 | } | ||
| 310 | 22680 | } | |
| 311 | |||
| 312 | //! \copydoc StaggeredCouplingManager::updateCouplingContext | ||
| 313 | template<std::size_t i, class LocalAssemblerI, std::enable_if_t<(i == stokesCellCenterIdx || i == stokesFaceIdx), int> = 0> | ||
| 314 | 183296 | void updateCouplingContext(Dune::index_constant<i> domainI, | |
| 315 | const LocalAssemblerI& localAssemblerI, | ||
| 316 | Dune::index_constant<darcyIdx> domainJ, | ||
| 317 | const std::size_t dofIdxGlobalJ, | ||
| 318 | const PrimaryVariables<darcyIdx>& priVarsJ, | ||
| 319 | int pvIdxJ) | ||
| 320 | { | ||
| 321 | 183296 | this->curSol(domainJ)[dofIdxGlobalJ] = priVarsJ; | |
| 322 | |||
| 323 |
2/2✓ Branch 0 taken 91648 times.
✓ Branch 1 taken 91648 times.
|
366592 | for (auto& data : stokesCouplingContext_) |
| 324 | { | ||
| 325 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 91648 times.
|
183296 | const auto darcyElemIdx = this->problem(darcyIdx).gridGeometry().elementMapper().index(data.element); |
| 326 | |||
| 327 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 91648 times.
|
183296 | if(darcyElemIdx != dofIdxGlobalJ) |
| 328 | ✗ | continue; | |
| 329 | |||
| 330 | 183296 | const auto darcyElemSol = elementSolution(data.element, this->curSol(darcyIdx), this->problem(darcyIdx).gridGeometry()); | |
| 331 | |||
| 332 |
2/2✓ Branch 0 taken 91648 times.
✓ Branch 1 taken 91648 times.
|
366592 | for(const auto& scv : scvs(data.fvGeometry)) |
| 333 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 91648 times.
|
183296 | data.volVars.update(darcyElemSol, this->problem(darcyIdx), data.element, scv); |
| 334 | } | ||
| 335 | 183296 | } | |
| 336 | |||
| 337 | // \} | ||
| 338 | |||
| 339 | /*! | ||
| 340 | * \brief Access the coupling data | ||
| 341 | */ | ||
| 342 | 1710850 | const auto& couplingData() const | |
| 343 | { | ||
| 344 | 1315724 | return *couplingData_; | |
| 345 | } | ||
| 346 | |||
| 347 | /*! | ||
| 348 | * \brief Access the coupling context needed for the Stokes domain | ||
| 349 | */ | ||
| 350 | 1605606 | const auto& stokesCouplingContext(const Element<stokesIdx>& element, const SubControlVolumeFace<stokesIdx>& scvf) const | |
| 351 | { | ||
| 352 |
2/4✓ Branch 0 taken 1605606 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1605606 times.
|
1605606 | if (stokesCouplingContext_.empty() || boundStokesElemIdx_ != scvf.insideScvIdx()) |
| 353 | ✗ | bindCouplingContext(stokesIdx, element); | |
| 354 | |||
| 355 |
2/4✓ Branch 0 taken 1605606 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1605606 times.
✗ Branch 3 not taken.
|
1605606 | for(const auto& context : stokesCouplingContext_) |
| 356 | { | ||
| 357 |
1/2✓ Branch 0 taken 1605606 times.
✗ Branch 1 not taken.
|
1605606 | if(scvf.index() == context.stokesScvfIdx) |
| 358 | 1605606 | return context; | |
| 359 | } | ||
| 360 | |||
| 361 | ✗ | DUNE_THROW(Dune::InvalidStateException, "No coupling context found at scvf " << scvf.center()); | |
| 362 | } | ||
| 363 | |||
| 364 | /*! | ||
| 365 | * \brief Access the coupling context needed for the Darcy domain | ||
| 366 | */ | ||
| 367 | 105244 | const auto& darcyCouplingContext(const Element<darcyIdx>& element, const SubControlVolumeFace<darcyIdx>& scvf) const | |
| 368 | { | ||
| 369 |
4/4✓ Branch 0 taken 105232 times.
✓ Branch 1 taken 12 times.
✓ Branch 2 taken 2808 times.
✓ Branch 3 taken 102424 times.
|
105244 | if (darcyCouplingContext_.empty() || boundDarcyElemIdx_ != scvf.insideScvIdx()) |
| 370 | 2820 | bindCouplingContext(darcyIdx, element); | |
| 371 | |||
| 372 |
2/4✓ Branch 0 taken 105244 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 105244 times.
✗ Branch 3 not taken.
|
105244 | for(const auto& context : darcyCouplingContext_) |
| 373 | { | ||
| 374 |
1/2✓ Branch 0 taken 105244 times.
✗ Branch 1 not taken.
|
105244 | if(scvf.index() == context.darcyScvfIdx) |
| 375 | 105244 | return context; | |
| 376 | } | ||
| 377 | |||
| 378 | ✗ | DUNE_THROW(Dune::InvalidStateException, "No coupling context found at scvf " << scvf.center()); | |
| 379 | } | ||
| 380 | |||
| 381 | /*! | ||
| 382 | * \brief The coupling stencils | ||
| 383 | */ | ||
| 384 | // \{ | ||
| 385 | |||
| 386 | using StaggeredCouplingManager<MDTraits>::couplingStencil; | ||
| 387 | |||
| 388 | /*! | ||
| 389 | * \brief The Stokes cell center coupling stencil w.r.t. Darcy DOFs. | ||
| 390 | * | ||
| 391 | * \param domainI The Stokes domain index. | ||
| 392 | * \param elementI The Sokes domain element. | ||
| 393 | * \param domainJ The Darcy domain index. | ||
| 394 | */ | ||
| 395 | 484096 | const CouplingStencil& couplingStencil(Dune::index_constant<stokesCellCenterIdx> domainI, | |
| 396 | const Element<stokesIdx>& elementI, | ||
| 397 | Dune::index_constant<darcyIdx> domainJ) const | ||
| 398 | { | ||
| 399 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 484096 times.
|
484096 | const auto eIdx = this->problem(domainI).gridGeometry().elementMapper().index(elementI); |
| 400 | 484096 | if(stokesCellCenterCouplingStencils_.count(eIdx)) | |
| 401 | 12096 | return stokesCellCenterCouplingStencils_.at(eIdx); | |
| 402 | else | ||
| 403 | 472000 | return emptyStencil_; | |
| 404 | } | ||
| 405 | |||
| 406 | /*! | ||
| 407 | * \brief The coupling stencil of domain I, i.e. which domain J DOFs | ||
| 408 | * the given domain I element's residual depends on. | ||
| 409 | * | ||
| 410 | * \param domainI the index of the domain in which the given element lives. | ||
| 411 | * \param elementI the coupled element of domainI | ||
| 412 | * \param domainJ the domain index of the coupled domain | ||
| 413 | */ | ||
| 414 | template<std::size_t i, std::size_t j> | ||
| 415 | const CouplingStencil& couplingStencil(Dune::index_constant<i> domainI, | ||
| 416 | const Element<i>& elementI, | ||
| 417 | Dune::index_constant<j> domainJ) const | ||
| 418 | { return emptyStencil_; } | ||
| 419 | |||
| 420 | /*! | ||
| 421 | * \brief Return the Stokes cell indices that influence the residual of | ||
| 422 | * an element in the Darcy domain. | ||
| 423 | * | ||
| 424 | * \param domainI The darcy domain index. | ||
| 425 | * \param elementI The element in the Darcy domain. | ||
| 426 | * \param domainJ the domain index of the Stokes domain | ||
| 427 | */ | ||
| 428 | 464448 | const CouplingStencil& couplingStencil(Dune::index_constant<darcyIdx> domainI, | |
| 429 | const Element<darcyIdx>& elementI, | ||
| 430 | Dune::index_constant<stokesCellCenterIdx> domainJ) const | ||
| 431 | { | ||
| 432 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 464448 times.
|
464448 | const auto eIdx = this->problem(domainI).gridGeometry().elementMapper().index(elementI); |
| 433 | 464448 | if(darcyToStokesCellCenterCouplingStencils_.count(eIdx)) | |
| 434 | 12096 | return darcyToStokesCellCenterCouplingStencils_.at(eIdx); | |
| 435 | else | ||
| 436 | 452352 | return emptyStencil_; | |
| 437 | } | ||
| 438 | |||
| 439 | /*! | ||
| 440 | * \brief Return the Stokes face indices that influence the residual of | ||
| 441 | * an element in the Darcy domain. | ||
| 442 | * | ||
| 443 | * \param domainI The darcy domain index. | ||
| 444 | * \param elementI The element in the Darcy domain. | ||
| 445 | * \param domainJ the domain index of the Stokes domain | ||
| 446 | */ | ||
| 447 | 464448 | const CouplingStencil& couplingStencil(Dune::index_constant<darcyIdx> domainI, | |
| 448 | const Element<darcyIdx>& elementI, | ||
| 449 | Dune::index_constant<stokesFaceIdx> domainJ) const | ||
| 450 | { | ||
| 451 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 464448 times.
|
464448 | const auto eIdx = this->problem(domainI).gridGeometry().elementMapper().index(elementI); |
| 452 | 464448 | if (darcyToStokesFaceCouplingStencils_.count(eIdx)) | |
| 453 | 12096 | return darcyToStokesFaceCouplingStencils_.at(eIdx); | |
| 454 | else | ||
| 455 | 452352 | return emptyStencil_; | |
| 456 | } | ||
| 457 | |||
| 458 | /*! | ||
| 459 | * \brief Return the dof indices of a subdomain that influence the residual of | ||
| 460 | * a sub-control volume face of the Stokes domain. | ||
| 461 | * | ||
| 462 | * \param domainI the index of the domain in which the given element lives. | ||
| 463 | * \param scvfI the coupled sub-control volume face of the Stokes domain | ||
| 464 | * \param domainJ the domain index of the coupled domain | ||
| 465 | */ | ||
| 466 | template<std::size_t i, std::size_t j> | ||
| 467 | const CouplingStencil& couplingStencil(Dune::index_constant<i> domainI, | ||
| 468 | const SubControlVolumeFace<stokesIdx>& scvfI, | ||
| 469 | Dune::index_constant<j> domainJ) const | ||
| 470 | { return emptyStencil_; } | ||
| 471 | |||
| 472 | /*! | ||
| 473 | * \brief The coupling stencil of a Stokes face w.r.t. Darcy DOFs. | ||
| 474 | * | ||
| 475 | * \param domainI the index of the Stokes domain | ||
| 476 | * \param scvfI the coupled subcontrolvolume face of the Stokes domain | ||
| 477 | * \param domainJ the index of the Darcy domain | ||
| 478 | */ | ||
| 479 | 1936384 | const CouplingStencil& couplingStencil(Dune::index_constant<stokesFaceIdx> domainI, | |
| 480 | const SubControlVolumeFace<stokesIdx>& scvfI, | ||
| 481 | Dune::index_constant<darcyIdx> domainJ) const | ||
| 482 | { | ||
| 483 | 1936384 | const auto faceDofIdx = scvfI.dofIndex(); | |
| 484 | 1936384 | if(stokesFaceCouplingStencils_.count(faceDofIdx)) | |
| 485 | 12096 | return stokesFaceCouplingStencils_.at(faceDofIdx); | |
| 486 | else | ||
| 487 | 1924288 | return emptyStencil_; | |
| 488 | } | ||
| 489 | |||
| 490 | // \} | ||
| 491 | |||
| 492 | /*! | ||
| 493 | * \brief There are no additional dof dependencies | ||
| 494 | */ | ||
| 495 | template<class IdType> | ||
| 496 | const std::vector<std::size_t>& getAdditionalDofDependencies(IdType id, std::size_t stokesElementIdx) const | ||
| 497 | { return emptyStencil_; } | ||
| 498 | |||
| 499 | /*! | ||
| 500 | * \brief There are no additional dof dependencies | ||
| 501 | */ | ||
| 502 | template<class IdType> | ||
| 503 | const std::vector<std::size_t>& getAdditionalDofDependenciesInverse(IdType id, std::size_t darcyElementIdx) const | ||
| 504 | { return emptyStencil_; } | ||
| 505 | |||
| 506 | /*! | ||
| 507 | * \brief Returns whether a given free flow scvf is coupled to the other domain | ||
| 508 | */ | ||
| 509 | 5213405 | bool isCoupledEntity(Dune::index_constant<stokesIdx>, const SubControlVolumeFace<stokesFaceIdx>& scvf) const | |
| 510 | { | ||
| 511 | 10426810 | return stokesFaceCouplingStencils_.count(scvf.dofIndex()); | |
| 512 | } | ||
| 513 | |||
| 514 | /*! | ||
| 515 | * \brief Returns whether a given free flow scvf is coupled to the other domain | ||
| 516 | */ | ||
| 517 | 467700 | bool isCoupledEntity(Dune::index_constant<darcyIdx>, const SubControlVolumeFace<darcyIdx>& scvf) const | |
| 518 | { | ||
| 519 |
9/10✓ Branch 0 taken 81164 times.
✓ Branch 1 taken 95642 times.
✓ Branch 2 taken 94444 times.
✓ Branch 3 taken 147562 times.
✓ Branch 4 taken 12880 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 11760 times.
✓ Branch 7 taken 11928 times.
✓ Branch 8 taken 3080 times.
✓ Branch 9 taken 9240 times.
|
467700 | return couplingMapper_.isCoupledDarcyScvf(scvf.index()); |
| 520 | } | ||
| 521 | |||
| 522 | protected: | ||
| 523 | |||
| 524 | //! Return a reference to an empty stencil | ||
| 525 | std::vector<std::size_t>& emptyStencil() | ||
| 526 | { return emptyStencil_; } | ||
| 527 | |||
| 528 | 3024 | void removeDuplicates_(std::vector<std::size_t>& stencil) | |
| 529 | { | ||
| 530 | 3024 | std::sort(stencil.begin(), stencil.end()); | |
| 531 | 3024 | stencil.erase(std::unique(stencil.begin(), stencil.end()), stencil.end()); | |
| 532 | 3024 | } | |
| 533 | |||
| 534 | private: | ||
| 535 | |||
| 536 | std::vector<bool> isCoupledDarcyDof_; | ||
| 537 | std::shared_ptr<CouplingData> couplingData_; | ||
| 538 | |||
| 539 | std::unordered_map<std::size_t, std::vector<std::size_t> > stokesCellCenterCouplingStencils_; | ||
| 540 | std::unordered_map<std::size_t, std::vector<std::size_t> > stokesFaceCouplingStencils_; | ||
| 541 | std::unordered_map<std::size_t, std::vector<std::size_t> > darcyToStokesCellCenterCouplingStencils_; | ||
| 542 | std::unordered_map<std::size_t, std::vector<std::size_t> > darcyToStokesFaceCouplingStencils_; | ||
| 543 | std::vector<std::size_t> emptyStencil_; | ||
| 544 | |||
| 545 | //////////////////////////////////////////////////////////////////////////// | ||
| 546 | //! The coupling context | ||
| 547 | //////////////////////////////////////////////////////////////////////////// | ||
| 548 | mutable std::vector<StationaryStokesCouplingContext> stokesCouplingContext_; | ||
| 549 | mutable std::vector<StationaryDarcyCouplingContext> darcyCouplingContext_; | ||
| 550 | |||
| 551 | mutable std::size_t boundStokesElemIdx_; | ||
| 552 | mutable std::size_t boundDarcyElemIdx_; | ||
| 553 | |||
| 554 | CouplingMapper couplingMapper_; | ||
| 555 | }; | ||
| 556 | |||
| 557 | } // end namespace Dumux | ||
| 558 | |||
| 559 | #endif | ||
| 560 |