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