| 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 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 |
6/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1102198 times.
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 3223028 times.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 2120830 times.
|
18887640 | 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 | 18010716 | void assembleMatrices(DataHandle& handle, IV& iv, const TensorFunc& getT, Scalar<IV> wijZeroThresh = 0.0) | |
| 70 | { | ||
| 71 | 108064296 | 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 9443820 times.
✗ Branch 1 not taken.
|
18010716 | if (iv.numUnknowns() > 0) |
| 75 | 18010716 | Helper::solveLocalSystem(this->fvGeometry(), handle, iv); | |
| 76 | 18010716 | } | |
| 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 | 18010716 | void assembleU(DataHandle& handle, const IV& iv, const GetU& getU) | |
| 88 | { | ||
| 89 | 18010716 | auto& u = handle.uj(); | |
| 90 | 18010716 | Helper::resizeVector(u, iv.numKnowns()); | |
| 91 | |||
| 92 | // put the cell unknowns first ... | ||
| 93 | 18010716 | typename IV::Traits::IndexSet::LocalIndexType i = 0; | |
| 94 |
4/4✓ Branch 0 taken 37002193 times.
✓ Branch 1 taken 9443820 times.
✓ Branch 2 taken 37002193 times.
✓ Branch 3 taken 9443820 times.
|
176609306 | for (; i < iv.numScvs(); i++) |
| 95 | 140587874 | u[i] = getU( this->elemVolVars()[iv.localScv(i).gridScvIndex()] ); | |
| 96 | |||
| 97 | // ... then put facet unknowns ... | ||
| 98 |
4/4✓ Branch 0 taken 2256944 times.
✓ Branch 1 taken 9443820 times.
✓ Branch 2 taken 2256944 times.
✓ Branch 3 taken 9443820 times.
|
80765128 | for (const auto& data : iv.interiorBoundaryData()) |
| 99 | { | ||
| 100 | 4361132 | const auto& scvf = this->fvGeometry().scvf(data.scvfIndex()); | |
| 101 |
1/2✓ Branch 1 taken 2256944 times.
✗ Branch 2 not taken.
|
8722264 | const auto element = this->fvGeometry().gridGeometry().element(scvf.insideScvIdx()); |
| 102 |
2/4✓ Branch 1 taken 2256944 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2256944 times.
✗ Branch 5 not taken.
|
8722264 | u[i++] = getU( this->problem().couplingManager().getLowDimVolVars(element, scvf) ); |
| 103 | } | ||
| 104 | |||
| 105 | // ... then put the Dirichlet unknowns | ||
| 106 |
4/4✓ Branch 0 taken 516080 times.
✓ Branch 1 taken 9443820 times.
✓ Branch 2 taken 516080 times.
✓ Branch 3 taken 9443820 times.
|
74078224 | for (const auto& data : iv.dirichletData()) |
| 107 | 1017680 | u[i++] = getU( this->elemVolVars()[data.volVarIndex()] ); | |
| 108 | 18010716 | } | |
| 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 |
5/8✓ 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.
✓ Branch 9 taken 3 times.
✗ Branch 10 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 |
6/6✓ Branch 0 taken 49896 times.
✓ Branch 1 taken 785484 times.
✓ Branch 2 taken 49896 times.
✓ Branch 3 taken 785484 times.
✓ Branch 4 taken 49896 times.
✓ Branch 5 taken 785484 times.
|
2506140 | const auto& gravity = this->problem().spatialParams().gravity(curGlobalScvf.ipGlobal()); |
| 151 | 835380 | const auto curIsInteriorBoundary = curLocalScvf.isOnInteriorBoundary(); | |
| 152 |
3/4✓ Branch 0 taken 49896 times.
✓ Branch 1 taken 785484 times.
✓ Branch 2 taken 49896 times.
✗ Branch 3 not taken.
|
835380 | const Scalar curXiFactor = curIsInteriorBoundary ? (curGlobalScvf.boundary() ? 1.0 : xi) : 1.0; |
| 153 | |||
| 154 | // get permeability tensor in "positive" sub volume | ||
| 155 |
2/2✓ Branch 0 taken 49896 times.
✓ Branch 1 taken 785484 times.
|
835380 | const auto& neighborScvIndices = curLocalScvf.neighboringLocalScvIndices(); |
| 156 |
6/6✓ Branch 0 taken 212388 times.
✓ Branch 1 taken 622992 times.
✓ Branch 2 taken 212388 times.
✓ Branch 3 taken 622992 times.
✓ Branch 4 taken 141592 times.
✓ Branch 5 taken 415328 times.
|
2227680 | const auto& posGlobalScv = this->fvGeometry().scv(iv.localScv(neighborScvIndices[0]).gridScvIndex()); |
| 157 | 835380 | const auto& posVolVars = this->elemVolVars()[posGlobalScv]; | |
| 158 | 1670760 | const auto alpha_inside = posVolVars.extrusionFactor()*vtmv(curGlobalScvf.unitOuterNormal(), | |
| 159 | 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 |
4/8✓ Branch 1 taken 278460 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 278460 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 272556 times.
✓ Branch 7 taken 5904 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
1664856 | OutsideAlphaStorage alpha_outside; alpha_outside.resize(numOutsideFaces); |
| 167 | 2506140 | 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 |
2/4✓ Branch 1 taken 261828 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 261828 times.
✗ Branch 5 not taken.
|
523656 | 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 | 831276 | const auto localDofIdx = curLocalScvf.localDofIndex(); | |
| 176 |
2/2✓ Branch 0 taken 49896 times.
✓ Branch 1 taken 781380 times.
|
831276 | const Scalar curOneMinusXi = curIsInteriorBoundary ? -(1.0 - xi) : 1.0; |
| 177 | |||
| 178 |
2/2✓ Branch 0 taken 817668 times.
✓ Branch 1 taken 13608 times.
|
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 |
4/4✓ Branch 0 taken 817668 times.
✓ Branch 1 taken 817668 times.
✓ Branch 2 taken 817668 times.
✓ Branch 3 taken 817668 times.
|
3270672 | 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 |
4/4✓ Branch 0 taken 209928 times.
✓ Branch 1 taken 607740 times.
✓ Branch 2 taken 209928 times.
✓ Branch 3 taken 607740 times.
|
1635336 | const auto& negGlobalScv = this->fvGeometry().scv(iv.localScv(negLocalScvIdx).gridScvIndex()); |
| 188 | 817668 | const auto& negVolVars = this->elemVolVars()[negGlobalScv]; | |
| 189 | 1090224 | 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.
|
2453004 | alpha_outside[idxInOutside] = negVolVars.extrusionFactor()*vtmv(flipScvf.unitOuterNormal(), |
| 193 | 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 | 1535544 | rho += getRho(negVolVars); | |
| 200 | |||
| 201 | 2180448 | 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.
|
99792 | const auto posElement = this->fvGeometry().gridGeometry().element(posGlobalScv.elementIndex()); |
| 208 |
2/4✓ Branch 1 taken 49896 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 49896 times.
✗ Branch 5 not taken.
|
99792 | const auto& facetVolVars = this->problem().couplingManager().getLowDimVolVars(posElement, curGlobalScvf); |
| 209 | 49896 | rho += getRho(facetVolVars); | |
| 210 | 49896 | rho /= 2.0; | |
| 211 | 99792 | const auto alphaFacet = posVolVars.extrusionFactor()*vtmv(curGlobalScvf.unitOuterNormal(), | |
| 212 | facetVolVars.permeability(), | ||
| 213 | gravity); | ||
| 214 | 99792 | deltaG[localDofIdx] += alphaFacet; | |
| 215 | } | ||
| 216 | else | ||
| 217 | 781380 | rho /= numOutsideFaces + 1; | |
| 218 | |||
| 219 | 831276 | deltaG[localDofIdx] -= curXiFactor*alpha_inside; | |
| 220 | 2493828 | 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 | 8208 | rho = getRho(this->elemVolVars()[curGlobalScvf.outsideScvIdx()]); | |
| 232 | |||
| 233 | // add "inside" & "outside" alphas to gravity containers | ||
| 234 |
4/4✓ Branch 0 taken 261828 times.
✓ Branch 1 taken 16632 times.
✓ Branch 2 taken 261828 times.
✓ Branch 3 taken 16632 times.
|
2227680 | 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 | 261828 | unsigned int i = 0; | |
| 239 |
4/4✓ Branch 0 taken 255924 times.
✓ Branch 1 taken 261828 times.
✓ Branch 2 taken 255924 times.
✓ Branch 3 taken 261828 times.
|
1297332 | for (const auto& alpha : alpha_outside) |
| 240 | 1023696 | outsideG[faceIdx][i++] = alpha*rho*Extrusion::area(this->fvGeometry(), curGlobalScvf); | |
| 241 | } | ||
| 242 | } | ||
| 243 | |||
| 244 | // add iv-wide contributions to gravity vectors | ||
| 245 | 409680 | handle.CA().umv(deltaG, g); | |
| 246 | if (isSurfaceGrid) | ||
| 247 | { | ||
| 248 | using FaceVector = typename IV::Traits::MatVecTraits::FaceVector; | ||
| 249 |
4/10✓ Branch 1 taken 68280 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 68280 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 68280 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 68280 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
|
273120 | FaceVector AG; |
| 250 |
1/2✓ Branch 1 taken 68280 times.
✗ Branch 2 not taken.
|
68280 | Helper::resizeVector(AG, iv.numUnknowns()); |
| 251 | 136560 | handle.A().mv(deltaG, AG); | |
| 252 | |||
| 253 | // compute gravitational accelerations | ||
| 254 |
4/4✓ Branch 0 taken 534384 times.
✓ Branch 1 taken 68280 times.
✓ Branch 2 taken 534384 times.
✓ Branch 3 taken 68280 times.
|
807504 | for (const auto& localFaceData : iv.localFaceData()) |
| 255 | { | ||
| 256 | // continue only for "outside" faces | ||
| 257 |
2/2✓ Branch 0 taken 255924 times.
✓ Branch 1 taken 278460 times.
|
534384 | if (!localFaceData.isOutsideFace()) |
| 258 | 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 | 767772 | 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 |
4/4✓ Branch 0 taken 511176 times.
✓ Branch 1 taken 672 times.
✓ Branch 2 taken 511176 times.
✓ Branch 3 taken 672 times.
|
1023696 | 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 | 2555880 | outsideG[localScvfIdx][idxInOutside] -= wijk[localDir]*AG[curLocalScvf.localDofIndex()]; | |
| 273 | } | ||
| 274 | } | ||
| 275 | } | ||
| 276 | 204840 | } | |
| 277 | |||
| 278 | private: | ||
| 279 | /*! | ||
| 280 | * \copydoc Dumux::MpfaOInteractionVolumeAssembler::assembleLocalMatrices_ | ||
| 281 | */ | ||
| 282 | template< class IV, class TensorFunc > | ||
| 283 | 18010716 | 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 |
5/8✓ Branch 0 taken 18 times.
✓ Branch 1 taken 9443802 times.
✓ Branch 3 taken 18 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 18 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 18 times.
✗ Branch 10 not taken.
|
18010716 | 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 4630588 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4630588 times.
|
9012956 | 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 |
2/6✓ Branch 1 taken 9443820 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 9443820 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
36021432 | std::vector< std::pair<LocalIndexType, LocalIndexType> > faceMarkers; |
| 307 |
1/2✓ Branch 1 taken 9443820 times.
✗ Branch 2 not taken.
|
18010716 | 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 9443820 times.
|
18010716 | 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 |
1/2✓ Branch 1 taken 9443820 times.
✗ Branch 2 not taken.
|
18010716 | Helper::resizeMatrix(A, iv.numUnknowns(), iv.numUnknowns()); A = 0.0; |
| 347 |
1/2✓ Branch 1 taken 9443820 times.
✗ Branch 2 not taken.
|
18010716 | Helper::resizeMatrix(B, iv.numUnknowns(), iv.numKnowns()); B = 0.0; |
| 348 |
1/2✓ Branch 1 taken 9443820 times.
✗ Branch 2 not taken.
|
18010716 | Helper::resizeMatrix(C, iv.numFaces(), iv.numUnknowns()); C = 0.0; |
| 349 |
1/2✓ Branch 1 taken 9443820 times.
✗ Branch 2 not taken.
|
18010716 | Helper::resizeMatrix(D, iv.numFaces(), iv.numKnowns()); D = 0.0; |
| 350 | |||
| 351 |
2/2✓ Branch 0 taken 39769901 times.
✓ Branch 1 taken 9443820 times.
|
93659909 | for (LocalIndexType faceIdx = 0; faceIdx < iv.numFaces(); ++faceIdx) |
| 352 | { | ||
| 353 | 75649193 | const auto& curLocalScvf = iv.localScvf(faceIdx); | |
| 354 | 75649193 | const auto& curGlobalScvf = this->fvGeometry().scvf(curLocalScvf.gridScvfIndex()); | |
| 355 | 75649193 | const auto curIsDirichlet = curLocalScvf.isDirichlet(); | |
| 356 | 75649193 | const auto curLocalDofIdx = curLocalScvf.localDofIndex(); | |
| 357 | 75649193 | const auto curIsInteriorBoundary = curLocalScvf.isOnInteriorBoundary(); | |
| 358 |
4/4✓ Branch 0 taken 4507084 times.
✓ Branch 1 taken 35262817 times.
✓ Branch 2 taken 4500280 times.
✓ Branch 3 taken 6804 times.
|
75649193 | 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 4507084 times.
✓ Branch 1 taken 35262817 times.
|
75649193 | const auto& neighborScvIndices = curLocalScvf.neighboringLocalScvIndices(); |
| 362 |
4/4✓ Branch 0 taken 10600654 times.
✓ Branch 1 taken 29169247 times.
✓ Branch 2 taken 5391236 times.
✓ Branch 3 taken 15014547 times.
|
113579380 | const auto& posLocalScv = iv.localScv(neighborScvIndices[0]); |
| 363 |
2/2✓ Branch 0 taken 10600654 times.
✓ Branch 1 taken 29169247 times.
|
75649193 | const auto& posGlobalScv = this->fvGeometry().scv(posLocalScv.gridScvIndex()); |
| 364 | 75649193 | const auto& posVolVars = this->elemVolVars()[posGlobalScv]; | |
| 365 |
1/2✓ Branch 1 taken 12489912 times.
✗ Branch 2 not taken.
|
113579380 | const auto& posElement = iv.element(neighborScvIndices[0]); |
| 366 |
1/2✓ Branch 1 taken 12489912 times.
✗ Branch 2 not taken.
|
75649193 | const auto tensor = getT(posVolVars); |
| 367 | |||
| 368 | // the omega factors of the "positive" sub volume | ||
| 369 |
3/6✓ Branch 1 taken 19364118 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 19364118 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 19364118 times.
✗ Branch 8 not taken.
|
189017392 | Helper::resizeVector(wijk[faceIdx], neighborScvIndices.size()); |
| 370 | 75649193 | wijk[faceIdx][0] = computeMpfaTransmissibility<EG>(this->fvGeometry(), posLocalScv, curGlobalScvf, tensor, posVolVars.extrusionFactor()); | |
| 371 | |||
| 372 | using std::abs; | ||
| 373 | 75649193 | bool isZeroWij = false; | |
| 374 | |||
| 375 | // go over the coordinate directions in the positive sub volume | ||
| 376 |
2/2✓ Branch 0 taken 79539802 times.
✓ Branch 1 taken 39769901 times.
|
226947579 | for (unsigned int localDir = 0; localDir < dim; localDir++) |
| 377 | { | ||
| 378 |
4/4✓ Branch 0 taken 39769901 times.
✓ Branch 1 taken 39769901 times.
✓ Branch 2 taken 39769901 times.
✓ Branch 3 taken 39769901 times.
|
302596772 | const auto& otherLocalScvf = iv.localScvf( posLocalScv.localScvfIndex(localDir) ); |
| 379 | 151298386 | const auto otherLocalDofIdx = otherLocalScvf.localDofIndex(); | |
| 380 | |||
| 381 |
2/2✓ Branch 0 taken 39769901 times.
✓ Branch 1 taken 39769901 times.
|
151298386 | if (otherLocalDofIdx == curLocalDofIdx) |
| 382 | { | ||
| 383 |
5/10✗ Branch 0 not taken.
✓ Branch 1 taken 39769901 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 39769901 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 39769901 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 39769901 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 20405783 times.
|
340526959 | 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 78750821 times.
✓ Branch 1 taken 788981 times.
|
151298386 | if (!otherLocalScvf.isDirichlet()) |
| 396 | { | ||
| 397 |
10/10✓ Branch 0 taken 78227937 times.
✓ Branch 1 taken 522884 times.
✓ Branch 2 taken 78227937 times.
✓ Branch 3 taken 522884 times.
✓ Branch 4 taken 78227937 times.
✓ Branch 5 taken 522884 times.
✓ Branch 6 taken 78227937 times.
✓ Branch 7 taken 522884 times.
✓ Branch 8 taken 40153713 times.
✓ Branch 9 taken 262224 times.
|
674104677 | C[faceIdx][otherLocalDofIdx] -= wijk[faceIdx][0][localDir]; |
| 398 |
2/2✓ Branch 0 taken 78227937 times.
✓ Branch 1 taken 522884 times.
|
149752445 | if (!curIsDirichlet) |
| 399 | 669499437 | 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 |
10/10✓ Branch 0 taken 266097 times.
✓ Branch 1 taken 522884 times.
✓ Branch 2 taken 266097 times.
✓ Branch 3 taken 522884 times.
✓ Branch 4 taken 266097 times.
✓ Branch 5 taken 522884 times.
✓ Branch 6 taken 266097 times.
✓ Branch 7 taken 522884 times.
✓ Branch 8 taken 133405 times.
✓ Branch 9 taken 262224 times.
|
6949241 | D[faceIdx][otherLocalDofIdx] -= wijk[faceIdx][0][localDir]; |
| 405 |
2/2✓ Branch 0 taken 266097 times.
✓ Branch 1 taken 522884 times.
|
1545941 | if (!curIsDirichlet) |
| 406 | 2344001 | B[curLocalDofIdx][otherLocalDofIdx] += curXiFactor*wijk[faceIdx][0][localDir]; | |
| 407 | } | ||
| 408 | |||
| 409 | // add entries related to pressures at the scv centers (dofs) | ||
| 410 | 151298386 | const auto posScvLocalDofIdx = posLocalScv.localDofIndex(); | |
| 411 |
10/10✓ Branch 0 taken 78494034 times.
✓ Branch 1 taken 1045768 times.
✓ Branch 2 taken 78494034 times.
✓ Branch 3 taken 1045768 times.
✓ Branch 4 taken 78494034 times.
✓ Branch 5 taken 1045768 times.
✓ Branch 6 taken 78494034 times.
✓ Branch 7 taken 1045768 times.
✓ Branch 8 taken 40287118 times.
✓ Branch 9 taken 524448 times.
|
681053918 | D[faceIdx][posScvLocalDofIdx] += wijk[faceIdx][0][localDir]; |
| 412 | |||
| 413 |
2/2✓ Branch 0 taken 78494034 times.
✓ Branch 1 taken 1045768 times.
|
151298386 | if (!curIsDirichlet) |
| 414 | 671843438 | 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 4500280 times.
✓ Branch 1 taken 35269621 times.
|
75649193 | if (curIsInteriorBoundary && !curIsDirichlet) |
| 419 | { | ||
| 420 |
2/4✓ Branch 1 taken 4500280 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4500280 times.
✗ Branch 5 not taken.
|
17417312 | const auto facetVolVars = this->problem().couplingManager().getLowDimVolVars(posElement, curGlobalScvf); |
| 421 | 8708656 | 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 | 8708656 | const auto wFacet = 2.0*Extrusion::area(this->fvGeometry(), curGlobalScvf)*posVolVars.extrusionFactor() | |
| 426 | 26125968 | *vtmv(curGlobalScvf.unitOuterNormal(), facetTensor, curGlobalScvf.unitOuterNormal()) | |
| 427 | 8708656 | / (dim < dimWorld ? sqrt(facetVolVars.extrusionFactor()) : facetVolVars.extrusionFactor()); | |
| 428 | |||
| 429 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 4500280 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4500280 times.
|
17417312 | A[curLocalDofIdx][curLocalDofIdx] -= wFacet; |
| 430 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 4500280 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4500280 times.
✓ Branch 4 taken 4500280 times.
✗ Branch 5 not taken.
|
17417312 | B[curLocalDofIdx][curLocalScvf.coupledFacetLocalDofIndex()] -= wFacet; |
| 431 | |||
| 432 | // check for zero transmissibilities (skip if inside has been zero already) | ||
| 433 |
3/6✓ Branch 0 taken 4500280 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 4500280 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 4500280 times.
|
8708656 | 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 38734765 times.
✓ Branch 1 taken 1035136 times.
✓ Branch 2 taken 38734765 times.
✗ Branch 3 not taken.
|
75649193 | 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 4500280 times.
✓ Branch 1 taken 34234485 times.
|
73647337 | const Scalar curOneMinusXi = curIsInteriorBoundary ? -(1.0 - xi) : 1.0; |
| 445 | |||
| 446 | // loop over all the outside neighbors of this face and add entries | ||
| 447 |
4/4✓ Branch 0 taken 38734765 times.
✓ Branch 1 taken 38734765 times.
✓ Branch 2 taken 38734765 times.
✓ Branch 3 taken 38734765 times.
|
147294674 | for (unsigned int idxInOutside = 0; idxInOutside < curGlobalScvf.numOutsideScvs(); ++idxInOutside) |
| 448 | { | ||
| 449 | 73647337 | const auto idxOnScvf = idxInOutside+1; | |
| 450 |
4/4✓ Branch 0 taken 9431860 times.
✓ Branch 1 taken 29302905 times.
✓ Branch 2 taken 9431860 times.
✓ Branch 3 taken 29302905 times.
|
147294674 | const auto& negLocalScv = iv.localScv( neighborScvIndices[idxOnScvf] ); |
| 451 |
2/2✓ Branch 0 taken 9431860 times.
✓ Branch 1 taken 29302905 times.
|
73647337 | const auto& negGlobalScv = this->fvGeometry().scv(negLocalScv.gridScvIndex()); |
| 452 | 73647337 | const auto& negVolVars = this->elemVolVars()[negGlobalScv]; | |
| 453 | 73647337 | const auto negTensor = getT(negVolVars); | |
| 454 | |||
| 455 | // On surface grids, use outside face for "negative" transmissibility calculation | ||
| 456 | 73647337 | const auto& scvf = dim < dimWorld ? this->fvGeometry().flipScvf(curGlobalScvf.index(), idxInOutside) | |
| 457 | : curGlobalScvf; | ||
| 458 | 73647337 | 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 | 36712598 | wijk[faceIdx][idxOnScvf] *= -1.0; | |
| 463 | |||
| 464 | // go over the coordinate directions in the negative sub volume | ||
| 465 |
2/2✓ Branch 0 taken 77469530 times.
✓ Branch 1 taken 38734765 times.
|
220942011 | for (int localDir = 0; localDir < dim; localDir++) |
| 466 | { | ||
| 467 |
2/2✓ Branch 0 taken 34234485 times.
✓ Branch 1 taken 43235045 times.
|
147294674 | const auto otherLocalScvfIdx = negLocalScv.localScvfIndex(localDir); |
| 468 |
2/2✓ Branch 0 taken 34234485 times.
✓ Branch 1 taken 43235045 times.
|
147294674 | const auto& otherLocalScvf = iv.localScvf(otherLocalScvfIdx); |
| 469 | 147294674 | const auto otherLocalDofIdx = otherLocalScvf.localDofIndex(); | |
| 470 | |||
| 471 | // check for zero transmissibilities (skip if inside has been zero already) | ||
| 472 |
3/4✓ Branch 0 taken 34234485 times.
✓ Branch 1 taken 43235045 times.
✓ Branch 2 taken 34234485 times.
✗ Branch 3 not taken.
|
147294674 | if (otherLocalDofIdx == curLocalDofIdx && !isZeroWij) |
| 473 |
5/10✗ Branch 0 not taken.
✓ Branch 1 taken 34234485 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 34234485 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 34234485 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 34234485 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 34234485 times.
|
324693405 | if (abs(wijk[faceIdx][idxOnScvf][localDir]) <= wijZeroThresh) |
| 474 | ✗ | faceMarkers.emplace_back( std::make_pair(curLocalDofIdx, faceIdx) ); | |
| 475 | |||
| 476 |
2/2✓ Branch 0 taken 77212743 times.
✓ Branch 1 taken 256787 times.
|
147294674 | if (!otherLocalScvf.isDirichlet()) |
| 477 | 880749882 | A[curLocalDofIdx][otherLocalDofIdx] += curOneMinusXi*wijk[faceIdx][idxOnScvf][localDir]; | |
| 478 | else | ||
| 479 | 3018162 | B[curLocalDofIdx][otherLocalDofIdx] -= curOneMinusXi*wijk[faceIdx][idxOnScvf][localDir]; | |
| 480 | |||
| 481 | // add entries to matrix B | ||
| 482 | 883768044 | B[curLocalDofIdx][negLocalScv.localDofIndex()] += curOneMinusXi*wijk[faceIdx][idxOnScvf][localDir]; | |
| 483 | } | ||
| 484 | } | ||
| 485 | } | ||
| 486 | } | ||
| 487 | } | ||
| 488 | 18010716 | } | |
| 489 | }; | ||
| 490 | |||
| 491 | } // end namespace | ||
| 492 | |||
| 493 | #endif | ||
| 494 |