| 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 PoroElastic | ||
| 10 | * \brief Coupling manager for porous medium flow problems coupled to a poro-mechanical problem | ||
| 11 | */ | ||
| 12 | |||
| 13 | #ifndef DUMUX_POROMECHANICS_COUPLING_MANAGER_HH | ||
| 14 | #define DUMUX_POROMECHANICS_COUPLING_MANAGER_HH | ||
| 15 | |||
| 16 | #include <algorithm> | ||
| 17 | #include <type_traits> | ||
| 18 | |||
| 19 | #include <dune/common/std/type_traits.hh> | ||
| 20 | #include <dumux/discretization/method.hh> | ||
| 21 | #include <dumux/common/properties.hh> | ||
| 22 | #include <dumux/multidomain/couplingmanager.hh> | ||
| 23 | |||
| 24 | namespace Dumux { | ||
| 25 | |||
| 26 | /*! | ||
| 27 | * \ingroup PoroElastic | ||
| 28 | * \brief Coupling manager for porous medium flow problems coupled to a poro-mechanical problem | ||
| 29 | * | ||
| 30 | * Coupling manager for porous medium flow problems coupled to a poro-mechanical | ||
| 31 | * problem both defined on the same grid. Coupling occurs via the change of porosity | ||
| 32 | * and permeability due to mechanical deformations and the influence of the pore | ||
| 33 | * pressure on the effective stresses acting on the porous medium. | ||
| 34 | * | ||
| 35 | * \tparam PMFlowId The porous medium flow domain id | ||
| 36 | * \tparam PoroMechId The poro-mechanical domain id | ||
| 37 | */ | ||
| 38 | template< class MDTraits, | ||
| 39 | std::size_t PMFlowId = 0, | ||
| 40 | std::size_t PoroMechId = PMFlowId+1 > | ||
| 41 | ✗ | class PoroMechanicsCouplingManager : public virtual CouplingManager< MDTraits > | |
| 42 | { | ||
| 43 | using ParentType = CouplingManager< MDTraits >; | ||
| 44 | |||
| 45 | // the sub-domain type tags | ||
| 46 | template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag; | ||
| 47 | |||
| 48 | // further types specific to the sub-problems | ||
| 49 | template<std::size_t id> using Scalar = GetPropType<SubDomainTypeTag<id>, Properties::Scalar>; | ||
| 50 | template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>; | ||
| 51 | template<std::size_t id> using GridVariables = GetPropType<SubDomainTypeTag<id>, Properties::GridVariables>; | ||
| 52 | template<std::size_t id> using PrimaryVariables = typename GridVariables<id>::PrimaryVariables; | ||
| 53 | template<std::size_t id> using GridVolumeVariables = typename GridVariables<id>::GridVolumeVariables; | ||
| 54 | template<std::size_t id> using ElementVolumeVariables = typename GridVolumeVariables<id>::LocalView; | ||
| 55 | template<std::size_t id> using VolumeVariables = typename GridVolumeVariables<id>::VolumeVariables; | ||
| 56 | template<std::size_t id> using GridGeometry = typename GridVariables<id>::GridGeometry; | ||
| 57 | template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView; | ||
| 58 | template<std::size_t id> using GridView = typename GridGeometry<id>::GridView; | ||
| 59 | template<std::size_t id> using GridIndexType = typename GridView<id>::IndexSet::IndexType; | ||
| 60 | template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity; | ||
| 61 | template<std::size_t id> using GlobalPosition = typename Element<id>::Geometry::GlobalCoordinate; | ||
| 62 | |||
| 63 | //! we assume that the two sub-problem operate on the same grid | ||
| 64 | static_assert(std::is_same< GridView<PMFlowId>, GridView<PoroMechId> >::value, | ||
| 65 | "The grid types of the two sub-problems have to be equal!"); | ||
| 66 | |||
| 67 | //! this coupling manager is for cc - box only | ||
| 68 | static_assert(GridGeometry<PoroMechId>::discMethod == DiscretizationMethods::box, | ||
| 69 | "Poro-mechanical problem must be discretized with the box scheme for this coupling manager!"); | ||
| 70 | |||
| 71 | static_assert(GridGeometry<PMFlowId>::discMethod == DiscretizationMethods::cctpfa || | ||
| 72 | GridGeometry<PMFlowId>::discMethod == DiscretizationMethods::ccmpfa, | ||
| 73 | "Porous medium flow problem must be discretized with a cell-centered scheme for this coupling manager!"); | ||
| 74 | |||
| 75 | //! this does not work for enabled grid volume variables caching (update of local view in context has no effect) | ||
| 76 | static_assert(!GetPropType<SubDomainTypeTag<PMFlowId>, Properties::GridVariables>::GridVolumeVariables::cachingEnabled, | ||
| 77 | "Poromechanics framework does not yet work for enabled grid volume variables caching"); | ||
| 78 | |||
| 79 | //! Types used for coupling stencils | ||
| 80 | template<std::size_t id> | ||
| 81 | using CouplingIndexType = typename std::conditional< id == PMFlowId, | ||
| 82 | GridIndexType<PoroMechId>, | ||
| 83 | GridIndexType<PMFlowId> >::type; | ||
| 84 | |||
| 85 | /*! | ||
| 86 | * \brief Porous medium flow domain data required for the residual calculation of an | ||
| 87 | * element of the poro-mechanidal domain. We store the data required to do an | ||
| 88 | * update of all primary/secondary variables at the dof of the element. | ||
| 89 | */ | ||
| 90 | ✗ | struct PoroMechanicsCouplingContext | |
| 91 | { | ||
| 92 | // We need unique ptrs because the local views have no default constructor | ||
| 93 | Element<PMFlowId> pmFlowElement; | ||
| 94 | std::unique_ptr< FVElementGeometry<PMFlowId> > pmFlowFvGeometry; | ||
| 95 | std::unique_ptr< ElementVolumeVariables<PMFlowId> > pmFlowElemVolVars; | ||
| 96 | }; | ||
| 97 | |||
| 98 | template<typename T> | ||
| 99 | using Storage = decltype(std::declval<T>().localResidual().evalStorage( | ||
| 100 | std::declval<T>().fvGeometry(), | ||
| 101 | std::declval<T>().curElemVolVars() | ||
| 102 | )); | ||
| 103 | |||
| 104 | template<typename LA> | ||
| 105 | static constexpr bool hasExperimentalEvalStorage = Dune::Std::is_detected_v<Storage, LA>; | ||
| 106 | |||
| 107 | public: | ||
| 108 | |||
| 109 | // export the domain ids | ||
| 110 | static constexpr auto pmFlowId = Dune::index_constant<PMFlowId>(); | ||
| 111 | static constexpr auto poroMechId = Dune::index_constant<PoroMechId>(); | ||
| 112 | |||
| 113 | /*! | ||
| 114 | * \brief The types used for coupling stencils. An element of the poro-mechanical | ||
| 115 | * domain always only couples to the single dof (because we use a cell-centered | ||
| 116 | * scheme in the porous medium flow domain) of this same element. | ||
| 117 | */ | ||
| 118 | template<std::size_t i, std::size_t j = (i == PMFlowId) ? PoroMechId : PMFlowId> | ||
| 119 | using CouplingStencilType = typename std::conditional< i == PMFlowId, | ||
| 120 | std::vector< CouplingIndexType<i> >, | ||
| 121 | std::array< CouplingIndexType<i>, 1> >::type; | ||
| 122 | |||
| 123 | //! the type of the solution vector | ||
| 124 | using SolutionVector = typename MDTraits::SolutionVector; | ||
| 125 | |||
| 126 | /*! | ||
| 127 | * \brief Initialize the coupling manager. | ||
| 128 | * | ||
| 129 | * \param pmFlowProblem The porous medium flow problem | ||
| 130 | * \param poroMechanicalProblem The poro-mechanical problem | ||
| 131 | * \param curSol The current solution | ||
| 132 | */ | ||
| 133 | 3 | void init(std::shared_ptr< Problem<PMFlowId> > pmFlowProblem, | |
| 134 | std::shared_ptr< Problem<PoroMechId> > poroMechanicalProblem, | ||
| 135 | const SolutionVector& curSol) | ||
| 136 | { | ||
| 137 | // set the sub problems | ||
| 138 |
2/4✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
|
3 | this->setSubProblem(pmFlowProblem, pmFlowId); |
| 139 |
2/4✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
|
3 | this->setSubProblem(poroMechanicalProblem, poroMechId); |
| 140 | |||
| 141 | // copy the solution vector | ||
| 142 | 3 | ParentType::updateSolution(curSol); | |
| 143 | // set up the coupling map pmfow -> poromechanics | ||
| 144 | 3 | initializeCouplingMap_(); | |
| 145 | 3 | } | |
| 146 | |||
| 147 | /*! | ||
| 148 | * \brief Return the coupling stencil for a given porous medium flow domain element | ||
| 149 | */ | ||
| 150 | ✗ | const CouplingStencilType<PMFlowId>& couplingStencil(Dune::index_constant<PMFlowId> pmFlowDomainId, | |
| 151 | const Element<PMFlowId>& element, | ||
| 152 | Dune::index_constant<PoroMechId> poroMechDomainId) const | ||
| 153 | { | ||
| 154 | ✗ | return pmFlowCouplingMap_[ this->problem(pmFlowId).gridGeometry().elementMapper().index(element) ]; | |
| 155 | } | ||
| 156 | |||
| 157 | /*! | ||
| 158 | * \brief Return the coupling element stencil for a given poro-mechanical domain element | ||
| 159 | */ | ||
| 160 | ✗ | const CouplingStencilType<PoroMechId> couplingStencil(Dune::index_constant<PoroMechId> poroMechDomainId, | |
| 161 | const Element<PoroMechId>& element, | ||
| 162 | Dune::index_constant<PMFlowId> pmFlowDomainId) const | ||
| 163 | { | ||
| 164 | ✗ | const auto eIdx = this->problem(pmFlowId).gridGeometry().elementMapper().index(element); | |
| 165 | ✗ | return CouplingStencilType<PoroMechId>{ {eIdx} }; | |
| 166 | } | ||
| 167 | |||
| 168 | //! Pull up the base class' default implementation | ||
| 169 | using ParentType::bindCouplingContext; | ||
| 170 | |||
| 171 | /*! | ||
| 172 | * \brief For the assembly of the element residual of an element of the poro-mechanics domain, | ||
| 173 | * we have to prepare the element variables of the porous medium flow domain. | ||
| 174 | */ | ||
| 175 | template< class Assembler > | ||
| 176 | 5036 | void bindCouplingContext(Dune::index_constant<PoroMechId> poroMechDomainId, | |
| 177 | const Element<PoroMechId>& element, | ||
| 178 | const Assembler& assembler) | ||
| 179 | { | ||
| 180 | // first reset the context | ||
| 181 |
2/2✓ Branch 0 taken 5033 times.
✓ Branch 1 taken 3 times.
|
5036 | poroMechCouplingContext_.pmFlowFvGeometry.reset(nullptr); |
| 182 |
2/2✓ Branch 0 taken 5033 times.
✓ Branch 1 taken 3 times.
|
5036 | poroMechCouplingContext_.pmFlowElemVolVars.reset(nullptr); |
| 183 | |||
| 184 | // prepare the fvGeometry and the element volume variables | ||
| 185 | // these quantities will be used later to obtain the effective pressure | ||
| 186 |
4/8✗ Branch 0 not taken.
✓ Branch 1 taken 5036 times.
✓ Branch 3 taken 5036 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 5036 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 5036 times.
✗ Branch 10 not taken.
|
15108 | const auto fvGeometry = localView( this->problem(pmFlowId).gridGeometry() ).bindElement(element); |
| 187 |
3/6✓ Branch 1 taken 5036 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5036 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 5036 times.
✗ Branch 8 not taken.
|
15108 | const auto elemVolVars = localView(assembler.gridVariables(Dune::index_constant<PMFlowId>()).curGridVolVars()).bindElement(element, |
| 188 | fvGeometry, | ||
| 189 | this->curSol(Dune::index_constant<PMFlowId>())); | ||
| 190 | |||
| 191 |
1/2✓ Branch 1 taken 5036 times.
✗ Branch 2 not taken.
|
5036 | poroMechCouplingContext_.pmFlowElement = element; |
| 192 |
3/6✓ Branch 1 taken 5036 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 5036 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 5036 times.
|
5036 | poroMechCouplingContext_.pmFlowFvGeometry = std::make_unique< FVElementGeometry<PMFlowId> >(fvGeometry); |
| 193 |
3/6✓ Branch 1 taken 5036 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 5036 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 5036 times.
|
5036 | poroMechCouplingContext_.pmFlowElemVolVars = std::make_unique< ElementVolumeVariables<PMFlowId> >(elemVolVars); |
| 194 | 5036 | } | |
| 195 | |||
| 196 | /*! | ||
| 197 | * \brief After deflection of the solution in the porous medium flow domain during | ||
| 198 | * element residual assembly in the poromechanics domain, we have to update | ||
| 199 | * the porous medium flow element variables of the coupling context | ||
| 200 | */ | ||
| 201 | template< class PoroMechLocalAssembler > | ||
| 202 | ✗ | void updateCouplingContext(Dune::index_constant<PoroMechId> poroMechDomainId, | |
| 203 | const PoroMechLocalAssembler& poroMechLocalAssembler, | ||
| 204 | Dune::index_constant<PMFlowId> pmFlowDomainId, | ||
| 205 | GridIndexType<PMFlowId> dofIdxGlobalJ, | ||
| 206 | const PrimaryVariables<PMFlowId>& priVarsJ, | ||
| 207 | unsigned int pvIdxJ) | ||
| 208 | { | ||
| 209 | // communicate the deflected pm flow domain primary variable | ||
| 210 | ✗ | ParentType::updateCouplingContext(poroMechDomainId, poroMechLocalAssembler, pmFlowDomainId, dofIdxGlobalJ, priVarsJ, pvIdxJ); | |
| 211 | |||
| 212 | // now, update the coupling context (i.e. elemVolVars) | ||
| 213 | ✗ | const auto& element = poroMechCouplingContext_.pmFlowElement; | |
| 214 | ✗ | const auto& fvGeometry = *poroMechCouplingContext_.pmFlowFvGeometry; | |
| 215 | ✗ | poroMechCouplingContext_.pmFlowElemVolVars->bindElement(element, fvGeometry, this->curSol(pmFlowDomainId)); | |
| 216 | ✗ | } | |
| 217 | |||
| 218 | /*! | ||
| 219 | * \brief After deflection of the solution in the poromechanics domain during | ||
| 220 | * element residual assembly in that same domain, we have to update the | ||
| 221 | * porous medium flow element variables of the coupling context because | ||
| 222 | * the porosity/permeability might depend on the mechanical deformation | ||
| 223 | */ | ||
| 224 | template< class PoroMechLocalAssembler > | ||
| 225 | ✗ | void updateCouplingContext(Dune::index_constant<PoroMechId> poroMechDomainIdI, | |
| 226 | const PoroMechLocalAssembler& poroMechLocalAssembler, | ||
| 227 | Dune::index_constant<PoroMechId> poroMechDomainIdJ, | ||
| 228 | GridIndexType<PoroMechId> dofIdxGlobalJ, | ||
| 229 | const PrimaryVariables<PoroMechId>& priVarsJ, | ||
| 230 | unsigned int pvIdxJ) | ||
| 231 | { | ||
| 232 | // communicate the deflected displacement | ||
| 233 | ✗ | ParentType::updateCouplingContext(poroMechDomainIdI, poroMechLocalAssembler, poroMechDomainIdJ, dofIdxGlobalJ, priVarsJ, pvIdxJ); | |
| 234 | |||
| 235 | // now, update the coupling context (i.e. elemVolVars) | ||
| 236 | ✗ | (*poroMechCouplingContext_.pmFlowElemVolVars).bindElement(poroMechCouplingContext_.pmFlowElement, | |
| 237 | ✗ | *poroMechCouplingContext_.pmFlowFvGeometry, | |
| 238 | this->curSol(Dune::index_constant<PMFlowId>())); | ||
| 239 | ✗ | } | |
| 240 | |||
| 241 | /*! | ||
| 242 | * \brief We need this overload to avoid ambiguity. However, for the porous medium flow | ||
| 243 | * domain weonly have to update the solution, which is done in the parent class. | ||
| 244 | */ | ||
| 245 | template< std::size_t j, class PMFlowLocalAssembler > | ||
| 246 | ✗ | void updateCouplingContext(Dune::index_constant<PMFlowId> pmFlowDomainId, | |
| 247 | const PMFlowLocalAssembler& pmFlowLocalAssembler, | ||
| 248 | Dune::index_constant<j> domainIdJ, | ||
| 249 | GridIndexType<j> dofIdxGlobalJ, | ||
| 250 | const PrimaryVariables<j>& priVarsJ, | ||
| 251 | unsigned int pvIdxJ) | ||
| 252 | { | ||
| 253 | // communicate the deflected displacement | ||
| 254 | 1134000 | ParentType::updateCouplingContext(pmFlowDomainId, pmFlowLocalAssembler, domainIdJ, dofIdxGlobalJ, priVarsJ, pvIdxJ); | |
| 255 | ✗ | } | |
| 256 | |||
| 257 | //! Pull up the base class' default implementation | ||
| 258 | using ParentType::updateCoupledVariables; | ||
| 259 | |||
| 260 | /*! | ||
| 261 | * \brief Update the porous medium flow domain volume variables and flux variables cache | ||
| 262 | * after the coupling context has been updated. This has to be done because the | ||
| 263 | * mechanical deformation enters the porosity/permeability relationships. | ||
| 264 | */ | ||
| 265 | template< class PMFlowLocalAssembler, class UpdatableFluxVarCache > | ||
| 266 | ✗ | void updateCoupledVariables(Dune::index_constant<PMFlowId> pmFlowDomainId, | |
| 267 | const PMFlowLocalAssembler& pmFlowLocalAssembler, | ||
| 268 | ElementVolumeVariables<PMFlowId>& elemVolVars, | ||
| 269 | UpdatableFluxVarCache& elemFluxVarsCache) | ||
| 270 | { | ||
| 271 | // update the element volume variables to obtain the updated porosity/permeability | ||
| 272 | ✗ | elemVolVars.bind(pmFlowLocalAssembler.element(), | |
| 273 | pmFlowLocalAssembler.fvGeometry(), | ||
| 274 | this->curSol(pmFlowDomainId)); | ||
| 275 | |||
| 276 | // update the transmissibilities subject to the new permeabilities | ||
| 277 | ✗ | elemFluxVarsCache.update(pmFlowLocalAssembler.element(), | |
| 278 | pmFlowLocalAssembler.fvGeometry(), | ||
| 279 | elemVolVars); | ||
| 280 | ✗ | } | |
| 281 | |||
| 282 | /*! | ||
| 283 | * \brief Update the poro-mechanics volume variables after the coupling context has been updated. | ||
| 284 | * This is necessary because the fluid density is stored in them and which potentially is | ||
| 285 | * solution-dependent. | ||
| 286 | */ | ||
| 287 | template< class PoroMechLocalAssembler, class UpdatableFluxVarCache > | ||
| 288 | ✗ | void updateCoupledVariables(Dune::index_constant<PoroMechId> poroMechDomainId, | |
| 289 | const PoroMechLocalAssembler& poroMechLocalAssembler, | ||
| 290 | ElementVolumeVariables<PoroMechId>& elemVolVars, | ||
| 291 | UpdatableFluxVarCache& elemFluxVarsCache) | ||
| 292 | { | ||
| 293 |
1/2✓ Branch 2 taken 3756 times.
✗ Branch 3 not taken.
|
10968 | elemVolVars.bind(poroMechLocalAssembler.element(), |
| 294 | poroMechLocalAssembler.fvGeometry(), | ||
| 295 | this->curSol(poroMechDomainId)); | ||
| 296 | ✗ | } | |
| 297 | |||
| 298 | /*! | ||
| 299 | * \brief Evaluates the coupling element residual of the porous medium flow domain with | ||
| 300 | * respect to the poro-mechanical domain. The deformation might has an effect on | ||
| 301 | * both the permeability as well as the porosity. Thus, we need to compute fluxes | ||
| 302 | * and the storage term here. | ||
| 303 | */ | ||
| 304 | template< class PMFlowLocalAssembler > | ||
| 305 | 369504 | auto evalCouplingResidual(Dune::index_constant<PMFlowId> pmFlowDomainId, | |
| 306 | const PMFlowLocalAssembler& pmFlowLocalAssembler, | ||
| 307 | Dune::index_constant<PoroMechId> poroMechDomainId, | ||
| 308 | GridIndexType<PoroMechId> dofIdxGlobalJ) | ||
| 309 | { | ||
| 310 | 2206944 | auto res = pmFlowLocalAssembler.localResidual().evalFluxAndSource(pmFlowLocalAssembler.element(), | |
| 311 | pmFlowLocalAssembler.fvGeometry(), | ||
| 312 | pmFlowLocalAssembler.curElemVolVars(), | ||
| 313 | pmFlowLocalAssembler.elemFluxVarsCache(), | ||
| 314 | pmFlowLocalAssembler.elemBcTypes()); | ||
| 315 | |||
| 316 | // If the residual instationary, evaluate storage | ||
| 317 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 10080 times.
|
369504 | if (!pmFlowLocalAssembler.localResidual().isStationary()) |
| 318 | { | ||
| 319 | if constexpr (hasExperimentalEvalStorage<PMFlowLocalAssembler>) | ||
| 320 | 359424 | res += pmFlowLocalAssembler.localResidual().evalStorage(pmFlowLocalAssembler.fvGeometry(), | |
| 321 | pmFlowLocalAssembler.curElemVolVars()); | ||
| 322 | else | ||
| 323 | ✗ | res += pmFlowLocalAssembler.localResidual().evalStorage(pmFlowLocalAssembler.element(), | |
| 324 | pmFlowLocalAssembler.fvGeometry(), | ||
| 325 | pmFlowLocalAssembler.prevElemVolVars(), | ||
| 326 | pmFlowLocalAssembler.curElemVolVars()); | ||
| 327 | } | ||
| 328 | |||
| 329 | 369504 | return res; | |
| 330 | } | ||
| 331 | |||
| 332 | /*! | ||
| 333 | * \brief Evaluates the coupling element residual of the poromechanical domain with | ||
| 334 | * respect to the porous medium flow domain. The pressure has an effect on the | ||
| 335 | * mechanical stresses as well as the body forces. Thus, we have to compute | ||
| 336 | * the fluxes as well as the source term here. | ||
| 337 | */ | ||
| 338 | template< class PoroMechLocalAssembler > | ||
| 339 | ✗ | auto evalCouplingResidual(Dune::index_constant<PoroMechId> poroMechDomainId, | |
| 340 | const PoroMechLocalAssembler& poroMechLocalAssembler, | ||
| 341 | Dune::index_constant<PMFlowId> pmFlowDomainId, | ||
| 342 | GridIndexType<PMFlowId> dofIdxGlobalJ) | ||
| 343 | { | ||
| 344 |
1/2✓ Branch 1 taken 3456 times.
✗ Branch 2 not taken.
|
10968 | return poroMechLocalAssembler.localResidual().evalFluxAndSource(poroMechLocalAssembler.element(), |
| 345 | poroMechLocalAssembler.fvGeometry(), | ||
| 346 | poroMechLocalAssembler.curElemVolVars(), | ||
| 347 | poroMechLocalAssembler.elemFluxVarsCache(), | ||
| 348 |
3/6✓ Branch 1 taken 3756 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3756 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3756 times.
✗ Branch 8 not taken.
|
47328 | poroMechLocalAssembler.elemBcTypes()); |
| 349 | } | ||
| 350 | |||
| 351 | //! Return the porous medium flow variables an element/scv of the poromech domain couples to | ||
| 352 | 1538960 | const VolumeVariables<PMFlowId>& getPMFlowVolVars(const Element<PoroMechId>& element) const | |
| 353 | { | ||
| 354 | //! If we do not yet have the queried object, build it first | ||
| 355 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1538960 times.
|
1538960 | const auto eIdx = this->problem(poroMechId).gridGeometry().elementMapper().index(element); |
| 356 | 3077920 | return (*poroMechCouplingContext_.pmFlowElemVolVars)[eIdx]; | |
| 357 | } | ||
| 358 | |||
| 359 | /*! | ||
| 360 | * \brief the solution vector of the subproblem | ||
| 361 | * \param domainIdx The domain index | ||
| 362 | * \note in case of numeric differentiation the solution vector always carries the deflected solution | ||
| 363 | */ | ||
| 364 | template<std::size_t i> | ||
| 365 | ✗ | const auto& curSol(Dune::index_constant<i> domainIdx) const | |
| 366 |
6/11✓ Branch 9 taken 3756 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 3756 times.
✗ Branch 13 not taken.
✓ Branch 19 taken 4736 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 300 times.
✓ Branch 22 taken 4736 times.
✗ Branch 23 not taken.
✓ Branch 24 taken 300 times.
✗ Branch 25 not taken.
|
10784368 | { return ParentType::curSol(domainIdx); } |
| 367 | |||
| 368 | |||
| 369 | private: | ||
| 370 | /*! | ||
| 371 | * \brief Initializes the pm flow domain coupling map. Since the elements | ||
| 372 | * of the poro-mechanical domain only couple to the same elements, we | ||
| 373 | * don't set up the maps here but return copy of the "stencil" always. | ||
| 374 | */ | ||
| 375 | 3 | void initializeCouplingMap_() | |
| 376 | { | ||
| 377 | // some references for convenience | ||
| 378 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
|
3 | const auto& pmFlowGridGeom = this->problem(pmFlowId).gridGeometry(); |
| 379 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
|
3 | const auto& poroMechGridGeom = this->problem(poroMechId).gridGeometry(); |
| 380 | |||
| 381 | // make sure the two grids are really the same. Note that if the two grids | ||
| 382 | // happen to have equal number of elements by chance, we don't detect this source of error. | ||
| 383 |
1/2✗ Branch 4 not taken.
✓ Branch 5 taken 3 times.
|
6 | if (pmFlowGridGeom.gridView().size(0) != poroMechGridGeom.gridView().size(0)) |
| 384 | ✗ | DUNE_THROW(Dune::InvalidStateException, "The two sub-problems are assumed to operate on the same mesh!"); | |
| 385 | |||
| 386 | 6 | pmFlowCouplingMap_.resize(pmFlowGridGeom.gridView().size(0)); | |
| 387 | static constexpr int dim = GridView<PMFlowId>::dimension; | ||
| 388 |
1/2✓ Branch 2 taken 231 times.
✗ Branch 3 not taken.
|
462 | for (const auto& element : elements(pmFlowGridGeom.gridView())) |
| 389 | { | ||
| 390 | 228 | const auto eIdx = pmFlowGridGeom.elementMapper().index(element); | |
| 391 | |||
| 392 | // firstly, the element couples to the nodal dofs in itself | ||
| 393 |
4/4✓ Branch 2 taken 1424 times.
✓ Branch 3 taken 228 times.
✓ Branch 4 taken 1424 times.
✓ Branch 5 taken 228 times.
|
1652 | for (int i = 0; i < element.geometry().corners(); ++i) |
| 394 | 4272 | pmFlowCouplingMap_[eIdx].push_back( poroMechGridGeom.vertexMapper().subIndex(element, i , dim) ); | |
| 395 | |||
| 396 | // the pm flow problem couples to the same elements as in its own stencil | ||
| 397 | // due to the dependency of the residual on all permeabilities in its stencil, | ||
| 398 | // which in turn depend on the mechanical deformations. | ||
| 399 | 456 | const auto& inverseConnectivity = pmFlowGridGeom.connectivityMap()[eIdx]; | |
| 400 |
4/4✓ Branch 0 taken 936 times.
✓ Branch 1 taken 228 times.
✓ Branch 2 taken 936 times.
✓ Branch 3 taken 228 times.
|
1620 | for (const auto& dataJ : inverseConnectivity) |
| 401 |
4/4✓ Branch 1 taken 6048 times.
✓ Branch 2 taken 936 times.
✓ Branch 3 taken 6048 times.
✓ Branch 4 taken 936 times.
|
6984 | for (int i = 0; i < element.geometry().corners(); ++i) |
| 402 | 18144 | pmFlowCouplingMap_[dataJ.globalJ].push_back( poroMechGridGeom.vertexMapper().subIndex(element, i , dim) ); | |
| 403 | } | ||
| 404 | |||
| 405 | // make stencils unique | ||
| 406 |
4/4✓ Branch 0 taken 228 times.
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 228 times.
✓ Branch 3 taken 3 times.
|
465 | for (auto& stencil : pmFlowCouplingMap_) |
| 407 | { | ||
| 408 | 684 | std::sort(stencil.begin(), stencil.end()); | |
| 409 | 1140 | stencil.erase(std::unique(stencil.begin(), stencil.end()), stencil.end()); | |
| 410 | } | ||
| 411 | 3 | } | |
| 412 | |||
| 413 | // Container for storing the coupling element stencils for the pm flow domain | ||
| 414 | std::vector< CouplingStencilType<PMFlowId> > pmFlowCouplingMap_; | ||
| 415 | |||
| 416 | // the coupling context of the poromechanics domain | ||
| 417 | PoroMechanicsCouplingContext poroMechCouplingContext_; | ||
| 418 | }; | ||
| 419 | |||
| 420 | } //end namespace Dumux | ||
| 421 | |||
| 422 | #endif | ||
| 423 |