| 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 FreeFlowPorousMediumCoupling | ||
| 10 | * \copydoc Dumux::FFMassPMCouplingManagerStaggeredCCTpfa | ||
| 11 | */ | ||
| 12 | |||
| 13 | #ifndef DUMUX_MULTIDOMAIN_BOUNDARY_FFPM_FFMASSPM_COUPLINGMANAGER_STAGGERED_TPFA_HH | ||
| 14 | #define DUMUX_MULTIDOMAIN_BOUNDARY_FFPM_FFMASSPM_COUPLINGMANAGER_STAGGERED_TPFA_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/discretization/staggered/elementsolution.hh> | ||
| 23 | #include <dumux/multidomain/couplingmanager.hh> | ||
| 24 | #include <dumux/multidomain/boundary/darcydarcy/couplingmapper.hh> | ||
| 25 | |||
| 26 | namespace Dumux { | ||
| 27 | |||
| 28 | /*! | ||
| 29 | * \ingroup FreeFlowPorousMediumCoupling | ||
| 30 | * \brief Coupling manager for Stokes and Darcy domains with equal dimension. | ||
| 31 | * Specialization for staggered-cctpfa coupling. | ||
| 32 | */ | ||
| 33 | template<class MDTraits> | ||
| 34 | class FFMassPMCouplingManagerStaggeredCCTpfa | ||
| 35 | : public CouplingManager<MDTraits> | ||
| 36 | { | ||
| 37 | using Scalar = typename MDTraits::Scalar; | ||
| 38 | using ParentType = CouplingManager<MDTraits>; | ||
| 39 | |||
| 40 | public: | ||
| 41 | static constexpr auto freeFlowMassIndex = typename MDTraits::template SubDomain<0>::Index(); | ||
| 42 | static constexpr auto porousMediumIndex = typename MDTraits::template SubDomain<1>::Index(); | ||
| 43 | |||
| 44 | using SolutionVectorStorage = typename ParentType::SolutionVectorStorage; | ||
| 45 | private: | ||
| 46 | |||
| 47 | // obtain the type tags of the sub problems | ||
| 48 | using FreeFlowMassTypeTag = typename MDTraits::template SubDomain<freeFlowMassIndex>::TypeTag; | ||
| 49 | using PorousMediumTypeTag = typename MDTraits::template SubDomain<porousMediumIndex>::TypeTag; | ||
| 50 | |||
| 51 | using CouplingStencils = std::unordered_map<std::size_t, std::vector<std::size_t> >; | ||
| 52 | using CouplingStencil = CouplingStencils::mapped_type; | ||
| 53 | |||
| 54 | // the sub domain type tags | ||
| 55 | template<std::size_t id> | ||
| 56 | using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag; | ||
| 57 | |||
| 58 | template<std::size_t id> using GridView = typename GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>::GridView; | ||
| 59 | template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>; | ||
| 60 | template<std::size_t id> using ElementVolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::LocalView; | ||
| 61 | template<std::size_t id> using GridVolumeVariables = GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>; | ||
| 62 | template<std::size_t id> using VolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::VolumeVariables; | ||
| 63 | template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>; | ||
| 64 | template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView; | ||
| 65 | template<std::size_t id> using GridVariables = GetPropType<SubDomainTypeTag<id>, Properties::GridVariables>; | ||
| 66 | template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity; | ||
| 67 | template<std::size_t id> using PrimaryVariables = GetPropType<SubDomainTypeTag<id>, Properties::PrimaryVariables>; | ||
| 68 | template<std::size_t id> using SubControlVolumeFace = typename FVElementGeometry<id>::SubControlVolumeFace; | ||
| 69 | template<std::size_t id> using SubControlVolume = typename FVElementGeometry<id>::SubControlVolume; | ||
| 70 | |||
| 71 | using VelocityVector = typename Element<freeFlowMassIndex>::Geometry::GlobalCoordinate; | ||
| 72 | |||
| 73 | 3920 | struct FreeFlowMassCouplingContext | |
| 74 | { | ||
| 75 | Element<porousMediumIndex> element; | ||
| 76 | VolumeVariables<porousMediumIndex> volVars; | ||
| 77 | FVElementGeometry<porousMediumIndex> fvGeometry; | ||
| 78 | std::size_t dofIdx; | ||
| 79 | std::size_t freeFlowMassScvfIdx; | ||
| 80 | std::size_t porousMediumScvfIdx; | ||
| 81 | mutable VelocityVector velocity; // velocity needs to be set externally, not available in this class | ||
| 82 | }; | ||
| 83 | |||
| 84 | struct PorousMediumCouplingContext | ||
| 85 | { | ||
| 86 | Element<freeFlowMassIndex> element; | ||
| 87 | VolumeVariables<freeFlowMassIndex> volVars; | ||
| 88 | FVElementGeometry<freeFlowMassIndex> fvGeometry; | ||
| 89 | std::size_t dofIdx; | ||
| 90 | std::size_t porousMediumScvfIdx; | ||
| 91 | std::size_t freeFlowMassScvfIdx; | ||
| 92 | mutable VelocityVector velocity; // velocity needs to be set externally, not available in this class | ||
| 93 | }; | ||
| 94 | |||
| 95 | using CouplingMapper = DarcyDarcyBoundaryCouplingMapper<MDTraits>; // TODO rename/generalize class | ||
| 96 | |||
| 97 | public: | ||
| 98 | |||
| 99 | /*! | ||
| 100 | * \brief Methods to be accessed by main | ||
| 101 | */ | ||
| 102 | // \{ | ||
| 103 | |||
| 104 | //! Initialize the coupling manager | ||
| 105 | 15 | void init(std::shared_ptr<Problem<freeFlowMassIndex>> freeFlowMassProblem, | |
| 106 | std::shared_ptr<Problem<porousMediumIndex>> darcyProblem, | ||
| 107 | const SolutionVectorStorage& curSol) | ||
| 108 | { | ||
| 109 | 30 | this->setSubProblems(std::make_tuple(freeFlowMassProblem, darcyProblem)); | |
| 110 | 15 | this->attachSolution(curSol); | |
| 111 | |||
| 112 | 15 | couplingMapper_.update(*this); | |
| 113 | 15 | } | |
| 114 | |||
| 115 | // \} | ||
| 116 | |||
| 117 | /*! | ||
| 118 | * \brief Methods to be accessed by the assembly | ||
| 119 | */ | ||
| 120 | // \{ | ||
| 121 | |||
| 122 | /*! | ||
| 123 | * \brief prepares all data and variables that are necessary to evaluate the residual (called from the local assembler) | ||
| 124 | */ | ||
| 125 | template<std::size_t i, class Assembler> | ||
| 126 | 772800 | void bindCouplingContext(Dune::index_constant<i> domainI, const Element<i>& element, const Assembler& assembler) const | |
| 127 | { | ||
| 128 | 772800 | bindCouplingContext_(domainI, element); | |
| 129 | } | ||
| 130 | |||
| 131 | /*! | ||
| 132 | * \brief Update the coupling context | ||
| 133 | */ | ||
| 134 | template<std::size_t i, std::size_t j, class LocalAssemblerI> | ||
| 135 | 3122560 | void updateCouplingContext(Dune::index_constant<i> domainI, | |
| 136 | const LocalAssemblerI& localAssemblerI, | ||
| 137 | Dune::index_constant<j> domainJ, | ||
| 138 | std::size_t dofIdxGlobalJ, | ||
| 139 | const PrimaryVariables<j>& priVarsJ, | ||
| 140 | int pvIdxJ) | ||
| 141 | { | ||
| 142 | 3122560 | this->curSol(domainJ)[dofIdxGlobalJ][pvIdxJ] = priVarsJ[pvIdxJ]; | |
| 143 | |||
| 144 | // we need to update all solution-depenent components of the coupling context | ||
| 145 | // the dof of domain J has been deflected | ||
| 146 | // if domainJ == freeFlowMassIndex: update volvars in the PorousMediumCouplingContext | ||
| 147 | // if domainJ == porousMediumIndex: update volvars in the FreeFlowMassCouplingContext | ||
| 148 | // as the update is symmetric we only need to write this once | ||
| 149 | 3122560 | auto& context = std::get<1-j>(couplingContext_); | |
| 150 |
2/2✓ Branch 0 taken 786880 times.
✓ Branch 1 taken 1561280 times.
|
4696320 | for (auto& c : context) |
| 151 | { | ||
| 152 |
2/2✓ Branch 0 taken 15826 times.
✓ Branch 1 taken 771054 times.
|
1573760 | if (c.dofIdx == dofIdxGlobalJ) |
| 153 | { | ||
| 154 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 15826 times.
|
31652 | const auto elemSol = elementSolution(c.element, this->curSol(domainJ), this->problem(domainJ).gridGeometry()); |
| 155 | 31652 | const auto& scv = *scvs(c.fvGeometry).begin(); | |
| 156 | 31652 | c.volVars.update(elemSol, this->problem(domainJ), c.element, scv); | |
| 157 | } | ||
| 158 | } | ||
| 159 | 3122560 | } | |
| 160 | |||
| 161 | // \} | ||
| 162 | |||
| 163 | /*! | ||
| 164 | * \brief Access the coupling context needed for the Stokes domain | ||
| 165 | */ | ||
| 166 | template<std::size_t i> | ||
| 167 | 146480 | const auto& couplingContext(Dune::index_constant<i> domainI, | |
| 168 | const FVElementGeometry<i>& fvGeometry, | ||
| 169 | const SubControlVolumeFace<i> scvf) const | ||
| 170 | { | ||
| 171 |
2/2✓ Branch 0 taken 73227 times.
✓ Branch 1 taken 13 times.
|
146480 | auto& contexts = std::get<i>(couplingContext_); |
| 172 | |||
| 173 |
4/4✓ Branch 0 taken 73227 times.
✓ Branch 1 taken 13 times.
✓ Branch 2 taken 2667 times.
✓ Branch 3 taken 70560 times.
|
146480 | if (contexts.empty() || couplingContextBoundForElement_[i] != scvf.insideScvIdx()) |
| 174 | 5360 | bindCouplingContext_(domainI, fvGeometry); | |
| 175 | |||
| 176 |
1/2✓ Branch 0 taken 73240 times.
✗ Branch 1 not taken.
|
146480 | for (const auto& context : contexts) |
| 177 | { | ||
| 178 |
1/2✓ Branch 0 taken 73240 times.
✗ Branch 1 not taken.
|
146480 | const auto expectedScvfIdx = domainI == freeFlowMassIndex ? context.freeFlowMassScvfIdx : context.porousMediumScvfIdx; |
| 179 |
1/2✓ Branch 0 taken 73240 times.
✗ Branch 1 not taken.
|
146480 | if (scvf.index() == expectedScvfIdx) |
| 180 | 146480 | return context; | |
| 181 | } | ||
| 182 | |||
| 183 | ✗ | DUNE_THROW(Dune::InvalidStateException, "No coupling context found at scvf " << scvf.center()); | |
| 184 | } | ||
| 185 | |||
| 186 | /*! | ||
| 187 | * \brief The coupling stencils | ||
| 188 | */ | ||
| 189 | // \{ | ||
| 190 | |||
| 191 | /*! | ||
| 192 | * \brief The Stokes cell center coupling stencil w.r.t. Darcy DOFs | ||
| 193 | */ | ||
| 194 | template<std::size_t i, std::size_t j> | ||
| 195 | 2088000 | const CouplingStencil& couplingStencil(Dune::index_constant<i> domainI, | |
| 196 | const Element<i>& element, | ||
| 197 | Dune::index_constant<j> domainJ) const | ||
| 198 | { | ||
| 199 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1044000 times.
|
2088000 | const auto eIdx = this->problem(domainI).gridGeometry().elementMapper().index(element); |
| 200 | 2088000 | return couplingMapper_.couplingStencil(domainI, eIdx, domainJ); | |
| 201 | } | ||
| 202 | |||
| 203 | // \} | ||
| 204 | |||
| 205 | /*! | ||
| 206 | * \brief Returns whether a given scvf is coupled to the other domain | ||
| 207 | */ | ||
| 208 | template<std::size_t i> | ||
| 209 | 557248 | bool isCoupled(Dune::index_constant<i> domainI, const SubControlVolumeFace<i>& scvf) const | |
| 210 |
10/10✓ Branch 1 taken 26200 times.
✓ Branch 2 taken 16808 times.
✓ Branch 4 taken 66304 times.
✓ Branch 5 taken 29176 times.
✓ Branch 7 taken 47040 times.
✓ Branch 8 taken 117900 times.
✓ Branch 10 taken 1568 times.
✓ Branch 11 taken 600 times.
✓ Branch 13 taken 30340 times.
✓ Branch 14 taken 11660 times.
|
515168 | { return couplingMapper_.isCoupled(domainI, scvf); } |
| 211 | |||
| 212 | private: | ||
| 213 | /*! | ||
| 214 | * \brief Returns whether a given scvf is coupled to the other domain | ||
| 215 | */ | ||
| 216 | template<std::size_t i> | ||
| 217 | 775480 | bool isCoupledElement_(Dune::index_constant<i> domainI, std::size_t eIdx) const | |
| 218 | 1550960 | { return couplingMapper_.isCoupledElement(domainI, eIdx); } | |
| 219 | |||
| 220 | //! Return the volume variables of domain i for a given element and scv | ||
| 221 | template<std::size_t i> | ||
| 222 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 10520 times.
|
21040 | VolumeVariables<i> volVars_(Dune::index_constant<i> domainI, |
| 223 | const Element<i>& element, | ||
| 224 | const SubControlVolume<i>& scv) const | ||
| 225 | { | ||
| 226 | 21040 | VolumeVariables<i> volVars; | |
| 227 | 21040 | const auto elemSol = elementSolution( | |
| 228 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 10520 times.
|
21040 | element, this->curSol(domainI), this->problem(domainI).gridGeometry() |
| 229 | ); | ||
| 230 | 21040 | volVars.update(elemSol, this->problem(domainI), element, scv); | |
| 231 | 21040 | return volVars; | |
| 232 | } | ||
| 233 | |||
| 234 | /*! | ||
| 235 | * \brief prepares all data and variables that are necessary to evaluate the residual of an Darcy element (i.e. Stokes information) | ||
| 236 | */ | ||
| 237 | template<std::size_t i> | ||
| 238 | 1545600 | void bindCouplingContext_(Dune::index_constant<i> domainI, const Element<i>& element) const | |
| 239 | { | ||
| 240 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 772800 times.
✓ Branch 3 taken 386400 times.
✗ Branch 4 not taken.
|
2318400 | const auto fvGeometry = localView(this->problem(domainI).gridGeometry()).bindElement(element); |
| 241 |
1/2✓ Branch 1 taken 386400 times.
✗ Branch 2 not taken.
|
1545600 | bindCouplingContext_(domainI, fvGeometry); |
| 242 | 1545600 | } | |
| 243 | |||
| 244 | /*! | ||
| 245 | * \brief prepares all data and variables that are necessary to evaluate the residual of an Darcy element (i.e. Stokes information) | ||
| 246 | */ | ||
| 247 | template<std::size_t i> | ||
| 248 | 1550960 | void bindCouplingContext_(Dune::index_constant<i> domainI, const FVElementGeometry<i>& fvGeometry) const | |
| 249 | { | ||
| 250 |
2/2✓ Branch 0 taken 10505 times.
✓ Branch 1 taken 764975 times.
|
1550960 | auto& context = std::get<domainI>(couplingContext_); |
| 251 | 1550960 | context.clear(); | |
| 252 | |||
| 253 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 775480 times.
|
1550960 | const auto eIdx = this->problem(domainI).gridGeometry().elementMapper().index(fvGeometry.element()); |
| 254 | |||
| 255 | // do nothing if the element is not coupled to the other domain | ||
| 256 |
2/2✓ Branch 1 taken 10520 times.
✓ Branch 2 taken 764960 times.
|
3101920 | if (!isCoupledElement_(domainI, eIdx)) |
| 257 | return; | ||
| 258 | |||
| 259 | 21040 | couplingContextBoundForElement_[domainI] = eIdx; | |
| 260 | |||
| 261 |
2/2✓ Branch 1 taken 42080 times.
✓ Branch 2 taken 10520 times.
|
105200 | for (const auto& scvf : scvfs(fvGeometry)) |
| 262 | { | ||
| 263 |
2/2✓ Branch 0 taken 10520 times.
✓ Branch 1 taken 31560 times.
|
84160 | if (isCoupled(domainI, scvf)) |
| 264 | { | ||
| 265 |
1/2✓ Branch 1 taken 6600 times.
✗ Branch 2 not taken.
|
21040 | const auto otherElementIdx = couplingMapper_.outsideElementIndex(domainI, scvf); |
| 266 | constexpr auto domainJ = Dune::index_constant<1-domainI>(); | ||
| 267 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 10520 times.
✓ Branch 3 taken 6600 times.
✗ Branch 4 not taken.
|
21040 | const auto& otherGridGeometry = this->problem(domainJ).gridGeometry(); |
| 268 |
2/3✓ Branch 1 taken 6600 times.
✓ Branch 2 taken 3920 times.
✗ Branch 3 not taken.
|
34240 | const auto& otherElement = otherGridGeometry.element(otherElementIdx); |
| 269 |
1/2✓ Branch 1 taken 3920 times.
✗ Branch 2 not taken.
|
28880 | auto otherFvGeometry = localView(otherGridGeometry).bindElement(otherElement); |
| 270 | |||
| 271 | // there is only one scv for TPFA | ||
| 272 |
2/4✓ Branch 1 taken 6600 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6600 times.
✗ Branch 5 not taken.
|
47440 | context.push_back({ |
| 273 | otherElement, | ||
| 274 |
1/2✓ Branch 1 taken 10520 times.
✗ Branch 2 not taken.
|
21040 | volVars_(domainJ, otherElement, *std::begin(scvs(otherFvGeometry))), |
| 275 | std::move(otherFvGeometry), | ||
| 276 | otherElementIdx, | ||
| 277 | 21040 | scvf.index(), | |
| 278 |
2/4✓ Branch 1 taken 10520 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6600 times.
✗ Branch 5 not taken.
|
21040 | couplingMapper_.flipScvfIndex(domainI, scvf), |
| 279 | VelocityVector{} | ||
| 280 | }); | ||
| 281 | 7840 | } | |
| 282 | } | ||
| 283 | 7840 | } | |
| 284 | |||
| 285 | mutable std::tuple<std::vector<FreeFlowMassCouplingContext>, std::vector<PorousMediumCouplingContext>> couplingContext_; | ||
| 286 | mutable std::array<std::size_t, 2> couplingContextBoundForElement_; | ||
| 287 | |||
| 288 | CouplingMapper couplingMapper_; | ||
| 289 | }; | ||
| 290 | |||
| 291 | } // end namespace Dumux | ||
| 292 | |||
| 293 | #endif | ||
| 294 |