Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 
36 #include <IncNavierStokesSolver/AdvectionTerms/LinearisedAdvection.h>
37 #include <StdRegions/StdSegExp.h>
38 
39 #include <MultiRegions/ContField3DHomogeneous1D.h>
40 #include <MultiRegions/ContField3DHomogeneous2D.h>
41 #include <MultiRegions/ContField1D.h>
42 #include <MultiRegions/ContField2D.h>
43 #include <MultiRegions/ContField3D.h>
44 #include <MultiRegions/DisContField1D.h>
45 #include <MultiRegions/DisContField2D.h>
46 #include <MultiRegions/DisContField3DHomogeneous1D.h>
47 #include <MultiRegions/DisContField3DHomogeneous2D.h>
48 
49 
50 namespace Nektar
51 {
52 string LinearisedAdvection::className
53  = SolverUtils::GetAdvectionFactory().RegisterCreatorFunction(
54  "Linearised",
55  LinearisedAdvection::create);
56 
57 /**
58  * Constructor. Creates ...
59  *
60  * \param
61  * \param
62  */
63 
64 LinearisedAdvection::LinearisedAdvection():
65  Advection()
66 {
67 }
68 
69 
70 void LinearisedAdvection::v_InitObject(
71  LibUtilities::SessionReaderSharedPtr pSession,
72  Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
73 {
74  Advection::v_InitObject(pSession, pFields);
75 
76  m_session = pSession;
77  m_boundaryConditions = MemoryManager<SpatialDomains::BoundaryConditions>
78  ::AllocateSharedPtr(m_session, pFields[0]->GetGraph());
79  m_spacedim = pFields[0]->GetGraph()->GetSpaceDimension();
80  m_expdim = pFields[0]->GetGraph()->GetMeshDimension();
81  m_CoeffState = MultiRegions::eLocal;
82 
83  //Setting parameters for homogeneous problems
84  m_HomoDirec = 0;
85  m_useFFT = false;
86  m_HomogeneousType = eNotHomogeneous;
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  {
100  m_HomogeneousType = eHomogeneous1D;
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  {
147  m_HomogeneousType = eHomogeneous2D;
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  {
158  m_HomogeneousType = eHomogeneous3D;
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  {
186  m_projectionType = MultiRegions::eGalerkin;
187  }
188  else if(ProjectStr == "DisContinuous")
189  {
190  m_projectionType = MultiRegions::eDiscontinuous;
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;
201  m_projectionType = MultiRegions::eGalerkin;
202  }
203 
204  int nvar = m_session->GetVariables().size();
205  m_baseflow = Array<OneD, Array<OneD, NekDouble> >(nvar);
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))
236  == LibUtilities::eFunctionTypeFile)
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  {
255  LibUtilities::EquationSharedPtr ifunc
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 
274 LinearisedAdvection::~LinearisedAdvection()
275 {
276 }
277 
278 
279 //Advection function
280 
281 void LinearisedAdvection::v_Advect(
282  const int nConvectiveFields,
283  const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
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)
321  == LibUtilities::eFunctionTypeFile)
322  {
323  for(int i=0; i<ndim;++i)
324  {
325  UpdateBase(m_slices,m_interp[i],m_baseflow[i],time,m_period);
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  grad1 = Array<OneD, NekDouble> (nPointsTot);
393  grad2 = Array<OneD, NekDouble> (nPointsTot);
394  grad_base_u1 = Array<OneD, NekDouble> (nPointsTot);
395  grad_base_v1 = Array<OneD, NekDouble> (nPointsTot);
396  grad_base_w1 = Array<OneD, NekDouble> (nPointsTot);
397 
398  grad_base_u2 = Array<OneD, NekDouble> (nPointsTot);
399  grad_base_v2 = Array<OneD, NekDouble> (nPointsTot);
400  grad_base_w2 = Array<OneD, NekDouble> (nPointsTot);
401 
402  fields[0]->PhysDeriv(m_baseflow[0], grad_base_u0, grad_base_u1,grad_base_u2);
403  fields[0]->PhysDeriv(m_baseflow[1], grad_base_v0, grad_base_v1,grad_base_v2);
404  fields[0]->PhysDeriv(m_baseflow[2], grad_base_w0, grad_base_w1, grad_base_w2);
405 
406  //HalfMode has W(x,y,t)=0
407  if(m_HalfMode || m_SingleMode)
408  {
409  for(int i=0; i<grad_base_u2.num_elements();++i)
410  {
411  grad_base_u2[i]=0;
412  grad_base_v2[i]=0;
413  grad_base_w2[i]=0;
414  }
415  }
416 
417  fields[0]->PhysDeriv(inarray[n], grad0, grad1, grad2);
418 
419  switch (n)
420  {
421  //x-equation
422  case 0:
423  if(m_homogen_dealiasing)
424  {
425  //U du'/dx
426  fields[0]->DealiasedProd(m_baseflow[0],grad0,grad0,m_CoeffState);
427 
428  //V du'/dy
429  fields[0]->DealiasedProd(m_baseflow[1],grad1,grad1,m_CoeffState);
430 
431  //W du'/dx
432  fields[0]->DealiasedProd(m_baseflow[2],grad2,grad2,m_CoeffState);
433 
434  // u' dU/dx
435  fields[0]->DealiasedProd(advVel[0],grad_base_u0,grad_base_u0,m_CoeffState);
436  // v' dU/dy
437  fields[0]->DealiasedProd(advVel[1],grad_base_u1,grad_base_u1,m_CoeffState);
438  // w' dU/dz
439  fields[0]->DealiasedProd(advVel[2],grad_base_u2,grad_base_u2,m_CoeffState);
440 
441  Vmath::Vadd(nPointsTot,grad0,1,grad1,1,outarray[n],1);
442  Vmath::Vadd(nPointsTot,grad2,1,outarray[n],1,outarray[n],1);
443  Vmath::Vadd(nPointsTot,grad_base_u0,1,outarray[n],1,outarray[n],1);
444  Vmath::Vadd(nPointsTot,grad_base_u1,1,outarray[n],1,outarray[n],1);
445  Vmath::Vadd(nPointsTot,grad_base_u2,1,outarray[n],1,outarray[n],1);
446  }
447  else
448  {
449  //Evaluate U du'/dx
450  Vmath::Vmul (nPointsTot,grad0,1,m_baseflow[0],1,outarray[n],1);
451  //Evaluate U du'/dx+ V du'/dy
452  Vmath::Vvtvp(nPointsTot,grad1,1,m_baseflow[1],1,outarray[n],1,outarray[n],1);
453  //Evaluate (U du'/dx+ V du'/dy)+u' dU/dx
454  Vmath::Vvtvp(nPointsTot,grad_base_u0,1,advVel[0],1,outarray[n],1,outarray[n],1);
455  //Evaluate (U du'/dx+ V du'/dy +u' dU/dx)+v' dU/dy
456  Vmath::Vvtvp(nPointsTot,grad_base_u1,1,advVel[1],1,outarray[n],1,outarray[n],1);
457  //Evaluate (U du'/dx+ V du'/dy +u' dU/dx +v' dU/dy) + W du'/dz
458  Vmath::Vvtvp(nPointsTot,grad2,1,m_baseflow[2],1,outarray[n],1,outarray[n],1);
459  //Evaluate (U du'/dx+ V du'/dy +u' dU/dx +v' dU/dy + W du'/dz)+ w' dU/dz
460  Vmath::Vvtvp(nPointsTot,grad_base_u2,1,advVel[2],1,outarray[n],1,outarray[n],1);
461  }
462  break;
463  //y-equation
464  case 1:
465  if(m_homogen_dealiasing)
466  {
467  //U dv'/dx
468  fields[0]->DealiasedProd(m_baseflow[0],grad0,grad0,m_CoeffState);
469  //V dv'/dy
470  fields[0]->DealiasedProd(m_baseflow[1],grad1,grad1,m_CoeffState);
471  //W dv'/dx
472  fields[0]->DealiasedProd(m_baseflow[2],grad2,grad2,m_CoeffState);
473  // u' dV/dx
474  fields[0]->DealiasedProd(advVel[0],grad_base_v0,grad_base_v0,m_CoeffState);
475  // v' dV/dy
476  fields[0]->DealiasedProd(advVel[1],grad_base_v1,grad_base_v1,m_CoeffState);
477  // w' dV/dz
478  fields[0]->DealiasedProd(advVel[2],grad_base_v2,grad_base_v2,m_CoeffState);
479 
480  Vmath::Vadd(nPointsTot,grad0,1,grad1,1,outarray[n],1);
481  Vmath::Vadd(nPointsTot,grad2,1,outarray[n],1,outarray[n],1);
482  Vmath::Vadd(nPointsTot,grad_base_v0,1,outarray[n],1,outarray[n],1);
483  Vmath::Vadd(nPointsTot,grad_base_v1,1,outarray[n],1,outarray[n],1);
484  Vmath::Vadd(nPointsTot,grad_base_v2,1,outarray[n],1,outarray[n],1);
485  }
486  else
487  {
488  //Evaluate U dv'/dx
489  Vmath::Vmul (nPointsTot,grad0,1,m_baseflow[0],1,outarray[n],1);
490  //Evaluate U dv'/dx+ V dv'/dy
491  Vmath::Vvtvp(nPointsTot,grad1,1,m_baseflow[1],1,outarray[n],1,outarray[n],1);
492  //Evaluate (U dv'/dx+ V dv'/dy)+u' dV/dx
493  Vmath::Vvtvp(nPointsTot,grad_base_v0,1,advVel[0],1,outarray[n],1,outarray[n],1);
494  //Evaluate (U du'/dx+ V du'/dy +u' dV/dx)+v' dV/dy
495  Vmath::Vvtvp(nPointsTot,grad_base_v1,1,advVel[1],1,outarray[n],1,outarray[n],1);
496  //Evaluate (U du'/dx+ V dv'/dy +u' dV/dx +v' dV/dy) + W du'/dz
497  Vmath::Vvtvp(nPointsTot,grad2,1,m_baseflow[2],1,outarray[n],1,outarray[n],1);
498  //Evaluate (U du'/dx+ V dv'/dy +u' dV/dx +v' dV/dy + W dv'/dz)+ w' dV/dz
499  Vmath::Vvtvp(nPointsTot,grad_base_v2,1,advVel[2],1,outarray[n],1,outarray[n],1);
500  }
501  break;
502 
503  //z-equation
504  case 2:
505  if(m_homogen_dealiasing)
506  {
507  //U dw'/dx
508  fields[0]->DealiasedProd(m_baseflow[0],grad0,grad0,m_CoeffState);
509  //V dw'/dy
510  fields[0]->DealiasedProd(m_baseflow[1],grad1,grad1,m_CoeffState);
511  //W dw'/dx
512  fields[0]->DealiasedProd(m_baseflow[2],grad2,grad2,m_CoeffState);
513  // u' dW/dx
514  fields[0]->DealiasedProd(advVel[0],grad_base_w0,grad_base_w0,m_CoeffState);
515  // v' dW/dy
516  fields[0]->DealiasedProd(advVel[1],grad_base_w1,grad_base_w1,m_CoeffState);
517  // w' dW/dz
518  fields[0]->DealiasedProd(advVel[2],grad_base_w2,grad_base_w2,m_CoeffState);
519 
520  Vmath::Vadd(nPointsTot,grad0,1,grad1,1,outarray[n],1);
521  Vmath::Vadd(nPointsTot,grad2,1,outarray[n],1,outarray[n],1);
522  Vmath::Vadd(nPointsTot,grad_base_w0,1,outarray[n],1,outarray[n],1);
523  Vmath::Vadd(nPointsTot,grad_base_w1,1,outarray[n],1,outarray[n],1);
524  Vmath::Vadd(nPointsTot,grad_base_w2,1,outarray[n],1,outarray[n],1);
525  }
526  else
527  {
528  //Evaluate U dw'/dx
529  Vmath::Vmul (nPointsTot,grad0,1,m_baseflow[0],1,outarray[n],1);
530  //Evaluate U dw'/dx+ V dw'/dx
531  Vmath::Vvtvp(nPointsTot,grad1,1,m_baseflow[1],1,outarray[n],1,outarray[n],1);
532  //Evaluate (U dw'/dx+ V dw'/dx)+u' dW/dx
533  Vmath::Vvtvp(nPointsTot,grad_base_w0,1,advVel[0],1,outarray[n],1,outarray[n],1);
534  //Evaluate (U dw'/dx+ V dw'/dx +w' dW/dx)+v' dW/dy
535  Vmath::Vvtvp(nPointsTot,grad_base_w1,1,advVel[1],1,outarray[n],1,outarray[n],1);
536  //Evaluate (U dw'/dx+ V dw'/dx +u' dW/dx +v' dW/dy) + W dw'/dz
537  Vmath::Vvtvp(nPointsTot,grad2,1,m_baseflow[2],1,outarray[n],1,outarray[n],1);
538  //Evaluate (U dw'/dx+ V dw'/dx +u' dW/dx +v' dW/dy + W dw'/dz)+ w' dW/dz
539  Vmath::Vvtvp(nPointsTot,grad_base_w2,1,advVel[2],1,outarray[n],1,outarray[n],1);
540  }
541  break;
542  }
543  break;
544  default:
545  ASSERTL0(false,"dimension unknown");
546  }
547 
548  Vmath::Neg(nqtot,outarray[n],1);
549  }
550 }
551 
552 void LinearisedAdvection::v_SetBaseFlow(
553  const Array<OneD, Array<OneD, NekDouble> > &inarray)
554 {
555  if (m_session->GetSolverInfo("EqType") == "UnsteadyNavierStokes")
556  {
557  // The SFD method is only applied to the velocity variables in
558  // incompressible
559  ASSERTL1(inarray.num_elements() == (m_baseflow.num_elements() - 1),
560  "Number of base flow variables does not match what is "
561  "expected.");
562  }
563  else
564  {
565  ASSERTL1(inarray.num_elements() == (m_baseflow.num_elements()),
566  "Number of base flow variables does not match what is expected.");
567  }
568 
569  int npts = inarray[0].num_elements();
570 
571  for (int i = 0; i < inarray.num_elements(); ++i)
572  {
573  ASSERTL1(npts == m_baseflow[i].num_elements(),
574  "Size of base flow array does not match expected.");
575  Vmath::Vcopy(npts, inarray[i], 1, m_baseflow[i], 1);
576  }
577 }
578 
579 
580 /**
581  * Import field from infile and load into \a m_fields. This routine will
582  * also perform a \a BwdTrans to ensure data is in both the physical and
583  * coefficient storage.
584  * @param infile Filename to read.
585  */
586 void LinearisedAdvection::ImportFldBase(std::string pInfile,
587  Array<OneD, MultiRegions::ExpListSharedPtr>& pFields, int slice)
588 {
589  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
590  std::vector<std::vector<NekDouble> > FieldData;
591 
592  int nqtot = m_baseflow[0].num_elements();
593  int nvar = m_session->GetVariables().size();
594  int s;
595  Array<OneD, NekDouble> tmp_coeff(pFields[0]->GetNcoeffs(), 0.0);
596 
597  //Get Homogeneous
598  LibUtilities::FieldIOSharedPtr fld =
599  MemoryManager<LibUtilities::FieldIO>::AllocateSharedPtr(
600  m_session->GetComm());
601  fld->Import(pInfile, FieldDef, FieldData);
602 
603 
604  if(m_session->DefinesSolverInfo("HOMOGENEOUS"))
605  {
606  std::string HomoStr = m_session->GetSolverInfo("HOMOGENEOUS");
607  }
608 
609  // copy FieldData into m_fields
610  for(int j = 0; j < nvar; ++j)
611  {
612  for(int i = 0; i < FieldDef.size(); ++i)
613  {
614  if((m_session->DefinesSolverInfo("HOMOGENEOUS") &&
615  (m_session->GetSolverInfo("HOMOGENEOUS")=="HOMOGENEOUS1D" ||
616  m_session->GetSolverInfo("HOMOGENEOUS")=="1D" ||
617  m_session->GetSolverInfo("HOMOGENEOUS")=="Homo1D")) &&
618  m_MultipleModes==false)
619  {
620  // w-component must be ignored and set to zero.
621  if (j != nvar - 2)
622  {
623  // p component (it is 4th variable of the 3D and corresponds 3nd variable of 2D)
624  s = (j == nvar - 1) ? 2 : j;
625 
626  //extraction of the 2D
627  pFields[j]->ExtractDataToCoeffs(
628  FieldDef[i],
629  FieldData[i],
630  FieldDef[i]->m_fields[s],
631  tmp_coeff);
632 
633  }
634 
635  //Put zero on higher modes
636  int ncplane = (pFields[0]->GetNcoeffs()) / m_npointsZ;
637 
638  if (m_npointsZ > 2)
639  {
640  Vmath::Zero(ncplane*(m_npointsZ-2),
641  &tmp_coeff[2*ncplane], 1);
642  }
643  }
644  // 2D cases and Homogeneous1D Base Flows
645  else
646  {
647  bool flag = FieldDef[i]->m_fields[j] ==
648  m_session->GetVariable(j);
649 
650  ASSERTL0(flag, (std::string("Order of ") + pInfile
651  + std::string(" data and that defined in "
652  "m_boundaryconditions differs")).c_str());
653 
654  pFields[j]->ExtractDataToCoeffs(FieldDef[i], FieldData[i],
655  FieldDef[i]->m_fields[j],
656  tmp_coeff);
657  }
658  }
659 
660  if(m_SingleMode || m_HalfMode)
661  {
662  //pFields[j]->SetWaveSpace(true);
663 
664  pFields[j]->GetPlane(0)->BwdTrans(tmp_coeff, m_baseflow[j]);
665 
666  if(m_SingleMode)
667  {
668  //copy the bwd into the second plane for single Mode Analysis
669  int ncplane=(pFields[0]->GetNpoints())/m_npointsZ;
670  Vmath::Vcopy(ncplane,&m_baseflow[j][0],1,&m_baseflow[j][ncplane],1);
671  }
672  }
673  else
674  {
675  pFields[j]->BwdTrans(tmp_coeff, m_baseflow[j]);
676 
677  }
678  }
679 
680  if(m_session->DefinesParameter("N_slices"))
681  {
682  int nConvectiveFields = pFields.num_elements()-1;
683 
684  for(int i=0; i<nConvectiveFields;++i)
685  {
686 
687  Vmath::Vcopy(nqtot, &m_baseflow[i][0], 1, &m_interp[i][slice*nqtot], 1);
688  }
689 
690  }
691 }
692 
693 
694 void LinearisedAdvection::UpdateBase( const NekDouble m_slices,
695  const Array<OneD, const NekDouble> &inarray,
696  Array<OneD, NekDouble> &outarray,
697  const NekDouble m_time,
698  const NekDouble m_period)
699 {
700  int npoints=m_baseflow[0].num_elements();
701 
702  NekDouble BetaT=2*M_PI*fmod (m_time, m_period) / m_period;
703  NekDouble phase;
704  Array<OneD, NekDouble> auxiliary(npoints);
705 
706  Vmath::Vcopy(npoints,&inarray[0],1,&outarray[0],1);
707  Vmath::Svtvp(npoints, cos(0.5*m_slices*BetaT),&inarray[npoints],1,&outarray[0],1,&outarray[0],1);
708 
709  for (int i = 2; i < m_slices; i += 2)
710  {
711  phase = (i>>1) * BetaT;
712 
713  Vmath::Svtvp(npoints, cos(phase),&inarray[i*npoints],1,&outarray[0],1,&outarray[0],1);
714  Vmath::Svtvp(npoints, sin(phase), &inarray[(i+1)*npoints], 1, &outarray[0], 1,&outarray[0],1);
715  }
716 
717 }
718 
719 DNekBlkMatSharedPtr LinearisedAdvection::GetFloquetBlockMatrix(FloquetMatType mattype, bool UseContCoeffs) const
720 {
721  DNekMatSharedPtr loc_mat;
722  DNekBlkMatSharedPtr BlkMatrix;
723  int n_exp = 0;
724 
725  n_exp = m_baseflow[0].num_elements(); // will operatore on m_phys
726 
727  Array<OneD,unsigned int> nrows(n_exp);
728  Array<OneD,unsigned int> ncols(n_exp);
729 
730  nrows = Array<OneD, unsigned int>(n_exp,m_slices);
731  ncols = Array<OneD, unsigned int>(n_exp,m_slices);
732 
733  MatrixStorage blkmatStorage = eDIAGONAL;
734  BlkMatrix = MemoryManager<DNekBlkMat>
735  ::AllocateSharedPtr(nrows,ncols,blkmatStorage);
736 
737 
738  const LibUtilities::PointsKey Pkey(m_slices,LibUtilities::eFourierEvenlySpaced);
739  const LibUtilities::BasisKey BK(LibUtilities::eFourier,m_slices,Pkey);
740  StdRegions::StdSegExp StdSeg(BK);
741 
742  StdRegions::StdMatrixKey matkey(StdRegions::eFwdTrans,
743  StdSeg.DetShapeType(),
744  StdSeg);
745 
746  loc_mat = StdSeg.GetStdMatrix(matkey);
747 
748  // set up array of block matrices.
749  for(int i = 0; i < n_exp; ++i)
750  {
751  BlkMatrix->SetBlock(i,i,loc_mat);
752  }
753 
754  return BlkMatrix;
755 }
756 
757 //Discrete Fourier Transform for Floquet analysis
758 void LinearisedAdvection::DFT(const string file,
759  Array<OneD, MultiRegions::ExpListSharedPtr>& pFields,
760  const NekDouble m_slices)
761 {
762  int npoints=m_baseflow[0].num_elements();
763 
764  //Convected fields
765  int ConvectedFields=m_baseflow.num_elements()-1;
766 
767  m_interp= Array<OneD, Array<OneD, NekDouble> > (ConvectedFields);
768  for(int i=0; i<ConvectedFields;++i)
769  {
770  m_interp[i]=Array<OneD,NekDouble>(npoints*m_slices);
771  }
772 
773  //Import the slides into auxiliary vector
774  //The base flow should be stored in the form filename_i.bse
775  for (int i=0; i< m_slices; ++i)
776  {
777  char chkout[16] = "";
778  sprintf(chkout, "%d", i);
779  ImportFldBase(file+"_"+chkout+".bse",pFields,i);
780  }
781 
782 
783  // Discrete Fourier Transform of the fields
784  for(int k=0; k< ConvectedFields;++k)
785  {
786 #ifdef NEKTAR_USING_FFTW
787 
788  //Discrete Fourier Transform using FFTW
789  Array<OneD, NekDouble> fft_in(npoints*m_slices);
790  Array<OneD, NekDouble> fft_out(npoints*m_slices);
791 
792  Array<OneD, NekDouble> m_tmpIN(m_slices);
793  Array<OneD, NekDouble> m_tmpOUT(m_slices);
794 
795  //Shuffle the data
796  for(int j= 0; j < m_slices; ++j)
797  {
798  Vmath::Vcopy(npoints,&m_interp[k][j*npoints],1,&(fft_in[j]),m_slices);
799  }
800 
801  m_FFT = LibUtilities::GetNektarFFTFactory().CreateInstance("NekFFTW", m_slices);
802 
803  //FFT Transform
804  for(int i=0; i<npoints; i++)
805  {
806  m_FFT->FFTFwdTrans(m_tmpIN =fft_in + i*m_slices, m_tmpOUT =fft_out + i*m_slices);
807 
808  }
809 
810  //Reshuffle data
811  for(int s = 0; s < m_slices; ++s)
812  {
813  Vmath::Vcopy(npoints,&fft_out[s],m_slices,&m_interp[k][s*npoints],1);
814 
815  }
816 
817  Vmath::Zero(fft_in.num_elements(),&fft_in[0],1);
818  Vmath::Zero(fft_out.num_elements(),&fft_out[0],1);
819 #else
820  //Discrete Fourier Transform using MVM
821  DNekBlkMatSharedPtr blkmat;
822  blkmat = GetFloquetBlockMatrix(eForwardsPhys);
823 
824  int nrows = blkmat->GetRows();
825  int ncols = blkmat->GetColumns();
826 
827  Array<OneD, NekDouble> sortedinarray(ncols);
828  Array<OneD, NekDouble> sortedoutarray(nrows);
829 
830  //Shuffle the data
831  for(int j= 0; j < m_slices; ++j)
832  {
833  Vmath::Vcopy(npoints,&m_interp[k][j*npoints],1,&(sortedinarray[j]),m_slices);
834  }
835 
836  // Create NekVectors from the given data arrays
837  NekVector<NekDouble> in (ncols,sortedinarray,eWrapper);
838  NekVector<NekDouble> out(nrows,sortedoutarray,eWrapper);
839 
840  // Perform matrix-vector multiply.
841  out = (*blkmat)*in;
842 
843  //Reshuffle data
844  for(int s = 0; s < m_slices; ++s)
845  {
846  Vmath::Vcopy(npoints,&sortedoutarray[s],m_slices,&m_interp[k][s*npoints],1);
847  }
848 
849  for(int r=0; r<sortedinarray.num_elements();++r)
850  {
851  sortedinarray[0]=0;
852  sortedoutarray[0]=0;
853  }
854 
855 #endif
856 
857  //scaling of the Fourier coefficients
858  NekDouble j=-1;
859  for (int i = 2; i < m_slices; i += 2)
860  {
861  Vmath::Smul(2*npoints,j,&m_interp[k][i*npoints],1,&m_interp[k][i*npoints],1);
862  j=-j;
863 
864  }
865 
866  }
867 
868 }
869 
870 } //end of namespace
871 
Generated on Wed Mar 4 2015 16:33:15 for Nektar++ by   doxygen 1.8.1.2