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 
37 #include <StdRegions/StdSegExp.h>
38 
48 
49 using namespace std;
50 
51 namespace Nektar
52 {
53 string LinearisedAdvection::className
55  "Linearised",
56  LinearisedAdvection::create);
57 
58 /**
59  * Constructor. Creates ...
60  *
61  * \param
62  * \param
63  */
64 
65 LinearisedAdvection::LinearisedAdvection():
66  Advection()
67 {
68 }
69 
70 
74 {
75  Advection::v_InitObject(pSession, pFields);
76 
77  m_session = pSession;
79  ::AllocateSharedPtr(m_session, pFields[0]->GetGraph());
80  m_spacedim = pFields[0]->GetGraph()->GetSpaceDimension();
81  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 
91  if(m_session->DefinesSolverInfo("HOMOGENEOUS"))
92  {
93  std::string HomoStr = m_session->GetSolverInfo("HOMOGENEOUS");
94  m_spacedim = 3;
95 
96  if((HomoStr == "HOMOGENEOUS1D")||(HomoStr == "Homogeneous1D")||
97  (HomoStr == "1D")||(HomoStr == "Homo1D"))
98  {
100  m_LhomZ = m_session->GetParameter("LZ");
101  m_HomoDirec = 1;
102 
103  ASSERTL0(m_session->DefinesSolverInfo("ModeType"),
104  "Need to specify ModeType as HalfMode,SingleMode or "
105  "MultipleModes");
106 
107  m_session->MatchSolverInfo("ModeType", "SingleMode",
108  m_singleMode, false);
109  m_session->MatchSolverInfo("ModeType", "HalfMode",
110  m_halfMode, false);
111  m_session->MatchSolverInfo("ModeType", "MultipleModes",
112  m_multipleModes, false);
113 
114  if(m_singleMode)
115  {
116  m_npointsZ = 2;
117  }
118  else if(m_halfMode)
119  {
120  m_npointsZ = 1;
121  }
122  else if(m_multipleModes)
123  {
124  m_npointsZ = m_session->GetParameter("HomModesZ");
125  }
126  }
127 
128  if((HomoStr == "HOMOGENEOUS2D")||(HomoStr == "Homogeneous2D")||
129  (HomoStr == "2D")||(HomoStr == "Homo2D"))
130  {
132  m_session->LoadParameter("HomModesY", m_npointsY);
133  m_session->LoadParameter("LY", m_LhomY);
134  m_session->LoadParameter("HomModesZ", m_npointsZ);
135  m_session->LoadParameter("LZ", m_LhomZ);
136  m_HomoDirec = 2;
137  }
138 
139  if((HomoStr == "HOMOGENEOUS3D")||(HomoStr == "Homogeneous3D")||
140  (HomoStr == "3D")||(HomoStr == "Homo3D"))
141  {
143  m_session->LoadParameter("HomModesX",m_npointsX);
144  m_session->LoadParameter("LX", m_LhomX );
145  m_session->LoadParameter("HomModesY",m_npointsY);
146  m_session->LoadParameter("LY", m_LhomY );
147  m_session->LoadParameter("HomModesZ",m_npointsZ);
148  m_session->LoadParameter("LZ", m_LhomZ );
149  m_HomoDirec = 3;
150  }
151 
152  if(m_session->DefinesSolverInfo("USEFFT"))
153  {
154  m_useFFT = true;
155  }
156  }
157  else
158  {
159  m_npointsZ = 1; // set to default value so can use to identify 2d or 3D (homogeneous) expansions
160  }
161 
162  int nvar = m_session->GetVariables().size();
164  for (int i = 0; i < nvar; ++i)
165  {
166  m_baseflow[i] = Array<OneD, NekDouble>(pFields[i]->GetTotPoints(), 0.0);
167  }
168 
169  ASSERTL0(m_session->DefinesFunction("BaseFlow"),
170  "Base flow must be defined for linearised forms.");
171  string file = m_session->GetFunctionFilename("BaseFlow", 0);
172 
173 
174  //Periodic base flows
175  if(m_session->DefinesParameter("N_slices"))
176  {
177  m_session->LoadParameter("N_slices",m_slices);
178  if(m_slices>1)
179  {
180  DFT(file,pFields,m_slices);
181  }
182  else
183  {
184  ASSERTL0(false, "Number of slices must be a positive number "
185  "greater than 1");
186  }
187  }
188  //Steady base-flow
189  else
190  {
191  m_slices=1;
192 
193  //BaseFlow from file
194  if (m_session->GetFunctionType("BaseFlow", m_session->GetVariable(0))
196  {
197  ImportFldBase(file,pFields,1);
198 
199  }
200  //analytic base flow
201  else
202  {
203  int nq = pFields[0]->GetNpoints();
204  Array<OneD,NekDouble> x0(nq);
205  Array<OneD,NekDouble> x1(nq);
206  Array<OneD,NekDouble> x2(nq);
207 
208  // get the coordinates (assuming all fields have the same
209  // discretisation)
210  pFields[0]->GetCoords(x0,x1,x2);
211 
212  for(unsigned int i = 0 ; i < pFields.num_elements(); i++)
213  {
215  = m_session->GetFunction("BaseFlow", i);
216 
217  ifunc->Evaluate(x0,x1,x2,m_baseflow[i]);
218  }
219  }
220  }
221 
222  if(m_session->DefinesParameter("period"))
223  {
224  m_period=m_session->GetParameter("period");
225  }
226  else
227  {
228  m_period=(m_session->GetParameter("TimeStep")*m_slices)/(m_slices-1.);
229  }
230 
231 }
232 
234 {
235 }
236 
237 
238 //Advection function
239 
241  const int nConvectiveFields,
243  const Array<OneD, Array<OneD, NekDouble> > &advVel,
244  const Array<OneD, Array<OneD, NekDouble> > &inarray,
245  Array<OneD, Array<OneD, NekDouble> > &outarray,
246  const NekDouble &time)
247 {
248  ASSERTL1(nConvectiveFields == inarray.num_elements(),
249  "Number of convective fields and Inarray are not compatible");
250 
251  int nPointsTot = fields[0]->GetNpoints();
252  int ndim = advVel.num_elements();
253 
255  = Array<OneD, NekDouble> (nPointsTot*nConvectiveFields);
256 
257  Array<OneD, NekDouble> grad0,grad1,grad2;
258 
259  // Evaluation of the gradient of each component of the base flow
260  // \nabla U
261  Array<OneD, NekDouble> grad_base_u0,grad_base_u1,grad_base_u2;
262 
263  // \nabla V
264  Array<OneD, NekDouble> grad_base_v0,grad_base_v1,grad_base_v2;
265 
266  // \nabla W
267  Array<OneD, NekDouble> grad_base_w0,grad_base_w1,grad_base_w2;
268 
269  grad0 = Array<OneD, NekDouble> (nPointsTot);
270  grad_base_u0 = Array<OneD, NekDouble> (nPointsTot);
271  grad_base_v0 = Array<OneD, NekDouble> (nPointsTot);
272  grad_base_w0 = Array<OneD, NekDouble> (nPointsTot);
273 
274  // Evaluation of the base flow for periodic cases
275  if (m_slices > 1)
276  {
277  ASSERTL0(m_session->GetFunctionType("BaseFlow", 0)
279  "Base flow should be a sequence of files.");
280 
281  for (int i = 0; i < ndim; ++i)
282  {
284  time, m_period);
285  }
286  }
287 
288 
289  //Evaluate the linearised advection term
290  switch(ndim)
291  {
292  case 1: // 1D
293  ASSERTL0(false,"Not set up for 1D");
294  break;
295  case 2: //2D
296  {
297  grad1 = Array<OneD, NekDouble> (nPointsTot);
298  grad_base_u1 = Array<OneD, NekDouble> (nPointsTot);
299  grad_base_v1 = Array<OneD, NekDouble> (nPointsTot);
300 
301  //Derivates of the base flow
302  fields[0]-> PhysDeriv(m_baseflow[0], grad_base_u0, grad_base_u1);
303  fields[0]-> PhysDeriv(m_baseflow[1], grad_base_v0, grad_base_v1);
304 
305  //x-equation
306  fields[0]->PhysDeriv(inarray[0],grad0,grad1);
307  // Evaluate U du'/dx
308  Vmath::Vmul (nPointsTot,grad0,1,m_baseflow[0],1,outarray[0],1);
309  //Evaluate U du'/dx+ V du'/dy
310  Vmath::Vvtvp(nPointsTot,grad1,1,m_baseflow[1],1,outarray[0],1,
311  outarray[0],1);
312  //Evaluate (U du'/dx+ V du'/dy)+u' dU/dx
313  Vmath::Vvtvp(nPointsTot,grad_base_u0,1,advVel[0],1,outarray[0],1,
314  outarray[0],1);
315  //Evaluate (U du'/dx+ V du'/dy +u' dU/dx)+v' dU/dy
316  Vmath::Vvtvp(nPointsTot,grad_base_u1,1,advVel[1],1,outarray[0],1,
317  outarray[0],1);
318  Vmath::Neg(nPointsTot,outarray[0],1);
319 
320  //y-equation
321  fields[0]->PhysDeriv(inarray[1],grad0,grad1);
322  // Evaluate U dv'/dx
323  Vmath::Vmul (nPointsTot,grad0,1,m_baseflow[0],1,outarray[1],1);
324  //Evaluate U dv'/dx+ V dv'/dy
325  Vmath::Vvtvp(nPointsTot,grad1,1,m_baseflow[1],1,outarray[1],1,
326  outarray[1],1);
327  //Evaluate (U dv'/dx+ V dv'/dy)+u' dV/dx
328  Vmath::Vvtvp(nPointsTot,grad_base_v0,1,advVel[0],1,outarray[1],1,
329  outarray[1],1);
330  //Evaluate (U dv'/dx+ V dv'/dy +u' dv/dx)+v' dV/dy
331  Vmath::Vvtvp(nPointsTot,grad_base_v1,1,advVel[1],1,outarray[1],1,
332  outarray[1],1);
333  Vmath::Neg(nPointsTot,outarray[1],1);
334  }
335  break;
336  case 3: //3D
337  {
338  grad1 = Array<OneD, NekDouble> (nPointsTot);
339  grad2 = Array<OneD, NekDouble> (nPointsTot);
340 
341  grad_base_u1 = Array<OneD, NekDouble> (nPointsTot);
342  grad_base_v1 = Array<OneD, NekDouble> (nPointsTot);
343  grad_base_w1 = Array<OneD, NekDouble> (nPointsTot,0.0);
344 
345  grad_base_u2 = Array<OneD, NekDouble> (nPointsTot,0.0);
346  grad_base_v2 = Array<OneD, NekDouble> (nPointsTot,0.0);
347  grad_base_w2 = Array<OneD, NekDouble> (nPointsTot,0.0);
348 
349  // note base flow should be moved to GetBaseFlow method
350  if(m_halfMode) // can assume W = 0 and d/dz == 0
351  {
352  // note base flow should be moved to GetBaseFlow method
353  fields[0]->PhysDeriv(m_baseflow[0],
354  grad_base_u0, grad_base_u1);
355  fields[0]->PhysDeriv(m_baseflow[1],
356  grad_base_v0, grad_base_v1);
357  }
358  else if(m_singleMode) // single mode where d/dz = 0
359  {
360  fields[0]->PhysDeriv(m_baseflow[0], grad_base_u0,
361  grad_base_u1);
362  fields[0]->PhysDeriv(m_baseflow[1], grad_base_v0,
363  grad_base_v1);
364  fields[0]->PhysDeriv(m_baseflow[2], grad_base_w0,
365  grad_base_w1);
366  }
367  else if(m_multipleModes)
368  {
369  // Differentiate base flow in physical space
370  bool oldwavespace = fields[0]->GetWaveSpace();
371  fields[0]->SetWaveSpace(false);
372  fields[0]->PhysDeriv(m_baseflow[0], grad_base_u0,
373  grad_base_u1, grad_base_u2);
374  fields[0]->PhysDeriv(m_baseflow[1], grad_base_v0,
375  grad_base_v1, grad_base_v2);
376  fields[0]->PhysDeriv(m_baseflow[2], grad_base_w0,
377  grad_base_w1, grad_base_w2);
378  fields[0]->SetWaveSpace(oldwavespace);
379  }
380  else
381  {
382  ASSERTL0(false, "ERROR: Must be one of half, single or multiple modes");
383  }
384 
385  //x-equation
386  //
387  // Could probably clean up independent looping if clean up
388  // base flow derivative evaluation
389  if(m_multipleModes)
390  {
391  fields[0]->PhysDeriv(inarray[0], grad0, grad1, grad2);
392  // transform gradients into physical Fourier space
393  fields[0]->HomogeneousBwdTrans(grad0, grad0);
394  fields[0]->HomogeneousBwdTrans(grad1, grad1);
395  fields[0]->HomogeneousBwdTrans(grad2, grad2);
396  }
397  else
398  {
399  if(m_halfMode) //W = 0 so no need for d/dz
400  {
401  fields[0]->PhysDeriv(inarray[0], grad0, grad1);
402  }
403  else // Single mode
404  {
405  fields[0]->PhysDeriv(inarray[0], grad0, grad1, grad2);
406  }
407  }
408 
409  //Evaluate: U du'/dx
410  Vmath::Vmul (nPointsTot, grad0, 1, m_baseflow[0], 1,
411  outarray[0], 1);
412  //Evaluate and add: V du'/dy
413  Vmath::Vvtvp(nPointsTot, grad1, 1, m_baseflow[1], 1,
414  outarray[0], 1, outarray[0], 1);
415  //Evaluate and add: u' dU/dx
416  Vmath::Vvtvp(nPointsTot,grad_base_u0,1,advVel[0],1,
417  outarray[0],1,outarray[0],1);
418  //Evaluate and add: v' dU/dy
419  Vmath::Vvtvp(nPointsTot,grad_base_u1,1,advVel[1],1,
420  outarray[0],1,outarray[0],1);
421 
422  if(!m_halfMode)
423  {
424  //Evaluate an add W du'/dz
425  Vmath::Vvtvp(nPointsTot, grad2, 1, m_baseflow[2],
426  1,outarray[0], 1, outarray[0], 1);
427  }
428 
429  if(m_multipleModes)
430  {
431  //Evaluate and add w' dU/dz
432  Vmath::Vvtvp(nPointsTot,grad_base_u2,1,
433  advVel[2],1,outarray[0],1,outarray[0],1);
434  fields[0]->HomogeneousFwdTrans(outarray[0],outarray[0]);
435  }
436  Vmath::Neg(nPointsTot,outarray[0],1);
437 
438  //y-equation------------
439  if(m_multipleModes)
440  {
441  fields[0]->PhysDeriv(inarray[1], grad0, grad1, grad2);
442  // transform gradients into physical fouier space
443  fields[0]->HomogeneousBwdTrans(grad0, grad0);
444  fields[0]->HomogeneousBwdTrans(grad1, grad1);
445  fields[0]->HomogeneousBwdTrans(grad2, grad2);
446  }
447  else
448  {
449  if(m_halfMode) //W = 0 so no need for d/dz
450  {
451  fields[0]->PhysDeriv(inarray[1], grad0, grad1);
452  }
453  else // Single mode
454  {
455  fields[0]->PhysDeriv(inarray[1], grad0, grad1, grad2);
456  }
457  }
458 
459  //Evaluate U dv'/dx
460  Vmath::Vmul (nPointsTot, grad0, 1, m_baseflow[0], 1,
461  outarray[1], 1);
462  //Evaluate V dv'/dy
463  Vmath::Vvtvp(nPointsTot, grad1, 1, m_baseflow[1], 1,
464  outarray[1], 1, outarray[1], 1);
465  //Evaluate u' dV/dx
466  Vmath::Vvtvp(nPointsTot,grad_base_v0,1,advVel[0],1
467  ,outarray[1],1,outarray[1],1);
468  //Evaluate v' dV/dy
469  Vmath::Vvtvp(nPointsTot,grad_base_v1,1,advVel[1],1,
470  outarray[1],1,outarray[1],1);
471 
472  if(!m_halfMode)
473  {
474  //Evaluate W du'/dz
475  Vmath::Vvtvp(nPointsTot, grad2, 1, m_baseflow[2], 1,
476  outarray[1], 1, outarray[1], 1);
477  }
478 
479  if(m_multipleModes)
480  {
481  //Evaluate w' dV/dz
482  Vmath::Vvtvp(nPointsTot,grad_base_v2,1,advVel[2],1,
483  outarray[1],1,outarray[1],1);
484  fields[0]->HomogeneousFwdTrans(outarray[1],outarray[1]);
485  }
486 
487  Vmath::Neg(nPointsTot,outarray[1],1);
488 
489  //z-equation ------------
490  if(m_multipleModes)
491  {
492  fields[0]->PhysDeriv(inarray[2], grad0, grad1, grad2);
493  // transform gradients into physical fouier space
494  fields[0]->HomogeneousBwdTrans(grad0, grad0);
495  fields[0]->HomogeneousBwdTrans(grad1, grad1);
496  fields[0]->HomogeneousBwdTrans(grad2, grad2);
497  }
498  else
499  {
500  if(m_halfMode) //W = 0 so no need for d/dz
501  {
502  fields[0]->PhysDeriv(inarray[2], grad0, grad1);
503  }
504  else // Single mode
505  {
506  fields[0]->PhysDeriv(inarray[2], grad0, grad1, grad2);
507  }
508  }
509 
510  //Evaluate U dw'/dx
511  Vmath::Vmul (nPointsTot, grad0, 1, m_baseflow[0], 1,
512  outarray[2], 1);
513  //Evaluate V dw'/dx
514  Vmath::Vvtvp(nPointsTot, grad1, 1, m_baseflow[1], 1,
515  outarray[2], 1, outarray[2], 1);
516 
517  if(!m_halfMode) // since if halfmode W = 0
518  {
519  //Evaluate u' dW/dx
520  Vmath::Vvtvp(nPointsTot,grad_base_w0,1,advVel[0],1,
521  outarray[2],1,outarray[2],1);
522  //Evaluate v' dW/dy
523  Vmath::Vvtvp(nPointsTot,grad_base_w1,1,advVel[1],1,
524  outarray[2],1,outarray[2],1);
525 
526  //Evaluate W dw'/dz
527  Vmath::Vvtvp(nPointsTot, grad2, 1, m_baseflow[2], 1,
528  outarray[2], 1, outarray[2], 1);
529  }
530 
531  if(m_multipleModes)
532  {
533  //Evaluate w' dW/dz
534  Vmath::Vvtvp(nPointsTot,grad_base_w2,1,advVel[2],1,
535  outarray[2],1,outarray[2],1);
536  fields[0]->HomogeneousFwdTrans(outarray[2],outarray[2]);
537  }
538  Vmath::Neg(nPointsTot,outarray[2],1);
539  }
540  }
541 }
542 
544  const Array<OneD, Array<OneD, NekDouble> > &inarray)
545 {
546  if (m_session->GetSolverInfo("EqType") == "UnsteadyNavierStokes")
547  {
548  // The SFD method is only applied to the velocity variables in
549  // incompressible
550  ASSERTL1(inarray.num_elements() == (m_baseflow.num_elements() - 1),
551  "Number of base flow variables does not match what is "
552  "expected.");
553  }
554  else
555  {
556  ASSERTL1(inarray.num_elements() == (m_baseflow.num_elements()),
557  "Number of base flow variables does not match what is expected.");
558  }
559 
560  int npts = inarray[0].num_elements();
561 
562  for (int i = 0; i < inarray.num_elements(); ++i)
563  {
564  ASSERTL1(npts == m_baseflow[i].num_elements(),
565  "Size of base flow array does not match expected.");
566  Vmath::Vcopy(npts, inarray[i], 1, m_baseflow[i], 1);
567  }
568 }
569 
570 
571 /**
572  * Import field from infile and load into \a m_fields. This routine will
573  * also perform a \a BwdTrans to ensure data is in both the physical and
574  * coefficient storage.
575  * @param pInFile Filename to read.
576  * @param pFields Array of expansion lists
577  */
579  std::string pInfile,
581  int pSlice)
582 {
583  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
584  std::vector<std::vector<NekDouble> > FieldData;
585 
586  int nqtot = m_baseflow[0].num_elements();
587  Array<OneD, NekDouble> tmp_coeff(pFields[0]->GetNcoeffs(), 0.0);
588 
589  int numexp = pFields[0]->GetExpSize();
590  Array<OneD,int> ElementGIDs(numexp);
591 
592  // Define list of global element ids
593  for(int i = 0; i < numexp; ++i)
594  {
595  ElementGIDs[i] = pFields[0]->GetExp(i)->GetGeom()->GetGlobalID();
596  }
597 
600  m_session->GetComm());
601  fld->Import(pInfile, FieldDef, FieldData,
603  ElementGIDs);
604 
605  int nSessionVar = m_session->GetVariables().size();
606  int nSessionConvVar = nSessionVar - 1;
607  int nFileVar = FieldDef[0]->m_fields.size();
608  int nFileConvVar = nFileVar - 1; // Ignore pressure
609  if (m_halfMode)
610  {
611  ASSERTL0(nFileVar == 3, "For half mode, expect 2D2C base flow.");
612  nFileConvVar = 2;
613  }
614 
615  for(int j = 0; j < nFileConvVar; ++j)
616  {
617  for(int i = 0; i < FieldDef.size(); ++i)
618  {
619  bool flag = FieldDef[i]->m_fields[j] ==
620  m_session->GetVariable(j);
621 
622  ASSERTL0(flag, (std::string("Order of ") + pInfile
623  + std::string(" data and that defined in "
624  "the session file differs")).c_str());
625 
626  pFields[j]->ExtractDataToCoeffs(
627  FieldDef[i],
628  FieldData[i],
629  FieldDef[i]->m_fields[j],
630  tmp_coeff);
631  }
632 
633  if(m_singleMode || m_halfMode)
634  {
635  pFields[j]->GetPlane(0)->BwdTrans(tmp_coeff, m_baseflow[j]);
636 
637  if(m_singleMode)
638  {
639  //copy the bwd trans into the second plane for single
640  //Mode Analysis
641  int ncplane=(pFields[0]->GetNpoints())/m_npointsZ;
642  Vmath::Vcopy(ncplane,&m_baseflow[j][0],1,&m_baseflow[j][ncplane],1);
643  }
644  }
645  else // fully 3D base flow - put in physical space.
646  {
647  bool oldwavespace = pFields[j]->GetWaveSpace();
648  pFields[j]->SetWaveSpace(false);
649  pFields[j]->BwdTrans(tmp_coeff, m_baseflow[j]);
650  pFields[j]->SetWaveSpace(oldwavespace);
651  }
652  }
653 
654  // Zero unused fields (e.g. w in a 2D2C base flow).
655  for (int j = nFileConvVar; j < nSessionConvVar; ++j) {
656  Vmath::Fill(nqtot, 0.0, m_baseflow[j], 1);
657  }
658 
659  // If time-periodic, put loaded data into the slice storage.
660  if(m_session->DefinesParameter("N_slices"))
661  {
662  for(int i = 0; i < nSessionConvVar; ++i)
663  {
664  Vmath::Vcopy(nqtot, &m_baseflow[i][0], 1, &m_interp[i][pSlice*nqtot], 1);
665  }
666  }
667 }
668 
669 
671  const NekDouble m_slices,
672  const Array<OneD, const NekDouble> &inarray,
673  Array<OneD, NekDouble> &outarray,
674  const NekDouble m_time,
675  const NekDouble m_period)
676 {
677  int npoints = m_baseflow[0].num_elements();
678  NekDouble BetaT = 2*M_PI*fmod (m_time, m_period) / m_period;
679  NekDouble phase;
680  Array<OneD, NekDouble> auxiliary(npoints);
681 
682  Vmath::Vcopy(npoints,&inarray[0],1,&outarray[0],1);
683  Vmath::Svtvp(npoints, cos(0.5*m_slices*BetaT),&inarray[npoints],1,&outarray[0],1,&outarray[0],1);
684 
685  for (int i = 2; i < m_slices; i += 2)
686  {
687  phase = (i>>1) * BetaT;
688 
689  Vmath::Svtvp(npoints, cos(phase),&inarray[i*npoints],1,&outarray[0],1,&outarray[0],1);
690  Vmath::Svtvp(npoints, -sin(phase), &inarray[(i+1)*npoints], 1, &outarray[0], 1,&outarray[0],1);
691  }
692 
693 }
694 
696 {
697  DNekMatSharedPtr loc_mat;
698  DNekBlkMatSharedPtr BlkMatrix;
699  int n_exp = 0;
700 
701  n_exp = m_baseflow[0].num_elements(); // will operatore on m_phys
702 
703  Array<OneD,unsigned int> nrows(n_exp);
704  Array<OneD,unsigned int> ncols(n_exp);
705 
706  nrows = Array<OneD, unsigned int>(n_exp,m_slices);
707  ncols = Array<OneD, unsigned int>(n_exp,m_slices);
708 
709  MatrixStorage blkmatStorage = eDIAGONAL;
710  BlkMatrix = MemoryManager<DNekBlkMat>
711  ::AllocateSharedPtr(nrows,ncols,blkmatStorage);
712 
713 
716  StdRegions::StdSegExp StdSeg(BK);
717 
719  StdSeg.DetShapeType(),
720  StdSeg);
721 
722  loc_mat = StdSeg.GetStdMatrix(matkey);
723 
724  // set up array of block matrices.
725  for(int i = 0; i < n_exp; ++i)
726  {
727  BlkMatrix->SetBlock(i,i,loc_mat);
728  }
729 
730  return BlkMatrix;
731 }
732 
733 //Discrete Fourier Transform for Floquet analysis
734 void LinearisedAdvection::DFT(const string file,
736  const NekDouble m_slices)
737 {
738  int ConvectedFields = m_baseflow.num_elements()-1;
739  int npoints = m_baseflow[0].num_elements();
740  m_interp = Array<OneD, Array<OneD, NekDouble> > (ConvectedFields);
741 
742  for (int i = 0; i < ConvectedFields; ++i)
743  {
744  m_interp[i] = Array<OneD,NekDouble>(npoints*m_slices, 0.0);
745  }
746 
747  // Import the slides into auxiliary vector
748  // The base flow should be stored in the form "filename_%d.ext"
749  // A subdirectory can also be included, such as "dir/filename_%d.ext"
750  size_t found = file.find("%d");
751  ASSERTL0(found != string::npos && file.find("%d", found+1) == string::npos,
752  "Since N_slices is specified, the filename provided for function "
753  "'BaseFlow' must include exactly one instance of the format "
754  "specifier '%d', to index the time-slices.");
755  char* buffer = new char[file.length() + 8];
756  for (int i = 0; i < m_slices; ++i)
757  {
758  sprintf(buffer, file.c_str(), i);
759  ImportFldBase(buffer,pFields,i);
760  }
761  delete[] buffer;
762 
763 
764  // Discrete Fourier Transform of the fields
765  for(int k=0; k< ConvectedFields;++k)
766  {
767 #ifdef NEKTAR_USING_FFTW
768 
769  //Discrete Fourier Transform using FFTW
770  Array<OneD, NekDouble> fft_in(npoints*m_slices);
771  Array<OneD, NekDouble> fft_out(npoints*m_slices);
772 
775 
776  //Shuffle the data
777  for(int j= 0; j < m_slices; ++j)
778  {
779  Vmath::Vcopy(npoints,&m_interp[k][j*npoints],1,&(fft_in[j]),m_slices);
780  }
781 
782  m_FFT = LibUtilities::GetNektarFFTFactory().CreateInstance("NekFFTW", m_slices);
783 
784  //FFT Transform
785  for(int i=0; i<npoints; i++)
786  {
787  m_FFT->FFTFwdTrans(m_tmpIN =fft_in + i*m_slices, m_tmpOUT =fft_out + i*m_slices);
788 
789  }
790 
791  //Reshuffle data
792  for(int s = 0; s < m_slices; ++s)
793  {
794  Vmath::Vcopy(npoints,&fft_out[s],m_slices,&m_interp[k][s*npoints],1);
795 
796  }
797 
798  Vmath::Zero(fft_in.num_elements(),&fft_in[0],1);
799  Vmath::Zero(fft_out.num_elements(),&fft_out[0],1);
800 #else
801  //Discrete Fourier Transform using MVM
802  DNekBlkMatSharedPtr blkmat;
804 
805  int nrows = blkmat->GetRows();
806  int ncols = blkmat->GetColumns();
807 
808  Array<OneD, NekDouble> sortedinarray(ncols);
809  Array<OneD, NekDouble> sortedoutarray(nrows);
810 
811  //Shuffle the data
812  for(int j= 0; j < m_slices; ++j)
813  {
814  Vmath::Vcopy(npoints,&m_interp[k][j*npoints],1,&(sortedinarray[j]),m_slices);
815  }
816 
817  // Create NekVectors from the given data arrays
818  NekVector<NekDouble> in (ncols,sortedinarray,eWrapper);
819  NekVector<NekDouble> out(nrows,sortedoutarray,eWrapper);
820 
821  // Perform matrix-vector multiply.
822  out = (*blkmat)*in;
823 
824  //Reshuffle data
825  for(int s = 0; s < m_slices; ++s)
826  {
827  Vmath::Vcopy(npoints,&sortedoutarray[s],m_slices,&m_interp[k][s*npoints],1);
828  }
829 
830  for(int r=0; r<sortedinarray.num_elements();++r)
831  {
832  sortedinarray[0]=0;
833  sortedoutarray[0]=0;
834  }
835 
836 #endif
837 
838  //scaling of the Fourier coefficients
839  NekDouble j=-1;
840  for (int i = 2; i < m_slices; i += 2)
841  {
842  Vmath::Smul(2*npoints,j,&m_interp[k][i*npoints],1,&m_interp[k][i*npoints],1);
843  j=-j;
844 
845  }
846 
847  }
848 
849 }
850 
851 } //end of namespace
852 
enum HomogeneousType m_HomogeneousType
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:470
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Array< OneD, NekDouble > m_tmpIN
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.
bool m_singleMode
flag to determine if use single mode or not
Array< OneD, Array< OneD, NekDouble > > m_baseflow
Storage for base flow.
LibUtilities::NektarFFTSharedPtr m_FFT
auxiliary variables
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:46
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:57
STL namespace.
int m_slices
number of slices
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:51
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:700
Array< OneD, Array< OneD, NekDouble > > m_interp
interpolation vector
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 DFT(const std::string file, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const NekDouble m_slices)
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.
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:225
double NekDouble
NekDouble m_period
period length
DNekBlkMatSharedPtr GetFloquetBlockMatrix(FloquetMatType mattype, bool UseContCoeffs=false) const
int m_npointsY
number of points in Y direction (if homogeneous)
boost::shared_ptr< Equation > EquationSharedPtr
bool m_multipleModes
flag to determine if use multiple mode or not
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
NekDouble m_LhomY
physical length in Y direction (if homogeneous)
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)
bool m_halfMode
flag to determine if use half mode or not
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:1047
static FieldMetaDataMap NullFieldMetaDataMap
Definition: FieldIO.h:54
bool m_useFFT
flag to determine if use or not the FFT for transformations
Describes the specification for a Basis.
Definition: Basis.h:50
virtual SOLVER_UTILS_EXPORT void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Initialises the advection object.
Definition: Advection.cpp:97
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
Defines a callback function which evaluates the flux vector.
Definition: Advection.h:69
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215