GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/discretization/cellcentered/mpfa/localassemblerhelper.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 57 61 93.4%
Functions: 50 58 86.2%
Branches: 81 128 63.3%
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 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 96460466 static void solveLocalSystem(const FVElementGeometry& fvGeometry,
71 DataHandle& handle,
72 IV& iv)
73 {
74
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 52461048 times.
 95046608 assert(iv.numUnknowns() > 0);
75
76 // T = C*(A^-1)*B + D
77 192920932 handle.A().invert();
78 289381398 handle.CA().rightmultiply(handle.A());
79 289381398 handle.T() += multiplyMatrices(handle.CA(), handle.AB());
80 289381398 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 9483420 auto& tijOut = handle.tijOutside();
90 9483420 tijOut.resize(iv.numFaces());
91 using LocalIndexType = typename IV::Traits::IndexSet::LocalIndexType;
92
2/2
✓ Branch 1 taken 21507824 times.
✓ Branch 2 taken 5101052 times.
 49346132 for (LocalIndexType fIdx = 0; fIdx < iv.numFaces(); ++fIdx)
93 {
94
4/4
✓ Branch 0 taken 1313786 times.
✓ Branch 1 taken 829920 times.
✓ Branch 2 taken 1313786 times.
✓ Branch 3 taken 829920 times.
 79725424 const auto& curGlobalScvf = fvGeometry.scvf(iv.localScvf(fIdx).gridScvfIndex());
95
2/2
✓ Branch 0 taken 20165880 times.
✓ Branch 1 taken 1341944 times.
 39862712 const auto numOutsideFaces = curGlobalScvf.boundary() ? 0 : curGlobalScvf.numOutsideScvs();
96 // resize each face entry to the right number of outside faces
97 79725424 tijOut[fIdx].resize(numOutsideFaces);
98 119588136 std::for_each(tijOut[fIdx].begin(),
99 79725424 tijOut[fIdx].end(),
100 76134600 [&](auto& v) { resizeVector(v, iv.numKnowns()); });
101 }
102
103 // compute outside transmissibilities
104
4/4
✓ Branch 0 taken 39509616 times.
✓ Branch 1 taken 5101052 times.
✓ Branch 2 taken 39509616 times.
✓ Branch 3 taken 5101052 times.
 111504144 for (const auto& localFaceData : iv.localFaceData())
105 {
106 // continue only for "outside" faces
107
2/2
✓ Branch 0 taken 18001792 times.
✓ Branch 1 taken 21507824 times.
 73570464 if (!localFaceData.isOutsideFace())
108 continue;
109
110 33707752 const auto scvIdx = localFaceData.ivLocalInsideScvIndex();
111 33707752 const auto scvfIdx = localFaceData.ivLocalScvfIndex();
112 33707752 const auto idxInOut = localFaceData.scvfLocalOutsideScvfIndex();
113
114 101123256 const auto& wijk = handle.omegas()[scvfIdx][idxInOut+1];
115 67415504 auto& tij = tijOut[scvfIdx][idxInOut];
116
117 33707752 tij = 0.0;
118
2/2
✓ Branch 0 taken 36003584 times.
✓ Branch 1 taken 18001792 times.
 101123256 for (unsigned int dir = 0; dir < dim; dir++)
119 {
120 // the scvf corresponding to this local direction in the scv
121
6/6
✓ Branch 0 taken 33166212 times.
✓ Branch 1 taken 127968 times.
✓ Branch 2 taken 35791136 times.
✓ Branch 3 taken 212448 times.
✓ Branch 4 taken 33166212 times.
✓ Branch 5 taken 127968 times.
 199537108 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 35791136 times.
✓ Branch 1 taken 212448 times.
 67415504 if (!scvf.isDirichlet())
125 {
126 268308512 auto tmp = handle.AB()[scvf.localDofIndex()];
127 5249848 tmp *= wijk[dir];
128
1/2
✓ Branch 0 taken 35791136 times.
✗ Branch 1 not taken.
 67077128 tij -= tmp;
129 }
130 else
131 1015128 tij[scvf.localDofIndex()] -= wijk[dir];
132
133 // add entry from the scv unknown
134 202246512 tij[scvIdx] += wijk[dir];
135 }
136 }
137 }
138 96460466 }
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 121 assembleFaceUnkowns(const DataHandle& handle, const IV& iv)
150 {
151
3/6
✓ Branch 1 taken 121 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 121 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 121 times.
✗ Branch 8 not taken.
 363 typename IV::Traits::MatVecTraits::FaceVector u;
152 242 resizeVector(u, iv.numFaces());
153
154 242 handle.AB().mv(handle.uj(), u);
155
156 // maybe add gravity terms
157
3/6
✓ Branch 0 taken 121 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 121 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 121 times.
✗ Branch 5 not taken.
 363 if (handle.deltaG().size() == iv.numUnknowns())
158 363 handle.A().umv(handle.deltaG(), u);
159
160 121 return u;
161 }
162
163 /*!
164 * \brief Assembles the solution gradients in the
165 * sub-control volumes within an interaction volume.
166 * \note This requires the data handle to be fully assembled already.
167 *
168 * \param handle The data handle in which the vector is stored
169 * \param iv The interaction volume
170 */
171 template< class DataHandle, class IV >
172 static std::vector< typename IV::Traits::LocalScvType::GlobalCoordinate >
173 121 assembleScvGradients(const DataHandle& handle, const IV& iv)
174 {
175
1/4
✓ Branch 1 taken 121 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 242 const auto u = assembleFaceUnkowns(handle, iv);
176
177 using LocalScv = typename IV::Traits::LocalScvType;
178 using Gradient = typename LocalScv::GlobalCoordinate;
179
180
3/6
✓ Branch 1 taken 121 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 121 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 121 times.
✗ Branch 8 not taken.
 363 std::vector<Gradient> result; result.reserve(iv.numScvs());
181
4/4
✓ Branch 0 taken 400 times.
✓ Branch 1 taken 121 times.
✓ Branch 2 taken 400 times.
✓ Branch 3 taken 121 times.
 642 for (unsigned int scvIdx = 0; scvIdx < iv.numScvs(); ++scvIdx)
182 {
183 400 const auto& scv = iv.localScv(scvIdx);
184
185 400 Gradient gradU(0.0);
186
2/2
✓ Branch 0 taken 800 times.
✓ Branch 1 taken 400 times.
 1200 for (unsigned int dir = 0; dir < LocalScv::myDimension; ++dir)
187 {
188 800 auto nu = scv.nu(dir);
189
190 // obtain face pressure
191
2/2
✓ Branch 1 taken 760 times.
✓ Branch 2 taken 40 times.
 800 const auto& scvf = iv.localScvf( scv.localScvfIndex(dir) );
192
2/2
✓ Branch 0 taken 760 times.
✓ Branch 1 taken 40 times.
 800 const auto faceU = !scvf.isDirichlet() ? u[scvf.localDofIndex()]
193 120 : handle.uj()[scvf.localDofIndex()];
194
195 2400 nu *= faceU - handle.uj()[scv.localDofIndex()];
196 800 gradU += nu;
197 }
198
199 400 gradU /= scv.detX();
200
1/2
✓ Branch 1 taken 400 times.
✗ Branch 2 not taken.
 400 result.emplace_back( std::move(gradU) );
201 }
202
203 121 return result;
204 }
205
206 //! resizes a matrix to the given sizes (specialization for dynamic matrix type)
207 template< class Matrix,
208 class size_type,
209 std::enable_if_t<matrixHasResizeFunction<Matrix>(), int> = 0 >
210 static void resizeMatrix(Matrix& M, size_type rows, size_type cols)
211 {
212
13/30
✓ Branch 1 taken 192 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 43015122 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 43015122 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 43015122 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 43015122 times.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✓ Branch 19 taken 6577894 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 6577894 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 6577894 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 6577894 times.
✗ Branch 29 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✓ Branch 34 taken 2868032 times.
✗ Branch 35 not taken.
✓ Branch 37 taken 2868032 times.
✗ Branch 38 not taken.
✓ Branch 40 taken 2868032 times.
✗ Branch 41 not taken.
✓ Branch 43 taken 2868032 times.
✗ Branch 44 not taken.
 209844384 M.resize(rows, cols);
213 }
214
215 //! resizes a matrix to the given sizes (specialization for static matrix type - do nothing)
216 template< class Matrix,
217 class size_type,
218 std::enable_if_t<!matrixHasResizeFunction<Matrix>(), int> = 0 >
219 ✗ static void resizeMatrix(Matrix& M, size_type rows, size_type cols)
220 ✗ {}
221
222 //! resizes a vector to the given size (specialization for dynamic matrix type)
223 template< class Vector,
224 class size_type,
225 std::enable_if_t<vectorHasResizeFunction<Vector>(), int> = 0 >
226 static void resizeVector(Vector& v, size_type size)
227 {
228
21/36
✓ Branch 1 taken 2209817 times.
✓ Branch 2 taken 36262452 times.
✓ Branch 3 taken 248220 times.
✓ Branch 4 taken 543752 times.
✓ Branch 5 taken 3221530 times.
✓ Branch 6 taken 53664 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 1580992 times.
✓ Branch 9 taken 83097992 times.
✓ Branch 10 taken 3014930 times.
✓ Branch 11 taken 4606478 times.
✓ Branch 12 taken 136558 times.
✓ Branch 13 taken 3051294 times.
✓ Branch 14 taken 543752 times.
✓ Branch 15 taken 2120832 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 26832 times.
✓ Branch 18 taken 261828 times.
✓ Branch 19 taken 2143706 times.
✓ Branch 20 taken 2267730 times.
✓ Branch 21 taken 5421956 times.
✗ Branch 22 not taken.
✓ Branch 25 taken 1098288 times.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✓ Branch 31 taken 4606476 times.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
✗ Branch 45 not taken.
✗ Branch 46 not taken.
 718622657 v.resize(size);
229 }
230
231 //! resizes a vector to the given size (specialization for static vector type - do nothing)
232 template< class Vector,
233 class size_type,
234 std::enable_if_t<!vectorHasResizeFunction<Vector>(), int> = 0 >
235 ✗ static void resizeVector(Vector& v, size_type rows)
236 ✗ {}
237 };
238
239 } // end namespace Dumux
240
241 #endif
242