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 = GetAdvectionTermFactory().RegisterCreatorFunction("Linearised", LinearisedAdvection::create);
53 
54  /**
55  * Constructor. Creates ...
56  *
57  * \param
58  * \param
59  */
60 
61  LinearisedAdvection::LinearisedAdvection(
62  const LibUtilities::SessionReaderSharedPtr& pSession,
63  const SpatialDomains::MeshGraphSharedPtr& pGraph):
64  AdvectionTerm(pSession, pGraph)
65  {
66  }
67 
68 
69  void LinearisedAdvection::v_InitObject()
70  {
71  AdvectionTerm::v_InitObject();
72 
73  m_boundaryConditions = MemoryManager<SpatialDomains::BoundaryConditions>
74  ::AllocateSharedPtr(m_session, m_graph);
75 
76  //Setting parameters for homogeneous problems
77  m_HomoDirec = 0;
78  m_useFFT = false;
79  m_HomogeneousType = eNotHomogeneous;
80  m_SingleMode = false;
81  m_HalfMode = false;
82  m_MultipleModes = false;
83 
84  if(m_session->DefinesSolverInfo("HOMOGENEOUS"))
85  {
86  std::string HomoStr = m_session->GetSolverInfo("HOMOGENEOUS");
87  m_spacedim = 3;
88 
89  if((HomoStr == "HOMOGENEOUS1D")||(HomoStr == "Homogeneous1D")||
90  (HomoStr == "1D")||(HomoStr == "Homo1D"))
91  {
92  m_HomogeneousType = eHomogeneous1D;
93  m_LhomZ = m_session->GetParameter("LZ");
94  m_HomoDirec = 1;
95 
96  if(m_session->DefinesSolverInfo("ModeType"))
97  {
98  m_session->MatchSolverInfo("ModeType","SingleMode",m_SingleMode,false);
99  m_session->MatchSolverInfo("ModeType","HalfMode",m_HalfMode,false);
100  m_session->MatchSolverInfo("ModeType","MultipleModes",m_MultipleModes,false);
101  }
102 
103  if(m_session->DefinesSolverInfo("ModeType"))
104  {
105  if(m_SingleMode)
106  {
107  m_npointsZ=2;
108 
109  }
110  else if(m_HalfMode)
111  {
112  m_npointsZ=1;
113 
114  }
115  else if(m_MultipleModes)
116  {
117  m_npointsZ = m_session->GetParameter("HomModesZ");
118  }
119  else
120  {
121  ASSERTL0(false, "SolverInfo ModeType not valid");
122 
123 
124  }
125  }
126  else
127  {
128  m_session->LoadParameter("HomModesZ",m_npointsZ);
129 
130  }
131 
132  }
133 
134  if((HomoStr == "HOMOGENEOUS2D")||(HomoStr == "Homogeneous2D")||
135  (HomoStr == "2D")||(HomoStr == "Homo2D"))
136  {
137  m_HomogeneousType = eHomogeneous2D;
138  m_session->LoadParameter("HomModesY", m_npointsY);
139  m_session->LoadParameter("LY", m_LhomY);
140  m_session->LoadParameter("HomModesZ", m_npointsZ);
141  m_session->LoadParameter("LZ", m_LhomZ);
142  m_HomoDirec = 2;
143  }
144 
145  if((HomoStr == "HOMOGENEOUS3D")||(HomoStr == "Homogeneous3D")||
146  (HomoStr == "3D")||(HomoStr == "Homo3D"))
147  {
148  m_HomogeneousType = eHomogeneous3D;
149  m_session->LoadParameter("HomModesX",m_npointsX);
150  m_session->LoadParameter("LX", m_LhomX );
151  m_session->LoadParameter("HomModesY",m_npointsY);
152  m_session->LoadParameter("LY", m_LhomY );
153  m_session->LoadParameter("HomModesZ",m_npointsZ);
154  m_session->LoadParameter("LZ", m_LhomZ );
155  m_HomoDirec = 3;
156  }
157 
158  if(m_session->DefinesSolverInfo("USEFFT"))
159  {
160  m_useFFT = true;
161  }
162  }
163  else
164  {
165  m_npointsZ = 1; // set to default value so can use to identify 2d or 3D (homogeneous) expansions
166  }
167 
168  if(m_session->DefinesSolverInfo("PROJECTION"))
169  {
170  std::string ProjectStr
171  = m_session->GetSolverInfo("PROJECTION");
172 
173  if((ProjectStr == "Continuous")||(ProjectStr == "Galerkin")||
174  (ProjectStr == "CONTINUOUS")||(ProjectStr == "GALERKIN"))
175  {
176  m_projectionType = MultiRegions::eGalerkin;
177  }
178  else if(ProjectStr == "DisContinuous")
179  {
180  m_projectionType = MultiRegions::eDiscontinuous;
181  }
182  else
183  {
184  ASSERTL0(false,"PROJECTION value not recognised");
185  }
186  }
187  else
188  {
189  cerr << "Projection type not specified in SOLVERINFO,"
190  "defaulting to continuous Galerkin" << endl;
191  m_projectionType = MultiRegions::eGalerkin;
192  }
193 
194  SetUpBaseFields(m_graph);
195  ASSERTL0(m_session->DefinesFunction("BaseFlow"),
196  "Base flow must be defined for linearised forms.");
197  string file = m_session->GetFunctionFilename("BaseFlow", 0);
198 
199 
200  //Periodic base flows
201  if(m_session->DefinesParameter("N_slices"))
202  {
203  m_session->LoadParameter("N_slices",m_slices);
204  if(m_slices>1)
205  {
206  DFT(file,m_slices);
207  }
208  else
209  {
210  ASSERTL0(false,"Number of slices must be a positive number greater than 1");
211  }
212  }
213  //Steady base-flow
214  else
215  {
216  m_slices=1;
217 
218  //BaseFlow from file
219  if (m_session->GetFunctionType("BaseFlow", m_session->GetVariable(0))
220  == LibUtilities::eFunctionTypeFile)
221  {
222  ImportFldBase(file,m_graph,1);
223 
224  }
225  //analytic base flow
226  else
227  {
228  int nq = m_base[0]->GetNpoints();
229  Array<OneD,NekDouble> x0(nq);
230  Array<OneD,NekDouble> x1(nq);
231  Array<OneD,NekDouble> x2(nq);
232 
233  // get the coordinates (assuming all fields have the same
234  // discretisation)
235  m_base[0]->GetCoords(x0,x1,x2);
236  for(unsigned int i = 0 ; i < m_base.num_elements(); i++)
237  {
238  LibUtilities::EquationSharedPtr ifunc
239  = m_session->GetFunction("BaseFlow", i);
240 
241  ifunc->Evaluate(x0,x1,x2,m_base[i]->UpdatePhys());
242 
243  m_base[i]->SetPhysState(true);
244  m_base[i]->FwdTrans_IterPerExp(m_base[i]->GetPhys(),
245  m_base[i]->UpdateCoeffs());
246  }
247  }
248  }
249  }
250 
251  LinearisedAdvection::~LinearisedAdvection()
252  {
253  }
254 
255  void LinearisedAdvection::SetUpBaseFields(SpatialDomains::MeshGraphSharedPtr &mesh)
256  {
257  int nvariables = m_session->GetVariables().size();
258  int i;
259  m_base = Array<OneD, MultiRegions::ExpListSharedPtr>(nvariables);
260  if (m_projectionType == MultiRegions::eGalerkin)
261  {
262  switch (m_expdim)
263  {
264  case 1:
265  {
266  if(m_HomogeneousType == eHomogeneous2D)
267  {
268  const LibUtilities::PointsKey PkeyY(m_npointsY,LibUtilities::eFourierEvenlySpaced);
269  const LibUtilities::BasisKey BkeyY(LibUtilities::eFourier,m_npointsY,PkeyY);
270  const LibUtilities::PointsKey PkeyZ(m_npointsZ,LibUtilities::eFourierEvenlySpaced);
271  const LibUtilities::BasisKey BkeyZ(LibUtilities::eFourier,m_npointsZ,PkeyZ);
272 
273  for(i = 0 ; i < m_base.num_elements(); i++)
274  {
275  m_base[i] = MemoryManager<MultiRegions::ContField3DHomogeneous2D>
276  ::AllocateSharedPtr(m_session,BkeyY,BkeyZ,m_LhomY,m_LhomZ,m_useFFT,m_homogen_dealiasing,m_graph,m_session->GetVariable(i));
277  }
278  }
279 
280  else {
281 
282  for(i = 0 ; i < m_base.num_elements(); i++)
283  {
284  m_base[i] = MemoryManager<MultiRegions::ContField1D>
285  ::AllocateSharedPtr(m_session,mesh,
286  m_session->GetVariable(i));
287  }
288  }
289 
290  }
291  break;
292  case 2:
293  {
294  if(m_HomogeneousType == eHomogeneous1D)
295  {
296  if(m_SingleMode)
297  {
298  const LibUtilities::PointsKey PkeyZ(m_npointsZ,LibUtilities::eFourierSingleModeSpaced);
299  const LibUtilities::BasisKey BkeyZ(LibUtilities::eFourier,m_npointsZ,PkeyZ);
300 
301  for(i = 0 ; i < m_base.num_elements(); i++)
302  {
303  m_base[i] = MemoryManager<MultiRegions::ContField3DHomogeneous1D>
304  ::AllocateSharedPtr(m_session,BkeyZ,m_LhomZ,m_useFFT,m_homogen_dealiasing,m_graph,m_session->GetVariable(i));
305  m_base[i]->SetWaveSpace(true);
306 
307 
308  }
309  }
310  else if(m_HalfMode)
311  {
312  //1 plane field (half mode expansion)
313  const LibUtilities::PointsKey PkeyZ(m_npointsZ,LibUtilities::eFourierSingleModeSpaced);
314  const LibUtilities::BasisKey BkeyZ(LibUtilities::eFourierHalfModeRe,m_npointsZ,PkeyZ);
315 
316  for(i = 0 ; i < m_base.num_elements(); i++)
317  {
318  m_base[i] = MemoryManager<MultiRegions::ContField3DHomogeneous1D>
319  ::AllocateSharedPtr(m_session,BkeyZ,m_LhomZ,m_useFFT,m_homogen_dealiasing,m_graph,m_session->GetVariable(i));
320  m_base[i]->SetWaveSpace(true);
321 
322  }
323 
324  }
325  else
326  {
327  const LibUtilities::PointsKey PkeyZ(m_npointsZ,LibUtilities::eFourierEvenlySpaced);
328  const LibUtilities::BasisKey BkeyZ(LibUtilities::eFourier,m_npointsZ,PkeyZ);
329 
330 
331  for(i = 0 ; i < m_base.num_elements(); i++)
332  {
333  m_base[i] = MemoryManager<MultiRegions::ContField3DHomogeneous1D>
334  ::AllocateSharedPtr(m_session,BkeyZ,m_LhomZ,m_useFFT,m_homogen_dealiasing,m_graph,m_session->GetVariable(i));
335  m_base[i]->SetWaveSpace(false);
336  }
337 
338  }
339  }
340  else
341  {
342  i = 0;
343  MultiRegions::ContField2DSharedPtr firstbase =
344  MemoryManager<MultiRegions::ContField2D>
345  ::AllocateSharedPtr(m_session,mesh,
346  m_session->GetVariable(i));
347  m_base[0]=firstbase;
348 
349  for(i = 1 ; i < m_base.num_elements(); i++)
350  {
351  m_base[i] = MemoryManager<MultiRegions::ContField2D>
352  ::AllocateSharedPtr(*firstbase,mesh,
353  m_session->GetVariable(i));
354  }
355  }
356  }
357  break;
358  case 3:
359  {
360  if(m_HomogeneousType == eHomogeneous3D)
361  {
362  ASSERTL0(false,"3D fully periodic problems not implemented yet");
363  }
364  else
365  {
366  MultiRegions::ContField3DSharedPtr firstbase =
367  MemoryManager<MultiRegions::ContField3D>
368  ::AllocateSharedPtr(m_session,mesh,
369  m_session->GetVariable(0));
370  m_base[0] = firstbase;
371 
372  for(i = 1 ; i < m_base.num_elements(); i++)
373  {
374  m_base[i] = MemoryManager<MultiRegions::ContField3D>
375  ::AllocateSharedPtr(*firstbase,mesh,
376  m_session->GetVariable(i));
377  }
378  }
379  }
380  break;
381  default:
382  ASSERTL0(false,"Expansion dimension not recognised");
383  break;
384  }
385  }
386  else
387  {
388  switch(m_expdim)
389  {
390  case 1:
391  {
392  if(m_HomogeneousType == eHomogeneous2D)
393  {
394  const LibUtilities::PointsKey PkeyY(m_npointsY,LibUtilities::eFourierEvenlySpaced);
395  const LibUtilities::BasisKey BkeyY(LibUtilities::eFourier,m_npointsY,PkeyY);
396  const LibUtilities::PointsKey PkeyZ(m_npointsZ,LibUtilities::eFourierEvenlySpaced);
397  const LibUtilities::BasisKey BkeyZ(LibUtilities::eFourier,m_npointsZ,PkeyZ);
398 
399  for(i = 0 ; i < m_base.num_elements(); i++)
400  {
401  m_base[i] = MemoryManager<MultiRegions::DisContField3DHomogeneous2D>
402  ::AllocateSharedPtr(m_session,BkeyY,BkeyZ,m_LhomY,m_LhomZ,m_useFFT,m_homogen_dealiasing,m_graph,m_session->GetVariable(i));
403  }
404  }
405  else
406  {
407  for(i = 0 ; i < m_base.num_elements(); i++)
408  {
409  m_base[i] = MemoryManager<MultiRegions
410  ::DisContField1D>::AllocateSharedPtr(m_session,mesh,
411  m_session->GetVariable(i));
412  }
413  }
414  break;
415  }
416  case 2:
417  {
418  if(m_HomogeneousType == eHomogeneous1D)
419  {
420 
421  const LibUtilities::PointsKey PkeyZ(m_npointsZ,LibUtilities::eFourierEvenlySpaced);
422  const LibUtilities::BasisKey BkeyZ(LibUtilities::eFourier,m_npointsZ,PkeyZ);
423 
424  for(i = 0 ; i < m_base.num_elements(); i++)
425  {
426  m_base[i] = MemoryManager<MultiRegions::DisContField3DHomogeneous1D>
427  ::AllocateSharedPtr(m_session,BkeyZ,m_LhomZ,m_useFFT,m_homogen_dealiasing,m_graph,m_session->GetVariable(i));
428  }
429 
430 
431 
432  }
433  else
434  {
435  for(i = 0 ; i < m_base.num_elements(); i++)
436  {
437  m_base[i] = MemoryManager<MultiRegions
438  ::DisContField2D>::AllocateSharedPtr(m_session, mesh,
439  m_session->GetVariable(i));
440  }
441  }
442  break;
443 
444  }
445  case 3:
446  ASSERTL0(false,"3 D not set up");
447  default:
448  ASSERTL0(false,"Expansion dimension not recognised");
449  break;
450  }
451  }
452 
453  }
454 
455  /**
456  * Import field from infile and load into \a m_fields. This routine will
457  * also perform a \a BwdTrans to ensure data is in both the physical and
458  * coefficient storage.
459  * @param infile Filename to read.
460  */
461  void LinearisedAdvection::ImportFldBase(std::string pInfile,
462  SpatialDomains::MeshGraphSharedPtr pGraph, int cnt)
463  {
464  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
465  std::vector<std::vector<NekDouble> > FieldData;
466  int nqtot = m_base[0]->GetTotPoints();
467 
468  //Get Homogeneous
469  LibUtilities::FieldIOSharedPtr fld =
470  MemoryManager<LibUtilities::FieldIO>::AllocateSharedPtr(
471  m_session->GetComm());
472  fld->Import(pInfile, FieldDef, FieldData);
473 
474  int nvar = m_session->GetVariables().size();
475  int s;
476 
477  if(m_session->DefinesSolverInfo("HOMOGENEOUS"))
478  {
479  std::string HomoStr = m_session->GetSolverInfo("HOMOGENEOUS");
480  }
481 
482  // copy FieldData into m_fields
483  for(int j = 0; j < nvar; ++j)
484  {
485  for(int i = 0; i < FieldDef.size(); ++i)
486  {
487  if((m_session->DefinesSolverInfo("HOMOGENEOUS") &&
488  (m_session->GetSolverInfo("HOMOGENEOUS")=="HOMOGENEOUS1D" ||
489  m_session->GetSolverInfo("HOMOGENEOUS")=="1D" ||
490  m_session->GetSolverInfo("HOMOGENEOUS")=="Homo1D")) &&
491  m_MultipleModes==false)
492  {
493  // w-component must be ignored and set to zero.
494  if (j != nvar - 2)
495  {
496  // p component (it is 4th variable of the 3D and corresponds 3nd variable of 2D)
497  s = (j == nvar - 1) ? 2 : j;
498 
499  //extraction of the 2D
500  m_base[j]->ExtractDataToCoeffs(
501  FieldDef[i],
502  FieldData[i],
503  FieldDef[i]->m_fields[s],
504  m_base[j]->UpdateCoeffs());
505 
506  }
507 
508  //Put zero on higher modes
509  int ncplane = (m_base[0]->GetNcoeffs()) / m_npointsZ;
510 
511  if (m_npointsZ > 2)
512  {
513  Vmath::Zero(ncplane*(m_npointsZ-2),
514  &m_base[j]->UpdateCoeffs()[2*ncplane], 1);
515  }
516  }
517  // 2D cases and Homogeneous1D Base Flows
518  else
519  {
520  bool flag = FieldDef[i]->m_fields[j] ==
521  m_session->GetVariable(j);
522 
523  ASSERTL1(flag, (std::string("Order of ") + pInfile
524  + std::string(" data and that defined in "
525  "m_boundaryconditions differs")).c_str());
526 
527  m_base[j]->ExtractDataToCoeffs(FieldDef[i], FieldData[i],
528  FieldDef[i]->m_fields[j],
529  m_base[j]->UpdateCoeffs());
530  }
531  }
532 
533  if(m_SingleMode || m_HalfMode)
534  {
535  m_base[j]->SetWaveSpace(true);
536 
537  m_base[j]->BwdTrans(m_base[j]->GetCoeffs(),
538  m_base[j]->UpdatePhys());
539 
540  if(m_SingleMode)
541  {
542  //copy the bwd into the second plane for single Mode Analysis
543  int ncplane=(m_base[0]->GetNpoints())/m_npointsZ;
544  Vmath::Vcopy(ncplane,&m_base[j]->GetPhys()[0],1,&m_base[j]->UpdatePhys()[ncplane],1);
545  }
546  }
547  else
548  {
549  m_base[j]->BwdTrans(m_base[j]->GetCoeffs(),
550  m_base[j]->UpdatePhys());
551 
552  }
553 
554  }
555 
556  //std::string outname ="BaseFlow.bse";
557  //WriteFldBase(outname);
558 
559  if(m_session->DefinesParameter("N_slices"))
560  {
561  m_nConvectiveFields = m_base.num_elements()-1;
562 
563  for(int i=0; i<m_nConvectiveFields;++i)
564  {
565 
566  Vmath::Vcopy(nqtot, &m_base[i]->GetPhys()[0], 1, &m_interp[i][cnt*nqtot], 1);
567  }
568 
569  }
570  }
571 
572 
573  //Evaluation of the advective terms
574  void LinearisedAdvection::v_ComputeAdvectionTerm(
575  Array<OneD, MultiRegions::ExpListSharedPtr > &pFields,
576  const Array<OneD, Array<OneD, NekDouble> > &pVelocity,
577  const Array<OneD, const NekDouble> &pU,
578  Array<OneD, NekDouble> &pOutarray,
579  int pVelocityComponent,
580  NekDouble m_time,
581  Array<OneD, NekDouble> &pWk)
582  {
583  int ndim = m_nConvectiveFields;
584  int nPointsTot = pFields[0]->GetNpoints();
585 
586  Array<OneD, NekDouble> grad0,grad1,grad2;
587 
588  //Evaluation of the gradiend of each component of the base flow
589  //\nabla U
590  Array<OneD, NekDouble> grad_base_u0,grad_base_u1,grad_base_u2;
591 
592  // \nabla V
593  Array<OneD, NekDouble> grad_base_v0,grad_base_v1,grad_base_v2;
594 
595  // \nabla W
596  Array<OneD, NekDouble> grad_base_w0,grad_base_w1,grad_base_w2;
597 
598  grad0 = Array<OneD, NekDouble> (nPointsTot);
599  grad_base_u0 = Array<OneD, NekDouble> (nPointsTot);
600  grad_base_v0 = Array<OneD, NekDouble> (nPointsTot);
601  grad_base_w0 = Array<OneD, NekDouble> (nPointsTot);
602 
603  //Evaluation of the base flow for periodic cases
604  //(it requires fld files)
605  if(m_slices>1)
606  {
607  if (m_session->GetFunctionType("BaseFlow", 0)
608  == LibUtilities::eFunctionTypeFile)
609  {
610  for(int i=0; i<m_nConvectiveFields;++i)
611  {
612  UpdateBase(m_slices,m_interp[i],m_base[i]->UpdatePhys(),m_time,m_period);
613  }
614  }
615  else
616  {
617  ASSERTL0(false, "Periodic Base flow requires filename_ files");
618  }
619  }
620 
621 
622  //Evaluate the linearised advection term
623  switch(ndim)
624  {
625  // 1D
626  case 1:
627  pFields[0]->PhysDeriv(pU,grad0);
628  pFields[0]->PhysDeriv(m_base[0]->GetPhys(),grad_base_u0);
629  //Evaluate U du'/dx
630  Vmath::Vmul(nPointsTot,grad0,1,m_base[0]->GetPhys(),1,pOutarray,1);
631  //Evaluate U du'/dx+ u' dU/dx
632  Vmath::Vvtvp(nPointsTot,grad_base_u0,1,pVelocity[0],1,pOutarray,1,pOutarray,1);
633  break;
634 
635  //2D
636  case 2:
637  grad1 = Array<OneD, NekDouble> (nPointsTot);
638  grad_base_u1 = Array<OneD, NekDouble> (nPointsTot);
639  grad_base_v1 = Array<OneD, NekDouble> (nPointsTot);
640 
641  pFields[0]->PhysDeriv(pU,grad0,grad1);
642 
643  //Derivates of the base flow
644  pFields[0]-> PhysDeriv(m_base[0]->GetPhys(), grad_base_u0, grad_base_u1);
645  pFields[0]-> PhysDeriv(m_base[1]->GetPhys(), grad_base_v0, grad_base_v1);
646 
647  //Since the components of the velocity are passed one by
648  //one, it is necessary to distinguish which term is
649  //consider
650  switch (pVelocityComponent)
651  {
652  //x-equation
653  case 0:
654  // Evaluate U du'/dx
655  Vmath::Vmul (nPointsTot,grad0,1,m_base[0]->GetPhys(),1,pOutarray,1);
656  //Evaluate U du'/dx+ V du'/dy
657  Vmath::Vvtvp(nPointsTot,grad1,1,m_base[1]->GetPhys(),1,pOutarray,1,pOutarray,1);
658  //Evaluate (U du'/dx+ V du'/dy)+u' dU/dx
659  Vmath::Vvtvp(nPointsTot,grad_base_u0,1,pVelocity[0],1,pOutarray,1,pOutarray,1);
660  //Evaluate (U du'/dx+ V du'/dy +u' dU/dx)+v' dU/dy
661  Vmath::Vvtvp(nPointsTot,grad_base_u1,1,pVelocity[1],1,pOutarray,1,pOutarray,1);
662  break;
663 
664  //y-equation
665  case 1:
666  // Evaluate U dv'/dx
667  Vmath::Vmul (nPointsTot,grad0,1,m_base[0]->GetPhys(),1,pOutarray,1);
668  //Evaluate U dv'/dx+ V dv'/dy
669  Vmath::Vvtvp(nPointsTot,grad1,1,m_base[1]->GetPhys(),1,pOutarray,1,pOutarray,1);
670  //Evaluate (U dv'/dx+ V dv'/dy)+u' dV/dx
671  Vmath::Vvtvp(nPointsTot,grad_base_v0,1,pVelocity[0],1,pOutarray,1,pOutarray,1);
672  //Evaluate (U dv'/dx+ V dv'/dy +u' dv/dx)+v' dV/dy
673  Vmath::Vvtvp(nPointsTot,grad_base_v1,1,pVelocity[1],1,pOutarray,1,pOutarray,1);
674  break;
675  }
676  break;
677 
678  //3D
679  case 3:
680  grad1 = Array<OneD, NekDouble> (nPointsTot);
681  grad2 = Array<OneD, NekDouble> (nPointsTot);
682  grad_base_u1 = Array<OneD, NekDouble> (nPointsTot);
683  grad_base_v1 = Array<OneD, NekDouble> (nPointsTot);
684  grad_base_w1 = Array<OneD, NekDouble> (nPointsTot);
685 
686  grad_base_u2 = Array<OneD, NekDouble> (nPointsTot);
687  grad_base_v2 = Array<OneD, NekDouble> (nPointsTot);
688  grad_base_w2 = Array<OneD, NekDouble> (nPointsTot);
689 
690  m_base[0]->PhysDeriv(m_base[0]->GetPhys(), grad_base_u0, grad_base_u1,grad_base_u2);
691  m_base[0]->PhysDeriv(m_base[1]->GetPhys(), grad_base_v0, grad_base_v1,grad_base_v2);
692  m_base[0]->PhysDeriv(m_base[2]->GetPhys(), grad_base_w0, grad_base_w1, grad_base_w2);
693 
694  //HalfMode has W(x,y,t)=0
695  if(m_HalfMode)
696  {
697  for(int i=0; i<grad_base_u2.num_elements();++i)
698  {
699  grad_base_u2[i]=0;
700  grad_base_v2[i]=0;
701  grad_base_w2[i]=0;
702  }
703  }
704 
705  pFields[0]->PhysDeriv(pU, grad0, grad1, grad2);
706 
707  switch (pVelocityComponent)
708  {
709  //x-equation
710  case 0:
711  if(m_homogen_dealiasing)
712  {
713  //U du'/dx
714  pFields[0]->DealiasedProd(m_base[0]->GetPhys(),grad0,grad0,m_CoeffState);
715 
716  //V du'/dy
717  pFields[0]->DealiasedProd(m_base[1]->GetPhys(),grad1,grad1,m_CoeffState);
718 
719  //W du'/dx
720  pFields[0]->DealiasedProd(m_base[2]->GetPhys(),grad2,grad2,m_CoeffState);
721 
722  // u' dU/dx
723  pFields[0]->DealiasedProd(pVelocity[0],grad_base_u0,grad_base_u0,m_CoeffState);
724  // v' dU/dy
725  pFields[0]->DealiasedProd(pVelocity[1],grad_base_u1,grad_base_u1,m_CoeffState);
726  // w' dU/dz
727  pFields[0]->DealiasedProd(pVelocity[2],grad_base_u2,grad_base_u2,m_CoeffState);
728 
729  Vmath::Vadd(nPointsTot,grad0,1,grad1,1,pOutarray,1);
730  Vmath::Vadd(nPointsTot,grad2,1,pOutarray,1,pOutarray,1);
731  Vmath::Vadd(nPointsTot,grad_base_u0,1,pOutarray,1,pOutarray,1);
732  Vmath::Vadd(nPointsTot,grad_base_u1,1,pOutarray,1,pOutarray,1);
733  Vmath::Vadd(nPointsTot,grad_base_u2,1,pOutarray,1,pOutarray,1);
734  }
735  else
736  {
737  //Evaluate U du'/dx
738  Vmath::Vmul (nPointsTot,grad0,1,m_base[0]->GetPhys(),1,pOutarray,1);
739  //Evaluate U du'/dx+ V du'/dy
740  Vmath::Vvtvp(nPointsTot,grad1,1,m_base[1]->GetPhys(),1,pOutarray,1,pOutarray,1);
741  //Evaluate (U du'/dx+ V du'/dy)+u' dU/dx
742  Vmath::Vvtvp(nPointsTot,grad_base_u0,1,pVelocity[0],1,pOutarray,1,pOutarray,1);
743  //Evaluate (U du'/dx+ V du'/dy +u' dU/dx)+v' dU/dy
744  Vmath::Vvtvp(nPointsTot,grad_base_u1,1,pVelocity[1],1,pOutarray,1,pOutarray,1);
745  //Evaluate (U du'/dx+ V du'/dy +u' dU/dx +v' dU/dy) + W du'/dz
746  Vmath::Vvtvp(nPointsTot,grad2,1,m_base[2]->GetPhys(),1,pOutarray,1,pOutarray,1);
747  //Evaluate (U du'/dx+ V du'/dy +u' dU/dx +v' dU/dy + W du'/dz)+ w' dU/dz
748  Vmath::Vvtvp(nPointsTot,grad_base_u2,1,pVelocity[2],1,pOutarray,1,pOutarray,1);
749  }
750  break;
751  //y-equation
752  case 1:
753  if(m_homogen_dealiasing)
754  {
755  //U dv'/dx
756  pFields[0]->DealiasedProd(m_base[0]->GetPhys(),grad0,grad0,m_CoeffState);
757  //V dv'/dy
758  pFields[0]->DealiasedProd(m_base[1]->GetPhys(),grad1,grad1,m_CoeffState);
759  //W dv'/dx
760  pFields[0]->DealiasedProd(m_base[2]->GetPhys(),grad2,grad2,m_CoeffState);
761  // u' dV/dx
762  pFields[0]->DealiasedProd(pVelocity[0],grad_base_v0,grad_base_v0,m_CoeffState);
763  // v' dV/dy
764  pFields[0]->DealiasedProd(pVelocity[1],grad_base_v1,grad_base_v1,m_CoeffState);
765  // w' dV/dz
766  pFields[0]->DealiasedProd(pVelocity[2],grad_base_v2,grad_base_v2,m_CoeffState);
767 
768  Vmath::Vadd(nPointsTot,grad0,1,grad1,1,pOutarray,1);
769  Vmath::Vadd(nPointsTot,grad2,1,pOutarray,1,pOutarray,1);
770  Vmath::Vadd(nPointsTot,grad_base_v0,1,pOutarray,1,pOutarray,1);
771  Vmath::Vadd(nPointsTot,grad_base_v1,1,pOutarray,1,pOutarray,1);
772  Vmath::Vadd(nPointsTot,grad_base_v2,1,pOutarray,1,pOutarray,1);
773  }
774  else
775  {
776  //Evaluate U dv'/dx
777  Vmath::Vmul (nPointsTot,grad0,1,m_base[0]->GetPhys(),1,pOutarray,1);
778  //Evaluate U dv'/dx+ V dv'/dy
779  Vmath::Vvtvp(nPointsTot,grad1,1,m_base[1]->GetPhys(),1,pOutarray,1,pOutarray,1);
780  //Evaluate (U dv'/dx+ V dv'/dy)+u' dV/dx
781  Vmath::Vvtvp(nPointsTot,grad_base_v0,1,pVelocity[0],1,pOutarray,1,pOutarray,1);
782  //Evaluate (U du'/dx+ V du'/dy +u' dV/dx)+v' dV/dy
783  Vmath::Vvtvp(nPointsTot,grad_base_v1,1,pVelocity[1],1,pOutarray,1,pOutarray,1);
784  //Evaluate (U du'/dx+ V dv'/dy +u' dV/dx +v' dV/dy) + W du'/dz
785  Vmath::Vvtvp(nPointsTot,grad2,1,m_base[2]->GetPhys(),1,pOutarray,1,pOutarray,1);
786  //Evaluate (U du'/dx+ V dv'/dy +u' dV/dx +v' dV/dy + W dv'/dz)+ w' dV/dz
787  Vmath::Vvtvp(nPointsTot,grad_base_v2,1,pVelocity[2],1,pOutarray,1,pOutarray,1);
788  }
789  break;
790 
791  //z-equation
792  case 2:
793  if(m_homogen_dealiasing)
794  {
795  //U dw'/dx
796  pFields[0]->DealiasedProd(m_base[0]->GetPhys(),grad0,grad0,m_CoeffState);
797  //V dw'/dy
798  pFields[0]->DealiasedProd(m_base[1]->GetPhys(),grad1,grad1,m_CoeffState);
799  //W dw'/dx
800  pFields[0]->DealiasedProd(m_base[2]->GetPhys(),grad2,grad2,m_CoeffState);
801  // u' dW/dx
802  pFields[0]->DealiasedProd(pVelocity[0],grad_base_w0,grad_base_w0,m_CoeffState);
803  // v' dW/dy
804  pFields[0]->DealiasedProd(pVelocity[1],grad_base_w1,grad_base_w1,m_CoeffState);
805  // w' dW/dz
806  pFields[0]->DealiasedProd(pVelocity[2],grad_base_w2,grad_base_w2,m_CoeffState);
807 
808  Vmath::Vadd(nPointsTot,grad0,1,grad1,1,pOutarray,1);
809  Vmath::Vadd(nPointsTot,grad2,1,pOutarray,1,pOutarray,1);
810  Vmath::Vadd(nPointsTot,grad_base_w0,1,pOutarray,1,pOutarray,1);
811  Vmath::Vadd(nPointsTot,grad_base_w1,1,pOutarray,1,pOutarray,1);
812  Vmath::Vadd(nPointsTot,grad_base_w2,1,pOutarray,1,pOutarray,1);
813  }
814  else
815  {
816  //Evaluate U dw'/dx
817  Vmath::Vmul (nPointsTot,grad0,1,m_base[0]->GetPhys(),1,pOutarray,1);
818  //Evaluate U dw'/dx+ V dw'/dx
819  Vmath::Vvtvp(nPointsTot,grad1,1,m_base[1]->GetPhys(),1,pOutarray,1,pOutarray,1);
820  //Evaluate (U dw'/dx+ V dw'/dx)+u' dW/dx
821  Vmath::Vvtvp(nPointsTot,grad_base_w0,1,pVelocity[0],1,pOutarray,1,pOutarray,1);
822  //Evaluate (U dw'/dx+ V dw'/dx +w' dW/dx)+v' dW/dy
823  Vmath::Vvtvp(nPointsTot,grad_base_w1,1,pVelocity[1],1,pOutarray,1,pOutarray,1);
824  //Evaluate (U dw'/dx+ V dw'/dx +u' dW/dx +v' dW/dy) + W dw'/dz
825  Vmath::Vvtvp(nPointsTot,grad2,1,m_base[2]->GetPhys(),1,pOutarray,1,pOutarray,1);
826  //Evaluate (U dw'/dx+ V dw'/dx +u' dW/dx +v' dW/dy + W dw'/dz)+ w' dW/dz
827  Vmath::Vvtvp(nPointsTot,grad_base_w2,1,pVelocity[2],1,pOutarray,1,pOutarray,1);
828  }
829  break;
830  }
831  break;
832  default:
833  ASSERTL0(false,"dimension unknown");
834  }
835  }
836 
837  void LinearisedAdvection::UpdateBase( const NekDouble m_slices,
838  Array<OneD, const NekDouble> &inarray,
839  Array<OneD, NekDouble> &outarray,
840  const NekDouble m_time,
841  const NekDouble m_period)
842  {
843 
844  int npoints=m_base[0]->GetTotPoints();
845 
846  NekDouble BetaT=2*M_PI*fmod (m_time, m_period) / m_period;
847  NekDouble phase;
848  Array<OneD, NekDouble> auxiliary(npoints);
849 
850  Vmath::Vcopy(npoints,&inarray[0],1,&outarray[0],1);
851  Vmath::Svtvp(npoints, cos(0.5*m_slices*BetaT),&inarray[npoints],1,&outarray[0],1,&outarray[0],1);
852 
853  for (int i = 2; i < m_slices; i += 2)
854  {
855  phase = (i>>1) * BetaT;
856 
857  Vmath::Svtvp(npoints, cos(phase),&inarray[i*npoints],1,&outarray[0],1,&outarray[0],1);
858  Vmath::Svtvp(npoints, sin(phase), &inarray[(i+1)*npoints], 1, &outarray[0], 1,&outarray[0],1);
859  }
860 
861  }
862 
863  void LinearisedAdvection::WriteFldBase(std::string &outname)
864  {
865 
866  Array<OneD, Array<OneD, NekDouble> > fieldcoeffs(m_base.num_elements());
867  Array<OneD, std::string> variables(m_base.num_elements());
868 
869  for(int i = 0; i < m_base.num_elements(); ++i)
870  {
871  fieldcoeffs[i] = m_base[i]->UpdateCoeffs();
872  variables[i] = m_boundaryConditions->GetVariable(i);
873  }
874  WriteFldBase(outname, m_base[0], fieldcoeffs, variables);
875 
876  }
877 
878 
879  /**
880  * Writes the field data to a file with the given filename.
881  * @param outname Filename to write to.
882  * @param field ExpList on which data is based
883  * @param fieldcoeffs An array of array of expansion coefficients
884  * @param variables An array of variable names
885  */
886  void LinearisedAdvection::WriteFldBase(std::string &outname, MultiRegions::ExpListSharedPtr &field, Array<OneD, Array<OneD, NekDouble> > &fieldcoeffs, Array<OneD, std::string> &variables)
887  {
888 
889  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef= field->GetFieldDefinitions();
890  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
891 
892  // copy Data into FieldData and set variable
893  for(int j = 0; j < fieldcoeffs.num_elements(); ++j)
894  {
895  for(int i = 0; i < FieldDef.size(); ++i)
896  {
897  // Could do a search here to find correct variable
898  FieldDef[i]->m_fields.push_back(variables[j]);
899  //cout<<"v="<<variables[j]<<endl;
900  field->AppendFieldData(FieldDef[i], FieldData[i], fieldcoeffs[j]);
901  }
902  }
903  LibUtilities::Write(outname,FieldDef,FieldData);
904  }
905 
906 
907  DNekBlkMatSharedPtr LinearisedAdvection::GetFloquetBlockMatrix(FloquetMatType mattype, bool UseContCoeffs) const
908  {
909  DNekMatSharedPtr loc_mat;
910  DNekBlkMatSharedPtr BlkMatrix;
911  int n_exp = 0;
912 
913  n_exp = m_base[0]->GetTotPoints(); // will operatore on m_phys
914 
915  Array<OneD,unsigned int> nrows(n_exp);
916  Array<OneD,unsigned int> ncols(n_exp);
917 
918  nrows = Array<OneD, unsigned int>(n_exp,m_slices);
919  ncols = Array<OneD, unsigned int>(n_exp,m_slices);
920 
921  MatrixStorage blkmatStorage = eDIAGONAL;
922  BlkMatrix = MemoryManager<DNekBlkMat>
923  ::AllocateSharedPtr(nrows,ncols,blkmatStorage);
924 
925 
926  const LibUtilities::PointsKey Pkey(m_slices,LibUtilities::eFourierEvenlySpaced);
927  const LibUtilities::BasisKey BK(LibUtilities::eFourier,m_slices,Pkey);
928  StdRegions::StdSegExp StdSeg(BK);
929 
930  StdRegions::StdMatrixKey matkey(StdRegions::eFwdTrans,
931  StdSeg.DetShapeType(),
932  StdSeg);
933 
934  loc_mat = StdSeg.GetStdMatrix(matkey);
935 
936  // set up array of block matrices.
937  for(int i = 0; i < n_exp; ++i)
938  {
939  BlkMatrix->SetBlock(i,i,loc_mat);
940  }
941 
942  return BlkMatrix;
943  }
944 
945  //Discrete Fourier Transform for Floquet analysis
946  void LinearisedAdvection::DFT(const string file, const NekDouble m_slices)
947  {
948  int npoints=m_base[0]->GetTotPoints();
949 
950  //Convected fields
951  int ConvectedFields=m_base.num_elements()-1;
952 
953  m_interp= Array<OneD, Array<OneD, NekDouble> > (ConvectedFields);
954  for(int i=0; i<ConvectedFields;++i)
955  {
956  m_interp[i]=Array<OneD,NekDouble>(npoints*m_slices);
957  }
958 
959  //Import the slides into auxiliary vector
960  //The base flow should be stored in the form filename_i.bse
961  for (int i=0; i< m_slices; ++i)
962  {
963  char chkout[16] = "";
964  sprintf(chkout, "%d", i);
965  ImportFldBase(file+"_"+chkout+".bse",m_graph,i);
966  }
967 
968 
969  // Discrete Fourier Transform of the fields
970  for(int k=0; k< ConvectedFields;++k)
971  {
972 #ifdef NEKTAR_USING_FFTW
973 
974  //Discrete Fourier Transform using FFTW
975  Array<OneD, NekDouble> fft_in(npoints*m_slices);
976  Array<OneD, NekDouble> fft_out(npoints*m_slices);
977 
978  Array<OneD, NekDouble> m_tmpIN(m_slices);
979  Array<OneD, NekDouble> m_tmpOUT(m_slices);
980 
981  //Shuffle the data
982  for(int j= 0; j < m_slices; ++j)
983  {
984  Vmath::Vcopy(npoints,&m_interp[k][j*npoints],1,&(fft_in[j]),m_slices);
985  }
986 
987  m_FFT = LibUtilities::GetNektarFFTFactory().CreateInstance("NekFFTW", m_slices);
988 
989  //FFT Transform
990  for(int i=0; i<npoints; i++)
991  {
992  m_FFT->FFTFwdTrans(m_tmpIN =fft_in + i*m_slices, m_tmpOUT =fft_out + i*m_slices);
993 
994  }
995 
996  //Reshuffle data
997  for(int s = 0; s < m_slices; ++s)
998  {
999  Vmath::Vcopy(npoints,&fft_out[s],m_slices,&m_interp[k][s*npoints],1);
1000 
1001  }
1002 
1003  Vmath::Zero(fft_in.num_elements(),&fft_in[0],1);
1004  Vmath::Zero(fft_out.num_elements(),&fft_out[0],1);
1005 #else
1006  //Discrete Fourier Transform using MVM
1007  DNekBlkMatSharedPtr blkmat;
1008  blkmat = GetFloquetBlockMatrix(eForwardsPhys);
1009 
1010  int nrows = blkmat->GetRows();
1011  int ncols = blkmat->GetColumns();
1012 
1013  Array<OneD, NekDouble> sortedinarray(ncols);
1014  Array<OneD, NekDouble> sortedoutarray(nrows);
1015 
1016  //Shuffle the data
1017  for(int j= 0; j < m_slices; ++j)
1018  {
1019  Vmath::Vcopy(npoints,&m_interp[k][j*npoints],1,&(sortedinarray[j]),m_slices);
1020  }
1021 
1022  // Create NekVectors from the given data arrays
1023  NekVector<NekDouble> in (ncols,sortedinarray,eWrapper);
1024  NekVector<NekDouble> out(nrows,sortedoutarray,eWrapper);
1025 
1026  // Perform matrix-vector multiply.
1027  out = (*blkmat)*in;
1028 
1029  //Reshuffle data
1030  for(int s = 0; s < m_slices; ++s)
1031  {
1032  Vmath::Vcopy(npoints,&sortedoutarray[s],m_slices,&m_interp[k][s*npoints],1);
1033  }
1034 
1035  for(int r=0; r<sortedinarray.num_elements();++r)
1036  {
1037  sortedinarray[0]=0;
1038  sortedoutarray[0]=0;
1039  }
1040 
1041 #endif
1042 
1043  //scaling of the Fourier coefficients
1044  NekDouble j=-1;
1045  for (int i = 2; i < m_slices; i += 2)
1046  {
1047  Vmath::Smul(2*npoints,j,&m_interp[k][i*npoints],1,&m_interp[k][i*npoints],1);
1048  j=-j;
1049 
1050  }
1051 
1052  }
1053 
1054  if(m_session->DefinesParameter("period"))
1055  {
1056  m_period=m_session->GetParameter("period");
1057  }
1058  else
1059  {
1060  m_period=(m_session->GetParameter("TimeStep")*m_slices)/(m_slices-1.);
1061  }
1062  }
1063 
1064 } //end of namespace
1065 
Generated on Sun Mar 15 2015 00:11:32 for Nektar++ by   doxygen 1.8.1.2