GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/freeflow/navierstokes/staggered/fluxoversurface.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 59 69 85.5%
Functions: 24 30 80.0%
Branches: 77 148 52.0%
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 NavierStokesModel
10 * \copydoc Dumux::FluxOverSurface
11 */
12 #ifndef DUMUX_FLUX_OVER_SURFACE_STAGGERED_HH
13 #define DUMUX_FLUX_OVER_SURFACE_STAGGERED_HH
14
15 #include <numeric>
16 #include <functional>
17 #include <type_traits>
18 #include <vector>
19
20 #include <dune/common/exceptions.hh>
21 #include <dune/common/fvector.hh>
22 #include <dune/geometry/type.hh>
23 #include <dune/geometry/multilineargeometry.hh>
24
25 #include <dumux/common/parameters.hh>
26 #include <dumux/geometry/makegeometry.hh>
27 #include <dumux/geometry/intersectspointgeometry.hh>
28 #include <dumux/discretization/extrusion.hh>
29
30 namespace Dumux {
31
32 /*!
33 * \ingroup NavierStokesModel
34 * \brief Class used to calculate fluxes over surfaces. This only works for the staggered grid discretization.
35 */
36 template<class GridVariables, class SolutionVector, class ModelTraits, class LocalResidual>
37
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 3 class FluxOverSurface
38 {
39 using Scalar = typename GridVariables::Scalar;
40 using GridGeometry = typename GridVariables::GridGeometry;
41 using FVElementGeometry = typename GridGeometry::LocalView;
42 using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
43 using Extrusion = Extrusion_t<GridGeometry>;
44 using GridView = typename GridGeometry::GridView;
45 using VolumeVariables = typename GridVariables::VolumeVariables;
46 using Element = typename GridView::template Codim<0>::Entity;
47
48 using CellCenterPrimaryVariables = std::decay_t<decltype(std::declval<SolutionVector>()[GridGeometry::cellCenterIdx()][0])>;
49
50 enum {
51 // Grid and world dimension
52 dim = GridView::dimension,
53 dimWorld = GridView::dimensionworld
54 };
55
56 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
57
58 static constexpr auto surfaceDim = dimWorld - 1;
59 using SurfaceGeometryType = Dune::MultiLinearGeometry< Scalar, surfaceDim, dimWorld >;
60
61 /*!
62 * \brief Auxiliary class that holds the surface-specific data.
63 */
64 template<int mydim, int coordDim, class Scalar=double>
65 class SurfaceData
66 {
67 using Geo = SurfaceGeometryType;
68 public:
69
70 ✗ SurfaceData() {}
71
72 SurfaceData(Geo&& geo)
73 {
74 values_.resize(1);
75 geometries_.push_back(std::move(geo));
76 }
77
78 SurfaceData(std::vector<Geo>&& geo)
79 {
80 values_.resize(geo.size());
81 geometries_ = std::forward<decltype(geo)>(geo);
82 }
83
84 void addValue(int surfaceIdx, const CellCenterPrimaryVariables& value)
85 {
86
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 100 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 100 times.
 5256 values_[surfaceIdx] += value;
87 }
88
89 void addSubSurface(Geo&& geo)
90 {
91 12 values_.emplace_back(0.0);
92
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 geometries_.push_back(std::move(geo));
93 }
94
95 auto& subSurfaces() const
96 {
97 return geometries_;
98 }
99
100 CellCenterPrimaryVariables& value(int surfaceIdx) const
101 {
102 return values_[surfaceIdx];
103 }
104
105 auto& values() const
106 {
107 320 return values_;
108 }
109
110 ✗ void printSurfaceBoundaries(int surfaceIdx) const
111 {
112 ✗ const auto& geometry = geometries_[surfaceIdx];
113 ✗ for(int i = 0; i < geometry.corners(); ++i)
114 ✗ std::cout << "(" << geometry.corner(i) << ") ";
115 ✗ }
116
117 void resetValues()
118 {
119 640 std::fill(values_.begin(), values_.end(), CellCenterPrimaryVariables(0.0));
120 }
121
122 private:
123
124 std::vector<Geo> geometries_;
125 std::vector<CellCenterPrimaryVariables> values_;
126 };
127
128 public:
129
130 using SurfaceList = std::vector<SurfaceGeometryType>;
131
132 /*!
133 * \brief The constructor
134 */
135 template<class Sol>
136 3 FluxOverSurface(const GridVariables& gridVariables,
137 const Sol& sol)
138 : gridVariables_(gridVariables),
139
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 sol_(sol)
140 {
141 static_assert(std::is_same<Sol, SolutionVector>::value, "Make sure that sol has the same type as SolutionVector."
142 "Use StaggeredVtkOutputModule<GridVariables, decltype(sol)> when calling the constructor.");
143
144
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 verbose_ = getParamFromGroup<bool>(problem_().paramGroup(), "FluxOverSurface.Verbose", false);
145 3 }
146
147 /*!
148 * \brief Add a collection of sub surfaces under a given name
149 *
150 * \param name The name of the surface
151 * \param surfaces The list of sub surfaces
152 */
153 void addSurface(const std::string& name, SurfaceList&& surfaces)
154 {
155 surfaces_[name] = SurfaceData<surfaceDim,dim>(std::forward<decltype(surfaces)>(surfaces));
156 }
157
158 /*!
159 * \brief Add a surface under a given name, specifying the surface's corner points.
160 * This is a specialization for 2D, therefore the surface is actually a line.
161 *
162 * \param name The name of the surface
163 * \param p0 The first corner
164 * \param p1 The second corner
165 */
166 6 void addSurface(const std::string& name, const GlobalPosition& p0, const GlobalPosition& p1)
167 {
168
6/16
✓ Branch 2 taken 6 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 6 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 6 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 6 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✓ Branch 14 taken 6 times.
✓ Branch 15 taken 6 times.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
 18 surfaces_[name].addSubSurface(makeSurface(std::vector<GlobalPosition>{p0, p1}));
169 6 }
170
171 /*!
172 * \brief Add a surface under a given name, specifying the surface's corner points.
173 * This is a specialization for 3D.
174 *
175 * \param name The name of the surface
176 * \param p0 The first corner
177 * \param p1 The second corner
178 * \param p2 The third corner
179 * \param p3 The fourth corner
180 */
181 void addSurface(const std::string& name,
182 const GlobalPosition& p0,
183 const GlobalPosition& p1,
184 const GlobalPosition& p2,
185 const GlobalPosition& p3)
186 {
187 surfaces_[name].addSubSurface(makeSurface(std::vector<GlobalPosition>{p0, p1, p2, p3}));
188 }
189
190 /*!
191 * \brief Creates a geometrical surface object for (2D).
192 *
193 * \param corners The vector storing the surface's corners
194 */
195 6 static SurfaceGeometryType makeSurface(const std::vector<Dune::FieldVector<Scalar, 2>>& corners)
196 {
197 6 return SurfaceGeometryType(Dune::GeometryTypes::line, corners);
198 }
199
200 /*!
201 * \brief Creates a geometrical surface object for (3D).
202 *
203 * \param corners The vector storing the surface's corners
204 */
205 static SurfaceGeometryType makeSurface(const std::vector<Dune::FieldVector<Scalar, 3>>& corners)
206 {
207 return makeDuneQuadrilaterial(corners);
208 }
209
210 /*!
211 * \brief Calculate the mass or mole fluxes over all surfaces
212 */
213 void calculateMassOrMoleFluxes()
214 {
215 2708 auto fluxType = [this](const auto& element,
216 const auto& fvGeometry,
217 const auto& elemVolVars,
218 const auto& elemFaceVars,
219 const auto& scvf,
220 2628 const auto& elemFluxVarsCache)
221 {
222 7884 LocalResidual localResidual(&problem_());
223
224
2/2
✓ Branch 0 taken 2578 times.
✓ Branch 1 taken 50 times.
 2628 if (scvf.boundary())
225 2578 return localResidual.computeBoundaryFluxForCellCenter(problem_(), element, fvGeometry, scvf, elemVolVars, elemFaceVars, /*elemBcTypes*/{}, elemFluxVarsCache);
226 else
227 100 return localResidual.computeFluxForCellCenter(problem_(), element, fvGeometry, elemVolVars, elemFaceVars, scvf, elemFluxVarsCache);
228 };
229
230
1/2
✓ Branch 1 taken 80 times.
✗ Branch 2 not taken.
 80 calculateFluxes(fluxType);
231 }
232
233 /*!
234 * \brief Calculate the volume fluxes over all surfaces.
235 */
236 void calculateVolumeFluxes()
237 {
238 auto fluxType = [](const auto& element,
239 const auto& fvGeometry,
240 const auto& elemVolVars,
241 const auto& elemFaceVars,
242 const auto& scvf,
243 2628 const auto& elemFluxVarsCache)
244 {
245 2628 CellCenterPrimaryVariables result(0.0);
246 5256 const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
247 5256 const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
248
1/2
✓ Branch 0 taken 2628 times.
✗ Branch 1 not taken.
 2628 const Scalar extrusionFactor = harmonicMean(insideVolVars.extrusionFactor(), outsideVolVars.extrusionFactor());
249 10412 result[0] = elemFaceVars[scvf].velocitySelf() * Extrusion::area(fvGeometry, scvf) * extrusionFactor * scvf.directionSign();
250 2628 return result;
251 };
252
253
1/2
✓ Branch 1 taken 80 times.
✗ Branch 2 not taken.
 80 calculateFluxes(fluxType);
254 }
255
256 /*!
257 * \brief Calculate the fluxes over all surfaces for a given flux type.
258 *
259 * \param fluxType The flux type. This can be a lambda of the following form:
260 * [](const auto& element,
261 const auto& fvGeometry,
262 const auto& elemVolVars,
263 const auto& elemFaceVars,
264 const auto& scvf,
265 const auto& elemFluxVarsCache)
266 { return ... ; }
267 */
268 template<class FluxType>
269 320 void calculateFluxes(const FluxType& fluxType)
270 {
271 // make sure to reset all the values of the surfaces, in case this method has been called already before
272
2/2
✓ Branch 0 taken 320 times.
✓ Branch 1 taken 160 times.
 1600 for(auto&& surface : surfaces_)
273 1280 surface.second.resetValues();
274
275 // make sure not to iterate over the same dofs twice
276
4/8
✓ Branch 1 taken 160 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 160 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 160 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 160 times.
✗ Branch 11 not taken.
 640 std::vector<bool> dofVisited(problem_().gridGeometry().numFaceDofs(), false);
277
278
1/2
✓ Branch 1 taken 160 times.
✗ Branch 2 not taken.
 320 auto fvGeometry = localView(problem_().gridGeometry());
279
2/4
✓ Branch 1 taken 160 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 160 times.
✗ Branch 5 not taken.
 960 auto elemVolVars = localView(gridVariables_.curGridVolVars());
280
2/4
✓ Branch 1 taken 160 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 160 times.
✗ Branch 5 not taken.
 640 auto elemFaceVars = localView(gridVariables_.curGridFaceVars());
281
2/4
✓ Branch 1 taken 160 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 160 times.
✗ Branch 5 not taken.
 640 auto elemFluxVarsCache = localView(gridVariables_.gridFluxVarsCache());
282
283
3/6
✓ Branch 1 taken 160 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 160 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 30384 times.
✗ Branch 7 not taken.
 121856 for(auto&& element : elements(problem_().gridGeometry().gridView()))
284 {
285 60448 fvGeometry.bind(element);
286
1/2
✓ Branch 1 taken 30224 times.
✗ Branch 2 not taken.
 60448 elemVolVars.bind(element, fvGeometry, sol_);
287 60448 elemFaceVars.bind(element, fvGeometry, sol_);
288 60448 elemFluxVarsCache.bind(element, fvGeometry, elemVolVars);
289
290
6/6
✓ Branch 0 taken 120896 times.
✓ Branch 1 taken 30224 times.
✓ Branch 2 taken 120896 times.
✓ Branch 3 taken 30224 times.
✓ Branch 4 taken 64540 times.
✓ Branch 5 taken 56356 times.
 362688 for(auto && scvf : scvfs(fvGeometry))
291 {
292
2/2
✓ Branch 0 taken 64540 times.
✓ Branch 1 taken 56356 times.
 241792 const auto dofIdx = scvf.dofIndex();
293 // do nothing of the dof was already visited
294
6/6
✓ Branch 0 taken 64540 times.
✓ Branch 1 taken 56356 times.
✓ Branch 2 taken 64540 times.
✓ Branch 3 taken 56356 times.
✓ Branch 4 taken 64540 times.
✓ Branch 5 taken 56356 times.
 725376 if(dofVisited[dofIdx])
295 continue;
296
297 258160 dofVisited[dofIdx] = true;
298
299 // iterate through all surfaces and check if the flux at the given position
300 // should be accounted for in the respective surface
301
2/2
✓ Branch 0 taken 129080 times.
✓ Branch 1 taken 64540 times.
 645400 for(auto&& surface : surfaces_)
302 {
303 const auto& subSurfaces = surface.second.subSurfaces();
304
305
4/4
✓ Branch 0 taken 129080 times.
✓ Branch 1 taken 129080 times.
✓ Branch 2 taken 129080 times.
✓ Branch 3 taken 129080 times.
 774480 for(int surfaceIdx = 0; surfaceIdx < subSurfaces.size(); ++surfaceIdx)
306 {
307
2/2
✓ Branch 3 taken 5256 times.
✓ Branch 4 taken 123824 times.
 774480 if(intersectsPointGeometry(scvf.center(), subSurfaces[surfaceIdx]))
308 {
309
1/2
✓ Branch 1 taken 2628 times.
✗ Branch 2 not taken.
 10512 const auto result = fluxType(element, fvGeometry, elemVolVars, elemFaceVars, scvf, elemFluxVarsCache);
310
311
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 200 times.
 10512 surface.second.addValue(surfaceIdx, result);
312
313
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5256 times.
 10512 if(verbose_)
314 {
315 ✗ std::cout << "Flux at face " << scvf.center() << " (" << surface.first << "): " << result;
316 ✗ std::cout << " (directionIndex: " << scvf.directionIndex() << "; surface boundaries: ";
317 ✗ surface.second.printSurfaceBoundaries(surfaceIdx);
318 ✗ std::cout << ", surfaceIdx " << surfaceIdx << ")" << std::endl;
319 }
320 }
321 }
322 }
323 }
324 }
325 320 }
326
327 /*!
328 * \brief Return the fluxes of the individual sub surface of a given name.
329 *
330 * \param name The name of the surface
331 */
332 auto& values(const std::string& name) const
333 {
334 960 return surfaces_.at(name).values();
335 }
336
337 /*!
338 * \brief Return the cumulative net fluxes of a surface of a given name.
339 *
340 * \param name The name of the surface
341 */
342 320 auto netFlux(const std::string& name) const
343 {
344 320 const auto& surfaceResults = values(name);
345 1276 return std::accumulate(surfaceResults.begin(), surfaceResults.end(), CellCenterPrimaryVariables(0.0));
346 }
347
348 private:
349
350
18/32
✓ Branch 1 taken 80 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 80 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 80 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 80 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 80 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 80 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 80 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 80 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 80 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 80 times.
✗ Branch 29 not taken.
✓ Branch 31 taken 80 times.
✗ Branch 32 not taken.
✓ Branch 34 taken 80 times.
✗ Branch 35 not taken.
✓ Branch 36 taken 2578 times.
✓ Branch 37 taken 50 times.
✓ Branch 38 taken 2578 times.
✓ Branch 39 taken 50 times.
✓ Branch 45 taken 3 times.
✗ Branch 46 not taken.
✓ Branch 48 taken 3 times.
✗ Branch 49 not taken.
 6222 const auto& problem_() const { return gridVariables_.curGridVolVars().problem(); }
351
352 std::map<std::string, SurfaceData<surfaceDim ,dim> > surfaces_;
353 const GridVariables& gridVariables_;
354 const SolutionVector& sol_;
355 bool verbose_;
356 };
357
358 } //end namespace
359
360 #endif
361