| 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 CCMpfaDiscretization | ||
| 10 | * \brief A class that contains helper functions as well as functionality | ||
| 11 | * which is common to different mpfa schemes and which solely | ||
| 12 | * operate on the interaction volume. | ||
| 13 | */ | ||
| 14 | #ifndef DUMUX_DISCRETIZATION_CC_MPFA_LOCAL_ASSEMBLER_HELPER_HH | ||
| 15 | #define DUMUX_DISCRETIZATION_CC_MPFA_LOCAL_ASSEMBLER_HELPER_HH | ||
| 16 | |||
| 17 | #include <algorithm> | ||
| 18 | #include <vector> | ||
| 19 | #include <cassert> | ||
| 20 | #include <utility> | ||
| 21 | #include <type_traits> | ||
| 22 | |||
| 23 | #include <dumux/common/typetraits/isvalid.hh> | ||
| 24 | |||
| 25 | namespace Dumux { | ||
| 26 | |||
| 27 | /*! | ||
| 28 | * \ingroup CCMpfaDiscretization | ||
| 29 | * \brief A class that contains helper functions as well as functionality | ||
| 30 | * which is common to different mpfa schemes and which solely | ||
| 31 | * operate on the interaction volume. | ||
| 32 | */ | ||
| 33 | class InteractionVolumeAssemblerHelper | ||
| 34 | { | ||
| 35 | // Helper structs to detect if matrix has resize function | ||
| 36 | struct HasMatrixResize | ||
| 37 | { | ||
| 38 | template<class M> | ||
| 39 | auto operator()(const M& m) -> decltype(std::declval<M>().resize(0, 0)) | ||
| 40 | {} | ||
| 41 | }; | ||
| 42 | |||
| 43 | // Helper structs to detect if vector has resize function | ||
| 44 | struct HasVectorResize | ||
| 45 | { | ||
| 46 | template<class V> | ||
| 47 | auto operator()(const V& v) -> decltype(std::declval<V>().resize(0)) | ||
| 48 | {} | ||
| 49 | }; | ||
| 50 | |||
| 51 | template<class Matrix> | ||
| 52 | static constexpr auto matrixHasResizeFunction() | ||
| 53 | { return decltype( isValid(HasMatrixResize())(std::declval<Matrix>()) )::value; } | ||
| 54 | |||
| 55 | template<class Vector> | ||
| 56 | static constexpr auto vectorHasResizeFunction() | ||
| 57 | { return decltype( isValid(HasVectorResize())(std::declval<Vector>()) )::value; } | ||
| 58 | |||
| 59 | public: | ||
| 60 | /*! | ||
| 61 | * \brief Solves a previously assembled iv-local system of equations | ||
| 62 | * and stores the resulting transmissibilities in the provided | ||
| 63 | * containers within the interaction volume data handle. | ||
| 64 | * | ||
| 65 | * \param fvGeometry The bound element finite volume geometry | ||
| 66 | * \param handle The data handle in which the matrices are stored | ||
| 67 | * \param iv The interaction volume | ||
| 68 | */ | ||
| 69 | template< class FVElementGeometry, class DataHandle, class IV > | ||
| 70 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 55823568 times.
|
103185506 | static void solveLocalSystem(const FVElementGeometry& fvGeometry, |
| 71 | DataHandle& handle, | ||
| 72 | IV& iv) | ||
| 73 | { | ||
| 74 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 55823568 times.
|
101771648 | assert(iv.numUnknowns() > 0); |
| 75 | |||
| 76 | // T = C*(A^-1)*B + D | ||
| 77 | 103185506 | handle.A().invert(); | |
| 78 | 103185506 | handle.CA().rightmultiply(handle.A()); | |
| 79 | 103185506 | handle.T() += multiplyMatrices(handle.CA(), handle.AB()); | |
| 80 | 103185506 | handle.AB().leftmultiply(handle.A()); | |
| 81 | |||
| 82 | // On surface grids, compute the "outside" transmissibilities | ||
| 83 | using GridView = typename IV::Traits::GridView; | ||
| 84 | static constexpr int dim = GridView::dimension; | ||
| 85 | static constexpr int dimWorld = GridView::dimensionworld; | ||
| 86 | if (dim < dimWorld) | ||
| 87 | { | ||
| 88 | // bring outside tij container to the right size | ||
| 89 | 15732300 | auto& tijOut = handle.tijOutside(); | |
| 90 | 15732300 | tijOut.resize(iv.numFaces()); | |
| 91 | using LocalIndexType = typename IV::Traits::IndexSet::LocalIndexType; | ||
| 92 |
2/2✓ Branch 1 taken 34612454 times.
✓ Branch 2 taken 8225492 times.
|
81804272 | for (LocalIndexType fIdx = 0; fIdx < iv.numFaces(); ++fIdx) |
| 93 | { | ||
| 94 |
3/3✓ Branch 1 taken 32436844 times.
✓ Branch 2 taken 861824 times.
✓ Branch 0 taken 1313786 times.
|
66071972 | const auto& curGlobalScvf = fvGeometry.scvf(iv.localScvf(fIdx).gridScvfIndex()); |
| 95 |
2/2✓ Branch 0 taken 32920710 times.
✓ Branch 1 taken 1691744 times.
|
66071972 | const auto numOutsideFaces = curGlobalScvf.boundary() ? 0 : curGlobalScvf.numOutsideScvs(); |
| 96 | // resize each face entry to the right number of outside faces | ||
| 97 | 66071972 | tijOut[fIdx].resize(numOutsideFaces); | |
| 98 | 132143944 | std::for_each(tijOut[fIdx].begin(), | |
| 99 | 66071972 | tijOut[fIdx].end(), | |
| 100 | 63576960 | [&](auto& v) { resizeVector(v, iv.numKnowns()); }); | |
| 101 | } | ||
| 102 | |||
| 103 | // compute outside transmissibilities | ||
| 104 |
4/4✓ Branch 0 taken 34612454 times.
✓ Branch 1 taken 29186482 times.
✓ Branch 2 taken 63798936 times.
✓ Branch 3 taken 8225492 times.
|
137881404 | for (const auto& localFaceData : iv.localFaceData()) |
| 105 | { | ||
| 106 | // continue only for "outside" faces | ||
| 107 |
2/2✓ Branch 0 taken 34612454 times.
✓ Branch 1 taken 29186482 times.
|
122149104 | if (!localFaceData.isOutsideFace()) |
| 108 | 66071972 | continue; | |
| 109 | |||
| 110 | 56077132 | const auto scvIdx = localFaceData.ivLocalInsideScvIndex(); | |
| 111 | 56077132 | const auto scvfIdx = localFaceData.ivLocalScvfIndex(); | |
| 112 | 56077132 | const auto idxInOut = localFaceData.scvfLocalOutsideScvfIndex(); | |
| 113 | |||
| 114 | 56077132 | const auto& wijk = handle.omegas()[scvfIdx][idxInOut+1]; | |
| 115 | 56077132 | auto& tij = tijOut[scvfIdx][idxInOut]; | |
| 116 | |||
| 117 | 56077132 | tij = 0.0; | |
| 118 |
2/2✓ Branch 0 taken 58372964 times.
✓ Branch 1 taken 29186482 times.
|
168231396 | for (unsigned int dir = 0; dir < dim; dir++) |
| 119 | { | ||
| 120 | // the scvf corresponding to this local direction in the scv | ||
| 121 |
3/3✓ Branch 0 taken 55446492 times.
✓ Branch 1 taken 2841992 times.
✓ Branch 2 taken 84480 times.
|
112154264 | const auto& scvf = iv.localScvf(iv.localScv(scvIdx).localScvfIndex(dir)); |
| 122 | |||
| 123 | // on interior faces the coefficients of the AB matrix come into play | ||
| 124 |
2/2✓ Branch 0 taken 58071416 times.
✓ Branch 1 taken 301548 times.
|
112154264 | if (!scvf.isDirichlet()) |
| 125 | { | ||
| 126 |
2/2✓ Branch 1 taken 228658916 times.
✓ Branch 2 taken 55446492 times.
|
670373364 | auto tmp = handle.AB()[scvf.localDofIndex()]; |
| 127 |
2/2✓ Branch 0 taken 241802764 times.
✓ Branch 1 taken 58071416 times.
|
577129372 | tmp *= wijk[dir]; |
| 128 |
1/2✓ Branch 0 taken 58071416 times.
✗ Branch 1 not taken.
|
111637688 | tij -= tmp; |
| 129 | 111637688 | } | |
| 130 | else | ||
| 131 | 516576 | tij[scvf.localDofIndex()] -= wijk[dir]; | |
| 132 | |||
| 133 | // add entry from the scv unknown | ||
| 134 | 112154264 | tij[scvIdx] += wijk[dir]; | |
| 135 | } | ||
| 136 | } | ||
| 137 | } | ||
| 138 | 103185506 | } | |
| 139 | |||
| 140 | /*! | ||
| 141 | * \brief Assembles the vector of face unknowns within an interaction volume. | ||
| 142 | * \note This requires the data handle to be fully assembled already. | ||
| 143 | * | ||
| 144 | * \param handle The data handle in which the vector is stored | ||
| 145 | * \param iv The interaction volume | ||
| 146 | */ | ||
| 147 | template< class DataHandle, class IV > | ||
| 148 | static typename IV::Traits::MatVecTraits::FaceVector | ||
| 149 |
1/2✓ Branch 1 taken 121 times.
✗ Branch 2 not taken.
|
121 | assembleFaceUnkowns(const DataHandle& handle, const IV& iv) |
| 150 | { | ||
| 151 |
1/2✓ Branch 1 taken 121 times.
✗ Branch 2 not taken.
|
121 | typename IV::Traits::MatVecTraits::FaceVector u; |
| 152 | |||
| 153 | // u is a vector, which will contain the values for non-Dirichlet faces (scvfs) | ||
| 154 | // i.e., with the length of number of unknowns in the current iv | ||
| 155 | 242 | resizeVector(u, iv.numUnknowns()); | |
| 156 | |||
| 157 | // handle.AB() is a NxM matrix containing geometrical weights to compute the face unknowns based on the known values | ||
| 158 | // with N: number of unknowns and M: number of knowns in the current iv. | ||
| 159 | // handle.uj() is a vector containing the known values, i.e., with the length of number of knowns in the current iv | ||
| 160 |
1/2✓ Branch 0 taken 121 times.
✗ Branch 1 not taken.
|
121 | handle.AB().mv(handle.uj(), u); |
| 161 | |||
| 162 | // maybe add gravity terms | ||
| 163 |
1/2✓ Branch 0 taken 121 times.
✗ Branch 1 not taken.
|
121 | if (handle.deltaG().size() == iv.numUnknowns()) |
| 164 | 121 | handle.A().umv(handle.deltaG(), u); | |
| 165 | |||
| 166 | 121 | return u; | |
| 167 | ✗ | } | |
| 168 | |||
| 169 | /*! | ||
| 170 | * \brief Assembles the solution gradients in the | ||
| 171 | * sub-control volumes within an interaction volume. | ||
| 172 | * \note This requires the data handle to be fully assembled already. | ||
| 173 | * | ||
| 174 | * \param handle The data handle in which the vector is stored | ||
| 175 | * \param iv The interaction volume | ||
| 176 | */ | ||
| 177 | template< class DataHandle, class IV > | ||
| 178 | static std::vector< typename IV::Traits::LocalScvType::GlobalCoordinate > | ||
| 179 | 121 | assembleScvGradients(const DataHandle& handle, const IV& iv) | |
| 180 | { | ||
| 181 | 121 | const auto u = assembleFaceUnkowns(handle, iv); | |
| 182 | |||
| 183 | using LocalScv = typename IV::Traits::LocalScvType; | ||
| 184 | using Gradient = typename LocalScv::GlobalCoordinate; | ||
| 185 | |||
| 186 |
1/2✓ Branch 1 taken 121 times.
✗ Branch 2 not taken.
|
121 | std::vector<Gradient> result; result.reserve(iv.numScvs()); |
| 187 |
2/2✓ Branch 0 taken 400 times.
✓ Branch 1 taken 121 times.
|
521 | for (unsigned int scvIdx = 0; scvIdx < iv.numScvs(); ++scvIdx) |
| 188 | { | ||
| 189 | 400 | const auto& scv = iv.localScv(scvIdx); | |
| 190 | |||
| 191 | 400 | Gradient gradU(0.0); | |
| 192 |
2/2✓ Branch 0 taken 800 times.
✓ Branch 1 taken 400 times.
|
1200 | for (unsigned int dir = 0; dir < LocalScv::myDimension; ++dir) |
| 193 | { | ||
| 194 | 800 | auto nu = scv.nu(dir); | |
| 195 | |||
| 196 | // obtain face pressure | ||
| 197 |
2/2✓ Branch 1 taken 760 times.
✓ Branch 2 taken 40 times.
|
800 | const auto& scvf = iv.localScvf( scv.localScvfIndex(dir) ); |
| 198 |
2/2✓ Branch 0 taken 760 times.
✓ Branch 1 taken 40 times.
|
800 | const auto faceU = !scvf.isDirichlet() ? u[scvf.localDofIndex()] |
| 199 | 40 | : handle.uj()[scvf.localDofIndex()]; | |
| 200 | |||
| 201 |
2/2✓ Branch 0 taken 1600 times.
✓ Branch 1 taken 800 times.
|
3200 | nu *= faceU - handle.uj()[scv.localDofIndex()]; |
| 202 | 800 | gradU += nu; | |
| 203 | } | ||
| 204 | |||
| 205 | 400 | gradU /= scv.detX(); | |
| 206 |
1/2✓ Branch 1 taken 400 times.
✗ Branch 2 not taken.
|
400 | result.emplace_back( std::move(gradU) ); |
| 207 | } | ||
| 208 | |||
| 209 | 121 | return result; | |
| 210 | 121 | } | |
| 211 | |||
| 212 | //! resizes a matrix to the given sizes (specialization for dynamic matrix type) | ||
| 213 | template< class Matrix, | ||
| 214 | class size_type, | ||
| 215 | std::enable_if_t<matrixHasResizeFunction<Matrix>(), int> = 0 > | ||
| 216 | 223294464 | static void resizeMatrix(Matrix& M, size_type rows, size_type cols) | |
| 217 | { | ||
| 218 |
13/30✓ Branch 1 taken 192 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 44135962 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 44135962 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 44135962 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 44135962 times.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✓ Branch 19 taken 7698734 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 7698734 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 7698734 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 7698734 times.
✗ Branch 29 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✓ Branch 34 taken 3988872 times.
✗ Branch 35 not taken.
✓ Branch 37 taken 3988872 times.
✗ Branch 38 not taken.
✓ Branch 40 taken 3988872 times.
✗ Branch 41 not taken.
✓ Branch 43 taken 3988872 times.
✗ Branch 44 not taken.
|
223294464 | M.resize(rows, cols); |
| 219 | 223294464 | } | |
| 220 | |||
| 221 | //! resizes a matrix to the given sizes (specialization for static matrix type - do nothing) | ||
| 222 | template< class Matrix, | ||
| 223 | class size_type, | ||
| 224 | std::enable_if_t<!matrixHasResizeFunction<Matrix>(), int> = 0 > | ||
| 225 | static void resizeMatrix(Matrix& M, size_type rows, size_type cols) | ||
| 226 | {} | ||
| 227 | |||
| 228 | //! resizes a vector to the given size (specialization for dynamic matrix type) | ||
| 229 | template< class Vector, | ||
| 230 | class size_type, | ||
| 231 | std::enable_if_t<vectorHasResizeFunction<Vector>(), int> = 0 > | ||
| 232 |
2/8✓ Branch 6 taken 13729184 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 1 taken 121 times.
✗ Branch 2 not taken.
|
390268096 | static void resizeVector(Vector& v, size_type size) |
| 233 | { | ||
| 234 |
13/25✓ Branch 1 taken 41536912 times.
✓ Branch 2 taken 2123242 times.
✓ Branch 4 taken 2407554 times.
✓ Branch 5 taken 2744220 times.
✓ Branch 7 taken 13744654 times.
✓ Branch 8 taken 826560 times.
✓ Branch 10 taken 2683520 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2143706 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 11242416 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 2139768 times.
✗ Branch 20 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✓ Branch 25 taken 8974686 times.
✗ Branch 26 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 21 not taken.
✓ Branch 12 taken 261828 times.
✓ Branch 3 taken 82088762 times.
✗ Branch 6 not taken.
✗ Branch 9 not taken.
|
563049748 | v.resize(size); |
| 235 | 172781268 | } | |
| 236 | |||
| 237 | //! resizes a vector to the given size (specialization for static vector type - do nothing) | ||
| 238 | template< class Vector, | ||
| 239 | class size_type, | ||
| 240 | std::enable_if_t<!vectorHasResizeFunction<Vector>(), int> = 0 > | ||
| 241 | static void resizeVector(Vector& v, size_type rows) | ||
| 242 | {} | ||
| 243 | }; | ||
| 244 | |||
| 245 | } // end namespace Dumux | ||
| 246 | |||
| 247 | #endif | ||
| 248 |