| 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-FileCopyrightInfo: 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 FreeFlowPoreNetworkCoupling | ||
| 10 | * \copydoc Dumux::FreeFlowMomentumPoreNetworkCouplingManager | ||
| 11 | */ | ||
| 12 | |||
| 13 | #ifndef DUMUX_MULTIDOMAIN_BOUNDARY_FFPM_FFMMOMENTUMPORENETWORK_COUPLINGMANAGER_HH | ||
| 14 | #define DUMUX_MULTIDOMAIN_BOUNDARY_FFPM_FFMMOMENTUMPORENETWORK_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/discretization/staggered/elementsolution.hh> | ||
| 23 | #include <dumux/multidomain/couplingmanager.hh> | ||
| 24 | #include <dumux/multidomain/boundary/freeflowporenetwork/couplingmapper.hh> | ||
| 25 | |||
| 26 | namespace Dumux { | ||
| 27 | |||
| 28 | /*! | ||
| 29 | * \ingroup FreeFlowPoreNetworkCoupling | ||
| 30 | * \brief Coupling manager for free-flow momentum and pore-network models. | ||
| 31 | */ | ||
| 32 | template<class MDTraits> | ||
| 33 | ✗ | class FreeFlowMomentumPoreNetworkCouplingManager | |
| 34 | : public CouplingManager<MDTraits> | ||
| 35 | { | ||
| 36 | using Scalar = typename MDTraits::Scalar; | ||
| 37 | using ParentType = CouplingManager<MDTraits>; | ||
| 38 | |||
| 39 | public: | ||
| 40 | static constexpr auto freeFlowMomentumIndex = typename MDTraits::template SubDomain<0>::Index(); | ||
| 41 | static constexpr auto poreNetworkIndex = typename MDTraits::template SubDomain<1>::Index(); | ||
| 42 | |||
| 43 | using SolutionVectorStorage = typename ParentType::SolutionVectorStorage; | ||
| 44 | private: | ||
| 45 | // obtain the type tags of the sub problems | ||
| 46 | using FreeFlowMomentumTypeTag = typename MDTraits::template SubDomain<freeFlowMomentumIndex>::TypeTag; | ||
| 47 | using PoreNetworkTypeTag = typename MDTraits::template SubDomain<poreNetworkIndex>::TypeTag; | ||
| 48 | |||
| 49 | using CouplingStencils = std::unordered_map<std::size_t, std::vector<std::size_t> >; | ||
| 50 | using CouplingStencil = CouplingStencils::mapped_type; | ||
| 51 | |||
| 52 | // the sub domain type tags | ||
| 53 | template<std::size_t id> | ||
| 54 | using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag; | ||
| 55 | |||
| 56 | template<std::size_t id> using GridView = typename GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>::GridView; | ||
| 57 | template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>; | ||
| 58 | template<std::size_t id> using ElementVolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::LocalView; | ||
| 59 | template<std::size_t id> using GridVolumeVariables = GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>; | ||
| 60 | template<std::size_t id> using VolumeVariables = typename GetPropType<SubDomainTypeTag<id>, Properties::GridVolumeVariables>::VolumeVariables; | ||
| 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 GridVariables = GetPropType<SubDomainTypeTag<id>, Properties::GridVariables>; | ||
| 64 | template<std::size_t id> using GridFluxVariablesCache = GetPropType<SubDomainTypeTag<id>, Properties::GridFluxVariablesCache>; | ||
| 65 | template<std::size_t id> using ElementFluxVariablesCache = typename GridFluxVariablesCache<id>::LocalView; | ||
| 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<freeFlowMomentumIndex>::Geometry::GlobalCoordinate; | ||
| 72 | |||
| 73 | struct FreeFlowMomentumCouplingContext | ||
| 74 | { | ||
| 75 | FVElementGeometry<poreNetworkIndex> fvGeometry; | ||
| 76 | ElementVolumeVariables<poreNetworkIndex> elemVolVars; | ||
| 77 | ElementFluxVariablesCache<poreNetworkIndex> elemFluxVarsCache; | ||
| 78 | std::size_t poreNetworkDofIdx; | ||
| 79 | }; | ||
| 80 | |||
| 81 | struct PoreNetworkCouplingContext | ||
| 82 | { | ||
| 83 | SubControlVolumeFace<freeFlowMomentumIndex> freeFlowMomentumScvf; | ||
| 84 | VelocityVector faceVelocity; | ||
| 85 | std::size_t freeFlowMomentumDofIdx; | ||
| 86 | }; | ||
| 87 | |||
| 88 | using CouplingMapper = StaggeredFreeFlowPoreNetworkCouplingMapper; | ||
| 89 | using GridVariablesTuple = typename MDTraits::template TupleOfSharedPtr<GridVariables>; | ||
| 90 | |||
| 91 | public: | ||
| 92 | |||
| 93 | /*! | ||
| 94 | * \brief Methods to be accessed by main | ||
| 95 | */ | ||
| 96 | // \{ | ||
| 97 | |||
| 98 | //! Initialize the coupling manager | ||
| 99 | 2 | void init(std::shared_ptr<Problem<freeFlowMomentumIndex>> freeFlowMomentumProblem, | |
| 100 | std::shared_ptr<Problem<poreNetworkIndex>> porousMediumProblem, | ||
| 101 | GridVariablesTuple&& gridVariables, | ||
| 102 | std::shared_ptr<CouplingMapper> couplingMapper, | ||
| 103 | SolutionVectorStorage& curSol) | ||
| 104 | { | ||
| 105 | 2 | couplingMapper_ = couplingMapper; | |
| 106 | 2 | gridVariables_ = gridVariables; | |
| 107 | 4 | this->setSubProblems(std::make_tuple(freeFlowMomentumProblem, porousMediumProblem)); | |
| 108 | 2 | this->attachSolution(curSol); | |
| 109 | 2 | } | |
| 110 | |||
| 111 | // \} | ||
| 112 | |||
| 113 | |||
| 114 | /*! | ||
| 115 | * \brief Methods to be accessed by the assembly | ||
| 116 | */ | ||
| 117 | // \{ | ||
| 118 | |||
| 119 | /*! | ||
| 120 | * \brief prepares all data and variables that are necessary to evaluate the residual (called from the local assembler) | ||
| 121 | */ | ||
| 122 | template<std::size_t i, class Assembler> | ||
| 123 | ✗ | void bindCouplingContext(Dune::index_constant<i> domainI, const Element<i>& element, const Assembler& assembler) const | |
| 124 | { | ||
| 125 |
0/18✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
|
966 | bindCouplingContext_(domainI, element); |
| 126 | ✗ | } | |
| 127 | |||
| 128 | /*! | ||
| 129 | * \brief prepares all data and variables that are necessary to evaluate the residual (called from the local assembler) | ||
| 130 | */ | ||
| 131 | template<std::size_t i> | ||
| 132 | void bindCouplingContext(Dune::index_constant<i> domainI, const Element<i>& element) const | ||
| 133 | { | ||
| 134 | bindCouplingContext_(domainI, element); | ||
| 135 | } | ||
| 136 | |||
| 137 | /*! | ||
| 138 | * \brief Update the coupling context | ||
| 139 | */ | ||
| 140 | template<std::size_t i, std::size_t j, class LocalAssemblerI> | ||
| 141 | 119748 | void updateCouplingContext(Dune::index_constant<i> domainI, | |
| 142 | const LocalAssemblerI& localAssemblerI, | ||
| 143 | Dune::index_constant<j> domainJ, | ||
| 144 | std::size_t dofIdxGlobalJ, | ||
| 145 | const PrimaryVariables<j>& priVarsJ, | ||
| 146 | int pvIdxJ) | ||
| 147 | { | ||
| 148 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 59716 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 59716 times.
|
239496 | this->curSol(domainJ)[dofIdxGlobalJ][pvIdxJ] = priVarsJ[pvIdxJ]; |
| 149 | |||
| 150 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 59716 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 59716 times.
|
240444 | const auto eIdx = localAssemblerI.fvGeometry().gridGeometry().elementMapper().index(localAssemblerI.fvGeometry().element()); |
| 151 | 239496 | if (!isCoupledElement_(domainI, eIdx)) | |
| 152 | return; | ||
| 153 | |||
| 154 | // we need to update all solution-depenent components of the coupling context | ||
| 155 | // the dof of domain J has been deflected | ||
| 156 | |||
| 157 | // update the faceVelocity in the PoreNetworkCouplingContext | ||
| 158 | if constexpr (domainJ == freeFlowMomentumIndex) | ||
| 159 | { | ||
| 160 | // we only need to update if we are assembling the porous medium domain | ||
| 161 | // since the freeflow domain will not use the velocity from the context | ||
| 162 | if constexpr (domainI == poreNetworkIndex) | ||
| 163 | { | ||
| 164 | 132 | auto& context = std::get<poreNetworkIndex>(couplingContext_); | |
| 165 |
4/4✓ Branch 0 taken 258 times.
✓ Branch 1 taken 66 times.
✓ Branch 2 taken 258 times.
✓ Branch 3 taken 66 times.
|
912 | for (auto& c : context) |
| 166 | { | ||
| 167 |
2/2✓ Branch 0 taken 66 times.
✓ Branch 1 taken 192 times.
|
516 | if (c.freeFlowMomentumDofIdx == dofIdxGlobalJ) |
| 168 | { | ||
| 169 |
2/4✓ Branch 0 taken 66 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 66 times.
|
132 | assert(c.freeFlowMomentumScvf.isFrontal() && c.freeFlowMomentumScvf.boundary()); |
| 170 | 132 | c.faceVelocity = faceVelocity(c.freeFlowMomentumScvf, c.freeFlowMomentumDofIdx); | |
| 171 | } | ||
| 172 | } | ||
| 173 | } | ||
| 174 | } | ||
| 175 | |||
| 176 | // update the elemVolVars and elemFluxVarsCache in the FreeFlowMomentumCouplingContext | ||
| 177 | else if constexpr (domainJ == poreNetworkIndex) | ||
| 178 | { | ||
| 179 |
2/3✗ Branch 5 not taken.
✓ Branch 6 taken 588 times.
✓ Branch 7 taken 36 times.
|
6312 | assert(couplingContextBoundForElement_[domainI] == localAssemblerI.fvGeometry().gridGeometry().elementMapper().index(localAssemblerI.fvGeometry().element())); |
| 180 | // there is only one context per coupled free-flow momentum dof | ||
| 181 | 1248 | auto& context = std::get<freeFlowMomentumIndex>(couplingContext_)[0]; | |
| 182 | 1248 | const auto& ggJ = context.fvGeometry.gridGeometry(); | |
| 183 | 1248 | const auto& element = context.fvGeometry.element(); | |
| 184 | 2496 | const auto elemSol = elementSolution(element, this->curSol(domainJ), ggJ); | |
| 185 | |||
| 186 |
2/2✓ Branch 0 taken 1248 times.
✓ Branch 1 taken 624 times.
|
4992 | for (const auto& scv : scvs(context.fvGeometry)) |
| 187 | { | ||
| 188 |
2/2✓ Branch 0 taken 608 times.
✓ Branch 1 taken 640 times.
|
2496 | if (scv.dofIndex() == dofIdxGlobalJ) |
| 189 | { | ||
| 190 | if constexpr (ElementVolumeVariables<poreNetworkIndex>::GridVolumeVariables::cachingEnabled) | ||
| 191 | gridVars_(poreNetworkIndex).curGridVolVars().volVars(scv).update(std::move(elemSol), this->problem(domainJ), element, scv); | ||
| 192 | else | ||
| 193 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 608 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 608 times.
|
2432 | context.elemVolVars[scv].update(std::move(elemSol), this->problem(domainJ), element, scv); |
| 194 | } | ||
| 195 | } | ||
| 196 | |||
| 197 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 624 times.
|
1248 | const auto& scvf = context.fvGeometry.scvf(0); |
| 198 | if constexpr (ElementFluxVariablesCache<poreNetworkIndex>::GridFluxVariablesCache::cachingEnabled) | ||
| 199 | { | ||
| 200 | const auto eIdx = ggJ.elementMapper().index(element); | ||
| 201 | gridVars_(poreNetworkIndex).gridFluxVarsCache().cache(eIdx, scvf.index()).update(this->problem(domainJ), element, context.fvGeometry, context.elemVolVars, scvf); | ||
| 202 | } | ||
| 203 | else | ||
| 204 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 624 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 624 times.
|
2496 | context.elemFluxVarsCache[scvf].update(this->problem(domainJ), element, context.fvGeometry, context.elemVolVars, scvf); |
| 205 | } | ||
| 206 | } | ||
| 207 | |||
| 208 | // \} | ||
| 209 | |||
| 210 | /*! | ||
| 211 | * \brief Access the coupling context needed for the Stokes domain | ||
| 212 | */ | ||
| 213 | ✗ | const auto& couplingContext(const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry, | |
| 214 | const SubControlVolumeFace<freeFlowMomentumIndex> scvf) const | ||
| 215 | { | ||
| 216 |
1/8✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 3606 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
3606 | auto& contexts = std::get<freeFlowMomentumIndex>(couplingContext_); |
| 217 | |||
| 218 |
5/40✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 3606 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 3606 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✓ Branch 15 taken 3606 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 3606 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 3606 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
|
7212 | if (contexts.empty() || couplingContextBoundForElement_[freeFlowMomentumIndex] != fvGeometry.elementIndex()) |
| 219 | ✗ | bindCouplingContext_(freeFlowMomentumIndex, fvGeometry); | |
| 220 | |||
| 221 | |||
| 222 |
0/2✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
3606 | return contexts[0]; |
| 223 | } | ||
| 224 | |||
| 225 | /*! | ||
| 226 | * \brief Access the coupling context needed for the PNM domain | ||
| 227 | */ | ||
| 228 | const auto& couplingContext(const FVElementGeometry<poreNetworkIndex>& fvGeometry, | ||
| 229 | const SubControlVolume<poreNetworkIndex> scv) const | ||
| 230 | { | ||
| 231 | auto& contexts = std::get<poreNetworkIndex>(couplingContext_); | ||
| 232 | |||
| 233 | const auto eIdx = fvGeometry.gridGeometry().elementMapper().index(fvGeometry.element()); | ||
| 234 | if (contexts.empty() || couplingContextBoundForElement_[poreNetworkIndex] != eIdx) | ||
| 235 | bindCouplingContext_(poreNetworkIndex, fvGeometry); | ||
| 236 | |||
| 237 | return contexts; | ||
| 238 | } | ||
| 239 | |||
| 240 | /*! | ||
| 241 | * \brief The coupling stencils | ||
| 242 | */ | ||
| 243 | // \{ | ||
| 244 | |||
| 245 | /*! | ||
| 246 | * \brief The Stokes cell center coupling stencil w.r.t. Darcy DOFs | ||
| 247 | */ | ||
| 248 | ✗ | const CouplingStencil& couplingStencil(Dune::index_constant<poreNetworkIndex> domainI, | |
| 249 | const Element<poreNetworkIndex>& element, | ||
| 250 | Dune::index_constant<freeFlowMomentumIndex> domainJ) const | ||
| 251 | { | ||
| 252 | ✗ | const auto eIdx = this->problem(domainI).gridGeometry().elementMapper().index(element); | |
| 253 | ✗ | return couplingMapper_->poreNetworkToFreeFlowMomentumCouplingStencil(eIdx); | |
| 254 | } | ||
| 255 | |||
| 256 | /*! | ||
| 257 | * \brief returns an iterable container of all indices of degrees of freedom of domain j | ||
| 258 | * that couple with / influence the residual of the given sub-control volume of domain i | ||
| 259 | * | ||
| 260 | * \param domainI the domain index of domain i | ||
| 261 | * \param elementI the coupled element of domain í | ||
| 262 | * \param scvI the sub-control volume of domain i | ||
| 263 | * \param domainJ the domain index of domain j | ||
| 264 | */ | ||
| 265 | ✗ | const CouplingStencil& couplingStencil(Dune::index_constant<freeFlowMomentumIndex> domainI, | |
| 266 | const Element<freeFlowMomentumIndex>& elementI, | ||
| 267 | const SubControlVolume<freeFlowMomentumIndex>& scvI, | ||
| 268 | Dune::index_constant<poreNetworkIndex> domainJ) const | ||
| 269 | { | ||
| 270 |
4/16✓ Branch 1 taken 3864 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3864 times.
✗ Branch 5 not taken.
✓ Branch 9 taken 1632 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 1632 times.
✗ Branch 13 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
|
10992 | return couplingMapper_->freeFlowMomentumToPoreNetworkCouplingStencil(scvI.dofIndex()); |
| 271 | } | ||
| 272 | |||
| 273 | using ParentType::evalCouplingResidual; | ||
| 274 | |||
| 275 | /*! | ||
| 276 | * \brief evaluate the coupling residual | ||
| 277 | * special interface for fcstaggered methods | ||
| 278 | */ | ||
| 279 | template<class LocalAssemblerI, std::size_t j> | ||
| 280 | ✗ | decltype(auto) evalCouplingResidual(Dune::index_constant<freeFlowMomentumIndex> domainI, | |
| 281 | const LocalAssemblerI& localAssemblerI, | ||
| 282 | const SubControlVolume<freeFlowMomentumIndex>& scvI, | ||
| 283 | Dune::index_constant<j> domainJ, | ||
| 284 | std::size_t dofIdxGlobalJ) const | ||
| 285 | { | ||
| 286 |
1/6✓ Branch 1 taken 330 times.
✗ Branch 2 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
660 | return localAssemblerI.evalLocalResidual(); |
| 287 | } | ||
| 288 | |||
| 289 | // \} | ||
| 290 | |||
| 291 | /*! | ||
| 292 | * \brief Returns whether a given scvf is coupled to the other domain | ||
| 293 | */ | ||
| 294 | bool isCoupledLateralScvf(Dune::index_constant<freeFlowMomentumIndex> domainI, const SubControlVolumeFace<freeFlowMomentumIndex>& scvf) const | ||
| 295 | { return couplingMapper_->isCoupledFreeFlowMomentumLateralScvf(scvf.index()); } | ||
| 296 | |||
| 297 | /*! | ||
| 298 | * \brief Returns whether a given scvf is coupled to the other domain | ||
| 299 | */ | ||
| 300 | ✗ | bool isCoupled(Dune::index_constant<freeFlowMomentumIndex> domainI, const SubControlVolumeFace<freeFlowMomentumIndex>& scvf) const | |
| 301 | { | ||
| 302 | ✗ | return couplingMapper_->isCoupledFreeFlowMomentumScvf(scvf.index()) || couplingMapper_->isCoupledFreeFlowMomentumLateralScvf(scvf.index()); | |
| 303 | } | ||
| 304 | |||
| 305 | /*! | ||
| 306 | * \brief If the boundary entity is on a coupling boundary | ||
| 307 | * \param domainI the domain index for which to compute the flux | ||
| 308 | * \param scv the sub control volume | ||
| 309 | */ | ||
| 310 | bool isCoupled(Dune::index_constant<poreNetworkIndex> domainI, | ||
| 311 | const SubControlVolume<poreNetworkIndex>& scv) const | ||
| 312 | { return couplingMapper_->isCoupledPoreNetworkDof(scv.dofIndex()); } | ||
| 313 | |||
| 314 | /*! | ||
| 315 | * \brief Returns the velocity at a given sub control volume face. | ||
| 316 | */ | ||
| 317 | 89 | auto faceVelocity(const SubControlVolumeFace<freeFlowMomentumIndex>& scvf, | |
| 318 | std::size_t freeFlowMomentumDofIdx) const | ||
| 319 | { | ||
| 320 | // create a unit normal vector oriented in positive coordinate direction | ||
| 321 | 89 | auto velocity = scvf.unitOuterNormal(); | |
| 322 | using std::abs; | ||
| 323 | 534 | std::for_each(velocity.begin(), velocity.end(), [](auto& v){ v = abs(v); }); | |
| 324 | |||
| 325 | // create the actual velocity vector | ||
| 326 | 267 | velocity *= this->curSol(freeFlowMomentumIndex)[freeFlowMomentumDofIdx]; | |
| 327 | |||
| 328 | 89 | return velocity; | |
| 329 | } | ||
| 330 | |||
| 331 | private: | ||
| 332 | |||
| 333 | /*! | ||
| 334 | * \brief Returns whether a given element is coupled to the other domain | ||
| 335 | */ | ||
| 336 | template<std::size_t i> | ||
| 337 | ✗ | bool isCoupledElement_(Dune::index_constant<i> domainI, std::size_t eIdx) const | |
| 338 | { | ||
| 339 | if constexpr (i == freeFlowMomentumIndex) | ||
| 340 |
2/2✓ Branch 4 taken 588 times.
✓ Branch 5 taken 72 times.
|
119432 | return couplingMapper_->isCoupledFreeFlowElement(eIdx); |
| 341 | else | ||
| 342 |
3/4✓ Branch 3 taken 36 times.
✓ Branch 4 taken 56 times.
✓ Branch 7 taken 66 times.
✗ Branch 8 not taken.
|
337 | return couplingMapper_->isCoupledPoreNetworkElement(eIdx); |
| 343 | } | ||
| 344 | |||
| 345 | /*! | ||
| 346 | * \brief prepares all data and variables that are necessary to evaluate the residual of an Darcy element (i.e. Stokes information) | ||
| 347 | */ | ||
| 348 | template<std::size_t i> | ||
| 349 | 46 | void bindCouplingContext_(Dune::index_constant<i> domainI, const Element<i>& element) const | |
| 350 | { | ||
| 351 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 23 times.
|
92 | const auto fvGeometry = localView(this->problem(domainI).gridGeometry()).bindElement(element); |
| 352 |
0/2✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
46 | bindCouplingContext_(domainI, fvGeometry); |
| 353 | 46 | } | |
| 354 | |||
| 355 | /*! | ||
| 356 | * \brief prepares all data and variables that are necessary to evaluate the residual of an free-flow momentum element | ||
| 357 | */ | ||
| 358 | 966 | void bindCouplingContext_(Dune::index_constant<freeFlowMomentumIndex> domainI, const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry) const | |
| 359 | { | ||
| 360 |
2/2✓ Branch 0 taken 956 times.
✓ Branch 1 taken 10 times.
|
966 | auto& context = std::get<domainI>(couplingContext_); |
| 361 |
2/2✓ Branch 0 taken 956 times.
✓ Branch 1 taken 10 times.
|
966 | const auto eIdx = fvGeometry.elementIndex(); |
| 362 | |||
| 363 | // do nothing if the element is already correctly bound | ||
| 364 |
6/8✓ Branch 0 taken 956 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 956 times.
✓ Branch 3 taken 10 times.
✓ Branch 4 taken 956 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 956 times.
✗ Branch 7 not taken.
|
1932 | if ((!context.empty() && couplingContextBoundForElement_[domainI] == eIdx)) |
| 365 | 915 | return; | |
| 366 | |||
| 367 | 966 | bool bindElement = false; | |
| 368 | std::size_t actuallyCoupledFreeFlowElementIndex; | ||
| 369 | |||
| 370 | // if the element is directly coupled to a lowDim dof, | ||
| 371 | // bind the element itself | ||
| 372 |
2/2✓ Branch 2 taken 33 times.
✓ Branch 3 taken 933 times.
|
1932 | if (couplingMapper_->isCoupledFreeFlowElement(eIdx)) |
| 373 | { | ||
| 374 | 33 | bindElement = true; | |
| 375 | 33 | actuallyCoupledFreeFlowElementIndex = eIdx; | |
| 376 | } | ||
| 377 | else | ||
| 378 | { | ||
| 379 | // if we assemble another element that is not directly coupled to the lowDim dof | ||
| 380 | // but shares an intersection (and hence, a dof) with a neighbor element that does, bind that neighbor | ||
| 381 |
3/4✗ Branch 0 not taken.
✓ Branch 1 taken 933 times.
✓ Branch 4 taken 933 times.
✓ Branch 5 taken 3732 times.
|
4665 | for (const auto& intersection : intersections(fvGeometry.gridGeometry().gridView(), fvGeometry.element())) |
| 382 | { | ||
| 383 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3732 times.
|
3732 | const auto dofIdx = fvGeometry.gridGeometry().intersectionMapper().globalIntersectionIndex(fvGeometry.element(), intersection.indexInInside()); |
| 384 |
6/6✓ Branch 0 taken 18 times.
✓ Branch 1 taken 3714 times.
✓ Branch 2 taken 18 times.
✓ Branch 3 taken 3714 times.
✓ Branch 4 taken 18 times.
✓ Branch 5 taken 3714 times.
|
11196 | if (couplingMapper_->isCoupledFreeFlowMomentumDof(dofIdx)) |
| 385 | { | ||
| 386 | 18 | bindElement = true; | |
| 387 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 18 times.
|
18 | actuallyCoupledFreeFlowElementIndex = fvGeometry.gridGeometry().elementMapper().index(intersection.outside()); |
| 388 | } | ||
| 389 | } | ||
| 390 | } | ||
| 391 | |||
| 392 | // do nothing if the element is not coupled to the other domain | ||
| 393 |
2/2✓ Branch 0 taken 18 times.
✓ Branch 1 taken 915 times.
|
933 | if (!bindElement) |
| 394 | return; | ||
| 395 | |||
| 396 |
2/2✓ Branch 0 taken 49 times.
✓ Branch 1 taken 2 times.
|
51 | context.clear(); |
| 397 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 51 times.
|
51 | couplingContextBoundForElement_[domainI] = eIdx; |
| 398 | |||
| 399 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 51 times.
|
51 | auto poreNetworkFVGeometry = localView(this->problem(poreNetworkIndex).gridGeometry()); |
| 400 |
2/6✓ Branch 2 taken 51 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 51 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
102 | auto poreNetworkElemVolVars = localView(gridVars_(poreNetworkIndex).curGridVolVars()); |
| 401 |
4/8✓ Branch 1 taken 51 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 51 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 51 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 51 times.
|
102 | auto poreNetworkElemFluxVarsCache = localView(gridVars_(poreNetworkIndex).gridFluxVarsCache()); |
| 402 | |||
| 403 |
3/6✓ Branch 1 taken 51 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 51 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 51 times.
✗ Branch 8 not taken.
|
153 | const auto poreNetworkElemIdx = couplingMapper_->freeFlowElementToPNMElementMap().at(actuallyCoupledFreeFlowElementIndex); |
| 404 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 51 times.
✓ Branch 3 taken 51 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 51 times.
✗ Branch 7 not taken.
|
51 | const auto& poreNetworkElement = this->problem(poreNetworkIndex).gridGeometry().element(poreNetworkElemIdx); |
| 405 | |||
| 406 |
1/2✓ Branch 1 taken 51 times.
✗ Branch 2 not taken.
|
51 | poreNetworkFVGeometry.bindElement(poreNetworkElement); |
| 407 |
2/4✓ Branch 1 taken 51 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 51 times.
✗ Branch 5 not taken.
|
102 | poreNetworkElemVolVars.bind(poreNetworkElement, poreNetworkFVGeometry, this->curSol(poreNetworkIndex)); |
| 408 |
1/2✓ Branch 1 taken 51 times.
✗ Branch 2 not taken.
|
51 | poreNetworkElemFluxVarsCache.bind(poreNetworkElement, poreNetworkFVGeometry, poreNetworkElemVolVars); |
| 409 | |||
| 410 | 51 | const std::size_t poreNetworkDofIdx = [&] | |
| 411 | { | ||
| 412 | 51 | std::size_t idx = 0; | |
| 413 | 51 | std::size_t counter = 0; | |
| 414 |
2/2✓ Branch 0 taken 102 times.
✓ Branch 1 taken 51 times.
|
204 | for (const auto& scv : scvs(poreNetworkFVGeometry)) |
| 415 | { | ||
| 416 |
4/4✓ Branch 0 taken 51 times.
✓ Branch 1 taken 51 times.
✓ Branch 2 taken 51 times.
✓ Branch 3 taken 51 times.
|
204 | if (couplingMapper_->isCoupledPoreNetworkDof(scv.dofIndex())) |
| 417 | { | ||
| 418 | 51 | idx = scv.dofIndex(); | |
| 419 | 51 | ++counter; | |
| 420 | } | ||
| 421 | } | ||
| 422 | |||
| 423 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 51 times.
|
51 | if (counter != 1) |
| 424 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Exactly one pore per throat needs to be coupled with the FF domain"); | |
| 425 | |||
| 426 | 51 | return idx; | |
| 427 |
3/4✓ Branch 1 taken 51 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 51 times.
✓ Branch 4 taken 51 times.
|
153 | }(); |
| 428 | |||
| 429 |
3/8✓ Branch 0 taken 51 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 51 times.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 51 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
102 | context.push_back({std::move(poreNetworkFVGeometry), |
| 430 |
1/2✓ Branch 1 taken 51 times.
✗ Branch 2 not taken.
|
51 | std::move(poreNetworkElemVolVars), |
| 431 |
1/2✓ Branch 1 taken 51 times.
✗ Branch 2 not taken.
|
51 | std::move(poreNetworkElemFluxVarsCache), |
| 432 | poreNetworkDofIdx} | ||
| 433 | ); | ||
| 434 | } | ||
| 435 | |||
| 436 | /*! | ||
| 437 | * \brief prepares all data and variables that are necessary to evaluate the residual of an pore-network model element | ||
| 438 | */ | ||
| 439 | 23 | void bindCouplingContext_(Dune::index_constant<poreNetworkIndex> domainI, const FVElementGeometry<poreNetworkIndex>& fvGeometry) const | |
| 440 | { | ||
| 441 | 23 | auto& context = std::get<domainI>(couplingContext_); | |
| 442 | 69 | const auto eIdx = fvGeometry.gridGeometry().elementMapper().index(fvGeometry.element()); | |
| 443 | |||
| 444 | // do nothing if the element is already bound or not coupled to the other domain | ||
| 445 |
10/10✓ Branch 0 taken 14 times.
✓ Branch 1 taken 9 times.
✓ Branch 2 taken 14 times.
✓ Branch 3 taken 9 times.
✓ Branch 4 taken 12 times.
✓ Branch 5 taken 2 times.
✓ Branch 6 taken 12 times.
✓ Branch 7 taken 2 times.
✓ Branch 9 taken 14 times.
✓ Branch 10 taken 7 times.
|
46 | if ((!context.empty() && couplingContextBoundForElement_[domainI] == eIdx) || !isCoupledElement_(domainI, eIdx)) |
| 446 | 16 | return; | |
| 447 | |||
| 448 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 2 times.
|
7 | context.clear(); |
| 449 | 7 | couplingContextBoundForElement_[domainI] = eIdx; | |
| 450 | |||
| 451 | 14 | const auto& stencil = couplingStencil(poreNetworkIndex, fvGeometry.element(), freeFlowMomentumIndex); | |
| 452 | 21 | const auto& freeFlowElements = couplingMapper_->pnmElementToFreeFlowElementsMap().at(eIdx); | |
| 453 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 7 times.
|
7 | auto ffFVGeometry = localView(this->problem(freeFlowMomentumIndex).gridGeometry()); |
| 454 | |||
| 455 |
4/4✓ Branch 0 taken 23 times.
✓ Branch 1 taken 7 times.
✓ Branch 2 taken 23 times.
✓ Branch 3 taken 7 times.
|
44 | for (const auto ffElementIdx : freeFlowElements) |
| 456 | { | ||
| 457 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 23 times.
✓ Branch 3 taken 23 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 23 times.
✗ Branch 7 not taken.
|
23 | const auto& ffElement = this->problem(freeFlowMomentumIndex).gridGeometry().element(ffElementIdx); |
| 458 | 23 | ffFVGeometry.bindElement(ffElement); | |
| 459 |
4/4✓ Branch 1 taken 92 times.
✓ Branch 2 taken 23 times.
✓ Branch 3 taken 92 times.
✓ Branch 4 taken 23 times.
|
115 | for (const auto& scv : scvs(ffFVGeometry)) |
| 460 | { | ||
| 461 |
6/6✓ Branch 0 taken 69 times.
✓ Branch 1 taken 23 times.
✓ Branch 2 taken 69 times.
✓ Branch 3 taken 23 times.
✓ Branch 4 taken 69 times.
✓ Branch 5 taken 23 times.
|
276 | if (couplingMapper_->isCoupledFreeFlowMomentumDof(scv.dofIndex())) |
| 462 | { | ||
| 463 |
16/16✓ Branch 0 taken 1 times.
✓ Branch 1 taken 14 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 13 times.
✓ Branch 4 taken 1 times.
✓ Branch 5 taken 12 times.
✓ Branch 6 taken 1 times.
✓ Branch 7 taken 11 times.
✓ Branch 8 taken 6 times.
✓ Branch 9 taken 48 times.
✓ Branch 10 taken 6 times.
✓ Branch 11 taken 42 times.
✓ Branch 12 taken 7 times.
✓ Branch 13 taken 46 times.
✓ Branch 17 taken 23 times.
✓ Branch 18 taken 46 times.
|
377 | if (std::any_of(stencil.begin(), stencil.end(), [&](const auto x){ return scv.dofIndex() == x; } )) |
| 464 | { | ||
| 465 | 23 | const auto& coupledScvf = ffFVGeometry.frontalScvfOnBoundary(scv); | |
| 466 |
1/2✓ Branch 2 taken 23 times.
✗ Branch 3 not taken.
|
23 | context.push_back({coupledScvf, |
| 467 | faceVelocity(coupledScvf, scv.dofIndex()), | ||
| 468 | scv.dofIndex()} | ||
| 469 | ); | ||
| 470 | } | ||
| 471 | } | ||
| 472 | } | ||
| 473 | } | ||
| 474 | } | ||
| 475 | |||
| 476 | /*! | ||
| 477 | * \brief Return a reference to the grid variables of a sub problem | ||
| 478 | * \param domainIdx The domain index | ||
| 479 | */ | ||
| 480 | template<std::size_t i> | ||
| 481 | 102 | const GridVariables<i>& gridVars_(Dune::index_constant<i> domainIdx) const | |
| 482 | { | ||
| 483 |
2/4✓ Branch 0 taken 102 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 102 times.
✗ Branch 3 not taken.
|
204 | if (std::get<i>(gridVariables_)) |
| 484 | 306 | return *std::get<i>(gridVariables_); | |
| 485 | else | ||
| 486 | ✗ | DUNE_THROW(Dune::InvalidStateException, "The gridVariables pointer was not set. Use setGridVariables() before calling this function"); | |
| 487 | } | ||
| 488 | |||
| 489 | /*! | ||
| 490 | * \brief Return a reference to the grid variables of a sub problem | ||
| 491 | * \param domainIdx The domain index | ||
| 492 | */ | ||
| 493 | template<std::size_t i> | ||
| 494 | GridVariables<i>& gridVars_(Dune::index_constant<i> domainIdx) | ||
| 495 | { | ||
| 496 | if (std::get<i>(gridVariables_)) | ||
| 497 | return *std::get<i>(gridVariables_); | ||
| 498 | else | ||
| 499 | DUNE_THROW(Dune::InvalidStateException, "The gridVariables pointer was not set. Use setGridVariables() before calling this function"); | ||
| 500 | } | ||
| 501 | |||
| 502 | //! A tuple of std::shared_ptrs to the grid variables of the sub problems | ||
| 503 | GridVariablesTuple gridVariables_; | ||
| 504 | |||
| 505 | mutable std::tuple<std::vector<FreeFlowMomentumCouplingContext>, std::vector<PoreNetworkCouplingContext>> couplingContext_; | ||
| 506 | mutable std::array<std::size_t, 2> couplingContextBoundForElement_; | ||
| 507 | |||
| 508 | std::shared_ptr<CouplingMapper> couplingMapper_; | ||
| 509 | }; | ||
| 510 | |||
| 511 | } // end namespace Dumux | ||
| 512 | |||
| 513 | #endif | ||
| 514 |