| 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 FacetCoupling | ||
| 10 | * \copydoc Dumux::MpfaOFacetCouplingInteractionVolumeAssembler | ||
| 11 | */ | ||
| 12 | #ifndef DUMUX_MULTIDOMAIN_FACET_CC_MPFA_O_LOCAL_ASSEMBLER_HH | ||
| 13 | #define DUMUX_MULTIDOMAIN_FACET_CC_MPFA_O_LOCAL_ASSEMBLER_HH | ||
| 14 | |||
| 15 | #include <vector> | ||
| 16 | |||
| 17 | #include <dune/common/exceptions.hh> | ||
| 18 | #include <dune/common/float_cmp.hh> | ||
| 19 | |||
| 20 | #include <dumux/common/math.hh> | ||
| 21 | #include <dumux/common/parameters.hh> | ||
| 22 | #include <dumux/discretization/extrusion.hh> | ||
| 23 | #include <dumux/discretization/cellcentered/mpfa/methods.hh> | ||
| 24 | #include <dumux/discretization/cellcentered/mpfa/localassemblerbase.hh> | ||
| 25 | #include <dumux/discretization/cellcentered/mpfa/localassemblerhelper.hh> | ||
| 26 | #include <dumux/discretization/cellcentered/mpfa/computetransmissibility.hh> | ||
| 27 | |||
| 28 | namespace Dumux { | ||
| 29 | |||
| 30 | /*! | ||
| 31 | * \ingroup FacetCoupling | ||
| 32 | * \brief Specialization of the interaction volume-local | ||
| 33 | * assembler class for the schemes using an mpfa-o type | ||
| 34 | * assembly in the context of facet coupling. | ||
| 35 | * | ||
| 36 | * \tparam P The problem type | ||
| 37 | * \tparam EG The element finite volume geometry | ||
| 38 | * \tparam EV The element volume variables type | ||
| 39 | */ | ||
| 40 | template< class P, class EG, class EV > | ||
| 41 | class MpfaOFacetCouplingInteractionVolumeAssembler | ||
| 42 | : public InteractionVolumeAssemblerBase< P, EG, EV > | ||
| 43 | { | ||
| 44 | using ParentType = InteractionVolumeAssemblerBase< P, EG, EV >; | ||
| 45 | using Helper = InteractionVolumeAssemblerHelper; | ||
| 46 | using Extrusion = Extrusion_t<typename EG::GridGeometry>; | ||
| 47 | |||
| 48 | template< class IV > | ||
| 49 | using Scalar = typename IV::Traits::MatVecTraits::FaceVector::value_type; | ||
| 50 | |||
| 51 | public: | ||
| 52 | //! Pull up constructor of the base class | ||
| 53 |
3/3✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1102200 times.
✓ Branch 2 taken 3162310 times.
|
12568260 | using ParentType::ParentType; |
| 54 | |||
| 55 | /*! | ||
| 56 | * \brief Assembles the matrices involved in the flux expressions | ||
| 57 | * and the local system of equations within an interaction volume. | ||
| 58 | * | ||
| 59 | * \param handle The data handle in which the matrices are stored | ||
| 60 | * \param iv The interaction volume | ||
| 61 | * \param getT Lambda to evaluate the scv-wise tensors | ||
| 62 | * \param wijZeroThresh Threshold below which transmissibilities are taken to be zero. | ||
| 63 | * On the basis of this threshold, trivial (0 = 0) rows in the | ||
| 64 | * A matrix are identified and modified accordingly in order to | ||
| 65 | * avoid ending up with singular matrices. This can occur when the | ||
| 66 | * tensor is zero in some cells. | ||
| 67 | */ | ||
| 68 | template< class DataHandle, class IV, class TensorFunc > | ||
| 69 | 24259596 | void assembleMatrices(DataHandle& handle, IV& iv, const TensorFunc& getT, Scalar<IV> wijZeroThresh = 0.0) | |
| 70 | { | ||
| 71 |
1/2✓ Branch 1 taken 12568260 times.
✗ Branch 2 not taken.
|
24259596 | assembleLocalMatrices_(handle.A(), handle.AB(), handle.CA(), handle.T(), handle.omegas(), iv, getT, wijZeroThresh); |
| 72 | |||
| 73 | // maybe solve the local system | ||
| 74 |
1/2✓ Branch 0 taken 12568260 times.
✗ Branch 1 not taken.
|
24259596 | if (iv.numUnknowns() > 0) |
| 75 | 48519192 | Helper::solveLocalSystem(this->fvGeometry(), handle, iv); | |
| 76 | 24259596 | } | |
| 77 | |||
| 78 | /*! | ||
| 79 | * \brief Assembles the vector of primary (cell) unknowns and (maybe) | ||
| 80 | * Dirichlet boundary conditions within an interaction volume. | ||
| 81 | * | ||
| 82 | * \param handle The data handle in which the vector is stored | ||
| 83 | * \param iv The interaction volume | ||
| 84 | * \param getU Lambda to obtain the desired cell/Dirichlet value from vol vars | ||
| 85 | */ | ||
| 86 | template< class DataHandle, class IV, class GetU > | ||
| 87 | 24259596 | void assembleU(DataHandle& handle, const IV& iv, const GetU& getU) | |
| 88 | { | ||
| 89 | 24259596 | auto& u = handle.uj(); | |
| 90 | 24259596 | Helper::resizeVector(u, iv.numKnowns()); | |
| 91 | |||
| 92 | // put the cell unknowns first ... | ||
| 93 | 24259596 | typename IV::Traits::IndexSet::LocalIndexType i = 0; | |
| 94 |
2/2✓ Branch 0 taken 49146853 times.
✓ Branch 1 taken 12568260 times.
|
118842853 | for (; i < iv.numScvs(); i++) |
| 95 | 94583257 | u[i] = getU( this->elemVolVars()[iv.localScv(i).gridScvIndex()] ); | |
| 96 | |||
| 97 | // ... then put facet unknowns ... | ||
| 98 |
2/2✓ Branch 0 taken 3042014 times.
✓ Branch 1 taken 12568260 times.
|
30190868 | for (const auto& data : iv.interiorBoundaryData()) |
| 99 | { | ||
| 100 | 5931272 | const auto& scvf = this->fvGeometry().scvf(data.scvfIndex()); | |
| 101 | 5931272 | const auto element = this->fvGeometry().gridGeometry().element(scvf.insideScvIdx()); | |
| 102 |
1/2✓ Branch 1 taken 3042014 times.
✗ Branch 2 not taken.
|
5931272 | u[i++] = getU( this->problem().couplingManager().getLowDimVolVars(element, scvf) ); |
| 103 | } | ||
| 104 | |||
| 105 | // ... then put the Dirichlet unknowns | ||
| 106 |
2/2✓ Branch 0 taken 697580 times.
✓ Branch 1 taken 12568260 times.
|
25640276 | for (const auto& data : iv.dirichletData()) |
| 107 | 1380680 | u[i++] = getU( this->elemVolVars()[data.volVarIndex()] ); | |
| 108 | 24259596 | } | |
| 109 | |||
| 110 | /*! | ||
| 111 | * \brief Assembles the gravitational flux contributions on the scvfs within an | ||
| 112 | * interaction volume. | ||
| 113 | * | ||
| 114 | * \param handle The data handle in which the vector is stored | ||
| 115 | * \param iv The interaction volume | ||
| 116 | * \param getRho Lambda to obtain the density from volume variables | ||
| 117 | */ | ||
| 118 | template< class DataHandle, class IV, class GetRho > | ||
| 119 | 204840 | void assembleGravity(DataHandle& handle, const IV& iv, const GetRho& getRho) | |
| 120 | { | ||
| 121 | using GridView = typename IV::Traits::GridView; | ||
| 122 | static constexpr int dim = GridView::dimension; | ||
| 123 | static constexpr int dimWorld = GridView::dimensionworld; | ||
| 124 | static constexpr bool isSurfaceGrid = dim < dimWorld; | ||
| 125 | |||
| 126 | // resize the gravity vectors | ||
| 127 | 204840 | auto& g = handle.g(); | |
| 128 | 204840 | auto& deltaG = handle.deltaG(); | |
| 129 | 204840 | auto& outsideG = handle.gOutside(); | |
| 130 | 204840 | Helper::resizeVector(g, iv.numFaces()); | |
| 131 | 204840 | Helper::resizeVector(deltaG, iv.numUnknowns()); | |
| 132 | if (isSurfaceGrid) | ||
| 133 | 68280 | Helper::resizeVector(outsideG, iv.numFaces()); | |
| 134 | |||
| 135 | //! For each face, we... | ||
| 136 | //! - arithmetically average the phase densities | ||
| 137 | //! - compute the term \f$ \alpha := \mathbf{A} \rho \ \mathbf{n}^T \mathbf{K} \mathbf{g} \f$ in each neighboring cell | ||
| 138 | //! - compute \f$ \alpha^* = \sum{\alpha_{outside, i}} - \alpha_{inside} \f$ | ||
| 139 | using Scalar = typename IV::Traits::MatVecTraits::TMatrix::value_type; | ||
| 140 | using LocalIndexType = typename IV::Traits::IndexSet::LocalIndexType; | ||
| 141 | |||
| 142 | // xi factor for coupling conditions | ||
| 143 |
4/6✓ Branch 0 taken 3 times.
✓ Branch 1 taken 204837 times.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 3 times.
✗ Branch 7 not taken.
|
204840 | static const Scalar xi = getParamFromGroup<Scalar>(this->problem().paramGroup(), "FacetCoupling.Xi", 1.0); |
| 144 | |||
| 145 |
2/2✓ Branch 0 taken 835380 times.
✓ Branch 1 taken 204840 times.
|
1040220 | for (LocalIndexType faceIdx = 0; faceIdx < iv.numFaces(); ++faceIdx) |
| 146 | { | ||
| 147 | // gravitational acceleration on this face | ||
| 148 | 835380 | const auto& curLocalScvf = iv.localScvf(faceIdx); | |
| 149 | 835380 | const auto& curGlobalScvf = this->fvGeometry().scvf(curLocalScvf.gridScvfIndex()); | |
| 150 |
2/2✓ Branch 0 taken 49896 times.
✓ Branch 1 taken 785484 times.
|
835380 | const auto& gravity = this->problem().spatialParams().gravity(curGlobalScvf.ipGlobal()); |
| 151 | 835380 | const auto curIsInteriorBoundary = curLocalScvf.isOnInteriorBoundary(); | |
| 152 |
5/6✓ Branch 0 taken 49896 times.
✓ Branch 1 taken 785484 times.
✓ Branch 2 taken 49896 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 49896 times.
✓ Branch 5 taken 785484 times.
|
835380 | const Scalar curXiFactor = curIsInteriorBoundary ? (curGlobalScvf.boundary() ? 1.0 : xi) : 1.0; |
| 153 | |||
| 154 | // get permeability tensor in "positive" sub volume | ||
| 155 | 835380 | const auto& neighborScvIndices = curLocalScvf.neighboringLocalScvIndices(); | |
| 156 |
2/2✓ Branch 0 taken 212388 times.
✓ Branch 1 taken 622992 times.
|
835380 | const auto& posGlobalScv = this->fvGeometry().scv(iv.localScv(neighborScvIndices[0]).gridScvIndex()); |
| 157 | 835380 | const auto& posVolVars = this->elemVolVars()[posGlobalScv]; | |
| 158 |
2/2✓ Branch 0 taken 817668 times.
✓ Branch 1 taken 17712 times.
|
1670760 | const auto alpha_inside = posVolVars.extrusionFactor()*vtmv(curGlobalScvf.unitOuterNormal(), |
| 159 | 835380 | posVolVars.permeability(), | |
| 160 | gravity); | ||
| 161 | |||
| 162 |
2/2✓ Branch 0 taken 817668 times.
✓ Branch 1 taken 17712 times.
|
835380 | const auto numOutsideFaces = !curGlobalScvf.boundary() ? curGlobalScvf.numOutsideScvs() : 0; |
| 163 | using OutsideAlphaStorage = std::conditional_t< isSurfaceGrid, | ||
| 164 | std::vector<Scalar>, | ||
| 165 | Dune::ReservedVector<Scalar, 1> >; | ||
| 166 |
1/2✓ Branch 1 taken 278460 times.
✗ Branch 2 not taken.
|
835380 | OutsideAlphaStorage alpha_outside; alpha_outside.resize(numOutsideFaces); |
| 167 | 835380 | std::fill(alpha_outside.begin(), alpha_outside.end(), 0.0); | |
| 168 | Scalar rho; | ||
| 169 | |||
| 170 |
2/2✓ Branch 0 taken 261828 times.
✓ Branch 1 taken 16632 times.
|
278460 | if (isSurfaceGrid && !curIsInteriorBoundary) |
| 171 |
1/2✓ Branch 1 taken 261828 times.
✗ Branch 2 not taken.
|
261828 | Helper::resizeVector(outsideG[faceIdx], numOutsideFaces); |
| 172 | |||
| 173 |
2/2✓ Branch 0 taken 831276 times.
✓ Branch 1 taken 4104 times.
|
835380 | if (!curLocalScvf.isDirichlet()) |
| 174 | { | ||
| 175 |
2/2✓ Branch 0 taken 49896 times.
✓ Branch 1 taken 781380 times.
|
831276 | const auto localDofIdx = curLocalScvf.localDofIndex(); |
| 176 |
4/4✓ Branch 0 taken 49896 times.
✓ Branch 1 taken 781380 times.
✓ Branch 2 taken 817668 times.
✓ Branch 3 taken 13608 times.
|
831276 | const Scalar curOneMinusXi = curIsInteriorBoundary ? -(1.0 - xi) : 1.0; |
| 177 | |||
| 178 | 831276 | rho = getRho(posVolVars); | |
| 179 |
2/2✓ Branch 0 taken 817668 times.
✓ Branch 1 taken 13608 times.
|
831276 | deltaG[localDofIdx] = 0.0; |
| 180 | |||
| 181 |
2/2✓ Branch 0 taken 817668 times.
✓ Branch 1 taken 13608 times.
|
831276 | if (!curGlobalScvf.boundary()) |
| 182 | { | ||
| 183 |
2/2✓ Branch 0 taken 817668 times.
✓ Branch 1 taken 817668 times.
|
1635336 | for (unsigned int idxInOutside = 0; idxInOutside < curGlobalScvf.numOutsideScvs(); ++idxInOutside) |
| 184 | { | ||
| 185 | // obtain outside tensor | ||
| 186 |
2/2✓ Branch 0 taken 209928 times.
✓ Branch 1 taken 607740 times.
|
817668 | const auto negLocalScvIdx = neighborScvIndices[idxInOutside+1]; |
| 187 |
2/2✓ Branch 0 taken 209928 times.
✓ Branch 1 taken 607740 times.
|
817668 | const auto& negGlobalScv = this->fvGeometry().scv(iv.localScv(negLocalScvIdx).gridScvIndex()); |
| 188 | 817668 | const auto& negVolVars = this->elemVolVars()[negGlobalScv]; | |
| 189 | 545112 | const auto& flipScvf = !isSurfaceGrid ? curGlobalScvf | |
| 190 | 272556 | : this->fvGeometry().flipScvf(curGlobalScvf.index(), idxInOutside); | |
| 191 | |||
| 192 |
2/2✓ Branch 0 taken 767772 times.
✓ Branch 1 taken 49896 times.
|
1635336 | alpha_outside[idxInOutside] = negVolVars.extrusionFactor()*vtmv(flipScvf.unitOuterNormal(), |
| 193 | 817668 | negVolVars.permeability(), | |
| 194 | gravity); | ||
| 195 | if (isSurfaceGrid) | ||
| 196 | 272556 | alpha_outside[idxInOutside] *= -1.0; | |
| 197 | |||
| 198 |
2/2✓ Branch 0 taken 767772 times.
✓ Branch 1 taken 49896 times.
|
817668 | if (!curIsInteriorBoundary) |
| 199 | 767772 | rho += getRho(negVolVars); | |
| 200 | |||
| 201 | 817668 | deltaG[localDofIdx] += curOneMinusXi*alpha_outside[idxInOutside]; | |
| 202 | } | ||
| 203 | } | ||
| 204 | |||
| 205 |
2/2✓ Branch 0 taken 49896 times.
✓ Branch 1 taken 781380 times.
|
831276 | if (curIsInteriorBoundary) |
| 206 | { | ||
| 207 |
1/2✓ Branch 1 taken 16632 times.
✗ Branch 2 not taken.
|
49896 | const auto posElement = this->fvGeometry().gridGeometry().element(posGlobalScv.elementIndex()); |
| 208 |
1/2✓ Branch 1 taken 49896 times.
✗ Branch 2 not taken.
|
49896 | const auto& facetVolVars = this->problem().couplingManager().getLowDimVolVars(posElement, curGlobalScvf); |
| 209 | 49896 | rho += getRho(facetVolVars); | |
| 210 | 49896 | rho /= 2.0; | |
| 211 | 49896 | const auto alphaFacet = posVolVars.extrusionFactor()*vtmv(curGlobalScvf.unitOuterNormal(), | |
| 212 | 49896 | facetVolVars.permeability(), | |
| 213 | gravity); | ||
| 214 | 49896 | deltaG[localDofIdx] += alphaFacet; | |
| 215 | 49896 | } | |
| 216 | else | ||
| 217 | 781380 | rho /= numOutsideFaces + 1; | |
| 218 | |||
| 219 | 831276 | deltaG[localDofIdx] -= curXiFactor*alpha_inside; | |
| 220 | 831276 | deltaG[localDofIdx] *= rho*Extrusion::area(this->fvGeometry(), curGlobalScvf); | |
| 221 | } | ||
| 222 | // use average density between facet and cell density on interior Dirichlet boundaries | ||
| 223 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 4104 times.
|
4104 | else if (curIsInteriorBoundary) |
| 224 | { | ||
| 225 | ✗ | const auto posElement = this->fvGeometry().gridGeometry().element(posGlobalScv.elementIndex()); | |
| 226 | ✗ | const auto& facetVolVars = this->problem().couplingManager().getLowDimVolVars(posElement, curGlobalScvf); | |
| 227 | ✗ | rho = 0.5*(getRho(facetVolVars) + getRho(posVolVars)); | |
| 228 | ✗ | } | |
| 229 | // use density resulting from Dirichlet BCs | ||
| 230 | else | ||
| 231 | 4104 | rho = getRho(this->elemVolVars()[curGlobalScvf.outsideScvIdx()]); | |
| 232 | |||
| 233 | // add "inside" & "outside" alphas to gravity containers | ||
| 234 |
2/2✓ Branch 0 taken 261828 times.
✓ Branch 1 taken 16632 times.
|
835380 | g[faceIdx] = alpha_inside*rho*Extrusion::area(this->fvGeometry(), curGlobalScvf); |
| 235 | |||
| 236 |
2/2✓ Branch 0 taken 261828 times.
✓ Branch 1 taken 16632 times.
|
278460 | if (isSurfaceGrid && !curIsInteriorBoundary) |
| 237 | { | ||
| 238 | unsigned int i = 0; | ||
| 239 |
2/2✓ Branch 0 taken 255924 times.
✓ Branch 1 taken 261828 times.
|
517752 | for (const auto& alpha : alpha_outside) |
| 240 | 255924 | outsideG[faceIdx][i++] = alpha*rho*Extrusion::area(this->fvGeometry(), curGlobalScvf); | |
| 241 | } | ||
| 242 | } | ||
| 243 | |||
| 244 | // add iv-wide contributions to gravity vectors | ||
| 245 |
1/2✓ Branch 2 taken 68280 times.
✗ Branch 3 not taken.
|
204840 | handle.CA().umv(deltaG, g); |
| 246 | if (isSurfaceGrid) | ||
| 247 | { | ||
| 248 | using FaceVector = typename IV::Traits::MatVecTraits::FaceVector; | ||
| 249 |
1/2✓ Branch 1 taken 68280 times.
✗ Branch 2 not taken.
|
68280 | FaceVector AG; |
| 250 |
1/2✓ Branch 1 taken 68280 times.
✗ Branch 2 not taken.
|
68280 | Helper::resizeVector(AG, iv.numUnknowns()); |
| 251 | 68280 | handle.A().mv(deltaG, AG); | |
| 252 | |||
| 253 | // compute gravitational accelerations | ||
| 254 |
4/4✓ Branch 0 taken 278460 times.
✓ Branch 1 taken 255924 times.
✓ Branch 2 taken 534384 times.
✓ Branch 3 taken 68280 times.
|
602664 | for (const auto& localFaceData : iv.localFaceData()) |
| 255 | { | ||
| 256 | // continue only for "outside" faces | ||
| 257 |
2/2✓ Branch 0 taken 278460 times.
✓ Branch 1 taken 255924 times.
|
534384 | if (!localFaceData.isOutsideFace()) |
| 258 | 278460 | continue; | |
| 259 | |||
| 260 | 255924 | const auto localScvIdx = localFaceData.ivLocalInsideScvIndex(); | |
| 261 | 255924 | const auto localScvfIdx = localFaceData.ivLocalScvfIndex(); | |
| 262 | 255924 | const auto idxInOutside = localFaceData.scvfLocalOutsideScvfIndex(); | |
| 263 | 255924 | const auto& posLocalScv = iv.localScv(localScvIdx); | |
| 264 | 255924 | const auto& wijk = handle.omegas()[localScvfIdx][idxInOutside+1]; | |
| 265 | |||
| 266 | // add contributions from all local directions | ||
| 267 |
2/2✓ Branch 0 taken 511848 times.
✓ Branch 1 taken 255924 times.
|
767772 | for (LocalIndexType localDir = 0; localDir < dim; localDir++) |
| 268 | { | ||
| 269 | // the scvf corresponding to this local direction in the scv | ||
| 270 |
2/2✓ Branch 0 taken 511176 times.
✓ Branch 1 taken 672 times.
|
511848 | const auto& curLocalScvf = iv.localScvf(posLocalScv.localScvfIndex(localDir)); |
| 271 |
2/2✓ Branch 0 taken 511176 times.
✓ Branch 1 taken 672 times.
|
511848 | if (!curLocalScvf.isDirichlet()) |
| 272 | 511176 | outsideG[localScvfIdx][idxInOutside] -= wijk[localDir]*AG[curLocalScvf.localDofIndex()]; | |
| 273 | } | ||
| 274 | } | ||
| 275 | 68280 | } | |
| 276 | 204840 | } | |
| 277 | |||
| 278 | private: | ||
| 279 | /*! | ||
| 280 | * \copydoc Dumux::MpfaOInteractionVolumeAssembler::assembleLocalMatrices_ | ||
| 281 | */ | ||
| 282 | template< class IV, class TensorFunc > | ||
| 283 | 24259596 | void assembleLocalMatrices_(typename IV::Traits::MatVecTraits::AMatrix& A, | |
| 284 | typename IV::Traits::MatVecTraits::BMatrix& B, | ||
| 285 | typename IV::Traits::MatVecTraits::CMatrix& C, | ||
| 286 | typename IV::Traits::MatVecTraits::DMatrix& D, | ||
| 287 | typename IV::Traits::MatVecTraits::OmegaStorage& wijk, | ||
| 288 | IV& iv, const TensorFunc& getT, | ||
| 289 | Scalar<IV> wijZeroThresh) | ||
| 290 | { | ||
| 291 | using Scalar = typename IV::Traits::MatVecTraits::TMatrix::value_type; | ||
| 292 | using LocalIndexType = typename IV::Traits::IndexSet::LocalIndexType; | ||
| 293 | static constexpr int dim = IV::Traits::GridView::dimension; | ||
| 294 | static constexpr int dimWorld = IV::Traits::GridView::dimensionworld; | ||
| 295 | |||
| 296 | // xi factor for coupling conditions | ||
| 297 |
4/6✓ Branch 0 taken 18 times.
✓ Branch 1 taken 12568242 times.
✓ Branch 3 taken 18 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 18 times.
✗ Branch 7 not taken.
|
24259596 | static const Scalar xi = getParamFromGroup<Scalar>(this->problem().paramGroup(), "FacetCoupling.Xi", 1.0); |
| 298 | |||
| 299 | // On surface grids only xi = 1.0 can be used, as the coupling condition | ||
| 300 | // for xi != 1.0 does not generalize for surface grids where the normal | ||
| 301 | // vectors of the inside/outside elements have different orientations. | ||
| 302 | if (dim < dimWorld) | ||
| 303 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 7755028 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 7755028 times.
|
15261836 | if (Dune::FloatCmp::ne(xi, 1.0, 1e-6)) |
| 304 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Xi != 1.0 cannot be used on surface grids"); | |
| 305 | |||
| 306 |
1/2✓ Branch 1 taken 12568260 times.
✗ Branch 2 not taken.
|
24259596 | std::vector< std::pair<LocalIndexType, LocalIndexType> > faceMarkers; |
| 307 |
2/4✓ Branch 1 taken 12568260 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 12568260 times.
|
48519192 | Helper::resizeVector(wijk, iv.numFaces()); |
| 308 | |||
| 309 | // if only Dirichlet faces are present in the iv, | ||
| 310 | // the matrices A, B & C are undefined and D = T | ||
| 311 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 12568260 times.
|
24259596 | if (iv.numUnknowns() == 0) |
| 312 | { | ||
| 313 | // resize & reset D matrix | ||
| 314 | ✗ | Helper::resizeMatrix(D, iv.numFaces(), iv.numKnowns()); D = 0.0; | |
| 315 | |||
| 316 | // Loop over all the faces, in this case these are all dirichlet boundaries | ||
| 317 | ✗ | for (LocalIndexType faceIdx = 0; faceIdx < iv.numFaces(); ++faceIdx) | |
| 318 | { | ||
| 319 | ✗ | const auto& curLocalScvf = iv.localScvf(faceIdx); | |
| 320 | ✗ | const auto& curGlobalScvf = this->fvGeometry().scvf(curLocalScvf.gridScvfIndex()); | |
| 321 | ✗ | const auto& neighborScvIndices = curLocalScvf.neighboringLocalScvIndices(); | |
| 322 | |||
| 323 | // get tensor in "positive" sub volume | ||
| 324 | ✗ | const auto& posLocalScv = iv.localScv(neighborScvIndices[0]); | |
| 325 | ✗ | const auto& posGlobalScv = this->fvGeometry().scv(posLocalScv.gridScvIndex()); | |
| 326 | ✗ | const auto& posVolVars = this->elemVolVars()[posGlobalScv]; | |
| 327 | ✗ | const auto tensor = getT(posVolVars); | |
| 328 | |||
| 329 | // the omega factors of the "positive" sub volume | ||
| 330 | ✗ | Helper::resizeVector(wijk[faceIdx], /*no outside scvs present*/1); | |
| 331 | ✗ | wijk[faceIdx][0] = computeMpfaTransmissibility<EG>(this->fvGeometry(), posLocalScv, curGlobalScvf, tensor, posVolVars.extrusionFactor()); | |
| 332 | |||
| 333 | ✗ | const auto posScvLocalDofIdx = posLocalScv.localDofIndex(); | |
| 334 | ✗ | for (LocalIndexType localDir = 0; localDir < dim; localDir++) | |
| 335 | { | ||
| 336 | ✗ | const auto& otherLocalScvf = iv.localScvf( posLocalScv.localScvfIndex(localDir) ); | |
| 337 | ✗ | const auto otherLocalDofIdx = otherLocalScvf.localDofIndex(); | |
| 338 | ✗ | D[faceIdx][otherLocalDofIdx] -= wijk[faceIdx][0][localDir]; | |
| 339 | ✗ | D[faceIdx][posScvLocalDofIdx] += wijk[faceIdx][0][localDir]; | |
| 340 | } | ||
| 341 | } | ||
| 342 | } | ||
| 343 | else | ||
| 344 | { | ||
| 345 | // resize & reset matrices | ||
| 346 |
2/4✓ Branch 1 taken 12568260 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 12568260 times.
✗ Branch 5 not taken.
|
72778788 | Helper::resizeMatrix(A, iv.numUnknowns(), iv.numUnknowns()); A = 0.0; |
| 347 |
2/4✓ Branch 1 taken 12568260 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 12568260 times.
✗ Branch 5 not taken.
|
72778788 | Helper::resizeMatrix(B, iv.numUnknowns(), iv.numKnowns()); B = 0.0; |
| 348 |
2/4✓ Branch 1 taken 12568260 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 12568260 times.
✗ Branch 5 not taken.
|
72778788 | Helper::resizeMatrix(C, iv.numFaces(), iv.numUnknowns()); C = 0.0; |
| 349 |
1/2✓ Branch 1 taken 12568260 times.
✗ Branch 2 not taken.
|
48519192 | Helper::resizeMatrix(D, iv.numFaces(), iv.numKnowns()); D = 0.0; |
| 350 | |||
| 351 |
2/2✓ Branch 0 taken 52874531 times.
✓ Branch 1 taken 12568260 times.
|
126118049 | for (LocalIndexType faceIdx = 0; faceIdx < iv.numFaces(); ++faceIdx) |
| 352 | { | ||
| 353 | 101858453 | const auto& curLocalScvf = iv.localScvf(faceIdx); | |
| 354 | 101858453 | const auto& curGlobalScvf = this->fvGeometry().scvf(curLocalScvf.gridScvfIndex()); | |
| 355 |
2/2✓ Branch 0 taken 6077224 times.
✓ Branch 1 taken 46797307 times.
|
101858453 | const auto curIsDirichlet = curLocalScvf.isDirichlet(); |
| 356 |
2/2✓ Branch 0 taken 6077224 times.
✓ Branch 1 taken 46797307 times.
|
101858453 | const auto curLocalDofIdx = curLocalScvf.localDofIndex(); |
| 357 | 101858453 | const auto curIsInteriorBoundary = curLocalScvf.isOnInteriorBoundary(); | |
| 358 |
6/6✓ Branch 0 taken 6077224 times.
✓ Branch 1 taken 46797307 times.
✓ Branch 2 taken 6070420 times.
✓ Branch 3 taken 6804 times.
✓ Branch 4 taken 6077224 times.
✓ Branch 5 taken 46797307 times.
|
101858453 | const Scalar curXiFactor = curIsInteriorBoundary ? (curGlobalScvf.boundary() ? 1.0 : xi) : 1.0; |
| 359 | |||
| 360 | // get tensor in "positive" sub volume | ||
| 361 |
2/2✓ Branch 0 taken 14140234 times.
✓ Branch 1 taken 38734297 times.
|
101858453 | const auto& neighborScvIndices = curLocalScvf.neighboringLocalScvIndices(); |
| 362 |
2/2✓ Branch 0 taken 14140234 times.
✓ Branch 1 taken 38734297 times.
|
101858453 | const auto& posLocalScv = iv.localScv(neighborScvIndices[0]); |
| 363 |
2/2✓ Branch 0 taken 14140234 times.
✓ Branch 1 taken 38734297 times.
|
101858453 | const auto& posGlobalScv = this->fvGeometry().scv(posLocalScv.gridScvIndex()); |
| 364 | 101858453 | const auto& posVolVars = this->elemVolVars()[posGlobalScv]; | |
| 365 |
1/2✓ Branch 1 taken 21226332 times.
✗ Branch 2 not taken.
|
101858453 | const auto& posElement = iv.element(neighborScvIndices[0]); |
| 366 |
1/2✓ Branch 1 taken 32468748 times.
✗ Branch 2 not taken.
|
101858453 | const auto tensor = getT(posVolVars); |
| 367 | |||
| 368 | // the omega factors of the "positive" sub volume | ||
| 369 |
1/2✓ Branch 1 taken 32468748 times.
✗ Branch 2 not taken.
|
101858453 | Helper::resizeVector(wijk[faceIdx], neighborScvIndices.size()); |
| 370 | 101858453 | wijk[faceIdx][0] = computeMpfaTransmissibility<EG>(this->fvGeometry(), posLocalScv, curGlobalScvf, tensor, posVolVars.extrusionFactor()); | |
| 371 | |||
| 372 | using std::abs; | ||
| 373 | 101858453 | bool isZeroWij = false; | |
| 374 | |||
| 375 | // go over the coordinate directions in the positive sub volume | ||
| 376 |
2/2✓ Branch 0 taken 105749062 times.
✓ Branch 1 taken 52874531 times.
|
305575359 | for (unsigned int localDir = 0; localDir < dim; localDir++) |
| 377 | { | ||
| 378 |
2/2✓ Branch 0 taken 52874531 times.
✓ Branch 1 taken 52874531 times.
|
203716906 | const auto& otherLocalScvf = iv.localScvf( posLocalScv.localScvfIndex(localDir) ); |
| 379 |
2/2✓ Branch 0 taken 52874531 times.
✓ Branch 1 taken 52874531 times.
|
203716906 | const auto otherLocalDofIdx = otherLocalScvf.localDofIndex(); |
| 380 | |||
| 381 |
2/2✓ Branch 0 taken 52874531 times.
✓ Branch 1 taken 52874531 times.
|
203716906 | if (otherLocalDofIdx == curLocalDofIdx) |
| 382 | { | ||
| 383 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 52874531 times.
|
101858453 | if (abs(wijk[faceIdx][0][localDir]) <= wijZeroThresh) |
| 384 | { | ||
| 385 | ✗ | if (!curIsDirichlet) | |
| 386 | { | ||
| 387 | ✗ | isZeroWij = true; | |
| 388 | ✗ | faceMarkers.emplace_back( std::make_pair(curLocalDofIdx, faceIdx) ); | |
| 389 | } | ||
| 390 | } | ||
| 391 | } | ||
| 392 | |||
| 393 | // if we are not on a Dirichlet face, add entries associated with unknown face pressures | ||
| 394 | // i.e. in matrix C and maybe A (if current face is not a Dirichlet face) | ||
| 395 |
2/2✓ Branch 0 taken 104686181 times.
✓ Branch 1 taken 1062881 times.
|
203716906 | if (!otherLocalScvf.isDirichlet()) |
| 396 | { | ||
| 397 |
2/2✓ Branch 0 taken 103981797 times.
✓ Branch 1 taken 704384 times.
|
201623165 | C[faceIdx][otherLocalDofIdx] -= wijk[faceIdx][0][localDir]; |
| 398 |
2/2✓ Branch 0 taken 103981797 times.
✓ Branch 1 taken 704384 times.
|
201623165 | if (!curIsDirichlet) |
| 399 | 200235681 | A[curLocalDofIdx][otherLocalDofIdx] -= curXiFactor*wijk[faceIdx][0][localDir]; | |
| 400 | } | ||
| 401 | // the current face is a Dirichlet face and creates entries in D & maybe B | ||
| 402 | else | ||
| 403 | { | ||
| 404 |
2/2✓ Branch 0 taken 358497 times.
✓ Branch 1 taken 704384 times.
|
2093741 | D[faceIdx][otherLocalDofIdx] -= wijk[faceIdx][0][localDir]; |
| 405 |
2/2✓ Branch 0 taken 358497 times.
✓ Branch 1 taken 704384 times.
|
2093741 | if (!curIsDirichlet) |
| 406 | 706257 | B[curLocalDofIdx][otherLocalDofIdx] += curXiFactor*wijk[faceIdx][0][localDir]; | |
| 407 | } | ||
| 408 | |||
| 409 | // add entries related to pressures at the scv centers (dofs) | ||
| 410 |
2/2✓ Branch 0 taken 104340294 times.
✓ Branch 1 taken 1408768 times.
|
203716906 | const auto posScvLocalDofIdx = posLocalScv.localDofIndex(); |
| 411 |
2/2✓ Branch 0 taken 104340294 times.
✓ Branch 1 taken 1408768 times.
|
203716906 | D[faceIdx][posScvLocalDofIdx] += wijk[faceIdx][0][localDir]; |
| 412 | |||
| 413 |
2/2✓ Branch 0 taken 104340294 times.
✓ Branch 1 taken 1408768 times.
|
203716906 | if (!curIsDirichlet) |
| 414 | 200941938 | B[curLocalDofIdx][posScvLocalDofIdx] -= curXiFactor*wijk[faceIdx][0][localDir]; | |
| 415 | } | ||
| 416 | |||
| 417 | // handle the entries related to the coupled facet element | ||
| 418 |
2/2✓ Branch 0 taken 6070420 times.
✓ Branch 1 taken 46804111 times.
|
101858453 | if (curIsInteriorBoundary && !curIsDirichlet) |
| 419 | { | ||
| 420 |
1/2✓ Branch 1 taken 6070420 times.
✗ Branch 2 not taken.
|
11848936 | const auto facetVolVars = this->problem().couplingManager().getLowDimVolVars(posElement, curGlobalScvf); |
| 421 | 11848936 | const auto facetTensor = getT(facetVolVars); | |
| 422 | |||
| 423 | // On surface grids we use the square root of the extrusion factor as approximation of the aperture | ||
| 424 | using std::sqrt; | ||
| 425 | 11848936 | const auto wFacet = 2.0*Extrusion::area(this->fvGeometry(), curGlobalScvf)*posVolVars.extrusionFactor() | |
| 426 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 6070420 times.
|
23697872 | *vtmv(curGlobalScvf.unitOuterNormal(), facetTensor, curGlobalScvf.unitOuterNormal()) |
| 427 | 11848936 | / (dim < dimWorld ? sqrt(facetVolVars.extrusionFactor()) : facetVolVars.extrusionFactor()); | |
| 428 | |||
| 429 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 6070420 times.
|
11848936 | A[curLocalDofIdx][curLocalDofIdx] -= wFacet; |
| 430 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 6070420 times.
✓ Branch 2 taken 6070420 times.
✗ Branch 3 not taken.
|
11848936 | B[curLocalDofIdx][curLocalScvf.coupledFacetLocalDofIndex()] -= wFacet; |
| 431 | |||
| 432 | // check for zero transmissibilities (skip if inside has been zero already) | ||
| 433 |
2/4✓ Branch 0 taken 6070420 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 6070420 times.
|
11848936 | if (!isZeroWij && abs(wFacet) <= wijZeroThresh) |
| 434 | { | ||
| 435 | ✗ | isZeroWij = true; | |
| 436 | ✗ | faceMarkers.emplace_back( std::make_pair(curLocalDofIdx, faceIdx) ); | |
| 437 | } | ||
| 438 | } | ||
| 439 | |||
| 440 | // If we are on an interior face (which isn't coupled via Dirichlet Bs), add values from negative sub volume | ||
| 441 |
3/4✓ Branch 0 taken 51489595 times.
✓ Branch 1 taken 1384936 times.
✓ Branch 2 taken 51489595 times.
✗ Branch 3 not taken.
|
101858453 | if (!curGlobalScvf.boundary() && !curIsDirichlet) |
| 442 | { | ||
| 443 | // the minus ensures the right sign of the coupling condition in the matrices | ||
| 444 |
2/2✓ Branch 0 taken 6070420 times.
✓ Branch 1 taken 45419175 times.
|
99156997 | const Scalar curOneMinusXi = curIsInteriorBoundary ? -(1.0 - xi) : 1.0; |
| 445 | |||
| 446 | // loop over all the outside neighbors of this face and add entries | ||
| 447 |
2/2✓ Branch 0 taken 51489595 times.
✓ Branch 1 taken 51489595 times.
|
198313994 | for (unsigned int idxInOutside = 0; idxInOutside < curGlobalScvf.numOutsideScvs(); ++idxInOutside) |
| 448 | { | ||
| 449 | 99156997 | const auto idxOnScvf = idxInOutside+1; | |
| 450 |
2/2✓ Branch 0 taken 12543100 times.
✓ Branch 1 taken 38946495 times.
|
99156997 | const auto& negLocalScv = iv.localScv( neighborScvIndices[idxOnScvf] ); |
| 451 |
2/2✓ Branch 0 taken 12543100 times.
✓ Branch 1 taken 38946495 times.
|
99156997 | const auto& negGlobalScv = this->fvGeometry().scv(negLocalScv.gridScvIndex()); |
| 452 | 99156997 | const auto& negVolVars = this->elemVolVars()[negGlobalScv]; | |
| 453 | 99156997 | const auto negTensor = getT(negVolVars); | |
| 454 | |||
| 455 | // On surface grids, use outside face for "negative" transmissibility calculation | ||
| 456 | 99156997 | const auto& scvf = dim < dimWorld ? this->fvGeometry().flipScvf(curGlobalScvf.index(), idxInOutside) | |
| 457 | : curGlobalScvf; | ||
| 458 | 99156997 | wijk[faceIdx][idxOnScvf] = computeMpfaTransmissibility<EG>(this->fvGeometry(), negLocalScv, scvf, negTensor, negVolVars.extrusionFactor()); | |
| 459 | |||
| 460 | // flip sign on surface grids (since we used the "outside" normal) | ||
| 461 | if (dim < dimWorld) | ||
| 462 | 62222258 | wijk[faceIdx][idxOnScvf] *= -1.0; | |
| 463 | |||
| 464 | // go over the coordinate directions in the negative sub volume | ||
| 465 |
2/2✓ Branch 0 taken 102979190 times.
✓ Branch 1 taken 51489595 times.
|
297470991 | for (int localDir = 0; localDir < dim; localDir++) |
| 466 | { | ||
| 467 |
2/2✓ Branch 0 taken 45419175 times.
✓ Branch 1 taken 57560015 times.
|
198313994 | const auto otherLocalScvfIdx = negLocalScv.localScvfIndex(localDir); |
| 468 | 198313994 | const auto& otherLocalScvf = iv.localScvf(otherLocalScvfIdx); | |
| 469 |
2/2✓ Branch 0 taken 45419175 times.
✓ Branch 1 taken 57560015 times.
|
198313994 | const auto otherLocalDofIdx = otherLocalScvf.localDofIndex(); |
| 470 | |||
| 471 | // check for zero transmissibilities (skip if inside has been zero already) | ||
| 472 |
3/4✓ Branch 0 taken 45419175 times.
✓ Branch 1 taken 57560015 times.
✓ Branch 2 taken 45419175 times.
✗ Branch 3 not taken.
|
198313994 | if (otherLocalDofIdx == curLocalDofIdx && !isZeroWij) |
| 473 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 45419175 times.
|
87308061 | if (abs(wijk[faceIdx][idxOnScvf][localDir]) <= wijZeroThresh) |
| 474 | ✗ | faceMarkers.emplace_back( std::make_pair(curLocalDofIdx, faceIdx) ); | |
| 475 | |||
| 476 |
2/2✓ Branch 0 taken 102633303 times.
✓ Branch 1 taken 345887 times.
|
198313994 | if (!otherLocalScvf.isDirichlet()) |
| 477 | 197632767 | A[curLocalDofIdx][otherLocalDofIdx] += curOneMinusXi*wijk[faceIdx][idxOnScvf][localDir]; | |
| 478 | else | ||
| 479 | 681227 | B[curLocalDofIdx][otherLocalDofIdx] -= curOneMinusXi*wijk[faceIdx][idxOnScvf][localDir]; | |
| 480 | |||
| 481 | // add entries to matrix B | ||
| 482 | 198313994 | B[curLocalDofIdx][negLocalScv.localDofIndex()] += curOneMinusXi*wijk[faceIdx][idxOnScvf][localDir]; | |
| 483 | } | ||
| 484 | } | ||
| 485 | } | ||
| 486 | } | ||
| 487 | } | ||
| 488 | 24259596 | } | |
| 489 | }; | ||
| 490 | |||
| 491 | } // end namespace | ||
| 492 | |||
| 493 | #endif | ||
| 494 |