GCC Code Coverage Report


Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/porenetwork/common/boundaryflux.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 72 76 94.7%
Functions: 27 30 90.0%
Branches: 155 267 58.1%

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 PoreNetworkModels
10 * \copydoc Dumux::PoreNetwork::BoundaryFlux
11 */
12 #ifndef DUMUX_PNM_BOUNDARYFLUX_HH
13 #define DUMUX_PNM_BOUNDARYFLUX_HH
14
15 #include <algorithm>
16 #include <numeric>
17 #include <vector>
18 #include <type_traits>
19 #include <unordered_map>
20 #include <string_view>
21 #include <iostream>
22
23 #include <dune/common/exceptions.hh>
24 #include <dumux/common/typetraits/problem.hh>
25 #include <dumux/discretization/cvfe/elementboundarytypes.hh>
26
27 namespace Dumux::PoreNetwork {
28
29 /*!
30 * \ingroup PoreNetworkModels
31 * \brief Class for the calculation of fluxes at the boundary
32 * of pore-network models.
33 */
34 template<class GridVariables, class LocalResidual, class SolutionVector>
35 class BoundaryFlux
36 {
37 using Problem = std::decay_t<decltype(std::declval<LocalResidual>().problem())>;
38 using GridGeometry = typename ProblemTraits<Problem>::GridGeometry;
39 using GridView = typename GridGeometry::GridView;
40 using Element = typename GridView::template Codim<0>::Entity;
41 using FVElementGeometry = typename GridGeometry::LocalView;
42 using SubControlVolume = typename FVElementGeometry::SubControlVolume;
43 using BoundaryTypes = typename ProblemTraits<Problem>::BoundaryTypes;
44 using ElementBoundaryTypes = CVFEElementBoundaryTypes<BoundaryTypes>;
45
46 using NumEqVector = typename SolutionVector::block_type;
47 static constexpr auto numEq = NumEqVector::dimension;
48
49 //! result struct that holds both the total flux and the flux per pore
50
5/9
✓ Branch 1 taken 12 times.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 10 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✓ Branch 11 taken 1 times.
✗ Branch 12 not taken.
59 struct Result
51 {
52 NumEqVector totalFlux;
53 std::unordered_map<std::size_t, NumEqVector> fluxPerPore;
54
55 //! make the total flux printable to e.g. std::cout
56 friend std::ostream& operator<< (std::ostream& stream, const Result& result)
57 {
58
8/16
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 1 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 1 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 1 times.
✗ Branch 23 not taken.
14 stream << result.totalFlux;
59 return stream;
60 }
61
62 //! allow to get the total flux for a given eqIdx
63 const auto& operator[] (int eqIdx) const
64
14/28
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 1 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 1 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 1 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 1 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 1 times.
✗ Branch 29 not taken.
✓ Branch 31 taken 1 times.
✗ Branch 32 not taken.
✓ Branch 37 taken 1 times.
✗ Branch 38 not taken.
✓ Branch 40 taken 1 times.
✗ Branch 41 not taken.
✓ Branch 43 taken 1 times.
✗ Branch 44 not taken.
37 { return totalFlux[eqIdx]; }
65
66 //! make the total flux assignable to NumEqVector through implicit conversion
67 operator NumEqVector() const
68 { return totalFlux; }
69 };
70
71 public:
72 // export the Scalar type
73 using Scalar = typename GridVariables::Scalar;
74
75 13 BoundaryFlux(const GridVariables& gridVariables,
76 const LocalResidual& localResidual,
77 const SolutionVector& sol)
78 : localResidual_(localResidual)
79 , gridVariables_(gridVariables)
80 , sol_(sol)
81
3/8
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 9 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
39 , isStationary_(localResidual.isStationary())
82 {
83
2/4
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9 times.
✗ Branch 5 not taken.
26 const auto numDofs = localResidual_.problem().gridGeometry().numDofs();
84
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
13 isConsidered_.resize(numDofs, false);
85
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
13 boundaryFluxes_.resize(numDofs);
86 13 }
87
88 /*!
89 * \brief Returns the cumulative flux in \f$\mathrm{[\frac{kg}{s}]}\f$, \f$\mathrm{[\frac{mole}{s}]}\f$ or \f$\mathrm{[\frac{J}{s}]}\f$ of several pores for a given list of pore labels to consider
90 *
91 * \param labels A list of pore labels which will be considered for the flux calculation
92 * \param verbose If set true, the fluxes at all individual SCVs are printed
93 */
94 template<class Label>
95 52 Result getFlux(const std::vector<Label>& labels, const bool verbose = false) const
96 {
97 // helper lambda to decide which scvs to consider for flux calculation
98 634558 auto restriction = [&] (const SubControlVolume& scv)
99 {
100 1269012 const Label poreLabel = localResidual_.problem().gridGeometry().poreLabel(scv.dofIndex());
101 1903518 return std::any_of(labels.begin(), labels.end(),
102
4/16
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 7234 times.
✓ Branch 7 taken 313550 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✓ Branch 14 taken 6421 times.
✓ Branch 15 taken 307301 times.
634506 [&](const Label l){ return l == poreLabel; });
103 };
104
105 162 std::fill(boundaryFluxes_.begin(), boundaryFluxes_.end(), NumEqVector(0.0));
106 156 std::fill(isConsidered_.begin(), isConsidered_.end(), false);
107
108 // sum up the fluxes
109
4/4
✓ Branch 3 taken 317253 times.
✓ Branch 4 taken 26 times.
✓ Branch 5 taken 317253 times.
✓ Branch 6 taken 26 times.
634662 for (const auto& element : elements(localResidual_.problem().gridGeometry().gridView()))
110 634506 computeBoundaryFlux_(element, restriction, verbose);
111
112 52 Result result;
113 162 result.totalFlux = std::accumulate(boundaryFluxes_.begin(), boundaryFluxes_.end(), NumEqVector(0.0));;
114
4/4
✓ Branch 0 taken 151128 times.
✓ Branch 1 taken 26 times.
✓ Branch 2 taken 151128 times.
✓ Branch 3 taken 26 times.
302308 for (int i = 0; i < isConsidered_.size(); ++i)
115 {
116
6/6
✓ Branch 0 taken 7117 times.
✓ Branch 1 taken 144011 times.
✓ Branch 2 taken 7117 times.
✓ Branch 3 taken 144011 times.
✓ Branch 4 taken 7117 times.
✓ Branch 5 taken 144011 times.
906768 if (isConsidered_[i])
117
2/4
✓ Branch 1 taken 7117 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7117 times.
✗ Branch 5 not taken.
28468 result.fluxPerPore[i] = boundaryFluxes_[i];
118 }
119
120 52 return result;
121 }
122
123 /*!
124 * \brief Returns the cumulative flux in \f$\mathrm{[\frac{kg}{s}]}\f$, \f$\mathrm{[\frac{mole}{s}]}\f$ or \f$\mathrm{[\frac{J}{s}]}\f$ of several pores at a given location on the boundary
125 *
126 * \param minMax Consider bBoxMin or bBoxMax by setting "min" or "max"
127 * \param coord x, y or z coordinate at which bBoxMin or bBoxMax is evaluated
128 * \param verbose If set true, the fluxes at all individual SCVs are printed
129 */
130 6 Result getFlux(std::string_view minMax, const int coord, const bool verbose = false) const
131 {
132
2/4
✗ Branch 1 not taken.
✓ Branch 2 taken 6 times.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
17 if (!(minMax == "min" || minMax == "max"))
133 DUNE_THROW(Dune::InvalidStateException,
134 "second argument must be either 'min' or 'max' (string) !");
135
136 6 const Scalar eps = 1e-6; //TODO
137 3110 auto onMeasuringBoundary = [&] (const Scalar pos)
138 {
139
9/10
✗ Branch 1 not taken.
✓ Branch 2 taken 776 times.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 1 times.
✓ Branch 5 taken 1 times.
✓ Branch 6 taken 1 times.
✓ Branch 7 taken 1 times.
✓ Branch 8 taken 1 times.
✓ Branch 9 taken 1 times.
✓ Branch 10 taken 1 times.
1558 return ( (minMax == "min" && pos < localResidual_.problem().gridGeometry().bBoxMin()[coord] + eps) ||
140
9/10
✓ Branch 1 taken 775 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 198 times.
✓ Branch 4 taken 576 times.
✓ Branch 5 taken 198 times.
✓ Branch 6 taken 576 times.
✓ Branch 7 taken 198 times.
✓ Branch 8 taken 576 times.
✓ Branch 9 taken 198 times.
✓ Branch 10 taken 576 times.
3872 (minMax == "max" && pos > localResidual_.problem().gridGeometry().bBoxMax()[coord] - eps) );
141 };
142
143 // helper lambda to decide which scvs to consider for flux calculation
144 3420 auto restriction = [&] (const SubControlVolume& scv)
145 {
146 1508 bool considerAllDirections = coord < 0 ? true : false;
147
148 //only consider SCVs on the boundary
149
8/8
✓ Branch 0 taken 776 times.
✓ Branch 1 taken 732 times.
✓ Branch 2 taken 776 times.
✓ Branch 3 taken 732 times.
✓ Branch 4 taken 776 times.
✓ Branch 5 taken 732 times.
✓ Branch 9 taken 577 times.
✓ Branch 10 taken 199 times.
4524 bool considerScv = localResidual_.problem().gridGeometry().dofOnBoundary(scv.dofIndex()) && onMeasuringBoundary(scv.dofPosition()[coord]);
150
151 //check whether a vertex lies on a boundary and also check whether this boundary shall be
152 // considered for the flux calculation
153
2/2
✓ Branch 0 taken 199 times.
✓ Branch 1 taken 1309 times.
1508 if (considerScv && !considerAllDirections)
154 {
155
2/2
✓ Branch 0 taken 198 times.
✓ Branch 1 taken 1 times.
199 const auto& pos = scv.dofPosition();
156
15/20
✓ Branch 0 taken 198 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 198 times.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 198 times.
✓ Branch 5 taken 1 times.
✓ Branch 6 taken 198 times.
✓ Branch 7 taken 1 times.
✓ Branch 8 taken 198 times.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 198 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 198 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 198 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 198 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 198 times.
995 if (!(pos[coord] < localResidual_.problem().gridGeometry().bBoxMin()[coord] + eps || pos[coord] > localResidual_.problem().gridGeometry().bBoxMax()[coord] -eps ))
157 considerScv = false;
158 }
159
160 1508 return considerScv;
161 };
162
163 18 std::fill(boundaryFluxes_.begin(), boundaryFluxes_.end(), NumEqVector(0.0));
164 18 std::fill(isConsidered_.begin(), isConsidered_.end(), false);
165
166 // sum up the fluxes
167
4/4
✓ Branch 3 taken 754 times.
✓ Branch 4 taken 6 times.
✓ Branch 5 taken 754 times.
✓ Branch 6 taken 6 times.
772 for (const auto& element : elements(localResidual_.problem().gridGeometry().gridView()))
168 754 computeBoundaryFlux_(element, restriction, verbose);
169
170 6 Result result;
171 18 result.totalFlux = std::accumulate(boundaryFluxes_.begin(), boundaryFluxes_.end(), NumEqVector(0.0));;
172
4/4
✓ Branch 0 taken 257 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 257 times.
✓ Branch 3 taken 6 times.
263 for (int i = 0; i < isConsidered_.size(); ++i)
173 {
174
6/6
✓ Branch 0 taken 27 times.
✓ Branch 1 taken 230 times.
✓ Branch 2 taken 27 times.
✓ Branch 3 taken 230 times.
✓ Branch 4 taken 27 times.
✓ Branch 5 taken 230 times.
771 if (isConsidered_[i])
175
2/4
✓ Branch 1 taken 27 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 27 times.
✗ Branch 5 not taken.
54 result.fluxPerPore[i] = boundaryFluxes_[i];
176 }
177
178 6 return result;
179 }
180
181 private:
182
183 /*!
184 * \brief Computes the cumulative flux in \f$\mathrm{[\frac{kg}{s}]}\f$, \f$\mathrm{[\frac{mole}{s}]}\f$ or \f$\mathrm{[\frac{J}{s}]}\f$ of an individual pore on the boundary
185 *
186 * \param element The element
187 * \param considerScv A lambda function to decide whether to consider a scv or not
188 * \param verbose If set true, the fluxes at all individual SCVs are printed
189 */
190 template<class RestrictingFunction>
191 635260 void computeBoundaryFlux_(const Element& element,
192 RestrictingFunction considerScv,
193 const bool verbose = false) const
194 {
195 // by default, all coordinate directions are considered for the definition of a boundary
196
197 // make sure FVElementGeometry and volume variables are bound to the element
198 1905780 const auto fvGeometry = localView(localResidual_.problem().gridGeometry()).bind(element);
199
4/10
✓ Branch 1 taken 318007 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 318007 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 318007 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 318007 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
1905780 const auto curElemVolVars = localView(gridVariables_.curGridVolVars()).bind(element, fvGeometry, sol_);
200
201
3/8
✗ Branch 0 not taken.
✓ Branch 1 taken 318007 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 318007 times.
✓ Branch 4 taken 318007 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
1905780 auto prevElemVolVars = localView(gridVariables_.prevGridVolVars());
202
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 318007 times.
635260 if (!isStationary_)
203 prevElemVolVars.bindElement(element, fvGeometry, sol_);
204
205
2/4
✓ Branch 1 taken 318007 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 318007 times.
1270520 ElementBoundaryTypes elemBcTypes;
206
1/2
✓ Branch 1 taken 318007 times.
✗ Branch 2 not taken.
635260 elemBcTypes.update(localResidual_.problem(), element, fvGeometry);
207
6/16
✓ Branch 1 taken 318007 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 318007 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 318007 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 318007 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 318007 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 318007 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
3176300 const auto elemFluxVarsCache = localView(gridVariables_.gridFluxVarsCache()).bindElement(element, fvGeometry, curElemVolVars);
208
1/2
✓ Branch 1 taken 318007 times.
✗ Branch 2 not taken.
635260 auto residual = localResidual_.evalFluxAndSource(element, fvGeometry, curElemVolVars, elemFluxVarsCache, elemBcTypes);
209
210
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 318007 times.
635260 if (!isStationary_)
211 residual += localResidual_.evalStorage(element, fvGeometry, prevElemVolVars, curElemVolVars);
212
213
2/2
✓ Branch 0 taken 636014 times.
✓ Branch 1 taken 318007 times.
1905780 for (auto&& scv : scvs(fvGeometry))
214 {
215 // compute the boundary flux using the local residual of the element's scv on the boundary
216
2/2
✓ Branch 1 taken 13854 times.
✓ Branch 2 taken 622160 times.
1270520 if (considerScv(scv))
217 {
218 55018 isConsidered_[scv.dofIndex()] = true;
219
220 // get the type of the boundary condition on the scv
221 27509 const auto& bcTypes = elemBcTypes.get(fvGeometry, scv);
222
223 27509 NumEqVector flux(0.0);
224
2/2
✓ Branch 0 taken 14365 times.
✓ Branch 1 taken 13854 times.
56040 for (int eqIdx = 0; eqIdx < NumEqVector::dimension; ++eqIdx)
225 {
226 // Check the type of the boundary condition.
227 // If BC is Dirichlet type, the flux is equal to the local residual of the element's scv on the boundary.
228 // Otherwise the flux is either zero or equal to a source term at the element's scv on the boundary.
229 // Technicaly, the PNM considers source terms instead of Neumann BCs.
230
2/2
✓ Branch 0 taken 762 times.
✓ Branch 1 taken 13603 times.
28531 if (!bcTypes.isDirichlet(eqIdx))
231 {
232
1/2
✓ Branch 1 taken 762 times.
✗ Branch 2 not taken.
1372 auto source = localResidual_.computeSource(localResidual_.problem(), element, fvGeometry, curElemVolVars, scv);
233 2744 source *= scv.volume() * curElemVolVars[scv].extrusionFactor();
234 1832 flux[eqIdx] = source[eqIdx];
235 }
236 else
237 57946 flux[eqIdx] = residual[scv.indexInElement()][eqIdx];
238 }
239
240 // The flux must be subtracted:
241 // On an inlet boundary, the flux part of the local residual will be positive, since all fluxes will leave the SCV towards to interior domain.
242 // For the domain itself, however, the sign has to be negative, since mass is entering the system.
243
4/4
✓ Branch 0 taken 194 times.
✓ Branch 1 taken 13149 times.
✓ Branch 2 taken 194 times.
✓ Branch 3 taken 13149 times.
55018 boundaryFluxes_[scv.dofIndex()] -= flux;
244
245
2/2
✓ Branch 0 taken 194 times.
✓ Branch 1 taken 13660 times.
27509 if (verbose)
246
4/8
✓ Branch 2 taken 194 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 194 times.
✗ Branch 7 not taken.
✓ Branch 10 taken 194 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 194 times.
✗ Branch 14 not taken.
776 std::cout << "SCV of element " << scv.elementIndex() << " at vertex " << scv.dofIndex() << " has flux: " << flux << std::endl;
247 }
248 }
249 635260 }
250
251 const LocalResidual localResidual_; // store a copy of the local residual
252 const GridVariables& gridVariables_;
253 const SolutionVector& sol_;
254 bool isStationary_;
255 mutable std::vector<bool> isConsidered_;
256 mutable std::vector<NumEqVector> boundaryFluxes_;
257 };
258
259 } // end Dumux::PoreNetwork
260
261 #endif
262