GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/freeflow/rans/oneeq/volumevariables.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 44 68 64.7%
Functions: 32 88 36.4%
Branches: 57 92 62.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 OneEqModel
10 *
11 * \copydoc Dumux::OneEqVolumeVariables
12 */
13 #ifndef DUMUX_ONEEQ_VOLUME_VARIABLES_HH
14 #define DUMUX_ONEEQ_VOLUME_VARIABLES_HH
15
16 #include <dune/common/math.hh>
17 #include <dumux/common/parameters.hh>
18 #include <dumux/freeflow/rans/volumevariables.hh>
19
20 namespace Dumux {
21
22 /*!
23 * \ingroup OneEqModel
24 * \brief Volume variables for the isothermal single-phase one-equation
25 * turbulence model by Spalart-Allmaras
26 */
27 template <class Traits, class NSVolumeVariables>
28
11/20
✓ Branch 3 taken 74676 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 25 times.
✓ Branch 7 taken 33032 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 25 times.
✓ Branch 11 taken 33032 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 25 times.
✗ Branch 15 not taken.
✓ Branch 18 taken 25 times.
✗ Branch 19 not taken.
✓ Branch 22 taken 25 times.
✗ Branch 23 not taken.
✓ Branch 26 taken 25 times.
✗ Branch 27 not taken.
✓ Branch 30 taken 25 times.
✗ Branch 31 not taken.
✓ Branch 34 taken 25 times.
✗ Branch 35 not taken.
 283260 class OneEqVolumeVariables
29 : public RANSVolumeVariables<Traits, NSVolumeVariables>
30 {
31 using RANSParentType = RANSVolumeVariables<Traits, NSVolumeVariables>;
32
33 using Scalar = typename Traits::PrimaryVariables::value_type;
34 using DimVector = Dune::FieldVector<Scalar, Traits::ModelTraits::dim()>;
35
36 public:
37 //! export the indices type
38 using Indices = typename Traits::ModelTraits::Indices;
39
40 /*!
41 * \brief Update all quantities for a given control volume
42 *
43 * \param elemSol A vector containing all primary variables connected to the element
44 * \param problem The object specifying the problem which ought to
45 * be simulated
46 * \param element An element which contains part of the control volume
47 * \param scv The sub-control volume
48 */
49 template<class ElementSolution, class Problem, class Element, class SubControlVolume>
50 10451052 void update(const ElementSolution &elemSol,
51 const Problem &problem,
52 const Element &element,
53 const SubControlVolume& scv)
54 {
55 10451052 RANSParentType::updateNavierStokesVolVars(elemSol, problem, element, scv);
56 10451052 updateRANSProperties(elemSol, problem, element, scv);
57 10451052 }
58
59 /*!
60 * \brief Update all turbulent quantities for a given control volume
61 *
62 * Wall and roughness related quantities are stored. Eddy viscosity is set.
63 *
64 * \param elemSol A vector containing all primary variables connected to the element
65 * \param problem The object specifying the problem which ought to be simulated
66 * \param element An element which contains part of the control volume
67 * \param scv The sub-control volume
68 */
69 template<class ElementSolution, class Problem, class Element, class SubControlVolume>
70 10451052 void updateRANSProperties(const ElementSolution &elemSol,
71 const Problem &problem,
72 const Element &element,
73 const SubControlVolume& scv)
74 {
75 10451052 RANSParentType::updateRANSProperties(elemSol, problem, element, scv);
76 10451052 viscosityTilde_ = elemSol[0][Indices::viscosityTildeIdx];
77 10451052 storedViscosityTilde_ = problem.storedViscosityTilde(RANSParentType::elementIdx());
78 10451052 storedViscosityTildeGradient_ = problem.storedViscosityTildeGradient(RANSParentType::elementIdx());
79 10451052 stressTensorScalarProduct_ = problem.stressTensorScalarProduct(RANSParentType::elementIdx());
80 10451052 vorticityTensorScalarProduct_ = problem.vorticityTensorScalarProduct(RANSParentType::elementIdx());
81
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 10450352 times.
 10451052 if (problem.useStoredEddyViscosity())
82 ✗ RANSParentType::setDynamicEddyViscosity_(problem.storedDynamicEddyViscosity(RANSParentType::elementIdx()));
83 else
84 20902104 RANSParentType::setDynamicEddyViscosity_(calculateEddyViscosity());
85 10451052 RANSParentType::calculateEddyDiffusivity(problem);
86 10451052 RANSParentType::calculateEddyThermalConductivity(problem);
87 10451052 }
88
89 /*!
90 * \brief Returns the dynamic eddy viscosity \f$\mathrm{[Pa s]}\f$.
91 */
92 Scalar calculateEddyViscosity()
93
4/8
✓ Branch 0 taken 33057 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 25 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 25 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 25 times.
✗ Branch 7 not taken.
 10483484 { return viscosityTilde() * fv1() * RANSParentType::density(); }
94
95 /*!
96 * \brief Returns the viscosity parameter \f$ m^2/s \f$
97 */
98 ✗ Scalar viscosityTilde() const
99 ✗ { return viscosityTilde_; }
100
101 /*!
102 * \brief Returns the viscosity parameter from the last iteration \f$ m^2/s \f$
103 */
104 Scalar storedViscosityTilde() const
105 { return storedViscosityTilde_; }
106
107 /*!
108 * \brief Returns the gradient of the viscosity parameter
109 */
110 ✗ DimVector storedViscosityTildeGradient() const
111 ✗ { return storedViscosityTildeGradient_; }
112
113 /*!
114 * \brief Returns the scalar product of the stress tensor
115 */
116 Scalar stressTensorScalarProduct() const
117 { return stressTensorScalarProduct_; }
118
119 /*!
120 * \brief Returns damping function for the eddy viscosity
121 */
122 Scalar fv1() const
123 {
124 153179728 return viscosityRatio() * viscosityRatio() * viscosityRatio()
125 153179728 / (viscosityRatio() * viscosityRatio() * viscosityRatio()
126 38294932 + cv1() * cv1() * cv1());
127 }
128
129 //! \brief Returns a model function
130 Scalar fv2() const
131 111245792 { return 1.0 - viscosityRatio() / (1.0 + viscosityRatio() * fv1()); }
132
133 //! \brief Returns a model function
134 ✗ Scalar ft2() const
135 {
136 using std::exp;
137 // the trip correction term is dropped according to Versteeg2009 and Wilcox2006
138 // return ct3() * exp(-ct4() * viscosityRatio() * viscosityRatio());
139 ✗ return 0.0;
140 }
141
142 //! \brief Returns a model function
143 3973064 Scalar fW() const
144 {
145 using std::pow;
146 using Dune::power;
147 7946128 return g() * pow(1.0 + power(cw3(), 6) / (power(g(), 6) + power(cw3(), 6)), 1.0/6.0);
148 }
149
150 //! \brief Returns a model function
151 7946128 Scalar g() const
152 {
153 using Dune::power;
154 15892256 return r() + cw2() * (power(r(), 6) - r());
155 }
156
157 //! \brief Returns a model function
158 23838384 Scalar r() const
159 {
160 using std::min;
161 71515152 return min(10.0,
162 23838384 viscosityTilde() / stressTensorScalarProductTilde()
163 23838384 / RANSParentType::karmanConstant() / RANSParentType::karmanConstant()
164
2/2
✓ Branch 0 taken 23821428 times.
✓ Branch 1 taken 16956 times.
 23838384 / RANSParentType::wallDistance() / RANSParentType::wallDistance());
165 }
166
167 //! \brief Returns the ratio of the kinematic viscosity and the viscosity parameter
168 Scalar viscosityRatio() const
169
35/56
✓ Branch 0 taken 150847 times.
✓ Branch 1 taken 27660626 times.
✓ Branch 2 taken 150847 times.
✓ Branch 3 taken 27660626 times.
✓ Branch 4 taken 150847 times.
✓ Branch 5 taken 27660626 times.
✓ Branch 6 taken 150847 times.
✓ Branch 7 taken 27660626 times.
✓ Branch 8 taken 150847 times.
✓ Branch 9 taken 27660626 times.
✓ Branch 10 taken 150847 times.
✓ Branch 11 taken 27660626 times.
✓ Branch 12 taken 150847 times.
✓ Branch 13 taken 27660626 times.
✓ Branch 14 taken 33057 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 33057 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 33057 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 33057 times.
✗ Branch 21 not taken.
✓ Branch 22 taken 33057 times.
✗ Branch 23 not taken.
✓ Branch 24 taken 33057 times.
✗ Branch 25 not taken.
✓ Branch 26 taken 33057 times.
✗ Branch 27 not taken.
✓ Branch 28 taken 25 times.
✗ Branch 29 not taken.
✓ Branch 30 taken 25 times.
✗ Branch 31 not taken.
✓ Branch 32 taken 25 times.
✗ Branch 33 not taken.
✓ Branch 34 taken 25 times.
✗ Branch 35 not taken.
✓ Branch 36 taken 25 times.
✗ Branch 37 not taken.
✓ Branch 38 taken 25 times.
✗ Branch 39 not taken.
✓ Branch 40 taken 25 times.
✗ Branch 41 not taken.
✓ Branch 42 taken 25 times.
✗ Branch 43 not taken.
✓ Branch 44 taken 25 times.
✗ Branch 45 not taken.
✓ Branch 46 taken 25 times.
✗ Branch 47 not taken.
✓ Branch 48 taken 25 times.
✗ Branch 49 not taken.
✓ Branch 50 taken 25 times.
✗ Branch 51 not taken.
✓ Branch 52 taken 25 times.
✗ Branch 53 not taken.
✓ Branch 54 taken 25 times.
✗ Branch 55 not taken.
 268064524 { return viscosityTilde() / RANSParentType::kinematicViscosity(); }
170
171 /*
172 * ! \brief Returns a modified version of the stress tensor scalar product
173 *
174 * According to <a href="https://turbmodels.larc.nasa.gov/spalart.html">NASA</a>
175 * this term should never be zero and different limiters might be used.
176 * The Implementation uses the one proposed in:
177 * Allmaras, S. R., Johnson, F. T., and Spalart, P. R.,
178 * "Modifications and Clarifications for the Implementation of the Spalart-Allmaras Turbulence Model," ICCFD7-1902
179 */
180 27811448 Scalar stressTensorScalarProductTilde() const
181 {
182 // original form
183 // return vorticityMagnitude()
184 // + viscosityTilde() * fv2()
185 // / RANSParentType::karmanConstant() / RANSParentType::karmanConstant()
186 // / RANSParentType::wallDistance() / RANSParentType::wallDistance();
187
188 // limiter form, literature source see above
189
2/2
✓ Branch 0 taken 150822 times.
✓ Branch 1 taken 27660626 times.
 27811448 Scalar sBar = viscosityTilde() * fv2()
190 27811448 / RANSParentType::karmanConstant() / RANSParentType::karmanConstant()
191 27811448 / RANSParentType::wallDistance() / RANSParentType::wallDistance();
192 55622896 return sBar < -c2() * vorticityMagnitude()
193 55622896 ? vorticityMagnitude()
194 603288 + (vorticityMagnitude() * (c2() * c2() * vorticityMagnitude() + c3() * sBar))
195 301644 / ((c3() - 2.0 * c2()) * vorticityMagnitude() - sBar)
196 55472074 : vorticityMagnitude() + sBar;
197 }
198
199 //! \brief Returns the magnitude of the vorticity
200 ✗ Scalar vorticityMagnitude() const
201 {
202 using std::sqrt;
203
2/2
✓ Branch 0 taken 150822 times.
✓ Branch 1 taken 27660626 times.
 27811448 return sqrt(2.0 * vorticityTensorScalarProduct_);
204 }
205
206 //! \brief Returns a model constant
207 ✗ Scalar c2() const
208 ✗ { return 0.7; }
209
210 //! \brief Returns a model constant
211 ✗ Scalar c3() const
212 ✗ { return 0.9; }
213
214 //! \brief Returns a model constant
215 ✗ Scalar sigma() const
216 ✗ { return 2.0/3.0; }
217
218 //! \brief Returns a model constant
219 ✗ Scalar cb1() const
220 ✗ { return 0.1355; }
221
222 //! \brief Returns a model constant
223 ✗ Scalar cb2() const
224 ✗ { return 0.622; }
225
226 //! \brief Returns a model constant
227 ✗ Scalar cv1() const
228 ✗ { return 7.1; }
229
230 //! \brief Returns a model constant
231 Scalar ct3() const
232 { return 1.2; }
233
234 //! \brief Returns a model constant
235 Scalar ct4() const
236 { return 0.5; }
237
238 //! \brief Returns a model constant
239 Scalar cw1() const
240 {
241 3973064 return cb1() / RANSParentType::karmanConstant() / RANSParentType::karmanConstant()
242 3973064 + (1.0 + cb2()) / sigma();
243 }
244
245 //! \brief Returns a model constant
246 ✗ Scalar cw2() const
247 ✗ { return 0.3; }
248
249 //! \brief Returns a model constant
250 ✗ Scalar cw3() const
251 ✗ { return 2.0; }
252
253 protected:
254 Scalar viscosityTilde_ = 0.0;
255 Scalar storedViscosityTilde_ = 0.0;
256 DimVector storedViscosityTildeGradient_ = DimVector(0.0);
257 Scalar stressTensorScalarProduct_;
258 Scalar vorticityTensorScalarProduct_;
259 };
260
261 } // end namespace Dumux
262
263 #endif
264