GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: dumux/dumux/assembly/jacobianpattern.hh
Date: 2025-04-12 19:19:20
Exec Total Coverage
Lines: 87 91 95.6%
Functions: 140 149 94.0%
Branches: 105 151 69.5%
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 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
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 621 Dune::MatrixIndexSet getJacobianPattern(const GridGeometry& gridGeometry)
29 {
30
1/2
✓ Branch 1 taken 457 times.
✗ Branch 2 not taken.
 621 const auto numDofs = gridGeometry.numDofs();
31 621 Dune::MatrixIndexSet pattern;
32
1/2
✓ Branch 1 taken 457 times.
✗ Branch 2 not taken.
 621 pattern.resize(numDofs, numDofs);
33
34 // matrix pattern for implicit Jacobians
35 if (isImplicit)
36 {
37
2/2
✓ Branch 0 taken 2175280 times.
✓ Branch 1 taken 447 times.
 3505077 for (unsigned int globalI = 0; globalI < numDofs; ++globalI)
38 {
39 3504474 pattern.add(globalI, globalI);
40
2/2
✓ Branch 0 taken 10610988 times.
✓ Branch 1 taken 2175280 times.
 21277960 for (const auto& dataJ : gridGeometry.connectivityMap()[globalI])
41
1/2
✓ Branch 1 taken 10610988 times.
✗ Branch 2 not taken.
 17773486 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 18092 times.
✓ Branch 1 taken 10 times.
 31202 for (unsigned int globalI = 0; globalI < numDofs; ++globalI)
49 31184 pattern.add(globalI, globalI);
50 }
51
52 621 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 204 auto getJacobianPattern(const GridGeometry& gridGeometry)
62 {
63 // resize the jacobian and the residual
64 204 const auto numDofs = gridGeometry.numDofs();
65 204 Dune::MatrixIndexSet pattern(numDofs, numDofs);
66
67
1/2
✓ Branch 1 taken 102 times.
✗ Branch 2 not taken.
 204 const auto& connectivityMap = gridGeometry.connectivityMap();
68
69 204 auto fvGeometry = localView(gridGeometry);
70 // evaluate the actual pattern
71
4/8
✓ Branch 1 taken 102 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 179280 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 14006 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 14000 times.
✗ Branch 10 not taken.
 1075308 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
3/5
✓ Branch 1 taken 7000 times.
✓ Branch 2 taken 82592 times.
✓ Branch 4 taken 7000 times.
✗ Branch 5 not taken.
✗ Branch 3 not taken.
 179184 const auto ccGlobalI = gridGeometry.elementMapper().index(element);
78 179184 pattern.add(ccGlobalI, ccGlobalI);
79
2/2
✓ Branch 0 taken 352542 times.
✓ Branch 1 taken 89592 times.
 884268 for (auto&& ccGlobalJ : connectivityMap(cellCenterIdx, cellCenterIdx, ccGlobalI))
80
1/2
✓ Branch 1 taken 352542 times.
✗ Branch 2 not taken.
 705084 pattern.add(ccGlobalI, ccGlobalJ);
81 }
82 else
83 {
84 static constexpr auto faceIdx = GridGeometry::faceIdx();
85
1/2
✓ Branch 1 taken 7000 times.
✗ Branch 2 not taken.
 179184 fvGeometry.bindElement(element);
86
87 // loop over sub control faces
88
3/4
✓ Branch 1 taken 358368 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 358368 times.
✓ Branch 4 taken 89592 times.
 895920 for (auto&& scvf : scvfs(fvGeometry))
89 {
90
1/2
✓ Branch 1 taken 358368 times.
✗ Branch 2 not taken.
 716736 const auto faceGlobalI = scvf.dofIndex();
91 716736 pattern.add(faceGlobalI, faceGlobalI);
92
2/2
✓ Branch 0 taken 2554872 times.
✓ Branch 1 taken 358368 times.
 5826480 for (auto&& faceGlobalJ : connectivityMap(faceIdx, faceIdx, scvf.index()))
93
1/2
✓ Branch 1 taken 2554872 times.
✗ Branch 2 not taken.
 5109744 pattern.add(faceGlobalI, faceGlobalJ);
94 }
95 }
96 }
97
98 204 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 21 Dune::MatrixIndexSet getFEJacobianPattern(const FEBasis& feBasis)
107 {
108
1/2
✓ Branch 1 taken 17 times.
✗ Branch 2 not taken.
 21 const auto numDofs = feBasis.size();
109
110 21 Dune::MatrixIndexSet pattern;
111
1/2
✓ Branch 1 taken 17 times.
✗ Branch 2 not taken.
 21 pattern.resize(numDofs, numDofs);
112
113
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 21 auto localView = feBasis.localView();
114 // matrix pattern for implicit Jacobians
115
5/7
✓ Branch 1 taken 16 times.
✓ Branch 2 taken 13 times.
✓ Branch 4 taken 77216 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 77200 times.
✗ Branch 8 not taken.
✓ Branch 3 taken 1 times.
 232573 for (const auto& element : elements(feBasis.gridView()))
116 {
117
1/2
✓ Branch 1 taken 77213 times.
✗ Branch 2 not taken.
 77526 localView.bind(element);
118
119 77526 const auto& finiteElement = localView.tree().finiteElement();
120 77526 const auto numLocalDofs = finiteElement.localBasis().size();
121
2/2
✓ Branch 0 taken 195139 times.
✓ Branch 1 taken 77213 times.
 274904 for (std::size_t i = 0; i < numLocalDofs; i++)
122 {
123 197378 const auto dofIdxI = localView.index(i);
124
2/2
✓ Branch 0 taken 678817 times.
✓ Branch 1 taken 195139 times.
 894112 for (std::size_t j = 0; j < numLocalDofs; j++)
125 {
126
1/2
✓ Branch 1 taken 678817 times.
✗ Branch 2 not taken.
 696734 const auto dofIdxJ = localView.index(j);
127
1/2
✓ Branch 1 taken 678817 times.
✗ Branch 2 not taken.
 696734 pattern.add(dofIdxI, dofIdxJ);
128 }
129 }
130 }
131
132 21 return pattern;
133 21 }
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 78 Dune::MatrixIndexSet getJacobianPattern(const GridGeometry& gridGeometry)
153 {
154 // resize the jacobian and the residual
155 78 const auto numDofs = gridGeometry.numDofs();
156 78 Dune::MatrixIndexSet pattern(numDofs, numDofs);
157
158
1/2
✓ Branch 1 taken 70 times.
✗ Branch 2 not taken.
 78 const auto& connectivityMap = gridGeometry.connectivityMap();
159 78 auto fvGeometry = localView(gridGeometry);
160
161 // set the pattern
162
11/15
✓ Branch 1 taken 74 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 371431 times.
✓ Branch 5 taken 100 times.
✓ Branch 7 taken 13591 times.
✓ Branch 8 taken 3 times.
✓ Branch 10 taken 3536 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 3536 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 3536 times.
✓ Branch 16 taken 6 times.
✓ Branch 6 taken 14279 times.
✓ Branch 9 taken 1935 times.
✗ Branch 3 not taken.
 1161901 for (const auto& element : elements(gridGeometry.gridView()))
163 {
164 399304 fvGeometry.bind(element);
165
5/5
✓ Branch 1 taken 1484312 times.
✓ Branch 2 taken 68302 times.
✓ Branch 3 taken 1484312 times.
✓ Branch 4 taken 439396 times.
✓ Branch 5 taken 16138 times.
 1939846 for (const auto& scv : scvs(fvGeometry))
166 {
167
1/2
✓ Branch 1 taken 1552614 times.
✗ Branch 2 not taken.
 1552614 const auto globalI = scv.dofIndex();
168 1552614 pattern.add(globalI, globalI);
169
170
2/2
✓ Branch 0 taken 10834006 times.
✓ Branch 1 taken 1552614 times.
 12386620 for (const auto& scvIdxJ : connectivityMap[scv.index()])
171 {
172
3/5
✗ Branch 0 not taken.
✓ Branch 1 taken 10751126 times.
✓ Branch 3 taken 506822 times.
✗ Branch 4 not taken.
✓ Branch 2 taken 82880 times.
 10834006 const auto globalJ = fvGeometry.scv(scvIdxJ).dofIndex();
173 10834006 pattern.add(globalI, globalJ);
174
175
2/2
✓ Branch 0 taken 334548 times.
✓ Branch 1 taken 10499458 times.
 10834006 if (gridGeometry.isPeriodic())
176 {
177
1/2
✓ Branch 0 taken 1650 times.
✗ Branch 1 not taken.
 334548 if (gridGeometry.dofOnPeriodicBoundary(globalI) && globalI != globalJ)
178 {
179
1/2
✓ Branch 1 taken 1650 times.
✗ Branch 2 not taken.
 1650 const auto globalIP = gridGeometry.periodicallyMappedDof(globalI);
180
1/2
✓ Branch 1 taken 1650 times.
✗ Branch 2 not taken.
 1650 pattern.add(globalIP, globalI);
181 1650 pattern.add(globalI, globalIP);
182
183
2/2
✓ Branch 0 taken 825 times.
✓ Branch 1 taken 825 times.
 1650 if (globalI > globalIP)
184 825 pattern.add(globalIP, globalJ);
185 }
186 }
187 }
188 }
189 }
190
191 78 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
1/3
✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 334 Dune::MatrixIndexSet getJacobianPattern(const GridGeometry& gridGeometry)
201 {
202 // resize the jacobian and the residual
203 334 const auto numDofs = gridGeometry.numDofs();
204
1/2
✓ Branch 2 taken 207 times.
✗ Branch 3 not taken.
 334 Dune::MatrixIndexSet pattern(numDofs, numDofs);
205
206 // matrix pattern for implicit Jacobians
207 if constexpr (isImplicit)
208 {
209 254 auto fvGeometry = localView(gridGeometry);
210
12/15
✓ Branch 1 taken 209 times.
✓ Branch 2 taken 43346 times.
✓ Branch 4 taken 458083 times.
✓ Branch 5 taken 1 times.
✓ Branch 7 taken 54970 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 119136 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 119136 times.
✓ Branch 13 taken 51 times.
✓ Branch 16 taken 22 times.
✗ Branch 17 not taken.
✓ Branch 3 taken 10056 times.
✓ Branch 6 taken 314505 times.
✓ Branch 9 taken 943 times.
 1963367 for (const auto& element : elements(gridGeometry.gridView()))
211 {
212
1/2
✓ Branch 1 taken 534838 times.
✗ Branch 2 not taken.
 905144 fvGeometry.bindElement(element);
213
2/2
✓ Branch 0 taken 2365716 times.
✓ Branch 1 taken 620406 times.
 4365700 for (const auto& scv : scvs(fvGeometry))
214 {
215 3460556 const auto globalI = scv.dofIndex();
216
3/4
✓ Branch 1 taken 9831894 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 9831894 times.
✓ Branch 4 taken 2365716 times.
 17744974 for (const auto& scvJ : scvs(fvGeometry))
217 {
218
1/2
✓ Branch 1 taken 9831894 times.
✗ Branch 2 not taken.
 14284418 const auto globalJ = scvJ.dofIndex();
219 14284418 pattern.add(globalI, globalJ);
220
221
2/2
✓ Branch 0 taken 160000 times.
✓ Branch 1 taken 8607986 times.
 14284418 if (gridGeometry.isPeriodic())
222 {
223
2/2
✓ Branch 0 taken 1200 times.
✓ Branch 1 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 1200 pattern.add(globalI, globalIP);
228
229
2/2
✓ Branch 0 taken 600 times.
✓ Branch 1 taken 600 times.
 1200 if (globalI > globalIP)
230 600 pattern.add(globalIP, globalJ);
231 }
232 }
233 }
234 }
235 }
236 206 }
237
238 // matrix pattern for explicit Jacobians -> diagonal matrix
239 else
240 {
241
2/2
✓ Branch 1 taken 6404 times.
✓ Branch 2 taken 4 times.
 7612 for (unsigned int globalI = 0; globalI < numDofs; ++globalI)
242 7606 pattern.add(globalI, globalI);
243 }
244
245 334 return pattern;
246 ✗ }
247
248
249 } // namespace Dumux
250
251 #endif
252