| 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 CCMpfaDiscretization | ||
| 10 | * \brief Defines the general interface of classes used for the assembly | ||
| 11 | * of the local systems of equations involved in the transmissibility | ||
| 12 | * computation in mpfa schemes. | ||
| 13 | */ | ||
| 14 | #ifndef DUMUX_DISCRETIZATION_CC_MPFA_LOCAL_ASSEMBLER_BASE_HH | ||
| 15 | #define DUMUX_DISCRETIZATION_CC_MPFA_LOCAL_ASSEMBLER_BASE_HH | ||
| 16 | |||
| 17 | #include <algorithm> | ||
| 18 | #include <vector> | ||
| 19 | #include <type_traits> | ||
| 20 | |||
| 21 | #include <dune/common/exceptions.hh> | ||
| 22 | #include <dune/common/reservedvector.hh> | ||
| 23 | |||
| 24 | #include <dumux/common/typetraits/isvalid.hh> | ||
| 25 | #include <dumux/discretization/extrusion.hh> | ||
| 26 | #include "localassemblerhelper.hh" | ||
| 27 | |||
| 28 | namespace Dumux { | ||
| 29 | |||
| 30 | /*! | ||
| 31 | * \ingroup CCMpfaDiscretization | ||
| 32 | * \brief Defines the general interface of the local assembler | ||
| 33 | * classes for the assembly of the interaction volume-local | ||
| 34 | * transmissibility matrix. Specializations have to be provided | ||
| 35 | * for the available interaction volume implementations. these | ||
| 36 | * should derive from this base class. | ||
| 37 | * | ||
| 38 | * \tparam P The problem type | ||
| 39 | * \tparam EG The element finite volume geometry | ||
| 40 | * \tparam EV The element volume variables type | ||
| 41 | */ | ||
| 42 | template< class P, class EG, class EV > | ||
| 43 | class InteractionVolumeAssemblerBase | ||
| 44 | : public InteractionVolumeAssemblerHelper | ||
| 45 | { | ||
| 46 | using Problem = P; | ||
| 47 | using FVElementGeometry = EG; | ||
| 48 | using ElementVolumeVariables = EV; | ||
| 49 | using Extrusion = Extrusion_t<typename EG::GridGeometry>; | ||
| 50 | using Helper = InteractionVolumeAssemblerHelper; | ||
| 51 | |||
| 52 | template< class IV > | ||
| 53 | using Scalar = typename IV::Traits::MatVecTraits::FaceVector::value_type; | ||
| 54 | |||
| 55 | public: | ||
| 56 | /*! | ||
| 57 | * \brief The constructor. | ||
| 58 | * Sets pointers to the objects required for a subsequent call to assemble(). | ||
| 59 | * | ||
| 60 | * \param problem The problem to be solved (boundary/initial conditions etc.) | ||
| 61 | * \param fvGeometry The local view on the finite volume grid geometry | ||
| 62 | * \param elemVolVars The local view on the primary/secondary variables | ||
| 63 | */ | ||
| 64 | 76855392 | InteractionVolumeAssemblerBase(const Problem& problem, | |
| 65 | const FVElementGeometry& fvGeometry, | ||
| 66 | const ElementVolumeVariables& elemVolVars) | ||
| 67 | { | ||
| 68 | 76855392 | problemPtr_ = &problem; | |
| 69 | 76855392 | fvGeometryPtr_ = &fvGeometry; | |
| 70 |
5/5✓ Branch 0 taken 398496 times.
✓ Branch 1 taken 48547560 times.
✓ Branch 2 taken 5272536 times.
✓ Branch 3 taken 12065201 times.
✓ Branch 4 taken 53662 times.
|
76855392 | elemVolVarsPtr_ = &elemVolVars; |
| 71 | } | ||
| 72 | |||
| 73 | // return functions to the local views & problem | ||
| 74 | ✗ | const Problem& problem() const { return *problemPtr_; } | |
| 75 | ✗ | const FVElementGeometry& fvGeometry() const { return *fvGeometryPtr_; } | |
| 76 | ✗ | const ElementVolumeVariables& elemVolVars() const { return *elemVolVarsPtr_; } | |
| 77 | |||
| 78 | /*! | ||
| 79 | * \brief Assembles the matrices involved in the flux | ||
| 80 | * expressions and the local system of equations | ||
| 81 | * within an mpfa interaction volume. | ||
| 82 | * | ||
| 83 | * \tparam DataHandle The data handle | ||
| 84 | * \tparam IV The interaction volume type implementation | ||
| 85 | * \tparam TensorFunc Lambda to obtain the tensor w.r.t. | ||
| 86 | * which the local system is to be solved | ||
| 87 | * | ||
| 88 | * \param handle The data handle in which the matrices are stored | ||
| 89 | * \param iv The interaction volume | ||
| 90 | * \param getT Lambda to evaluate the scv-wise tensors | ||
| 91 | * \param wijZeroThresh the zero threshold wij | ||
| 92 | */ | ||
| 93 | template< class DataHandle, class IV, class TensorFunc > | ||
| 94 | void assembleMatrices(DataHandle& handle, IV& iv, const TensorFunc& getT, Scalar<IV> wijZeroThresh = 0.0) | ||
| 95 | { | ||
| 96 | DUNE_THROW(Dune::NotImplemented, "Implementation does not provide an assembleMatrices() function"); | ||
| 97 | } | ||
| 98 | |||
| 99 | /*! | ||
| 100 | * \brief Assembles the vector of primary (cell) unknowns and (maybe) | ||
| 101 | * Dirichlet boundary conditions within an interaction volume. | ||
| 102 | * | ||
| 103 | * \tparam IV The interaction volume type implementation | ||
| 104 | * \tparam GetU Lambda to obtain the cell unknowns from grid indices | ||
| 105 | * | ||
| 106 | * \param handle The data handle in which the vector is stored | ||
| 107 | * \param iv The interaction volume | ||
| 108 | * \param getU Lambda to obtain the desired cell/Dirichlet value from vol vars | ||
| 109 | */ | ||
| 110 | template< class DataHandle, class IV, class GetU > | ||
| 111 | void assembleU(DataHandle& handle, const IV& iv, const GetU& getU) | ||
| 112 | { | ||
| 113 | DUNE_THROW(Dune::NotImplemented, "Implementation does not provide an assemble() function for the cell/Dirichlet unknowns"); | ||
| 114 | } | ||
| 115 | |||
| 116 | /*! | ||
| 117 | * \brief Assembles the gravitational flux contributions on the scvfs within an | ||
| 118 | * interaction volume. | ||
| 119 | * | ||
| 120 | * \param handle The data handle in which the vector is stored | ||
| 121 | * \param iv The interaction volume | ||
| 122 | * \param getRho Lambda to obtain the density from volume variables | ||
| 123 | */ | ||
| 124 | template< class DataHandle, class IV, class GetRho > | ||
| 125 | 47073922 | void assembleGravity(DataHandle& handle, const IV& iv, const GetRho& getRho) | |
| 126 | { | ||
| 127 | using GridView = typename IV::Traits::GridView; | ||
| 128 | static constexpr int dim = GridView::dimension; | ||
| 129 | static constexpr int dimWorld = GridView::dimensionworld; | ||
| 130 | static constexpr bool isSurfaceGrid = dim < dimWorld; | ||
| 131 | |||
| 132 | // resize the gravity vectors | ||
| 133 | 47073922 | auto& g = handle.g(); | |
| 134 | 47073922 | auto& deltaG = handle.deltaG(); | |
| 135 | 47073922 | auto& outsideG = handle.gOutside(); | |
| 136 | 47073922 | Helper::resizeVector(g, iv.numFaces()); | |
| 137 | 47073922 | Helper::resizeVector(deltaG, iv.numUnknowns()); | |
| 138 | if (isSurfaceGrid) | ||
| 139 | 13660904 | Helper::resizeVector(outsideG, iv.numFaces()); | |
| 140 | |||
| 141 | //! For each face, we... | ||
| 142 | //! - arithmetically average the phase densities | ||
| 143 | //! - compute the term \f$ \alpha := \mathbf{A} \rho \ \mathbf{n}^T \mathbf{K} \mathbf{g} \f$ in each neighboring cell | ||
| 144 | //! - compute \f$ \alpha^* = \sum{\alpha_{outside, i}} - \alpha_{inside} \f$ | ||
| 145 | using Scalar = typename IV::Traits::MatVecTraits::TMatrix::value_type; | ||
| 146 | using LocalIndexType = typename IV::Traits::IndexSet::LocalIndexType; | ||
| 147 | |||
| 148 |
3/3✓ Branch 0 taken 131235760 times.
✓ Branch 1 taken 114971294 times.
✓ Branch 2 taken 13660904 times.
|
262437668 | for (LocalIndexType faceIdx = 0; faceIdx < iv.numFaces(); ++faceIdx) |
| 149 | { | ||
| 150 | // gravitational acceleration on this face | ||
| 151 | 215363746 | const auto& curLocalScvf = iv.localScvf(faceIdx); | |
| 152 | 215363746 | const auto& curGlobalScvf = fvGeometry().scvf(curLocalScvf.gridScvfIndex()); | |
| 153 |
6/6✓ Branch 0 taken 3766636 times.
✓ Branch 1 taken 10754072 times.
✓ Branch 2 taken 3766636 times.
✓ Branch 3 taken 10754072 times.
✓ Branch 4 taken 3766636 times.
✓ Branch 5 taken 10754072 times.
|
646091238 | const auto& gravity = problem().spatialParams().gravity(curGlobalScvf.ipGlobal()); |
| 154 | |||
| 155 | // get permeability tensor in "positive" sub volume | ||
| 156 | 215363746 | const auto& neighborScvIndices = curLocalScvf.neighboringLocalScvIndices(); | |
| 157 |
6/6✓ Branch 0 taken 3766636 times.
✓ Branch 1 taken 10754072 times.
✓ Branch 2 taken 3766636 times.
✓ Branch 3 taken 10754072 times.
✓ Branch 4 taken 3766636 times.
✓ Branch 5 taken 10754072 times.
|
564008012 | const auto& posGlobalScv = fvGeometry().scv(iv.localScv(neighborScvIndices[0]).gridScvIndex()); |
| 158 | 215363746 | const auto& posVolVars = elemVolVars()[posGlobalScv]; | |
| 159 | 430727492 | const auto alpha_inside = posVolVars.extrusionFactor()*vtmv(curGlobalScvf.unitOuterNormal(), | |
| 160 | posVolVars.permeability(), | ||
| 161 | gravity); | ||
| 162 | |||
| 163 |
2/2✓ Branch 0 taken 208155190 times.
✓ Branch 1 taken 5163796 times.
|
215363746 | const auto numOutsideFaces = !curGlobalScvf.boundary() ? curGlobalScvf.numOutsideScvs() : 0; |
| 164 | using OutsideAlphaStorage = std::conditional_t< isSurfaceGrid, | ||
| 165 | std::vector<Scalar>, | ||
| 166 | Dune::ReservedVector<Scalar, 1> >; | ||
| 167 |
4/8✓ Branch 1 taken 82083226 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 82083226 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 79292974 times.
✓ Branch 7 taken 2790252 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
458823172 | OutsideAlphaStorage alpha_outside; alpha_outside.resize(numOutsideFaces); |
| 168 | 646091238 | std::fill(alpha_outside.begin(), alpha_outside.end(), 0.0); | |
| 169 | Scalar rho; | ||
| 170 | |||
| 171 | if (isSurfaceGrid) | ||
| 172 |
2/4✓ Branch 1 taken 82083226 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 82083226 times.
✗ Branch 5 not taken.
|
164166452 | Helper::resizeVector(outsideG[faceIdx], numOutsideFaces); |
| 173 | |||
| 174 |
2/2✓ Branch 0 taken 212397602 times.
✓ Branch 1 taken 921384 times.
|
215363746 | if (!curLocalScvf.isDirichlet()) |
| 175 | { | ||
| 176 | 214421034 | const auto localDofIdx = curLocalScvf.localDofIndex(); | |
| 177 | |||
| 178 |
2/2✓ Branch 0 taken 208155190 times.
✓ Branch 1 taken 4242412 times.
|
214421034 | rho = getRho(posVolVars); |
| 179 |
2/2✓ Branch 0 taken 208155190 times.
✓ Branch 1 taken 4242412 times.
|
214421034 | deltaG[localDofIdx] = 0.0; |
| 180 | |||
| 181 |
2/2✓ Branch 0 taken 208155190 times.
✓ Branch 1 taken 4242412 times.
|
214421034 | if (!curGlobalScvf.boundary()) |
| 182 | { | ||
| 183 |
4/4✓ Branch 0 taken 209412694 times.
✓ Branch 1 taken 208155190 times.
✓ Branch 2 taken 209412694 times.
✓ Branch 3 taken 208155190 times.
|
842885496 | for (unsigned int idxInOutside = 0; idxInOutside < curGlobalScvf.numOutsideScvs(); ++idxInOutside) |
| 184 | { | ||
| 185 | // obtain outside tensor | ||
| 186 |
2/2✓ Branch 0 taken 3609148 times.
✓ Branch 1 taken 10593856 times.
|
211350126 | const auto negLocalScvIdx = neighborScvIndices[idxInOutside+1]; |
| 187 |
4/4✓ Branch 0 taken 3609148 times.
✓ Branch 1 taken 10593856 times.
✓ Branch 2 taken 3609148 times.
✓ Branch 3 taken 10593856 times.
|
422700252 | const auto& negGlobalScv = fvGeometry().scv(iv.localScv(negLocalScvIdx).gridScvIndex()); |
| 188 | 211350126 | const auto& negVolVars = elemVolVars()[negGlobalScv]; | |
| 189 | 291900604 | const auto& flipScvf = !isSurfaceGrid ? curGlobalScvf | |
| 190 | 80550478 | : fvGeometry().flipScvf(curGlobalScvf.index(), idxInOutside); | |
| 191 | |||
| 192 | 634050378 | alpha_outside[idxInOutside] = negVolVars.extrusionFactor()*vtmv(flipScvf.unitOuterNormal(), | |
| 193 | negVolVars.permeability(), | ||
| 194 | gravity); | ||
| 195 | if (isSurfaceGrid) | ||
| 196 | 80550478 | alpha_outside[idxInOutside] *= -1.0; | |
| 197 | |||
| 198 | 211350126 | rho += getRho(negVolVars); | |
| 199 | 553499900 | deltaG[localDofIdx] += alpha_outside[idxInOutside]; | |
| 200 | } | ||
| 201 | } | ||
| 202 | |||
| 203 | 214421034 | rho /= numOutsideFaces + 1; | |
| 204 | 214421034 | deltaG[localDofIdx] -= alpha_inside; | |
| 205 | 643263102 | deltaG[localDofIdx] *= rho*Extrusion::area(fvGeometry(), curGlobalScvf); | |
| 206 | } | ||
| 207 | // use density resulting from Dirichlet BCs | ||
| 208 | else | ||
| 209 | 1885424 | rho = getRho(elemVolVars()[curGlobalScvf.outsideScvIdx()]); | |
| 210 | |||
| 211 | // add "inside" & "outside" alphas to gravity containers | ||
| 212 | 564008012 | g[faceIdx] = alpha_inside*rho*Extrusion::area(fvGeometry(), curGlobalScvf); | |
| 213 | |||
| 214 | if (isSurfaceGrid) | ||
| 215 | { | ||
| 216 | 82083226 | unsigned int i = 0; | |
| 217 |
4/4✓ Branch 0 taken 80550478 times.
✓ Branch 1 taken 82083226 times.
✓ Branch 2 taken 80550478 times.
✓ Branch 3 taken 82083226 times.
|
407350634 | for (const auto& alpha : alpha_outside) |
| 218 | 322201912 | outsideG[faceIdx][i++] = alpha*rho*Extrusion::area(fvGeometry(), curGlobalScvf); | |
| 219 | } | ||
| 220 | } | ||
| 221 | |||
| 222 | // add iv-wide contributions to gravity vectors | ||
| 223 | 94147844 | handle.CA().umv(deltaG, g); | |
| 224 | if (isSurfaceGrid) | ||
| 225 | { | ||
| 226 | using FaceVector = typename IV::Traits::MatVecTraits::FaceVector; | ||
| 227 |
4/10✓ Branch 1 taken 13660904 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13660904 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 13660904 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 13660904 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
|
54643616 | FaceVector AG; |
| 228 |
2/4✓ Branch 1 taken 13660904 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13660904 times.
✗ Branch 5 not taken.
|
27321808 | Helper::resizeVector(AG, iv.numUnknowns()); |
| 229 | 27321808 | handle.A().mv(deltaG, AG); | |
| 230 | |||
| 231 | // compute gravitational accelerations | ||
| 232 |
4/4✓ Branch 0 taken 162633704 times.
✓ Branch 1 taken 13660904 times.
✓ Branch 2 taken 162633704 times.
✓ Branch 3 taken 13660904 times.
|
217277320 | for (const auto& localFaceData : iv.localFaceData()) |
| 233 | { | ||
| 234 | // continue only for "outside" faces | ||
| 235 |
2/2✓ Branch 0 taken 80550478 times.
✓ Branch 1 taken 82083226 times.
|
162633704 | if (!localFaceData.isOutsideFace()) |
| 236 | continue; | ||
| 237 | |||
| 238 | 80550478 | const auto localScvIdx = localFaceData.ivLocalInsideScvIndex(); | |
| 239 | 80550478 | const auto localScvfIdx = localFaceData.ivLocalScvfIndex(); | |
| 240 | 80550478 | const auto idxInOutside = localFaceData.scvfLocalOutsideScvfIndex(); | |
| 241 | 161100956 | const auto& posLocalScv = iv.localScv(localScvIdx); | |
| 242 | 241651434 | const auto& wijk = handle.omegas()[localScvfIdx][idxInOutside+1]; | |
| 243 | |||
| 244 | // add contributions from all local directions | ||
| 245 |
2/2✓ Branch 0 taken 161100956 times.
✓ Branch 1 taken 80550478 times.
|
241651434 | for (LocalIndexType localDir = 0; localDir < dim; localDir++) |
| 246 | { | ||
| 247 | // the scvf corresponding to this local direction in the scv | ||
| 248 |
2/2✓ Branch 1 taken 160815176 times.
✓ Branch 2 taken 285780 times.
|
161100956 | const auto& curLocalScvf = iv.localScvf(posLocalScv.localScvfIndex(localDir)); |
| 249 |
2/2✓ Branch 0 taken 160815176 times.
✓ Branch 1 taken 285780 times.
|
161100956 | if (!curLocalScvf.isDirichlet()) |
| 250 | 804075880 | outsideG[localScvfIdx][idxInOutside] -= wijk[localDir]*AG[curLocalScvf.localDofIndex()]; | |
| 251 | } | ||
| 252 | } | ||
| 253 | } | ||
| 254 | 47073922 | } | |
| 255 | |||
| 256 | private: | ||
| 257 | // pointers to the data required for assembly | ||
| 258 | const Problem* problemPtr_; | ||
| 259 | const FVElementGeometry* fvGeometryPtr_; | ||
| 260 | const ElementVolumeVariables* elemVolVarsPtr_; | ||
| 261 | }; | ||
| 262 | |||
| 263 | } // end namespace Dumux | ||
| 264 | |||
| 265 | #endif | ||
| 266 |