GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/freeflow/navierstokes/momentum/cvfe/flux.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 55 67 82.1%
Functions: 36 108 33.3%
Branches: 46 54 85.2%
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::NavierStokesFluxVariablesImpl
11 */
12 #ifndef DUMUX_NAVIERSTOKES_MOMENTUM_CVFE_FLUXVARIABLES_HH
13 #define DUMUX_NAVIERSTOKES_MOMENTUM_CVFE_FLUXVARIABLES_HH
14
15 #include <dune/common/fmatrix.hh>
16
17 #include <dumux/common/math.hh>
18 #include <dumux/common/exceptions.hh>
19 #include <dumux/common/parameters.hh>
20
21 #include <dumux/discretization/extrusion.hh>
22
23 namespace Dumux {
24
25 /*!
26 * \ingroup NavierStokesModel
27 * \brief Context for computing fluxes
28 *
29 * \tparam Problem the problem type to solve
30 * \tparam FVElementGeometry the element geometry type
31 * \tparam ElementVolumeVariables the element volume variables type
32 * \tparam ElementFluxVariablesCache the element flux variables cache type
33 */
34 template<class Problem,
35 class FVElementGeometry,
36 class ElementVolumeVariables,
37 class ElementFluxVariablesCache>
38 class NavierStokesMomentumFluxContext
39 {
40 using Element = typename FVElementGeometry::Element;
41 using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
42 public:
43
44 //! Initialize the flux variables storing some temporary pointers
45 42520552 NavierStokesMomentumFluxContext(
46 const Problem& problem,
47 const FVElementGeometry& fvGeometry,
48 const ElementVolumeVariables& elemVolVars,
49 const ElementFluxVariablesCache& elemFluxVarsCache,
50 const SubControlVolumeFace& scvf
51 )
52 : problem_(problem)
53 , fvGeometry_(fvGeometry)
54 , elemVolVars_(elemVolVars)
55 , elemFluxVarsCache_(elemFluxVarsCache)
56 42520552 , scvf_(scvf)
57 {}
58
59 ✗ const Problem& problem() const
60 ✗ { return problem_; }
61
62 ✗ const Element& element() const
63 228678288 { return fvGeometry_.element(); }
64
65 ✗ const SubControlVolumeFace& scvFace() const
66 ✗ { return scvf_; }
67
68 ✗ const FVElementGeometry& fvGeometry() const
69 ✗ { return fvGeometry_; }
70
71 ✗ const ElementVolumeVariables& elemVolVars() const
72 ✗ { return elemVolVars_; }
73
74 ✗ const ElementFluxVariablesCache& elemFluxVarsCache() const
75 ✗ { return elemFluxVarsCache_; }
76
77 private:
78 const Problem& problem_;
79 const FVElementGeometry& fvGeometry_;
80 const ElementVolumeVariables& elemVolVars_;
81 const ElementFluxVariablesCache& elemFluxVarsCache_;
82 const SubControlVolumeFace& scvf_;
83 };
84
85 /*!
86 * \ingroup NavierStokesModel
87 * \brief The flux variables class for the Navier-Stokes model using control-volume finite element schemes
88 */
89 template<class GridGeometry, class NumEqVector>
90 class NavierStokesMomentumFluxCVFE
91 {
92 using GridView = typename GridGeometry::GridView;
93 using Element = typename GridView::template Codim<0>::Entity;
94 using Scalar = typename NumEqVector::value_type;
95
96 using Extrusion = Extrusion_t<GridGeometry>;
97
98 static constexpr int dim = GridView::dimension;
99 static constexpr int dimWorld = GridView::dimensionworld;
100
101 using Tensor = Dune::FieldMatrix<Scalar, dim, dimWorld>;
102 static_assert(NumEqVector::dimension == dimWorld, "Wrong dimension of velocity vector");
103
104 public:
105 /*!
106 * \brief Returns the diffusive momentum flux due to viscous forces
107 */
108 template<class Context>
109 42520552 NumEqVector advectiveMomentumFlux(const Context& context) const
110 {
111
2/2
✓ Branch 0 taken 13222512 times.
✓ Branch 1 taken 29298040 times.
 42520552 if (!context.problem().enableInertiaTerms())
112 23314672 return NumEqVector(0.0);
113
114 29298040 const auto& fvGeometry = context.fvGeometry();
115 29298040 const auto& elemVolVars = context.elemVolVars();
116 29298040 const auto& scvf = context.scvFace();
117 29298040 const auto& fluxVarCache = context.elemFluxVarsCache()[scvf];
118 29298040 const auto& shapeValues = fluxVarCache.shapeValues();
119
120 // interpolate velocity at scvf
121 29298040 NumEqVector v(0.0);
122
4/4
✓ Branch 0 taken 123910760 times.
✓ Branch 1 taken 29298040 times.
✓ Branch 2 taken 123910760 times.
✓ Branch 3 taken 29298040 times.
 306417600 for (const auto& scv : scvs(fvGeometry))
123 619553800 v.axpy(shapeValues[scv.indexInElement()][0], elemVolVars[scv].velocity());
124
125 // get density from the problem
126 58596080 const Scalar density = context.problem().density(context.element(), context.fvGeometry(), scvf);
127
128 29298040 const auto vn = v*scvf.unitOuterNormal();
129
4/4
✓ Branch 0 taken 13517626 times.
✓ Branch 1 taken 15780414 times.
✓ Branch 2 taken 13517626 times.
✓ Branch 3 taken 15780414 times.
 58596080 const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
130
4/4
✓ Branch 0 taken 13517626 times.
✓ Branch 1 taken 15780414 times.
✓ Branch 2 taken 13517626 times.
✓ Branch 3 taken 15780414 times.
 58596080 const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
131
2/2
✓ Branch 0 taken 13517626 times.
✓ Branch 1 taken 15780414 times.
 29298040 const auto upwindVelocity = vn > 0 ? insideVolVars.velocity() : outsideVolVars.velocity();
132
2/2
✓ Branch 0 taken 13517626 times.
✓ Branch 1 taken 15780414 times.
 29298040 const auto downwindVelocity = vn > 0 ? outsideVolVars.velocity() : insideVolVars.velocity();
133
5/8
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 29298037 times.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 3 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 3 times.
✗ Branch 10 not taken.
 29298040 static const auto upwindWeight = getParamFromGroup<Scalar>(context.problem().paramGroup(), "Flux.UpwindWeight");
134 117192160 const auto advectiveTermIntegrand = density*vn * (upwindWeight * upwindVelocity + (1.0-upwindWeight)*downwindVelocity);
135
136 117192160 return advectiveTermIntegrand * Extrusion::area(fvGeometry, scvf) * insideVolVars.extrusionFactor();
137 }
138
139 /*!
140 * \brief Returns the diffusive momentum flux due to viscous forces
141 */
142 template<class Context>
143 42520552 NumEqVector diffusiveMomentumFlux(const Context& context) const
144 {
145 42520552 const auto& element = context.element();
146 42520552 const auto& fvGeometry = context.fvGeometry();
147 42520552 const auto& elemVolVars = context.elemVolVars();
148 42520552 const auto& scvf = context.scvFace();
149 42520552 const auto& fluxVarCache = context.elemFluxVarsCache()[scvf];
150
151 // interpolate velocity gradient at scvf
152 42520552 Tensor gradV(0.0);
153
4/4
✓ Branch 0 taken 186782768 times.
✓ Branch 1 taken 42520552 times.
✓ Branch 2 taken 186782768 times.
✓ Branch 3 taken 42520552 times.
 271823872 for (const auto& scv : scvs(fvGeometry))
154 {
155 186782768 const auto& volVars = elemVolVars[scv];
156
2/2
✓ Branch 0 taken 390397544 times.
✓ Branch 1 taken 186782768 times.
 577180312 for (int dir = 0; dir < dim; ++dir)
157 1951987720 gradV[dir].axpy(volVars.velocity(dir), fluxVarCache.gradN(scv.indexInElement()));
158 }
159
160 // get viscosity from the problem
161
2/2
✓ Branch 0 taken 34 times.
✓ Branch 1 taken 6334366 times.
 42520552 const auto mu = context.problem().effectiveViscosity(element, fvGeometry, scvf);
162
163
4/4
✓ Branch 0 taken 43 times.
✓ Branch 1 taken 42520509 times.
✓ Branch 3 taken 27 times.
✓ Branch 4 taken 16 times.
 42520579 static const bool enableUnsymmetrizedVelocityGradient
164
2/4
✓ Branch 1 taken 27 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 27 times.
✗ Branch 5 not taken.
 81 = getParamFromGroup<bool>(context.problem().paramGroup(), "FreeFlow.EnableUnsymmetrizedVelocityGradient", false);
165
166 // compute -mu*gradV*n*dA
167
2/2
✓ Branch 0 taken 40347752 times.
✓ Branch 1 taken 2172800 times.
 42520552 NumEqVector diffusiveFlux = enableUnsymmetrizedVelocityGradient ?
168 43478104 mv(gradV, scvf.unitOuterNormal())
169 6768000 : mv(gradV + getTransposed(gradV),scvf.unitOuterNormal());
170
171 42520552 diffusiveFlux *= -mu;
172
173
6/8
✓ Branch 0 taken 30 times.
✓ Branch 1 taken 42520522 times.
✓ Branch 3 taken 27 times.
✓ Branch 4 taken 3 times.
✓ Branch 6 taken 27 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 27 times.
✗ Branch 10 not taken.
 42520552 static const bool enableDilatationTerm = getParamFromGroup<bool>(context.problem().paramGroup(), "FreeFlow.EnableDilatationTerm", false);
174
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 42520552 times.
 42520552 if (enableDilatationTerm)
175 ✗ diffusiveFlux += 2.0/3.0 * mu * trace(gradV) * scvf.unitOuterNormal();
176
177 170082208 diffusiveFlux *= Extrusion::area(fvGeometry, scvf) * elemVolVars[scvf.insideScvIdx()].extrusionFactor();
178 42520552 return diffusiveFlux;
179 }
180
181 template<class Context>
182 42520552 NumEqVector pressureContribution(const Context& context) const
183 {
184 42520552 const auto& element = context.element();
185 42520552 const auto& fvGeometry = context.fvGeometry();
186 42520552 const auto& elemVolVars = context.elemVolVars();
187 42520552 const auto& scvf = context.scvFace();
188
189 // The pressure force needs to take the extruded scvf area into account
190 42520552 const auto pressure = context.problem().pressure(element, fvGeometry, scvf);
191
192 // The pressure contribution calculated above might have a much larger numerical value compared to the viscous or inertial forces.
193 // This may lead to numerical inaccuracies due to loss of significance (cancellation) for the final residual value.
194 // In the end, we are only interested in a pressure gradient between the two relevant faces so we can
195 // subtract a constant reference value from the actual pressure contribution.
196 42520552 const auto referencePressure = context.problem().referencePressure();
197
198 42520552 NumEqVector pn(scvf.unitOuterNormal());
199 170082208 pn *= (pressure-referencePressure)*Extrusion::area(fvGeometry, scvf)*elemVolVars[scvf.insideScvIdx()].extrusionFactor();
200
201 42520552 return pn;
202 }
203 };
204
205 } // end namespace Dumux
206
207 #endif
208