Nektar++
LinearisedAdvection.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File LinearisedAdvection.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: Evaluation of the linearised advective term
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
37 #include <StdRegions/StdSegExp.h>
38 
48 
49 
50 namespace Nektar
51 {
54  "Linearised",
56 
57 /**
58  * Constructor. Creates ...
59  *
60  * \param
61  * \param
62  */
63 
65  Advection()
66 {
67 }
68 
69 
73 {
74  Advection::v_InitObject(pSession, pFields);
75 
76  m_session = pSession;
78  ::AllocateSharedPtr(m_session, pFields[0]->GetGraph());
79  m_spacedim = pFields[0]->GetGraph()->GetSpaceDimension();
80  m_expdim = pFields[0]->GetGraph()->GetMeshDimension();
82 
83  //Setting parameters for homogeneous problems
84  m_HomoDirec = 0;
85  m_useFFT = false;
87  m_SingleMode = false;
88  m_HalfMode = false;
89  m_MultipleModes = false;
90  m_homogen_dealiasing = false;
91 
92  if(m_session->DefinesSolverInfo("HOMOGENEOUS"))
93  {
94  std::string HomoStr = m_session->GetSolverInfo("HOMOGENEOUS");
95  m_spacedim = 3;
96 
97  if((HomoStr == "HOMOGENEOUS1D")||(HomoStr == "Homogeneous1D")||
98  (HomoStr == "1D")||(HomoStr == "Homo1D"))
99  {
101  m_LhomZ = m_session->GetParameter("LZ");
102  m_HomoDirec = 1;
103 
104  m_homogen_dealiasing = pSession->DefinesSolverInfo("DEALIASING");
105 
106  if(m_session->DefinesSolverInfo("ModeType"))
107  {
108  m_session->MatchSolverInfo("ModeType","SingleMode",m_SingleMode,false);
109  m_session->MatchSolverInfo("ModeType","HalfMode",m_HalfMode,false);
110  m_session->MatchSolverInfo("ModeType","MultipleModes",m_MultipleModes,false);
111  }
112 
113  if(m_session->DefinesSolverInfo("ModeType"))
114  {
115  if(m_SingleMode)
116  {
117  m_npointsZ=2;
118 
119  }
120  else if(m_HalfMode)
121  {
122  m_npointsZ=1;
123 
124  }
125  else if(m_MultipleModes)
126  {
127  m_npointsZ = m_session->GetParameter("HomModesZ");
128  }
129  else
130  {
131  ASSERTL0(false, "SolverInfo ModeType not valid");
132 
133 
134  }
135  }
136  else
137  {
138  m_session->LoadParameter("HomModesZ",m_npointsZ);
139 
140  }
141 
142  }
143 
144  if((HomoStr == "HOMOGENEOUS2D")||(HomoStr == "Homogeneous2D")||
145  (HomoStr == "2D")||(HomoStr == "Homo2D"))
146  {
148  m_session->LoadParameter("HomModesY", m_npointsY);
149  m_session->LoadParameter("LY", m_LhomY);
150  m_session->LoadParameter("HomModesZ", m_npointsZ);
151  m_session->LoadParameter("LZ", m_LhomZ);
152  m_HomoDirec = 2;
153  }
154 
155  if((HomoStr == "HOMOGENEOUS3D")||(HomoStr == "Homogeneous3D")||
156  (HomoStr == "3D")||(HomoStr == "Homo3D"))
157  {
159  m_session->LoadParameter("HomModesX",m_npointsX);
160  m_session->LoadParameter("LX", m_LhomX );
161  m_session->LoadParameter("HomModesY",m_npointsY);
162  m_session->LoadParameter("LY", m_LhomY );
163  m_session->LoadParameter("HomModesZ",m_npointsZ);
164  m_session->LoadParameter("LZ", m_LhomZ );
165  m_HomoDirec = 3;
166  }
167 
168  if(m_session->DefinesSolverInfo("USEFFT"))
169  {
170  m_useFFT = true;
171  }
172  }
173  else
174  {
175  m_npointsZ = 1; // set to default value so can use to identify 2d or 3D (homogeneous) expansions
176  }
177 
178  if(m_session->DefinesSolverInfo("PROJECTION"))
179  {
180  std::string ProjectStr
181  = m_session->GetSolverInfo("PROJECTION");
182 
183  if((ProjectStr == "Continuous")||(ProjectStr == "Galerkin")||
184  (ProjectStr == "CONTINUOUS")||(ProjectStr == "GALERKIN"))
185  {
187  }
188  else if(ProjectStr == "DisContinuous")
189  {
191  }
192  else
193  {
194  ASSERTL0(false,"PROJECTION value not recognised");
195  }
196  }
197  else
198  {
199  cerr << "Projection type not specified in SOLVERINFO,"
200  "defaulting to continuous Galerkin" << endl;
202  }
203 
204  int nvar = m_session->GetVariables().size();
206  for (int i = 0; i < nvar; ++i)
207  {
208  m_baseflow[i] = Array<OneD, NekDouble>(pFields[i]->GetTotPoints(), 0.0);
209  }
210 
211  ASSERTL0(m_session->DefinesFunction("BaseFlow"),
212  "Base flow must be defined for linearised forms.");
213  string file = m_session->GetFunctionFilename("BaseFlow", 0);
214 
215 
216  //Periodic base flows
217  if(m_session->DefinesParameter("N_slices"))
218  {
219  m_session->LoadParameter("N_slices",m_slices);
220  if(m_slices>1)
221  {
222  DFT(file,pFields,m_slices);
223  }
224  else
225  {
226  ASSERTL0(false,"Number of slices must be a positive number greater than 1");
227  }
228  }
229  //Steady base-flow
230  else
231  {
232  m_slices=1;
233 
234  //BaseFlow from file
235  if (m_session->GetFunctionType("BaseFlow", m_session->GetVariable(0))
237  {
238  ImportFldBase(file,pFields,1);
239 
240  }
241  //analytic base flow
242  else
243  {
244  int nq = pFields[0]->GetNpoints();
245  Array<OneD,NekDouble> x0(nq);
246  Array<OneD,NekDouble> x1(nq);
247  Array<OneD,NekDouble> x2(nq);
248 
249  // get the coordinates (assuming all fields have the same
250  // discretisation)
251  pFields[0]->GetCoords(x0,x1,x2);
252 
253  for(unsigned int i = 0 ; i < pFields.num_elements(); i++)
254  {
256  = m_session->GetFunction("BaseFlow", i);
257 
258  ifunc->Evaluate(x0,x1,x2,m_baseflow[i]);
259  }
260  }
261  }
262 
263  if(m_session->DefinesParameter("period"))
264  {
265  m_period=m_session->GetParameter("period");
266  }
267  else
268  {
269  m_period=(m_session->GetParameter("TimeStep")*m_slices)/(m_slices-1.);
270  }
271 
272 }
273 
275 {
276 }
277 
278 
279 //Advection function
280 
282  const int nConvectiveFields,
284  const Array<OneD, Array<OneD, NekDouble> > &advVel,
285  const Array<OneD, Array<OneD, NekDouble> > &inarray,
286  Array<OneD, Array<OneD, NekDouble> > &outarray,
287  const NekDouble &time)
288 {
289  int nqtot = fields[0]->GetTotPoints();
290  ASSERTL1(nConvectiveFields == inarray.num_elements(),"Number of convective fields and Inarray are not compatible");
291 
292  Array<OneD, NekDouble > Deriv = Array<OneD, NekDouble> (nqtot*nConvectiveFields);
293 
294  for(int n = 0; n < nConvectiveFields; ++n)
295  {
296  int ndim = advVel.num_elements();
297  int nPointsTot = fields[0]->GetNpoints();
298 
299  Array<OneD, NekDouble> grad0,grad1,grad2;
300 
301  //Evaluation of the gradiend of each component of the base flow
302  //\nabla U
303  Array<OneD, NekDouble> grad_base_u0,grad_base_u1,grad_base_u2;
304 
305  // \nabla V
306  Array<OneD, NekDouble> grad_base_v0,grad_base_v1,grad_base_v2;
307 
308  // \nabla W
309  Array<OneD, NekDouble> grad_base_w0,grad_base_w1,grad_base_w2;
310 
311  grad0 = Array<OneD, NekDouble> (nPointsTot);
312  grad_base_u0 = Array<OneD, NekDouble> (nPointsTot);
313  grad_base_v0 = Array<OneD, NekDouble> (nPointsTot);
314  grad_base_w0 = Array<OneD, NekDouble> (nPointsTot);
315 
316  //Evaluation of the base flow for periodic cases
317  //(it requires fld files)
318  if(m_slices>1)
319  {
320  if (m_session->GetFunctionType("BaseFlow", 0)
322  {
323  for(int i=0; i<ndim;++i)
324  {
326  }
327  }
328  else
329  {
330  ASSERTL0(false, "Periodic Base flow requires filename_ files");
331  }
332  }
333 
334 
335  //Evaluate the linearised advection term
336  switch(ndim)
337  {
338  // 1D
339  case 1:
340  fields[0]->PhysDeriv(inarray[n],grad0);
341  fields[0]->PhysDeriv(m_baseflow[0],grad_base_u0);
342  //Evaluate U du'/dx
343  Vmath::Vmul(nPointsTot,grad0,1,m_baseflow[0],1,outarray[n],1);
344  //Evaluate U du'/dx+ u' dU/dx
345  Vmath::Vvtvp(nPointsTot,grad_base_u0,1,advVel[0],1,outarray[n],1,outarray[n],1);
346  break;
347 
348  //2D
349  case 2:
350  grad1 = Array<OneD, NekDouble> (nPointsTot);
351  grad_base_u1 = Array<OneD, NekDouble> (nPointsTot);
352  grad_base_v1 = Array<OneD, NekDouble> (nPointsTot);
353  fields[0]->PhysDeriv(inarray[n],grad0,grad1);
354 
355  //Derivates of the base flow
356  fields[0]-> PhysDeriv(m_baseflow[0], grad_base_u0, grad_base_u1);
357  fields[0]-> PhysDeriv(m_baseflow[1], grad_base_v0, grad_base_v1);
358 
359  //Since the components of the velocity are passed one by
360  //one, it is necessary to distinguish which term is
361  //consider
362  switch (n)
363  {
364  //x-equation
365  case 0:
366  // Evaluate U du'/dx
367  Vmath::Vmul (nPointsTot,grad0,1,m_baseflow[0],1,outarray[n],1);
368  //Evaluate U du'/dx+ V du'/dy
369  Vmath::Vvtvp(nPointsTot,grad1,1,m_baseflow[1],1,outarray[n],1,outarray[n],1);
370  //Evaluate (U du'/dx+ V du'/dy)+u' dU/dx
371  Vmath::Vvtvp(nPointsTot,grad_base_u0,1,advVel[0],1,outarray[n],1,outarray[n],1);
372  //Evaluate (U du'/dx+ V du'/dy +u' dU/dx)+v' dU/dy
373  Vmath::Vvtvp(nPointsTot,grad_base_u1,1,advVel[1],1,outarray[n],1,outarray[n],1);
374  break;
375 
376  //y-equation
377  case 1:
378  // Evaluate U dv'/dx
379  Vmath::Vmul (nPointsTot,grad0,1,m_baseflow[0],1,outarray[n],1);
380  //Evaluate U dv'/dx+ V dv'/dy
381  Vmath::Vvtvp(nPointsTot,grad1,1,m_baseflow[1],1,outarray[n],1,outarray[n],1);
382  //Evaluate (U dv'/dx+ V dv'/dy)+u' dV/dx
383  Vmath::Vvtvp(nPointsTot,grad_base_v0,1,advVel[0],1,outarray[n],1,outarray[n],1);
384  //Evaluate (U dv'/dx+ V dv'/dy +u' dv/dx)+v' dV/dy
385  Vmath::Vvtvp(nPointsTot,grad_base_v1,1,advVel[1],1,outarray[n],1,outarray[n],1);
386  break;
387  }
388  break;
389 
390  //3D
391  case 3:
392  {
393  grad1 = Array<OneD, NekDouble> (nPointsTot);
394  grad2 = Array<OneD, NekDouble> (nPointsTot);
395  grad_base_u1 = Array<OneD, NekDouble> (nPointsTot);
396  grad_base_v1 = Array<OneD, NekDouble> (nPointsTot);
397  grad_base_w1 = Array<OneD, NekDouble> (nPointsTot);
398 
399  grad_base_u2 = Array<OneD, NekDouble> (nPointsTot);
400  grad_base_v2 = Array<OneD, NekDouble> (nPointsTot);
401  grad_base_w2 = Array<OneD, NekDouble> (nPointsTot);
402 
403  // Differentiate base flow in physical space
404  bool oldwavespace = fields[0]->GetWaveSpace();
405  fields[0]->SetWaveSpace(false);
406  fields[0]->PhysDeriv(m_baseflow[0],
407  grad_base_u0, grad_base_u1, grad_base_u2);
408  fields[0]->PhysDeriv(m_baseflow[1],
409  grad_base_v0, grad_base_v1, grad_base_v2);
410  fields[0]->PhysDeriv(m_baseflow[2],
411  grad_base_w0, grad_base_w1, grad_base_w2);
412  fields[0]->SetWaveSpace(oldwavespace);
413 
414  //HalfMode has W(x,y,t)=0
415  if(m_HalfMode || m_SingleMode)
416  {
417  for(int i=0; i<grad_base_u2.num_elements();++i)
418  {
419  grad_base_u2[i]=0;
420  grad_base_v2[i]=0;
421  grad_base_w2[i]=0;
422  }
423  }
424 
425  fields[0]->PhysDeriv(inarray[n], grad0, grad1, grad2);
426  if (oldwavespace)
427  {
428  fields[0]->HomogeneousBwdTrans(grad0, grad0);
429  fields[0]->HomogeneousBwdTrans(grad1, grad1);
430  fields[0]->HomogeneousBwdTrans(grad2, grad2);
431  }
432 
433  switch (n)
434  {
435  //x-equation
436  case 0:
438  {
439  //U du'/dx
440  fields[0]->DealiasedProd(m_baseflow[0],grad0,grad0,m_CoeffState);
441 
442  //V du'/dy
443  fields[0]->DealiasedProd(m_baseflow[1],grad1,grad1,m_CoeffState);
444 
445  //W du'/dx
446  fields[0]->DealiasedProd(m_baseflow[2],grad2,grad2,m_CoeffState);
447 
448  // u' dU/dx
449  fields[0]->DealiasedProd(advVel[0],grad_base_u0,grad_base_u0,m_CoeffState);
450  // v' dU/dy
451  fields[0]->DealiasedProd(advVel[1],grad_base_u1,grad_base_u1,m_CoeffState);
452  // w' dU/dz
453  fields[0]->DealiasedProd(advVel[2],grad_base_u2,grad_base_u2,m_CoeffState);
454 
455  Vmath::Vadd(nPointsTot,grad0,1,grad1,1,outarray[n],1);
456  Vmath::Vadd(nPointsTot,grad2,1,outarray[n],1,outarray[n],1);
457  Vmath::Vadd(nPointsTot,grad_base_u0,1,outarray[n],1,outarray[n],1);
458  Vmath::Vadd(nPointsTot,grad_base_u1,1,outarray[n],1,outarray[n],1);
459  Vmath::Vadd(nPointsTot,grad_base_u2,1,outarray[n],1,outarray[n],1);
460  }
461  else
462  {
463  //Evaluate U du'/dx
464  Vmath::Vmul (nPointsTot, grad0, 1, m_baseflow[0], 1,
465  outarray[n], 1);
466  //Evaluate U du'/dx+ V du'/dy
467  Vmath::Vvtvp(nPointsTot, grad1, 1, m_baseflow[1], 1,
468  outarray[n], 1, outarray[n], 1);
469  //Evaluate (U du'/dx+ V du'/dy)+u' dU/dx
470  Vmath::Vvtvp(nPointsTot,grad_base_u0,1,advVel[0],1,outarray[n],1,outarray[n],1);
471  //Evaluate (U du'/dx+ V du'/dy +u' dU/dx)+v' dU/dy
472  Vmath::Vvtvp(nPointsTot,grad_base_u1,1,advVel[1],1,outarray[n],1,outarray[n],1);
473  //Evaluate (U du'/dx+ V du'/dy +u' dU/dx +v' dU/dy) + W du'/dz
474  Vmath::Vvtvp(nPointsTot, grad2, 1, m_baseflow[2], 1,
475  outarray[n], 1, outarray[n], 1);
476  //Evaluate (U du'/dx+ V du'/dy +u' dU/dx +v' dU/dy + W du'/dz)+ w' dU/dz
477  Vmath::Vvtvp(nPointsTot,grad_base_u2,1,advVel[2],1,outarray[n],1,outarray[n],1);
478  }
479  break;
480  //y-equation
481  case 1:
483  {
484  //U dv'/dx
485  fields[0]->DealiasedProd(m_baseflow[0],grad0,grad0,m_CoeffState);
486  //V dv'/dy
487  fields[0]->DealiasedProd(m_baseflow[1],grad1,grad1,m_CoeffState);
488  //W dv'/dx
489  fields[0]->DealiasedProd(m_baseflow[2],grad2,grad2,m_CoeffState);
490  // u' dV/dx
491  fields[0]->DealiasedProd(advVel[0],grad_base_v0,grad_base_v0,m_CoeffState);
492  // v' dV/dy
493  fields[0]->DealiasedProd(advVel[1],grad_base_v1,grad_base_v1,m_CoeffState);
494  // w' dV/dz
495  fields[0]->DealiasedProd(advVel[2],grad_base_v2,grad_base_v2,m_CoeffState);
496 
497  Vmath::Vadd(nPointsTot,grad0,1,grad1,1,outarray[n],1);
498  Vmath::Vadd(nPointsTot,grad2,1,outarray[n],1,outarray[n],1);
499  Vmath::Vadd(nPointsTot,grad_base_v0,1,outarray[n],1,outarray[n],1);
500  Vmath::Vadd(nPointsTot,grad_base_v1,1,outarray[n],1,outarray[n],1);
501  Vmath::Vadd(nPointsTot,grad_base_v2,1,outarray[n],1,outarray[n],1);
502  }
503  else
504  {
505  //Evaluate U dv'/dx
506  Vmath::Vmul (nPointsTot, grad0, 1, m_baseflow[0], 1,
507  outarray[n], 1);
508  //Evaluate U dv'/dx+ V dv'/dy
509  Vmath::Vvtvp(nPointsTot, grad1, 1, m_baseflow[1], 1,
510  outarray[n], 1, outarray[n], 1);
511  //Evaluate (U dv'/dx+ V dv'/dy)+u' dV/dx
512  Vmath::Vvtvp(nPointsTot,grad_base_v0,1,advVel[0],1,outarray[n],1,outarray[n],1);
513  //Evaluate (U du'/dx+ V du'/dy +u' dV/dx)+v' dV/dy
514  Vmath::Vvtvp(nPointsTot,grad_base_v1,1,advVel[1],1,outarray[n],1,outarray[n],1);
515  //Evaluate (U du'/dx+ V dv'/dy +u' dV/dx +v' dV/dy) + W du'/dz
516  Vmath::Vvtvp(nPointsTot, grad2, 1, m_baseflow[2], 1,
517  outarray[n], 1, outarray[n], 1);
518  //Evaluate (U du'/dx+ V dv'/dy +u' dV/dx +v' dV/dy + W dv'/dz)+ w' dV/dz
519  Vmath::Vvtvp(nPointsTot,grad_base_v2,1,advVel[2],1,outarray[n],1,outarray[n],1);
520  }
521  break;
522 
523  //z-equation
524  case 2:
526  {
527  //U dw'/dx
528  fields[0]->DealiasedProd(m_baseflow[0],grad0,grad0,m_CoeffState);
529  //V dw'/dy
530  fields[0]->DealiasedProd(m_baseflow[1],grad1,grad1,m_CoeffState);
531  //W dw'/dx
532  fields[0]->DealiasedProd(m_baseflow[2],grad2,grad2,m_CoeffState);
533  // u' dW/dx
534  fields[0]->DealiasedProd(advVel[0],grad_base_w0,grad_base_w0,m_CoeffState);
535  // v' dW/dy
536  fields[0]->DealiasedProd(advVel[1],grad_base_w1,grad_base_w1,m_CoeffState);
537  // w' dW/dz
538  fields[0]->DealiasedProd(advVel[2],grad_base_w2,grad_base_w2,m_CoeffState);
539 
540  Vmath::Vadd(nPointsTot,grad0,1,grad1,1,outarray[n],1);
541  Vmath::Vadd(nPointsTot,grad2,1,outarray[n],1,outarray[n],1);
542  Vmath::Vadd(nPointsTot,grad_base_w0,1,outarray[n],1,outarray[n],1);
543  Vmath::Vadd(nPointsTot,grad_base_w1,1,outarray[n],1,outarray[n],1);
544  Vmath::Vadd(nPointsTot,grad_base_w2,1,outarray[n],1,outarray[n],1);
545  }
546  else
547  {
548  //Evaluate U dw'/dx
549  Vmath::Vmul (nPointsTot, grad0, 1, m_baseflow[0], 1,
550  outarray[n], 1);
551  //Evaluate U dw'/dx+ V dw'/dx
552  Vmath::Vvtvp(nPointsTot, grad1, 1, m_baseflow[1], 1,
553  outarray[n], 1, outarray[n], 1);
554  //Evaluate (U dw'/dx+ V dw'/dx)+u' dW/dx
555  Vmath::Vvtvp(nPointsTot,grad_base_w0,1,advVel[0],1,outarray[n],1,outarray[n],1);
556  //Evaluate (U dw'/dx+ V dw'/dx +w' dW/dx)+v' dW/dy
557  Vmath::Vvtvp(nPointsTot,grad_base_w1,1,advVel[1],1,outarray[n],1,outarray[n],1);
558  //Evaluate (U dw'/dx+ V dw'/dx +u' dW/dx +v' dW/dy) + W dw'/dz
559  Vmath::Vvtvp(nPointsTot, grad2, 1, m_baseflow[2], 1,
560  outarray[n], 1, outarray[n], 1);
561  //Evaluate (U dw'/dx+ V dw'/dx +u' dW/dx +v' dW/dy + W dw'/dz)+ w' dW/dz
562  Vmath::Vvtvp(nPointsTot,grad_base_w2,1,advVel[2],1,outarray[n],1,outarray[n],1);
563  }
564  break;
565  }
566 
567  if (oldwavespace)
568  {
569  fields[0]->HomogeneousFwdTrans(outarray[n],outarray[n]);
570  }
571  break;
572  }
573  default:
574  ASSERTL0(false,"dimension unknown");
575  }
576 
577  Vmath::Neg(nqtot,outarray[n],1);
578  }
579 }
580 
582  const Array<OneD, Array<OneD, NekDouble> > &inarray)
583 {
584  if (m_session->GetSolverInfo("EqType") == "UnsteadyNavierStokes")
585  {
586  // The SFD method is only applied to the velocity variables in
587  // incompressible
588  ASSERTL1(inarray.num_elements() == (m_baseflow.num_elements() - 1),
589  "Number of base flow variables does not match what is "
590  "expected.");
591  }
592  else
593  {
594  ASSERTL1(inarray.num_elements() == (m_baseflow.num_elements()),
595  "Number of base flow variables does not match what is expected.");
596  }
597 
598  int npts = inarray[0].num_elements();
599 
600  for (int i = 0; i < inarray.num_elements(); ++i)
601  {
602  ASSERTL1(npts == m_baseflow[i].num_elements(),
603  "Size of base flow array does not match expected.");
604  Vmath::Vcopy(npts, inarray[i], 1, m_baseflow[i], 1);
605  }
606 }
607 
608 
609 /**
610  * Import field from infile and load into \a m_fields. This routine will
611  * also perform a \a BwdTrans to ensure data is in both the physical and
612  * coefficient storage.
613  * @param infile Filename to read.
614  */
615 void LinearisedAdvection::ImportFldBase(std::string pInfile,
617 {
618  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
619  std::vector<std::vector<NekDouble> > FieldData;
620 
621  int nqtot = m_baseflow[0].num_elements();
622  int nvar = m_session->GetVariables().size();
623  int s;
624  Array<OneD, NekDouble> tmp_coeff(pFields[0]->GetNcoeffs(), 0.0);
625 
626  //Get Homogeneous
629  m_session->GetComm());
630  fld->Import(pInfile, FieldDef, FieldData);
631 
632 
633  if(m_session->DefinesSolverInfo("HOMOGENEOUS"))
634  {
635  std::string HomoStr = m_session->GetSolverInfo("HOMOGENEOUS");
636  }
637 
638  // copy FieldData into m_fields
639  for(int j = 0; j < nvar; ++j)
640  {
641  for(int i = 0; i < FieldDef.size(); ++i)
642  {
643  if((m_session->DefinesSolverInfo("HOMOGENEOUS") &&
644  (m_session->GetSolverInfo("HOMOGENEOUS")=="HOMOGENEOUS1D" ||
645  m_session->GetSolverInfo("HOMOGENEOUS")=="1D" ||
646  m_session->GetSolverInfo("HOMOGENEOUS")=="Homo1D")) &&
647  m_MultipleModes==false)
648  {
649  // w-component must be ignored and set to zero.
650  if (j != nvar - 2)
651  {
652  // p component (it is 4th variable of the 3D and corresponds 3nd variable of 2D)
653  s = (j == nvar - 1) ? 2 : j;
654 
655  //extraction of the 2D
656  pFields[j]->ExtractDataToCoeffs(
657  FieldDef[i],
658  FieldData[i],
659  FieldDef[i]->m_fields[s],
660  tmp_coeff);
661 
662  }
663 
664  //Put zero on higher modes
665  int ncplane = (pFields[0]->GetNcoeffs()) / m_npointsZ;
666 
667  if (m_npointsZ > 2)
668  {
669  Vmath::Zero(ncplane*(m_npointsZ-2),
670  &tmp_coeff[2*ncplane], 1);
671  }
672  }
673  // 2D cases and Homogeneous1D Base Flows
674  else
675  {
676  bool flag = FieldDef[i]->m_fields[j] ==
677  m_session->GetVariable(j);
678 
679  ASSERTL0(flag, (std::string("Order of ") + pInfile
680  + std::string(" data and that defined in "
681  "m_boundaryconditions differs")).c_str());
682 
683  pFields[j]->ExtractDataToCoeffs(FieldDef[i], FieldData[i],
684  FieldDef[i]->m_fields[j],
685  tmp_coeff);
686  }
687  }
688 
689  if(m_SingleMode || m_HalfMode)
690  {
691  //pFields[j]->SetWaveSpace(true);
692 
693  pFields[j]->GetPlane(0)->BwdTrans(tmp_coeff, m_baseflow[j]);
694 
695  if(m_SingleMode)
696  {
697  //copy the bwd into the second plane for single Mode Analysis
698  int ncplane=(pFields[0]->GetNpoints())/m_npointsZ;
699  Vmath::Vcopy(ncplane,&m_baseflow[j][0],1,&m_baseflow[j][ncplane],1);
700  }
701  }
702  else
703  {
704  bool oldwavespace = pFields[j]->GetWaveSpace();
705  pFields[j]->SetWaveSpace(false);
706  pFields[j]->BwdTrans(tmp_coeff, m_baseflow[j]);
707  pFields[j]->SetWaveSpace(oldwavespace);
708  }
709  }
710 
711  if(m_session->DefinesParameter("N_slices"))
712  {
713  int nConvectiveFields = pFields.num_elements()-1;
714 
715  for(int i=0; i<nConvectiveFields;++i)
716  {
717 
718  Vmath::Vcopy(nqtot, &m_baseflow[i][0], 1, &m_interp[i][slice*nqtot], 1);
719  }
720 
721  }
722 }
723 
724 
726  const Array<OneD, const NekDouble> &inarray,
727  Array<OneD, NekDouble> &outarray,
728  const NekDouble m_time,
729  const NekDouble m_period)
730 {
731  int npoints=m_baseflow[0].num_elements();
732 
733  NekDouble BetaT=2*M_PI*fmod (m_time, m_period) / m_period;
734  NekDouble phase;
735  Array<OneD, NekDouble> auxiliary(npoints);
736 
737  Vmath::Vcopy(npoints,&inarray[0],1,&outarray[0],1);
738  Vmath::Svtvp(npoints, cos(0.5*m_slices*BetaT),&inarray[npoints],1,&outarray[0],1,&outarray[0],1);
739 
740  for (int i = 2; i < m_slices; i += 2)
741  {
742  phase = (i>>1) * BetaT;
743 
744  Vmath::Svtvp(npoints, cos(phase),&inarray[i*npoints],1,&outarray[0],1,&outarray[0],1);
745  Vmath::Svtvp(npoints, sin(phase), &inarray[(i+1)*npoints], 1, &outarray[0], 1,&outarray[0],1);
746  }
747 
748 }
749 
751 {
752  DNekMatSharedPtr loc_mat;
753  DNekBlkMatSharedPtr BlkMatrix;
754  int n_exp = 0;
755 
756  n_exp = m_baseflow[0].num_elements(); // will operatore on m_phys
757 
758  Array<OneD,unsigned int> nrows(n_exp);
759  Array<OneD,unsigned int> ncols(n_exp);
760 
761  nrows = Array<OneD, unsigned int>(n_exp,m_slices);
762  ncols = Array<OneD, unsigned int>(n_exp,m_slices);
763 
764  MatrixStorage blkmatStorage = eDIAGONAL;
765  BlkMatrix = MemoryManager<DNekBlkMat>
766  ::AllocateSharedPtr(nrows,ncols,blkmatStorage);
767 
768 
771  StdRegions::StdSegExp StdSeg(BK);
772 
774  StdSeg.DetShapeType(),
775  StdSeg);
776 
777  loc_mat = StdSeg.GetStdMatrix(matkey);
778 
779  // set up array of block matrices.
780  for(int i = 0; i < n_exp; ++i)
781  {
782  BlkMatrix->SetBlock(i,i,loc_mat);
783  }
784 
785  return BlkMatrix;
786 }
787 
788 //Discrete Fourier Transform for Floquet analysis
789 void LinearisedAdvection::DFT(const string file,
791  const NekDouble m_slices)
792 {
793  int npoints=m_baseflow[0].num_elements();
794 
795  //Convected fields
796  int ConvectedFields=m_baseflow.num_elements()-1;
797 
798  m_interp= Array<OneD, Array<OneD, NekDouble> > (ConvectedFields);
799  for(int i=0; i<ConvectedFields;++i)
800  {
802  }
803 
804  // Import the slides into auxiliary vector
805  // The base flow should be stored in the form "filename_%d.ext"
806  // A subdirectory can also be included, such as "dir/filename_%d.ext"
807  size_t found = file.find("%d");
808  ASSERTL0(found != string::npos && file.find("%d", found+1) == string::npos,
809  "Since N_slices is specified, the filename provided for function "
810  "'BaseFlow' must include exactly one instance of the format "
811  "specifier '%d', to index the time-slices.");
812  char* buffer = new char[file.length() + 8];
813  for (int i = 0; i < m_slices; ++i)
814  {
815  sprintf(buffer, file.c_str(), i);
816  ImportFldBase(buffer,pFields,i);
817  }
818  delete[] buffer;
819 
820 
821  // Discrete Fourier Transform of the fields
822  for(int k=0; k< ConvectedFields;++k)
823  {
824 #ifdef NEKTAR_USING_FFTW
825 
826  //Discrete Fourier Transform using FFTW
827  Array<OneD, NekDouble> fft_in(npoints*m_slices);
828  Array<OneD, NekDouble> fft_out(npoints*m_slices);
829 
832 
833  //Shuffle the data
834  for(int j= 0; j < m_slices; ++j)
835  {
836  Vmath::Vcopy(npoints,&m_interp[k][j*npoints],1,&(fft_in[j]),m_slices);
837  }
838 
839  m_FFT = LibUtilities::GetNektarFFTFactory().CreateInstance("NekFFTW", m_slices);
840 
841  //FFT Transform
842  for(int i=0; i<npoints; i++)
843  {
844  m_FFT->FFTFwdTrans(m_tmpIN =fft_in + i*m_slices, m_tmpOUT =fft_out + i*m_slices);
845 
846  }
847 
848  //Reshuffle data
849  for(int s = 0; s < m_slices; ++s)
850  {
851  Vmath::Vcopy(npoints,&fft_out[s],m_slices,&m_interp[k][s*npoints],1);
852 
853  }
854 
855  Vmath::Zero(fft_in.num_elements(),&fft_in[0],1);
856  Vmath::Zero(fft_out.num_elements(),&fft_out[0],1);
857 #else
858  //Discrete Fourier Transform using MVM
859  DNekBlkMatSharedPtr blkmat;
861 
862  int nrows = blkmat->GetRows();
863  int ncols = blkmat->GetColumns();
864 
865  Array<OneD, NekDouble> sortedinarray(ncols);
866  Array<OneD, NekDouble> sortedoutarray(nrows);
867 
868  //Shuffle the data
869  for(int j= 0; j < m_slices; ++j)
870  {
871  Vmath::Vcopy(npoints,&m_interp[k][j*npoints],1,&(sortedinarray[j]),m_slices);
872  }
873 
874  // Create NekVectors from the given data arrays
875  NekVector<NekDouble> in (ncols,sortedinarray,eWrapper);
876  NekVector<NekDouble> out(nrows,sortedoutarray,eWrapper);
877 
878  // Perform matrix-vector multiply.
879  out = (*blkmat)*in;
880 
881  //Reshuffle data
882  for(int s = 0; s < m_slices; ++s)
883  {
884  Vmath::Vcopy(npoints,&sortedoutarray[s],m_slices,&m_interp[k][s*npoints],1);
885  }
886 
887  for(int r=0; r<sortedinarray.num_elements();++r)
888  {
889  sortedinarray[0]=0;
890  sortedoutarray[0]=0;
891  }
892 
893 #endif
894 
895  //scaling of the Fourier coefficients
896  NekDouble j=-1;
897  for (int i = 2; i < m_slices; i += 2)
898  {
899  Vmath::Smul(2*npoints,j,&m_interp[k][i*npoints],1,&m_interp[k][i*npoints],1);
900  j=-j;
901 
902  }
903 
904  }
905 
906 }
907 
908 } //end of namespace
909 
enum HomogeneousType m_HomogeneousType
void DFT(const string file, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const NekDouble m_slices)
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:454
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Array< OneD, NekDouble > m_tmpIN
Local coefficients.
int m_HomoDirec
number of homogenous directions
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:162
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
Array< OneD, Array< OneD, NekDouble > > m_baseflow
Storage for base flow.
MultiRegions::ProjectionType m_projectionType
LibUtilities::NektarFFTSharedPtr m_FFT
void ImportFldBase(std::string pInfile, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, int slice)
Import Base flow.
void Svtvp(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:471
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:428
array buffer
Definition: GsLib.hpp:56
Array< OneD, NekDouble > m_tmpOUT
Fourier Expansion .
Definition: BasisType.h:52
boost::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:70
NektarFFTFactory & GetNektarFFTFactory()
Definition: NektarFFT.cpp:69
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:50
virtual void v_SetBaseFlow(const Array< OneD, Array< OneD, NekDouble > > &inarray)
Overrides the base flow used during linearised advection.
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:684
Array< OneD, Array< OneD, NekDouble > > m_interp
Class representing a segment element in reference space.
Definition: StdSegExp.h:54
void UpdateBase(const NekDouble m_slices, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const NekDouble m_time, const NekDouble m_period)
1D Evenly-spaced points using Fourier Fit
Definition: PointsType.h:64
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:199
virtual void v_Advect(const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &advVel, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble &time)
Advects a vector field.
bool m_HalfMode
flag to determine if use half mode or not
int m_npointsZ
number of points in Z direction (if homogeneous)
static std::string npts
Definition: InputFld.cpp:43
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Definition: Advection.cpp:46
Defines a specification for a set of points.
Definition: Points.h:58
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:382
boost::shared_ptr< FieldIO > FieldIOSharedPtr
Definition: FieldIO.h:236
double NekDouble
DNekBlkMatSharedPtr GetFloquetBlockMatrix(FloquetMatType mattype, bool UseContCoeffs=false) const
MultiRegions::CoeffState m_CoeffState
int m_npointsY
number of points in Y direction (if homogeneous)
boost::shared_ptr< Equation > EquationSharedPtr
bool m_SingleMode
flag to determine if use single mode or not
static SolverUtils::AdvectionSharedPtr create(std::string)
Creates an instance of this class.
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
NekDouble m_LhomY
physical length in Y direction (if homogeneous)
bool m_MultipleModes
flag to determine if use multiple mode or not
LibUtilities::SessionReaderSharedPtr m_session
boost::shared_ptr< DNekBlkMat > DNekBlkMatSharedPtr
Definition: NekTypeDefs.hpp:72
int m_npointsX
number of points in X direction (if homogeneous)
NekDouble m_LhomX
physical length in X direction (if homogeneous)
static std::string className
Name of class.
virtual void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Initialises the advection object.
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
NekDouble m_LhomZ
physical length in Z direction (if homogeneous)
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:191
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1038
bool m_useFFT
flag to determine if use or not the FFT for transformations
Describes the specification for a Basis.
Definition: Basis.h:50
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:285
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