Nektar++
FilterReynoldsStresses.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: FilterReynoldsStresses.cpp
4//
5// For more information, please see: http://www.nektar.info
6//
7// The MIT License
8//
9// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10// Department of Aeronautics, Imperial College London (UK), and Scientific
11// Computing and Imaging Institute, University of Utah (USA).
12//
13// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description: Append Reynolds stresses to the average fields
32//
33///////////////////////////////////////////////////////////////////////////////
34
36
37namespace Nektar::SolverUtils
38{
39
41 GetFilterFactory().RegisterCreatorFunction("ReynoldsStresses",
43
44/**
45 * @class FilterReynoldsStresses
46 *
47 * @brief Append Reynolds stresses to the average fields
48 *
49 * This class appends the average fields with the Reynolds stresses of the form
50 * \f$ \overline{u' v'} \f$.
51 *
52 * For the default case, this is achieved by calculating
53 * \f$ C_{n} = \Sigma_{i=1}^{n} (u_i - \bar{u}_n)(v_i - \bar{v}_n)\f$
54 * using the recursive relation:
55 *
56 * \f[ C_{n} = C_{n-1} + \frac{n}{n-1} (u_n - \bar{u}_n)(v_n - \bar{v}_n) \f]
57 *
58 * The FilterSampler base class then divides the result by n, leading
59 * to the Reynolds stress.
60 *
61 * It is also possible to perform the averages using an exponential moving
62 * average, in which case either the moving average parameter \f$ \alpha \f$
63 * or the time constant \f$ \tau \f$ must be prescribed.
64 */
67 const std::shared_ptr<SolverUtils::EquationSystem> &pEquation,
68 const std::map<std::string, std::string> &pParams)
69 : FilterFieldConvert(pSession, pEquation, pParams)
70{
71 // Load sampling frequency
72 auto it = pParams.find("SampleFrequency");
73 if (it == pParams.end())
74 {
76 }
77 else
78 {
79 LibUtilities::Equation equ(m_session->GetInterpreter(), it->second);
80 m_sampleFrequency = round(equ.Evaluate());
81 }
82
83 // Check if should use moving average
84 it = pParams.find("MovingAverage");
85 if (it == pParams.end())
86 {
87 m_movAvg = false;
88 }
89 else
90 {
91 std::string sOption = it->second.c_str();
92 m_movAvg = (boost::iequals(sOption, "true")) ||
93 (boost::iequals(sOption, "yes"));
94 }
95
96 if (!m_movAvg)
97 {
98 return;
99 }
100
101 // Load alpha parameter for moving average
102 it = pParams.find("alpha");
103 if (it == pParams.end())
104 {
105 it = pParams.find("tau");
106 if (it == pParams.end())
107 {
108 ASSERTL0(false, "MovingAverage needs either alpha or tau.");
109 }
110 else
111 {
112 // Load time constant
113 LibUtilities::Equation equ(m_session->GetInterpreter(), it->second);
114 NekDouble tau = equ.Evaluate();
115 // Load delta T between samples
116 NekDouble dT;
117 m_session->LoadParameter("TimeStep", dT);
118 dT = dT * m_sampleFrequency;
119 // Calculate alpha
120 m_alpha = dT / (tau + dT);
121 }
122 }
123 else
124 {
125 LibUtilities::Equation equ(m_session->GetInterpreter(), it->second);
126 m_alpha = equ.Evaluate();
127 // Check if tau was also defined
128 it = pParams.find("tau");
129 if (it != pParams.end())
130 {
131 ASSERTL0(false,
132 "Cannot define both alpha and tau in MovingAverage.");
133 }
134 }
135 // Check bounds of m_alpha
136 ASSERTL0(m_alpha > 0 && m_alpha < 1, "Alpha out of bounds.");
137}
138
141 const NekDouble &time)
142{
143 size_t dim = pFields.size() - 1;
144 size_t nExtraFields = (dim + 1) * dim / 2;
145 size_t origFields = pFields.size();
146 size_t nqtot = pFields[0]->GetTotPoints();
147
148 // Allocate storage
149 m_fields.resize(origFields + nExtraFields);
150 m_delta.resize(dim);
151
152 for (size_t n = 0; n < m_fields.size(); ++n)
153 {
154 m_fields[n] = Array<OneD, NekDouble>(nqtot, 0.0);
155 }
156 for (size_t n = 0; n < m_delta.size(); ++n)
157 {
158 m_delta[n] = Array<OneD, NekDouble>(nqtot, 0.0);
159 }
160
161 // Initialise output arrays
163
164 // Update m_fields if using restart file
165 if (m_numSamples)
166 {
167 for (size_t j = 0; j < m_fields.size(); ++j)
168 {
169 pFields[0]->BwdTrans(m_outFields[j], m_fields[j]);
170 if (pFields[0]->GetWaveSpace())
171 {
172 pFields[0]->HomogeneousBwdTrans(nqtot, m_fields[j],
173 m_fields[j]);
174 }
175 }
176 }
177}
178
181{
182 size_t dim = pFields.size() - 1;
183 size_t origFields = pFields.size();
184
185 // Fill name of variables
186 for (size_t n = 0; n < origFields; ++n)
187 {
188 m_variables.push_back(pFields[n]->GetSession()->GetVariable(n));
189 }
190 for (int i = 0; i < dim; ++i)
191 {
192 for (int j = i; j < dim; ++j)
193 {
194 std::string var = pFields[i]->GetSession()->GetVariable(i) +
195 pFields[j]->GetSession()->GetVariable(j);
196 m_variables.push_back(var);
197 }
198 }
199}
200
203 [[maybe_unused]] std::vector<Array<OneD, NekDouble>> &fieldcoeffs,
204 [[maybe_unused]] const NekDouble &time)
205{
206 size_t i, j, n;
207 size_t nq = pFields[0]->GetTotPoints();
208 size_t dim = pFields.size() - 1;
209 bool waveSpace = pFields[0]->GetWaveSpace();
210 NekDouble nSamples = (NekDouble)m_numSamples;
211
212 // For moving average, take first sample as initial vector
213 NekDouble alpha = m_alpha;
214 if (m_numSamples == 1)
215 {
216 alpha = 1.0;
217 }
218
219 // Define auxiliary constants for averages
220 NekDouble facOld, facAvg, facStress, facDelta;
221 if (m_movAvg)
222 {
223 facOld = 1.0 - alpha;
224 facAvg = alpha;
225 facStress = alpha;
226 facDelta = 1.0;
227 }
228 else
229 {
230 facOld = 1.0;
231 facAvg = 1.0;
232 facStress = nSamples / (nSamples - 1);
233 facDelta = 1.0 / nSamples;
234 }
235
238
239 // Update original velocities in phys space and calculate (\bar{u} - u_n)
240 for (n = 0; n < dim; ++n)
241 {
242 if (waveSpace)
243 {
244 pFields[n]->HomogeneousBwdTrans(nq, pFields[n]->GetPhys(), vel);
245 }
246 else
247 {
248 vel = pFields[n]->GetPhys();
249 }
250 Vmath::Svtsvtp(nq, facAvg, vel, 1, facOld, m_fields[n], 1, m_fields[n],
251 1);
252 Vmath::Svtvm(nq, facDelta, m_fields[n], 1, vel, 1, m_delta[n], 1);
253 }
254 // Update pressure (directly to outFields)
255 Vmath::Svtsvtp(m_outFields[dim].size(), facAvg, pFields[dim]->GetCoeffs(),
256 1, facOld, m_outFields[dim], 1, m_outFields[dim], 1);
257
258 // Ignore Reynolds stress for first sample (its contribution is zero)
259 if (m_numSamples == 1)
260 {
261 return;
262 }
263
264 // Calculate C_{n} = facOld * C_{n-1} + facStress * deltaI * deltaJ
265 for (i = 0, n = dim + 1; i < dim; ++i)
266 {
267 for (j = i; j < dim; ++j, ++n)
268 {
269 Vmath::Vmul(nq, m_delta[i], 1, m_delta[j], 1, tmp, 1);
270 Vmath::Svtsvtp(nq, facStress, tmp, 1, facOld, m_fields[n], 1,
271 m_fields[n], 1);
272 }
273 }
274}
275
278 [[maybe_unused]] const NekDouble &time)
279{
280 size_t dim = pFields.size() - 1;
281
282 m_fieldMetaData["NumberOfFieldDumps"] = std::to_string(m_numSamples);
283
284 // Set wavespace to false, as calculations were performed in physical space
285 bool waveSpace = pFields[0]->GetWaveSpace();
286 pFields[0]->SetWaveSpace(false);
287
288 // Forward transform and put into m_outFields (except pressure)
289 for (size_t i = 0; i < m_fields.size(); ++i)
290 {
291 if (i != dim)
292 {
293 pFields[0]->FwdTransLocalElmt(m_fields[i], m_outFields[i]);
294 }
295 }
296
297 // Restore waveSpace
298 pFields[0]->SetWaveSpace(waveSpace);
299}
300
302{
303 if (m_movAvg)
304 {
305 return 1.0;
306 }
307 else
308 {
309 return 1.0 / m_numSamples;
310 }
311}
312
313} // namespace Nektar::SolverUtils
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
LibUtilities::FieldMetaDataMap m_fieldMetaData
std::vector< Array< OneD, NekDouble > > m_outFields
SOLVER_UTILS_EXPORT void v_Initialise(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time) override
LibUtilities::SessionReaderSharedPtr m_session
Definition: Filter.h:93
FilterReynoldsStresses(const LibUtilities::SessionReaderSharedPtr &pSession, const std::shared_ptr< SolverUtils::EquationSystem > &pEquation, const std::map< std::string, std::string > &pParams)
static std::string className
Name of the class.
std::vector< Array< OneD, NekDouble > > m_delta
void v_Initialise(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time) override
void v_FillVariablesName(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields) override
static FilterSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const std::shared_ptr< SolverUtils::EquationSystem > &pEquation, const std::map< std::string, std::string > &pParams)
Creates an instance of this class.
std::vector< Array< OneD, NekDouble > > m_fields
void v_PrepareOutput(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time) override
void v_ProcessSample(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, std::vector< Array< OneD, NekDouble > > &fieldcoeffs, const NekDouble &time) override
std::shared_ptr< SessionReader > SessionReaderSharedPtr
FilterFactory & GetFilterFactory()
double NekDouble
void Svtsvtp(int n, const T alpha, const T *x, int incx, const T beta, const T *y, int incy, T *z, int incz)
Svtsvtp (scalar times vector plus scalar times vector):
Definition: Vmath.hpp:473
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.hpp:72
void Svtvm(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Svtvm (scalar times vector minus vector): z = alpha*x - y.
Definition: Vmath.hpp:424