GCC Code Coverage Report


Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/flux/cctpfa/fickslaw.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 62 73 84.9%
Functions: 119 349 34.1%
Branches: 50 78 64.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 CCTpfaFlux
10 * \brief Fick's law for cell-centered finite volume schemes with two-point flux approximation
11 */
12 #ifndef DUMUX_DISCRETIZATION_CC_TPFA_FICKS_LAW_HH
13 #define DUMUX_DISCRETIZATION_CC_TPFA_FICKS_LAW_HH
14
15 #include <dune/common/fvector.hh>
16
17 #include <dumux/common/parameters.hh>
18 #include <dumux/common/properties.hh>
19
20 #include <dumux/discretization/method.hh>
21 #include <dumux/discretization/extrusion.hh>
22 #include <dumux/discretization/cellcentered/tpfa/computetransmissibility.hh>
23 #include <dumux/flux/fickiandiffusioncoefficients.hh>
24
25 #include <dumux/flux/referencesystemformulation.hh>
26
27 namespace Dumux {
28
29 // forward declaration
30 template<class TypeTag, class DiscretizationMethod, ReferenceSystemFormulation referenceSystem>
31 class FicksLawImplementation;
32
33 /*!
34 * \ingroup CCTpfaFlux
35 * \brief Fick's law for cell-centered finite volume schemes with two-point flux approximation
36 */
37 template <class TypeTag, ReferenceSystemFormulation referenceSystem>
38 class FicksLawImplementation<TypeTag, DiscretizationMethods::CCTpfa, referenceSystem>
39 {
40 using Implementation = FicksLawImplementation<TypeTag, DiscretizationMethods::CCTpfa, referenceSystem>;
41 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
42 using Problem = GetPropType<TypeTag, Properties::Problem>;
43 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
44 using FVElementGeometry = typename GridGeometry::LocalView;
45 using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
46 using Extrusion = Extrusion_t<GridGeometry>;
47 using GridView = typename GridGeometry::GridView;
48 using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
49 using VolumeVariables = typename ElementVolumeVariables::VolumeVariables;
50 using Element = typename GridView::template Codim<0>::Entity;
51 using GridFluxVariablesCache = GetPropType<TypeTag, Properties::GridFluxVariablesCache>;
52 using ElementFluxVariablesCache = typename GridFluxVariablesCache::LocalView;
53 using FluxVariablesCache = typename GridFluxVariablesCache::FluxVariablesCache;
54 using BalanceEqOpts = GetPropType<TypeTag, Properties::BalanceEqOpts>;
55
56 using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
57 static const int dim = GridView::dimension;
58 static const int dimWorld = GridView::dimensionworld;
59 static const int numPhases = ModelTraits::numFluidPhases();
60 static const int numComponents = ModelTraits::numFluidComponents();
61
62 using ComponentFluxVector = Dune::FieldVector<Scalar, numComponents>;
63
64 //! Class that fills the cache corresponding to tpfa Fick's Law
65 class TpfaFicksLawCacheFiller
66 {
67 public:
68 //! Function to fill a TpfaFicksLawCache of a given scvf
69 //! This interface has to be met by any diffusion-related cache filler class
70 template<class FluxVariablesCacheFiller>
71 static void fill(FluxVariablesCache& scvfFluxVarsCache,
72 unsigned int phaseIdx, unsigned int compIdx,
73 const Problem& problem,
74 const Element& element,
75 const FVElementGeometry& fvGeometry,
76 const ElementVolumeVariables& elemVolVars,
77 const SubControlVolumeFace& scvf,
78 const FluxVariablesCacheFiller& fluxVarsCacheFiller)
79 {
80 477869064 scvfFluxVarsCache.updateDiffusion(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx, compIdx);
81 }
82 };
83
84 //! Class that caches the transmissibility
85 class TpfaFicksLawCache
86 {
87 public:
88 using Filler = TpfaFicksLawCacheFiller;
89
90 void updateDiffusion(const Problem& problem,
91 const Element& element,
92 const FVElementGeometry& fvGeometry,
93 const ElementVolumeVariables& elemVolVars,
94 const SubControlVolumeFace &scvf,
95 const unsigned int phaseIdx,
96 const unsigned int compIdx)
97 {
98 236336076 tij_[phaseIdx][compIdx] = Implementation::calculateTransmissibility(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx, compIdx);
99 }
100
101 const Scalar& diffusionTij(unsigned int phaseIdx, unsigned int compIdx) const
102
2/3
✓ Branch 0 taken 22398504 times.
✓ Branch 1 taken 4009824 times.
✗ Branch 2 not taken.
213786500 { return tij_[phaseIdx][compIdx]; }
103
104 private:
105 std::array< std::array<Scalar, numComponents>, numPhases> tij_;
106 };
107
108 public:
109 using DiscretizationMethod = DiscretizationMethods::CCTpfa;
110 //! state the discretization method this implementation belongs to
111 static constexpr DiscretizationMethod discMethod{};
112
113 //! Return the reference system
114 static constexpr ReferenceSystemFormulation referenceSystemFormulation()
115 { return referenceSystem; }
116
117 //! state the type for the corresponding cache and its filler
118 using Cache = TpfaFicksLawCache;
119
120 using DiffusionCoefficientsContainer = FickianDiffusionCoefficients<Scalar, numPhases, numComponents>;
121
122 /*!
123 * \brief Returns the diffusive fluxes of all components within
124 * a fluid phase across the given sub-control volume face.
125 * The computed fluxes are given in mole/s or kg/s, depending
126 * on the template parameter ReferenceSystemFormulation.
127 *
128 * \note This overload allows to explicitly specify the inside and outside volume variables
129 * which can be useful to evaluate diffusive fluxes at boundaries with given outside values.
130 * This only works if scvf.numOutsideScv() == 1.
131 */
132 static ComponentFluxVector flux(const Problem& problem,
133 const Element& element,
134 const FVElementGeometry& fvGeometry,
135 const VolumeVariables& insideVolVars,
136 const VolumeVariables& outsideVolVars,
137 const SubControlVolumeFace& scvf,
138 const int phaseIdx,
139 const ElementFluxVariablesCache& elemFluxVarsCache)
140 {
141 if constexpr (isMixedDimensional_)
142 if (scvf.numOutsideScv() != 1)
143 DUNE_THROW(Dune::Exception, "This flux overload requires scvf.numOutsideScv() == 1");
144
145 // helper lambda to get the outside mole or mass fraction
146 const auto getOutsideX = [&](const Scalar xInside, const Scalar tij, const int phaseIdx, const int compIdx)
147 {
148 return massOrMoleFraction(outsideVolVars, referenceSystem, phaseIdx, compIdx);
149 };
150
151 // helper lambda to get the averaged density at the scvf
152 const auto getRho = [&](const int phaseIdx, const Scalar rhoInside, const Scalar rhoOutside)
153 {
154 return 0.5*(rhoInside + rhoOutside);
155 };
156
157 return flux_(element, fvGeometry, insideVolVars, outsideVolVars, elemFluxVarsCache, scvf, phaseIdx, getOutsideX, getRho);
158 }
159
160
161 /*!
162 * \brief Returns the diffusive fluxes of all components within
163 * a fluid phase across the given sub-control volume face.
164 * The computed fluxes are given in mole/s or kg/s, depending
165 * on the template parameter ReferenceSystemFormulation.
166 */
167 14598246 static ComponentFluxVector flux(const Problem& problem,
168 const Element& element,
169 const FVElementGeometry& fvGeometry,
170 const ElementVolumeVariables& elemVolVars,
171 const SubControlVolumeFace& scvf,
172 const int phaseIdx,
173 const ElementFluxVariablesCache& elemFluxVarsCache)
174 {
175 // get inside/outside volume variables
176 29196492 const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
177 29196492 const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
178
179 // helper lambda to get the outside mole or mass fraction
180 17565312 const auto getOutsideX = [&](const Scalar xInside, const Scalar tij, const int phaseIdx, const int compIdx)
181 {
182
0/3
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
6926266 const Scalar massOrMoleFractionOutside = massOrMoleFraction(outsideVolVars, referenceSystem, phaseIdx, compIdx);
183 if constexpr (isMixedDimensional_)
184 {
185
6/8
✓ Branch 0 taken 82452 times.
✓ Branch 1 taken 2867610 times.
✓ Branch 2 taken 82452 times.
✓ Branch 3 taken 2867610 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 17004 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 17004 times.
5934132 return scvf.numOutsideScvs() == 1 ? massOrMoleFractionOutside
186 82452 : branchingFacetX_(problem, element, fvGeometry, elemVolVars,
187 82452 elemFluxVarsCache, scvf, xInside, tij, phaseIdx, compIdx);
188 }
189 else
190 return massOrMoleFractionOutside;
191 };
192
193 // helper lambda to get the averaged density at the scvf
194 5899136 const auto getRho = [&](const int phaseIdx, const Scalar rhoInside, const Scalar rhoOutside)
195 {
196 if constexpr (isMixedDimensional_)
197 {
198
3/12
✓ Branch 0 taken 2867610 times.
✓ Branch 1 taken 82452 times.
✓ Branch 2 taken 17004 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
2967066 return scvf.numOutsideScvs() == 1 ? 0.5*(rhoInside + rhoOutside)
199 82452 : branchingFacetDensity_(elemVolVars, scvf, phaseIdx, rhoInside);
200 }
201 else
202 214177346 return 0.5*(rhoInside + rhoOutside);
203 };
204
205
0/2
✗ Branch 1 not taken.
✗ Branch 2 not taken.
14598246 return flux_(element, fvGeometry, insideVolVars, outsideVolVars, elemFluxVarsCache, scvf, phaseIdx, getOutsideX, getRho);
206 }
207
208 //! compute diffusive transmissibilities
209 246694810 static Scalar calculateTransmissibility(const Problem& problem,
210 const Element& element,
211 const FVElementGeometry& fvGeometry,
212 const ElementVolumeVariables& elemVolVars,
213 const SubControlVolumeFace& scvf,
214 const int phaseIdx, const int compIdx)
215 {
216
217
218
2/2
✓ Branch 0 taken 66957314 times.
✓ Branch 1 taken 61221132 times.
246694810 const auto insideScvIdx = scvf.insideScvIdx();
219
2/2
✓ Branch 0 taken 66957314 times.
✓ Branch 1 taken 61221132 times.
246694810 const auto& insideScv = fvGeometry.scv(insideScvIdx);
220 246694810 const auto& insideVolVars = elemVolVars[insideScvIdx];
221 340573642 const auto getDiffCoeff = [&](const auto& vv)
222 {
223 using FluidSystem = typename ElementVolumeVariables::VolumeVariables::FluidSystem;
224 if constexpr (FluidSystem::isTracerFluidSystem())
225 116583472 return vv.effectiveDiffusionCoefficient(0, 0, compIdx);
226 else
227
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 15646472 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
72425552 return vv.effectiveDiffusionCoefficient(phaseIdx, FluidSystem::getMainComponent(phaseIdx), compIdx);
228 };
229
230 246694810 const auto insideDiffCoeff = getDiffCoeff(insideVolVars);
231
232 246694810 const Scalar ti = computeTpfaTransmissibility(fvGeometry, scvf, insideScv, insideDiffCoeff, insideVolVars.extrusionFactor());
233
234 // for the boundary (dirichlet) or at branching points we only need ti
235 Scalar tij;
236
6/6
✓ Branch 0 taken 234172717 times.
✓ Branch 1 taken 7282403 times.
✓ Branch 2 taken 439 times.
✓ Branch 3 taken 234172278 times.
✓ Branch 4 taken 439 times.
✓ Branch 5 taken 234172278 times.
246694810 if (scvf.boundary() || scvf.numOutsideScvs() > 1)
237 14757900 tij = Extrusion::area(fvGeometry, scvf)*ti;
238
239 // otherwise we compute a tpfa harmonic mean
240 else
241 {
242
2/2
✓ Branch 0 taken 61221132 times.
✓ Branch 1 taken 61221132 times.
239315860 const auto outsideScvIdx = scvf.outsideScvIdx();
243
2/2
✓ Branch 0 taken 61221132 times.
✓ Branch 1 taken 61221132 times.
239315860 const auto& outsideScv = fvGeometry.scv(outsideScvIdx);
244 239315860 const auto& outsideVolVars = elemVolVars[outsideScvIdx];
245 239315860 const auto outsideDiffCoeff = getDiffCoeff(outsideVolVars);
246
247 Scalar tj;
248 if constexpr (dim == dimWorld)
249 // assume the normal vector from outside is anti parallel so we save flipping a vector
250 234166616 tj = -1.0*computeTpfaTransmissibility(fvGeometry, scvf, outsideScv, outsideDiffCoeff, outsideVolVars.extrusionFactor());
251 else
252 5161888 tj = computeTpfaTransmissibility(fvGeometry, fvGeometry.flipScvf(scvf.index()),
253 outsideScv,
254 outsideDiffCoeff,
255 outsideVolVars.extrusionFactor());
256
257 // check if we are dividing by zero!
258
2/2
✓ Branch 0 taken 132579476 times.
✓ Branch 1 taken 101592802 times.
239315860 if (ti*tj <= 0.0)
259 tij = 0;
260 else
261 275446116 tij = Extrusion::area(fvGeometry, scvf)*(ti * tj)/(ti + tj);
262 }
263
264 246694810 return tij;
265 }
266
267 private:
268 template<class OutsideFractionHelper, class DensityHelper>
269 220944836 static ComponentFluxVector flux_(const Element& element,
270 const FVElementGeometry& fvGeometry,
271 const VolumeVariables& insideVolVars,
272 const VolumeVariables& outsideVolVars,
273 const ElementFluxVariablesCache& elemFluxVarsCache,
274 const SubControlVolumeFace& scvf,
275 const int phaseIdx,
276 const OutsideFractionHelper& getOutsideX,
277 const DensityHelper& getRho)
278 {
279 220944836 ComponentFluxVector componentFlux(0.0);
280
2/2
✓ Branch 0 taken 418203992 times.
✓ Branch 1 taken 207455502 times.
661085604 for (int compIdx = 0; compIdx < numComponents; compIdx++)
281 {
282 using FluidSystem = typename VolumeVariables::FluidSystem;
283 if constexpr (!FluidSystem::isTracerFluidSystem())
284
4/4
✓ Branch 0 taken 202282786 times.
✓ Branch 1 taken 208899714 times.
✓ Branch 2 taken 7475904 times.
✓ Branch 3 taken 3737952 times.
439291240 if (compIdx == FluidSystem::getMainComponent(phaseIdx))
285 continue;
286
287 // diffusion tensors are always solution dependent
288
3/4
✓ Branch 0 taken 20950878 times.
✓ Branch 1 taken 5457450 times.
✓ Branch 2 taken 20950878 times.
✗ Branch 3 not taken.
288366860 const Scalar tij = elemFluxVarsCache[scvf].diffusionTij(phaseIdx, compIdx);
289
290 // the inside and outside mass/mole fractions fractions
291
1/2
✓ Branch 0 taken 26408328 times.
✗ Branch 1 not taken.
230633746 const Scalar xInside = massOrMoleFraction(insideVolVars, referenceSystem, phaseIdx, compIdx);
292
1/2
✓ Branch 0 taken 7967424 times.
✗ Branch 1 not taken.
230633746 const Scalar xOutside = getOutsideX(xInside, tij, phaseIdx, compIdx);
293
294
2/2
✓ Branch 0 taken 10852038 times.
✓ Branch 1 taken 82452 times.
230633746 const Scalar rhoInside = massOrMolarDensity(insideVolVars, referenceSystem, phaseIdx);
295
2/2
✓ Branch 0 taken 10852038 times.
✓ Branch 1 taken 82452 times.
230633746 const Scalar rhoOutside = massOrMolarDensity(outsideVolVars, referenceSystem, phaseIdx);
296
297
2/2
✓ Branch 0 taken 2884614 times.
✓ Branch 1 taken 82452 times.
455368356 const Scalar rho = getRho(phaseIdx, rhoInside, rhoOutside);
298
299 232613346 componentFlux[compIdx] = rho*tij*(xInside - xOutside);
300 if constexpr (!FluidSystem::isTracerFluidSystem())
301 218570362 if (BalanceEqOpts::mainComponentIsBalanced(phaseIdx))
302 654618396 componentFlux[FluidSystem::getMainComponent(phaseIdx)] -= componentFlux[compIdx];
303 }
304
305 220944836 return componentFlux;
306 }
307
308 //! compute the mole/mass fraction at branching facets for network grids
309 82452 static Scalar branchingFacetX_(const Problem& problem,
310 const Element& element,
311 const FVElementGeometry& fvGeometry,
312 const ElementVolumeVariables& elemVolVars,
313 const ElementFluxVariablesCache& elemFluxVarsCache,
314 const SubControlVolumeFace& scvf,
315 const Scalar insideX, const Scalar insideTi,
316 const int phaseIdx, const int compIdx)
317 {
318 82452 Scalar sumTi(insideTi);
319 82452 Scalar sumXTi(insideTi*insideX);
320
321
4/4
✓ Branch 0 taken 208944 times.
✓ Branch 1 taken 82452 times.
✓ Branch 2 taken 208944 times.
✓ Branch 3 taken 82452 times.
500340 for (unsigned int i = 0; i < scvf.numOutsideScvs(); ++i)
322 {
323 208944 const auto outsideScvIdx = scvf.outsideScvIdx(i);
324 208944 const auto& outsideVolVars = elemVolVars[outsideScvIdx];
325 208944 const Scalar massOrMoleFractionOutside = massOrMoleFraction(outsideVolVars, referenceSystem, phaseIdx, compIdx);
326 208944 const auto& flippedScvf = fvGeometry.flipScvf(scvf.index(), i);
327
328 417888 const Scalar outsideTi = elemFluxVarsCache[flippedScvf].diffusionTij(phaseIdx, compIdx);
329 208944 sumTi += outsideTi;
330 208944 sumXTi += outsideTi*massOrMoleFractionOutside;
331 }
332
333
1/2
✓ Branch 0 taken 82452 times.
✗ Branch 1 not taken.
82452 return sumTi > 0 ? sumXTi/sumTi : 0;
334 }
335
336 //! compute the density at branching facets for network grids as arithmetic mean
337 static Scalar branchingFacetDensity_(const ElementVolumeVariables& elemVolVars,
338 const SubControlVolumeFace& scvf,
339 const int phaseIdx,
340 const Scalar insideRho)
341 {
342 Scalar rho(insideRho);
343 for (unsigned int i = 0; i < scvf.numOutsideScvs(); ++i)
344 {
345 const auto outsideScvIdx = scvf.outsideScvIdx(i);
346 const auto& outsideVolVars = elemVolVars[outsideScvIdx];
347 const Scalar rhoOutside = massOrMolarDensity(outsideVolVars, referenceSystem, phaseIdx);
348 rho += rhoOutside;
349 }
350 return rho/(scvf.numOutsideScvs()+1);
351 }
352
353 static constexpr bool isMixedDimensional_ = static_cast<int>(GridView::dimension) < static_cast<int>(GridView::dimensionworld);
354 };
355
356 } // end namespace Dumux
357
358 #endif
359