Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator 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  : FilterSampler(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  // Fill name of variables
143  for (int n = 0; n < origFields; ++n)
144  {
145  m_variables.push_back(pFields[n]->GetSession()->GetVariable(n));
146  }
147  if (dim == 2)
148  {
149  m_variables.push_back("uu");
150  m_variables.push_back("uv");
151  m_variables.push_back("vv");
152  }
153  else if (dim == 3)
154  {
155  m_variables.push_back("uu");
156  m_variables.push_back("uv");
157  m_variables.push_back("uw");
158  m_variables.push_back("vv");
159  m_variables.push_back("vw");
160  m_variables.push_back("ww");
161  }
162  else
163  {
164  ASSERTL0(false, "Unsupported dimension");
165  }
166 
167  // Allocate storage
168  m_fields.resize(origFields + nExtraFields);
169  m_delta.resize(dim);
170 
171  for (int n = 0; n < m_fields.size(); ++n)
172  {
173  m_fields[n] = Array<OneD, NekDouble>(pFields[0]->GetTotPoints(), 0.0);
174  }
175  for (int n = 0; n < m_delta.size(); ++n)
176  {
177  m_delta[n] = Array<OneD, NekDouble>(pFields[0]->GetTotPoints(), 0.0);
178  }
179 
180  // Initialise output arrays
181  FilterSampler::v_Initialise(pFields, time);
182 }
183 
186  const NekDouble &time)
187 {
188  int i, j, n;
189  int nq = pFields[0]->GetTotPoints();
190  int dim = pFields.num_elements() - 1;
191  bool waveSpace = pFields[0]->GetWaveSpace();
192  NekDouble nSamples = (NekDouble)m_numSamples;
193 
194  // For moving average, take first sample as initial vector
195  NekDouble alpha = m_alpha;
196  if (m_numSamples == 1)
197  {
198  alpha = 1.0;
199  }
200 
201  // Define auxiliary constants for averages
202  NekDouble facOld, facAvg, facStress, facDelta;
203  if (m_movAvg)
204  {
205  facOld = 1.0 - alpha;
206  facAvg = alpha;
207  facStress = alpha;
208  facDelta = 1.0;
209  }
210  else
211  {
212  facOld = 1.0;
213  facAvg = 1.0;
214  facStress = nSamples / (nSamples - 1);
215  facDelta = 1.0 / nSamples;
216  }
217 
218  Array<OneD, NekDouble> vel(nq);
219  Array<OneD, NekDouble> tmp(nq);
220 
221  // Update original velocities in phys space and calculate (\bar{u} - u_n)
222  for (n = 0; n < dim; ++n)
223  {
224  if (waveSpace)
225  {
226  pFields[n]->HomogeneousBwdTrans(pFields[n]->GetPhys(), vel);
227  }
228  else
229  {
230  vel = pFields[n]->GetPhys();
231  }
233  nq, facAvg, vel, 1, facOld, m_fields[n], 1, m_fields[n], 1);
234  Vmath::Svtvm(nq, facDelta, m_fields[n], 1, vel, 1, m_delta[n], 1);
235  }
236  // Update pressure (directly to outFields)
237  Vmath::Svtsvtp(m_outFields[dim].num_elements(),
238  facAvg,
239  pFields[dim]->GetCoeffs(),
240  1,
241  facOld,
242  m_outFields[dim],
243  1,
244  m_outFields[dim],
245  1);
246 
247  // Ignore Reynolds stress for first sample (its contribution is zero)
248  if (m_numSamples == 1)
249  {
250  return;
251  }
252 
253  // Calculate C_{n} = facOld * C_{n-1} + facStress * deltaI * deltaJ
254  for (i = 0, n = dim + 1; i < dim; ++i)
255  {
256  for (j = i; j < dim; ++j, ++n)
257  {
258  Vmath::Vmul(nq, m_delta[i], 1, m_delta[j], 1, tmp, 1);
260  nq, facStress, tmp, 1, facOld, m_fields[n], 1, m_fields[n], 1);
261  }
262  }
263 }
264 
267  const NekDouble &time)
268 {
269  int dim = pFields.num_elements() - 1;
270 
271  if (m_movAvg)
272  {
273  m_scale = 1.0;
274  }
275  else
276  {
277  m_scale = 1.0 / m_numSamples;
278  }
279 
280  // Set wavespace to false, as calculations were performed in physical space
281  bool waveSpace = pFields[0]->GetWaveSpace();
282  pFields[0]->SetWaveSpace(false);
283 
284  // Forward transform and put into m_outFields (except pressure)
285  for (int i = 0; i < m_fields.size(); ++i)
286  {
287  if (i != dim)
288  {
289  pFields[0]->FwdTrans_IterPerExp(m_fields[i], m_outFields[i]);
290  }
291  }
292 
293  // Restore waveSpace
294  pFields[0]->SetWaveSpace(waveSpace);
295 }
296 
298 {
299  return true;
300 }
301 }
302 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
std::vector< Array< OneD, NekDouble > > m_outFields
Definition: FilterSampler.h:90
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
std::vector< std::string > m_variables
Definition: FilterSampler.h:91
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:504
virtual SOLVER_UTILS_EXPORT void v_Initialise(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
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:577
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:169
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215