GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/assembly/jacobianpattern.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 85 85 100.0%
Functions: 93 97 95.9%
Branches: 118 170 69.4%
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 Assembly
10 * \brief Helper function to generate Jacobian pattern for different discretization methods
11 */
12 #ifndef DUMUX_JACOBIAN_PATTERN_HH
13 #define DUMUX_JACOBIAN_PATTERN_HH
14
15 #include <type_traits>
16 #include <dune/istl/matrixindexset.hh>
17 #include <dumux/discretization/method.hh>
18
19 namespace Dumux {
20
21 /*!
22 * \ingroup Assembly
23 * \brief Helper function to generate Jacobian pattern for cell-centered methods
24 */
25 template<bool isImplicit, class GridGeometry,
26 typename std::enable_if_t<( (GridGeometry::discMethod == DiscretizationMethods::cctpfa)
27 || (GridGeometry::discMethod == DiscretizationMethods::ccmpfa) ), int> = 0>
28 442 Dune::MatrixIndexSet getJacobianPattern(const GridGeometry& gridGeometry)
29 {
30
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 442 const auto numDofs = gridGeometry.numDofs();
31
1/2
✓ Branch 1 taken 282 times.
✗ Branch 2 not taken.
 442 Dune::MatrixIndexSet pattern;
32
1/2
✓ Branch 1 taken 282 times.
✗ Branch 2 not taken.
 442 pattern.resize(numDofs, numDofs);
33
34 // matrix pattern for implicit Jacobians
35 if (isImplicit)
36 {
37
2/2
✓ Branch 0 taken 1039913 times.
✓ Branch 1 taken 280 times.
 2359889 for (unsigned int globalI = 0; globalI < numDofs; ++globalI)
38 {
39
2/4
✓ Branch 1 taken 1039913 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1039913 times.
✗ Branch 5 not taken.
 4718914 pattern.add(globalI, globalI);
40
4/4
✓ Branch 0 taken 4249278 times.
✓ Branch 1 taken 1039913 times.
✓ Branch 2 taken 4249278 times.
✓ Branch 3 taken 1039913 times.
 34457643 for (const auto& dataJ : gridGeometry.connectivityMap()[globalI])
41
1/2
✓ Branch 1 taken 4249278 times.
✗ Branch 2 not taken.
 11362634 pattern.add(dataJ.globalJ, globalI);
42 }
43 }
44
45 // matrix pattern for explicit Jacobians -> diagonal matrix
46 else
47 {
48
2/2
✓ Branch 0 taken 5000 times.
✓ Branch 1 taken 2 times.
 18102 for (unsigned int globalI = 0; globalI < numDofs; ++globalI)
49
1/2
✓ Branch 1 taken 5000 times.
✗ Branch 2 not taken.
 18092 pattern.add(globalI, globalI);
50 }
51
52 442 return pattern;
53 }
54
55 /*!
56 * \ingroup Assembly
57 * \brief Helper function to generate Jacobian pattern for the staggered method
58 */
59 template<bool isImplicit, class GridGeometry,
60 typename std::enable_if_t<( (GridGeometry::discMethod == DiscretizationMethods::staggered) ), int> = 0>
61 100 auto getJacobianPattern(const GridGeometry& gridGeometry)
62 {
63 // resize the jacobian and the residual
64 100 const auto numDofs = gridGeometry.numDofs();
65 100 Dune::MatrixIndexSet pattern(numDofs, numDofs);
66
67 100 const auto& connectivityMap = gridGeometry.connectivityMap();
68
69 100 auto fvGeometry = localView(gridGeometry);
70 // evaluate the actual pattern
71 550552 for (const auto& element : elements(gridGeometry.gridView()))
72 {
73 if(gridGeometry.isCellCenter())
74 {
75 // the global index of the element at hand
76 static constexpr auto cellCenterIdx = GridGeometry::cellCenterIdx();
77 178784 const auto ccGlobalI = gridGeometry.elementMapper().index(element);
78 89392 pattern.add(ccGlobalI, ccGlobalI);
79 1061172 for (auto&& ccGlobalJ : connectivityMap(cellCenterIdx, cellCenterIdx, ccGlobalI))
80 351802 pattern.add(ccGlobalI, ccGlobalJ);
81 }
82 else
83 {
84 static constexpr auto faceIdx = GridGeometry::faceIdx();
85 89392 fvGeometry.bindElement(element);
86
87 // loop over sub control faces
88 536352 for (auto&& scvf : scvfs(fvGeometry))
89 {
90 357568 const auto faceGlobalI = scvf.dofIndex();
91 357568 pattern.add(faceGlobalI, faceGlobalI);
92 6529296 for (auto&& faceGlobalJ : connectivityMap(faceIdx, faceIdx, scvf.index()))
93 2549512 pattern.add(faceGlobalI, faceGlobalJ);
94 }
95 }
96 }
97
98 100 return pattern;
99 }
100
101 /*!
102 * \ingroup Assembly
103 * \brief Helper function to generate Jacobian pattern for finite element scheme
104 */
105 template<class FEBasis>
106 17 Dune::MatrixIndexSet getFEJacobianPattern(const FEBasis& feBasis)
107 {
108 17 const auto numDofs = feBasis.size();
109
110
1/2
✓ Branch 1 taken 13 times.
✗ Branch 2 not taken.
 17 Dune::MatrixIndexSet pattern;
111
1/2
✓ Branch 1 taken 13 times.
✗ Branch 2 not taken.
 17 pattern.resize(numDofs, numDofs);
112
113
1/2
✓ Branch 1 taken 13 times.
✗ Branch 2 not taken.
 18 auto localView = feBasis.localView();
114 // matrix pattern for implicit Jacobians
115
4/8
✓ Branch 1 taken 13 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 76913 times.
✗ Branch 7 not taken.
✓ Branch 10 taken 76900 times.
✗ Branch 11 not taken.
 154463 for (const auto& element : elements(feBasis.gridView()))
116 {
117
1/2
✓ Branch 1 taken 76900 times.
✗ Branch 2 not taken.
 77213 localView.bind(element);
118
119 77213 const auto& finiteElement = localView.tree().finiteElement();
120
2/4
✓ Branch 1 taken 76900 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 76900 times.
✗ Branch 5 not taken.
 77213 const auto numLocalDofs = finiteElement.localBasis().size();
121
2/2
✓ Branch 0 taken 192900 times.
✓ Branch 1 taken 76900 times.
 272352 for (std::size_t i = 0; i < numLocalDofs; i++)
122 {
123 195139 const auto dofIdxI = localView.index(i);
124
2/2
✓ Branch 0 taken 660900 times.
✓ Branch 1 taken 192900 times.
 873956 for (std::size_t j = 0; j < numLocalDofs; j++)
125 {
126
1/2
✓ Branch 1 taken 660900 times.
✗ Branch 2 not taken.
 678817 const auto dofIdxJ = localView.index(j);
127
3/6
✓ Branch 1 taken 660900 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 660900 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 660900 times.
✗ Branch 8 not taken.
 2036451 pattern.add(dofIdxI, dofIdxJ);
128 }
129 }
130 }
131
132 17 return pattern;
133 }
134
135 /*!
136 * \ingroup Assembly
137 * \brief Helper function to generate Jacobian pattern for finite element scheme
138 * \note This interface is for compatibility with the other schemes. The pattern
139 * in fem is the same independent of the time discretization scheme.
140 */
141 template<bool isImplicit, class GridGeometry,
142 typename std::enable_if_t<(GridGeometry::discMethod == DiscretizationMethods::fem), int> = 0>
143 Dune::MatrixIndexSet getJacobianPattern(const GridGeometry& gridGeometry)
144 { return getFEJacobianPattern(gridGeometry.feBasis()); }
145
146 /*!
147 * \ingroup Assembly
148 * \brief Helper function to generate Jacobian pattern for the face-centered methods
149 */
150 template<bool isImplicit, class GridGeometry,
151 typename std::enable_if_t<(GridGeometry::discMethod == DiscretizationMethods::fcstaggered), int> = 0>
152 65 Dune::MatrixIndexSet getJacobianPattern(const GridGeometry& gridGeometry)
153 {
154 // resize the jacobian and the residual
155 65 const auto numDofs = gridGeometry.numDofs();
156 65 Dune::MatrixIndexSet pattern(numDofs, numDofs);
157
158
1/2
✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
 65 const auto& connectivityMap = gridGeometry.connectivityMap();
159
1/2
✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
 65 auto fvGeometry = localView(gridGeometry);
160
161 // set the pattern
162
14/19
✓ Branch 1 taken 64 times.
✓ Branch 2 taken 101 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 164 times.
✓ Branch 5 taken 1 times.
✓ Branch 6 taken 348096 times.
✓ Branch 7 taken 6 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 1250 times.
✓ Branch 10 taken 998 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 3536 times.
✓ Branch 13 taken 6 times.
✓ Branch 14 taken 3536 times.
✓ Branch 15 taken 6 times.
✓ Branch 17 taken 3536 times.
✗ Branch 18 not taken.
✓ Branch 20 taken 3536 times.
✗ Branch 21 not taken.
 699945 for (const auto& element : elements(gridGeometry.gridView()))
163 {
164
1/2
✓ Branch 1 taken 5778 times.
✗ Branch 2 not taken.
 351674 fvGeometry.bind(element);
165
5/5
✓ Branch 0 taken 12072 times.
✓ Branch 1 taken 1399442 times.
✓ Branch 2 taken 348656 times.
✓ Branch 3 taken 1396424 times.
✓ Branch 4 taken 348656 times.
 1760170 for (const auto& scv : scvs(fvGeometry))
166 {
167 1408496 const auto globalI = scv.dofIndex();
168
1/2
✓ Branch 1 taken 1408496 times.
✗ Branch 2 not taken.
 1408496 pattern.add(globalI, globalI);
169
170
4/4
✓ Branch 0 taken 9820560 times.
✓ Branch 1 taken 1408496 times.
✓ Branch 2 taken 9820560 times.
✓ Branch 3 taken 1408496 times.
 15454544 for (const auto& scvIdxJ : connectivityMap[scv.index()])
171 {
172
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9737680 times.
 9820560 const auto globalJ = fvGeometry.scv(scvIdxJ).dofIndex();
173
1/2
✓ Branch 1 taken 9820560 times.
✗ Branch 2 not taken.
 9820560 pattern.add(globalI, globalJ);
174
175
2/2
✓ Branch 0 taken 2720 times.
✓ Branch 1 taken 9817840 times.
 9820560 if (gridGeometry.isPeriodic())
176 {
177
3/4
✓ Branch 1 taken 136 times.
✓ Branch 2 taken 2584 times.
✓ Branch 3 taken 136 times.
✗ Branch 4 not taken.
 2720 if (gridGeometry.dofOnPeriodicBoundary(globalI) && globalI != globalJ)
178 {
179
1/2
✓ Branch 1 taken 136 times.
✗ Branch 2 not taken.
 136 const auto globalIP = gridGeometry.periodicallyMappedDof(globalI);
180
1/2
✓ Branch 1 taken 136 times.
✗ Branch 2 not taken.
 136 pattern.add(globalIP, globalI);
181
1/2
✓ Branch 1 taken 136 times.
✗ Branch 2 not taken.
 136 pattern.add(globalI, globalIP);
182
183
2/2
✓ Branch 0 taken 68 times.
✓ Branch 1 taken 68 times.
 136 if (globalI > globalIP)
184
1/2
✓ Branch 1 taken 68 times.
✗ Branch 2 not taken.
 68 pattern.add(globalIP, globalJ);
185 }
186 }
187 }
188 }
189 }
190
191 65 return pattern;
192 }
193
194 /*!
195 * \ingroup Assembly
196 * \brief Compute the Jacobian matrix sparsity pattern for control volume finite element schemes
197 */
198 template<bool isImplicit, class GridGeometry,
199 typename std::enable_if_t<DiscretizationMethods::isCVFE<typename GridGeometry::DiscretizationMethod>, int> = 0>
200 256 Dune::MatrixIndexSet getJacobianPattern(const GridGeometry& gridGeometry)
201 {
202 // resize the jacobian and the residual
203
1/3
✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 256 const auto numDofs = gridGeometry.numDofs();
204 256 Dune::MatrixIndexSet pattern(numDofs, numDofs);
205
206 // matrix pattern for implicit Jacobians
207 if constexpr (isImplicit)
208 {
209
1/4
✓ Branch 1 taken 162 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 301 auto fvGeometry = localView(gridGeometry);
210
18/21
✓ Branch 1 taken 164 times.
✓ Branch 2 taken 52439 times.
✓ Branch 3 taken 28 times.
✓ Branch 4 taken 52602 times.
✓ Branch 5 taken 29 times.
✓ Branch 6 taken 203147 times.
✓ Branch 7 taken 35503 times.
✓ Branch 8 taken 75345 times.
✓ Branch 9 taken 36 times.
✓ Branch 10 taken 56421 times.
✓ Branch 11 taken 75345 times.
✓ Branch 12 taken 50386 times.
✓ Branch 13 taken 27 times.
✓ Branch 14 taken 125731 times.
✓ Branch 15 taken 27 times.
✓ Branch 17 taken 50386 times.
✗ Branch 18 not taken.
✓ Branch 20 taken 50386 times.
✗ Branch 21 not taken.
✓ Branch 25 taken 19 times.
✗ Branch 26 not taken.
 769607 for (const auto& element : elements(gridGeometry.gridView()))
211 {
212
1/2
✓ Branch 1 taken 227565 times.
✗ Branch 2 not taken.
 503380 fvGeometry.bindElement(element);
213
4/4
✓ Branch 0 taken 1180837 times.
✓ Branch 1 taken 305870 times.
✓ Branch 2 taken 1107607 times.
✓ Branch 3 taken 269255 times.
 2954806 for (const auto& scv : scvs(fvGeometry))
214 {
215 1984661 const auto globalI = scv.dofIndex();
216
4/4
✓ Branch 0 taken 5100847 times.
✓ Branch 1 taken 1180837 times.
✓ Branch 2 taken 4954387 times.
✓ Branch 3 taken 1107607 times.
 13449749 for (const auto& scvJ : scvs(fvGeometry))
217 {
218 8562995 const auto globalJ = scvJ.dofIndex();
219
1/2
✓ Branch 1 taken 5100847 times.
✗ Branch 2 not taken.
 8562995 pattern.add(globalI, globalJ);
220
221
2/2
✓ Branch 0 taken 160000 times.
✓ Branch 1 taken 4658059 times.
 8562995 if (gridGeometry.isPeriodic())
222 {
223
4/4
✓ Branch 1 taken 1600 times.
✓ Branch 2 taken 158400 times.
✓ Branch 3 taken 1200 times.
✓ Branch 4 taken 400 times.
 160000 if (gridGeometry.dofOnPeriodicBoundary(globalI) && globalI != globalJ)
224 {
225
1/2
✓ Branch 1 taken 1200 times.
✗ Branch 2 not taken.
 1200 const auto globalIP = gridGeometry.periodicallyMappedDof(globalI);
226
1/2
✓ Branch 1 taken 1200 times.
✗ Branch 2 not taken.
 1200 pattern.add(globalIP, globalI);
227
1/2
✓ Branch 1 taken 1200 times.
✗ Branch 2 not taken.
 1200 pattern.add(globalI, globalIP);
228
229
2/2
✓ Branch 0 taken 600 times.
✓ Branch 1 taken 600 times.
 1200 if (globalI > globalIP)
230
1/2
✓ Branch 1 taken 600 times.
✗ Branch 2 not taken.
 600 pattern.add(globalIP, globalJ);
231 }
232 }
233 }
234 }
235 }
236 }
237
238 // matrix pattern for explicit Jacobians -> diagonal matrix
239 else
240 {
241
2/2
✓ Branch 1 taken 5202 times.
✓ Branch 2 taken 2 times.
 6408 for (unsigned int globalI = 0; globalI < numDofs; ++globalI)
242
1/2
✓ Branch 1 taken 5202 times.
✗ Branch 2 not taken.
 6404 pattern.add(globalI, globalI);
243 }
244
245 256 return pattern;
246 }
247
248
249 } // namespace Dumux
250
251 #endif
252