| 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 DarcyDarcyCoupling | ||
| 10 | * \brief Coupling manager for equal-dimension boundary coupling | ||
| 11 | */ | ||
| 12 | |||
| 13 | #ifndef DUMUX_MULTIDOMAIN_DARCYDARCY_BOUNDARY_COUPLINGMANAGER_HH | ||
| 14 | #define DUMUX_MULTIDOMAIN_DARCYDARCY_BOUNDARY_COUPLINGMANAGER_HH | ||
| 15 | |||
| 16 | #include <iostream> | ||
| 17 | #include <vector> | ||
| 18 | |||
| 19 | #include <dune/common/indices.hh> | ||
| 20 | #include <dune/common/exceptions.hh> | ||
| 21 | |||
| 22 | #include <dumux/common/properties.hh> | ||
| 23 | #include <dumux/common/math.hh> | ||
| 24 | #include <dumux/common/numeqvector.hh> | ||
| 25 | #include <dumux/discretization/elementsolution.hh> | ||
| 26 | #include <dumux/discretization/method.hh> | ||
| 27 | #include <dumux/discretization/cellcentered/tpfa/computetransmissibility.hh> | ||
| 28 | #include <dumux/multidomain/couplingmanager.hh> | ||
| 29 | #include <dumux/multidomain/boundary/darcydarcy/couplingmapper.hh> | ||
| 30 | |||
| 31 | namespace Dumux { | ||
| 32 | |||
| 33 | /*! | ||
| 34 | * \ingroup DarcyDarcyCoupling | ||
| 35 | * \brief Coupling manager for equal-dimension boundary coupling of darcy models | ||
| 36 | */ | ||
| 37 | template<class MDTraits> | ||
| 38 |
0/2✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
3 | class DarcyDarcyBoundaryCouplingManager |
| 39 | : public CouplingManager<MDTraits> | ||
| 40 | { | ||
| 41 | using ParentType = CouplingManager<MDTraits>; | ||
| 42 | |||
| 43 | using Scalar = typename MDTraits::Scalar; | ||
| 44 | using SolutionVector = typename MDTraits::SolutionVector; | ||
| 45 | |||
| 46 | template<std::size_t i> using SubDomainTypeTag = typename MDTraits::template SubDomain<i>::TypeTag; | ||
| 47 | template<std::size_t i> using Problem = GetPropType<SubDomainTypeTag<i>, Properties::Problem>; | ||
| 48 | template<std::size_t i> using PrimaryVariables = GetPropType<SubDomainTypeTag<i>, Properties::PrimaryVariables>; | ||
| 49 | template<std::size_t i> using NumEqVector = Dumux::NumEqVector<PrimaryVariables<i>>; | ||
| 50 | template<std::size_t i> using ElementVolumeVariables = typename GetPropType<SubDomainTypeTag<i>, Properties::GridVolumeVariables>::LocalView; | ||
| 51 | template<std::size_t i> using VolumeVariables = typename GetPropType<SubDomainTypeTag<i>, Properties::GridVolumeVariables>::VolumeVariables; | ||
| 52 | template<std::size_t i> using GridGeometry = typename MDTraits::template SubDomain<i>::GridGeometry; | ||
| 53 | template<std::size_t i> using FVElementGeometry = typename GridGeometry<i>::LocalView; | ||
| 54 | template<std::size_t i> using SubControlVolumeFace = typename GridGeometry<i>::SubControlVolumeFace; | ||
| 55 | template<std::size_t i> using SubControlVolume = typename GridGeometry<i>::SubControlVolume; | ||
| 56 | template<std::size_t i> using GridView = typename GridGeometry<i>::GridView; | ||
| 57 | template<std::size_t i> using Element = typename GridView<i>::template Codim<0>::Entity; | ||
| 58 | |||
| 59 | template<std::size_t i> | ||
| 60 | static constexpr auto domainIdx() | ||
| 61 | { return typename MDTraits::template SubDomain<i>::Index{}; } | ||
| 62 | |||
| 63 | template<std::size_t i> | ||
| 64 | static constexpr bool isCCTpfa() | ||
| 65 | { return GridGeometry<i>::discMethod == DiscretizationMethods::cctpfa; } | ||
| 66 | |||
| 67 | using CouplingStencil = std::vector<std::size_t>; | ||
| 68 | public: | ||
| 69 | //! export traits | ||
| 70 | using MultiDomainTraits = MDTraits; | ||
| 71 | //! export stencil types | ||
| 72 | using CouplingStencils = std::unordered_map<std::size_t, CouplingStencil>; | ||
| 73 | |||
| 74 | /*! | ||
| 75 | * \brief Methods to be accessed by main | ||
| 76 | */ | ||
| 77 | // \{ | ||
| 78 | |||
| 79 | 3 | void init(std::shared_ptr<Problem<domainIdx<0>()>> problem0, | |
| 80 | std::shared_ptr<Problem<domainIdx<1>()>> problem1, | ||
| 81 | const SolutionVector& curSol) | ||
| 82 | { | ||
| 83 | static_assert(isCCTpfa<0>() && isCCTpfa<1>(), "Only cctpfa implemented!"); | ||
| 84 | |||
| 85 | 6 | this->setSubProblems(std::make_tuple(problem0, problem1)); | |
| 86 | 3 | this->updateSolution(curSol); | |
| 87 | 3 | couplingMapper_.update(*this); | |
| 88 | 3 | } | |
| 89 | |||
| 90 | // \} | ||
| 91 | |||
| 92 | /*! | ||
| 93 | * \brief Methods to be accessed by the assembly | ||
| 94 | */ | ||
| 95 | // \{ | ||
| 96 | |||
| 97 | /*! | ||
| 98 | * \brief returns an iterable container of all indices of degrees of freedom of domain j | ||
| 99 | * that couple with / influence the element residual of the given element of domain i | ||
| 100 | * | ||
| 101 | * \param domainI the domain index of domain i | ||
| 102 | * \param element the coupled element of domain í | ||
| 103 | * \param domainJ the domain index of domain j | ||
| 104 | */ | ||
| 105 | template<std::size_t i, std::size_t j> | ||
| 106 | ✗ | const CouplingStencil& couplingStencil(Dune::index_constant<i> domainI, | |
| 107 | const Element<i>& element, | ||
| 108 | Dune::index_constant<j> domainJ) const | ||
| 109 | { | ||
| 110 | static_assert(i != j, "A domain cannot be coupled to itself!"); | ||
| 111 | ✗ | const auto eIdx = this->problem(domainI).gridGeometry().elementMapper().index(element); | |
| 112 | ✗ | return couplingMapper_.couplingStencil(domainI, eIdx, domainJ); | |
| 113 | } | ||
| 114 | |||
| 115 | // \} | ||
| 116 | |||
| 117 | /*! | ||
| 118 | * \brief If the boundary entity is on a coupling boundary | ||
| 119 | * \param domainI the domain index for which to compute the flux | ||
| 120 | * \param scvf the sub control volume face | ||
| 121 | */ | ||
| 122 | template<std::size_t i> | ||
| 123 | ✗ | bool isCoupled(Dune::index_constant<i> domainI, | |
| 124 | const SubControlVolumeFace<i>& scvf) const | ||
| 125 | { | ||
| 126 |
4/6✓ Branch 1 taken 105878 times.
✓ Branch 2 taken 744 times.
✓ Branch 4 taken 1014 times.
✓ Branch 5 taken 264 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
107900 | return couplingMapper_.isCoupled(domainI, scvf); |
| 127 | } | ||
| 128 | |||
| 129 | /*! | ||
| 130 | * \brief Evaluate the boundary conditions for a coupled Neumann boundary segment. | ||
| 131 | * | ||
| 132 | * This is the method for the case where the Neumann condition is | ||
| 133 | * potentially solution dependent | ||
| 134 | * | ||
| 135 | * \param domainI the domain index for which to compute the flux | ||
| 136 | * \param element The finite element | ||
| 137 | * \param fvGeometry The finite-volume geometry | ||
| 138 | * \param elemVolVars All volume variables for the element | ||
| 139 | * \param scvf The sub control volume face | ||
| 140 | * \param phaseIdx the phase for which to compute the flux | ||
| 141 | * | ||
| 142 | * Negative values mean influx. | ||
| 143 | * E.g. for the mass balance that would the mass flux in \f$ [ kg / (m^2 \cdot s)] \f$. | ||
| 144 | */ | ||
| 145 | template<std::size_t i> | ||
| 146 | 201152 | Scalar advectiveFluxCoupling(Dune::index_constant<i> domainI, | |
| 147 | const Element<i>& element, | ||
| 148 | const FVElementGeometry<i>& fvGeometry, | ||
| 149 | const ElementVolumeVariables<i>& elemVolVars, | ||
| 150 | const SubControlVolumeFace<i>& scvf, | ||
| 151 | int phaseIdx = 0) const | ||
| 152 | { | ||
| 153 | 201152 | Scalar flux = 0.0; | |
| 154 |
6/10✓ Branch 0 taken 6 times.
✓ Branch 1 taken 100570 times.
✓ Branch 3 taken 6 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 6 times.
✓ Branch 8 taken 6 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 6 times.
✗ Branch 12 not taken.
|
201152 | static const bool enableGravity = getParamFromGroup<bool>(this->problem(domainI).paramGroup(), "Problem.EnableGravity"); |
| 155 | constexpr auto otherDomainIdx = domainIdx<1-i>(); | ||
| 156 | |||
| 157 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 100576 times.
|
201152 | const auto& outsideElement = this->problem(otherDomainIdx).gridGeometry().element(couplingMapper_.outsideElementIndex(domainI, scvf)); |
| 158 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 100576 times.
✓ Branch 3 taken 100576 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 100576 times.
✗ Branch 7 not taken.
|
201152 | auto fvGeometryOutside = localView(this->problem(otherDomainIdx).gridGeometry()); |
| 159 |
1/2✓ Branch 1 taken 100576 times.
✗ Branch 2 not taken.
|
201152 | fvGeometryOutside.bindElement(outsideElement); |
| 160 | |||
| 161 |
1/2✓ Branch 1 taken 100576 times.
✗ Branch 2 not taken.
|
201152 | const auto& flipScvf = fvGeometryOutside.scvf(couplingMapper_.flipScvfIndex(domainI, scvf)); |
| 162 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 100576 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 100576 times.
|
402304 | const auto& outsideScv = fvGeometryOutside.scv(flipScvf.insideScvIdx()); |
| 163 |
1/2✓ Branch 1 taken 100576 times.
✗ Branch 2 not taken.
|
201152 | const auto outsideVolVars = volVars(otherDomainIdx, outsideElement, outsideScv); |
| 164 | |||
| 165 |
2/4✓ Branch 0 taken 100576 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 100576 times.
✗ Branch 3 not taken.
|
402304 | const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx()); |
| 166 | 201152 | const auto& insideVolVars = elemVolVars[insideScv]; | |
| 167 | |||
| 168 | using Extrusion = typename GridGeometry<i>::Extrusion; | ||
| 169 | 201152 | const auto ti = computeTpfaTransmissibility(fvGeometry, scvf, insideScv, insideVolVars.permeability(), insideVolVars.extrusionFactor()); | |
| 170 | 201152 | const auto tj = computeTpfaTransmissibility(fvGeometryOutside, flipScvf, outsideScv, outsideVolVars.permeability(), outsideVolVars.extrusionFactor()); | |
| 171 | 201152 | Scalar tij = 0.0; | |
| 172 |
1/2✓ Branch 0 taken 100576 times.
✗ Branch 1 not taken.
|
201152 | if (ti*tj > 0.0) |
| 173 | 402304 | tij = Extrusion::area(fvGeometry, scvf)*(ti * tj)/(ti + tj); | |
| 174 | |||
| 175 |
2/2✓ Branch 0 taken 864 times.
✓ Branch 1 taken 99712 times.
|
201152 | if (enableGravity) |
| 176 | { | ||
| 177 | // do averaging for the density over all neighboring elements | ||
| 178 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 864 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 864 times.
|
3456 | const auto rho = (insideVolVars.density(phaseIdx) + outsideVolVars.density(phaseIdx))*0.5; |
| 179 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 864 times.
|
1728 | const auto& g = this->problem(domainI).spatialParams().gravity(scvf.ipGlobal()); |
| 180 | |||
| 181 | //! compute alpha := n^T*K*g | ||
| 182 | 3456 | const auto alpha_inside = vtmv(scvf.unitOuterNormal(), insideVolVars.permeability(), g)*insideVolVars.extrusionFactor(); | |
| 183 | 3456 | const auto alpha_outside = vtmv(scvf.unitOuterNormal(), outsideVolVars.permeability(), g)*outsideVolVars.extrusionFactor(); | |
| 184 | |||
| 185 | // Obtain inside and outside pressures | ||
| 186 | 1728 | const auto pInside = insideVolVars.pressure(0); | |
| 187 | 1728 | const auto pOutside = outsideVolVars.pressure(0); | |
| 188 | |||
| 189 | 3456 | flux = tij*(pInside - pOutside) + rho*Extrusion::area(fvGeometry, scvf)*alpha_inside - rho*tij/tj*(alpha_inside - alpha_outside); | |
| 190 | |||
| 191 | } | ||
| 192 | else | ||
| 193 | { | ||
| 194 | // Obtain inside and outside pressures | ||
| 195 | 199424 | const auto pInside = insideVolVars.pressure(phaseIdx); | |
| 196 | 199424 | const auto pOutside = outsideVolVars.pressure(phaseIdx); | |
| 197 | |||
| 198 | // return flux | ||
| 199 | 199424 | flux = tij*(pInside - pOutside); | |
| 200 | } | ||
| 201 | |||
| 202 | // upwind scheme | ||
| 203 |
4/6✓ Branch 0 taken 6 times.
✓ Branch 1 taken 100570 times.
✓ Branch 3 taken 6 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 6 times.
✗ Branch 7 not taken.
|
201152 | static const Scalar upwindWeight = getParam<Scalar>("Flux.UpwindWeight"); |
| 204 | 402304 | auto upwindTerm = [phaseIdx](const auto& volVars){ return volVars.density(phaseIdx)*volVars.mobility(phaseIdx); }; | |
| 205 |
4/4✓ Branch 0 taken 51140 times.
✓ Branch 1 taken 49436 times.
✓ Branch 2 taken 51140 times.
✓ Branch 3 taken 49436 times.
|
402304 | if (std::signbit(flux)) // if sign of flux is negative |
| 206 | 102280 | flux *= (upwindWeight*upwindTerm(outsideVolVars) | |
| 207 | 204560 | + (1.0 - upwindWeight)*upwindTerm(insideVolVars)); | |
| 208 | else | ||
| 209 | 98872 | flux *= (upwindWeight*upwindTerm(insideVolVars) | |
| 210 | 197744 | + (1.0 - upwindWeight)*upwindTerm(outsideVolVars)); | |
| 211 | |||
| 212 | // make this a area-specific flux | ||
| 213 | 402304 | flux /= Extrusion::area(fvGeometry, scvf)*insideVolVars.extrusionFactor(); | |
| 214 | 201152 | return flux; | |
| 215 | } | ||
| 216 | |||
| 217 | //! Return the volume variables of domain i for a given element and scv | ||
| 218 | template<std::size_t i> | ||
| 219 | 100576 | VolumeVariables<i> volVars(Dune::index_constant<i> domainI, | |
| 220 | const Element<i>& element, | ||
| 221 | const SubControlVolume<i>& scv) const | ||
| 222 | { | ||
| 223 | 100576 | VolumeVariables<i> volVars; | |
| 224 | 100576 | const auto elemSol = elementSolution(element, this->curSol(domainI), this->problem(domainI).gridGeometry()); | |
| 225 | 100816 | volVars.update(elemSol, this->problem(domainI), element, scv); | |
| 226 | 100576 | return volVars; | |
| 227 | } | ||
| 228 | |||
| 229 | private: | ||
| 230 | DarcyDarcyBoundaryCouplingMapper<MDTraits> couplingMapper_; | ||
| 231 | }; | ||
| 232 | |||
| 233 | } // end namespace Dumux | ||
| 234 | |||
| 235 | #endif | ||
| 236 |