Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties 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 
39 using namespace std;
40 
41 namespace Nektar
42 {
43 string LinearisedAdvection::className
45  "Linearised",
46  LinearisedAdvection::create);
47 
48 /**
49  * Constructor. Creates ...
50  *
51  * \param
52  * \param
53  */
54 
55 LinearisedAdvection::LinearisedAdvection():
56  Advection()
57 {
58 }
59 
60 
64 {
65  Advection::v_InitObject(pSession, pFields);
66 
67  m_session = pSession;
69  ::AllocateSharedPtr(m_session, pFields[0]->GetGraph());
70  m_spacedim = pFields[0]->GetGraph()->GetSpaceDimension();
71  m_expdim = pFields[0]->GetGraph()->GetMeshDimension();
72 
73  //Setting parameters for homogeneous problems
74  m_HomoDirec = 0;
75  m_useFFT = false;
77  m_singleMode = false;
78  m_halfMode = false;
79  m_multipleModes = false;
80 
81  if(m_session->DefinesSolverInfo("HOMOGENEOUS"))
82  {
83  std::string HomoStr = m_session->GetSolverInfo("HOMOGENEOUS");
84  m_spacedim = 3;
85 
86  if((HomoStr == "HOMOGENEOUS1D")||(HomoStr == "Homogeneous1D")||
87  (HomoStr == "1D")||(HomoStr == "Homo1D"))
88  {
90  m_LhomZ = m_session->GetParameter("LZ");
91  m_HomoDirec = 1;
92 
93  ASSERTL0(m_session->DefinesSolverInfo("ModeType"),
94  "Need to specify ModeType as HalfMode,SingleMode or "
95  "MultipleModes");
96 
97  m_session->MatchSolverInfo("ModeType", "SingleMode",
98  m_singleMode, false);
99  m_session->MatchSolverInfo("ModeType", "HalfMode",
100  m_halfMode, false);
101  m_session->MatchSolverInfo("ModeType", "MultipleModes",
102  m_multipleModes, false);
103 
104  if(m_singleMode)
105  {
106  m_npointsZ = 2;
107  }
108  else if(m_halfMode)
109  {
110  m_npointsZ = 1;
111  }
112  else if(m_multipleModes)
113  {
114  m_npointsZ = m_session->GetParameter("HomModesZ");
115  }
116  }
117 
118  if((HomoStr == "HOMOGENEOUS2D")||(HomoStr == "Homogeneous2D")||
119  (HomoStr == "2D")||(HomoStr == "Homo2D"))
120  {
122  m_session->LoadParameter("HomModesY", m_npointsY);
123  m_session->LoadParameter("LY", m_LhomY);
124  m_session->LoadParameter("HomModesZ", m_npointsZ);
125  m_session->LoadParameter("LZ", m_LhomZ);
126  m_HomoDirec = 2;
127  }
128 
129  if((HomoStr == "HOMOGENEOUS3D")||(HomoStr == "Homogeneous3D")||
130  (HomoStr == "3D")||(HomoStr == "Homo3D"))
131  {
133  m_session->LoadParameter("HomModesX",m_npointsX);
134  m_session->LoadParameter("LX", m_LhomX );
135  m_session->LoadParameter("HomModesY",m_npointsY);
136  m_session->LoadParameter("LY", m_LhomY );
137  m_session->LoadParameter("HomModesZ",m_npointsZ);
138  m_session->LoadParameter("LZ", m_LhomZ );
139  m_HomoDirec = 3;
140  }
141 
142  if(m_session->DefinesSolverInfo("USEFFT"))
143  {
144  m_useFFT = true;
145  }
146  }
147  else
148  {
149  m_npointsZ = 1; // set to default value so can use to identify 2d or 3D (homogeneous) expansions
150  }
151 
152  int nvar = m_session->GetVariables().size();
154  for (int i = 0; i < nvar; ++i)
155  {
156  m_baseflow[i] = Array<OneD, NekDouble>(pFields[i]->GetTotPoints(), 0.0);
157  }
158 
159  int nBaseDerivs = (m_halfMode || m_singleMode) ? 2 : m_spacedim;
160  m_gradBase = Array<OneD, Array<OneD, NekDouble> >(nvar*nBaseDerivs);
161  for (int i = 0; i < nvar; ++i)
162  {
163  for (int j = 0; j < nBaseDerivs; ++j)
164  {
165  m_gradBase[i*nBaseDerivs + j ] = Array<OneD, NekDouble>
166  (pFields[i]->GetTotPoints(), 0.0);
167  }
168  }
169 
170  ASSERTL0(m_session->DefinesFunction("BaseFlow"),
171  "Base flow must be defined for linearised forms.");
172  string file = m_session->GetFunctionFilename("BaseFlow", 0);
173 
174 
175  //Periodic base flows
176  if(m_session->DefinesParameter("N_slices"))
177  {
178  m_session->LoadParameter("N_slices",m_slices);
179  if(m_slices>1)
180  {
181  ASSERTL0(m_session->GetFunctionType("BaseFlow", 0)
183  "Base flow should be a sequence of files.");
184  DFT(file,pFields,m_slices);
185  }
186  else
187  {
188  ASSERTL0(false, "Number of slices must be a positive number "
189  "greater than 1");
190  }
191  }
192  //Steady base-flow
193  else
194  {
195  m_slices=1;
196 
197  //BaseFlow from file
198  if (m_session->GetFunctionType("BaseFlow", m_session->GetVariable(0))
200  {
201  ImportFldBase(file,pFields,1);
202 
203  }
204  //analytic base flow
205  else
206  {
207  int nq = pFields[0]->GetNpoints();
208  Array<OneD,NekDouble> x0(nq);
209  Array<OneD,NekDouble> x1(nq);
210  Array<OneD,NekDouble> x2(nq);
211 
212  // get the coordinates (assuming all fields have the same
213  // discretisation)
214  pFields[0]->GetCoords(x0,x1,x2);
215 
216  for(unsigned int i = 0 ; i < pFields.num_elements(); i++)
217  {
219  = m_session->GetFunction("BaseFlow", i);
220 
221  ifunc->Evaluate(x0,x1,x2,m_baseflow[i]);
222  }
223  }
224  }
225 
226  for (int i = 0; i < nvar; ++i)
227  {
228  UpdateGradBase(i, pFields[i]);
229  }
230 
231  if(m_session->DefinesParameter("period"))
232  {
233  m_period=m_session->GetParameter("period");
234  }
235  else
236  {
237  m_period=(m_session->GetParameter("TimeStep")*m_slices)/(m_slices-1.);
238  }
239 
240 }
241 
243 {
244 }
245 
246 
247 //Advection function
248 
250  const int nConvectiveFields,
252  const Array<OneD, Array<OneD, NekDouble> > &advVel,
253  const Array<OneD, Array<OneD, NekDouble> > &inarray,
254  Array<OneD, Array<OneD, NekDouble> > &outarray,
255  const NekDouble &time,
256  const Array<OneD, Array<OneD, NekDouble> > &pFwd,
257  const Array<OneD, Array<OneD, NekDouble> > &pBwd)
258 {
259  ASSERTL1(nConvectiveFields == inarray.num_elements(),
260  "Number of convective fields and Inarray are not compatible");
261 
262  int nPointsTot = fields[0]->GetNpoints();
263  int ndim = advVel.num_elements();
264  int nBaseDerivs = (m_halfMode || m_singleMode) ? 2 : m_spacedim;
265  int nDerivs = (m_halfMode) ? 2 : m_spacedim;
266 
267  Array<OneD, Array<OneD, NekDouble> > velocity(ndim);
268  for(int i = 0; i < ndim; ++i)
269  {
270  if(fields[i]->GetWaveSpace() && !m_singleMode && !m_halfMode)
271  {
272  velocity[i] = Array<OneD, NekDouble>(nPointsTot,0.0);
273  fields[i]->HomogeneousBwdTrans(advVel[i],velocity[i]);
274  }
275  else
276  {
277  velocity[i] = advVel[i];
278  }
279  }
280 
281  Array<OneD, Array<OneD, NekDouble> > grad (nDerivs);
282  for( int i = 0; i < nDerivs; ++i)
283  {
284  grad[i] = Array<OneD, NekDouble> (nPointsTot);
285  }
286 
287  // Evaluation of the base flow for periodic cases
288  if (m_slices > 1)
289  {
290  for (int i = 0; i < ndim; ++i)
291  {
293  time, m_period);
294  UpdateGradBase(i, fields[i]);
295  }
296  }
297 
298  //Evaluate the linearised advection term
299  for( int i = 0; i < ndim; ++i)
300  {
301  // Calculate gradient
302  switch(nDerivs)
303  {
304  case 1:
305  {
306  fields[i]->PhysDeriv(inarray[i], grad[0]);
307  }
308  break;
309  case 2:
310  {
311  fields[i]->PhysDeriv(inarray[i], grad[0], grad[1]);
312  }
313  break;
314  case 3:
315  {
316  fields[i]->PhysDeriv(inarray[i], grad[0], grad[1], grad[2]);
317  if(m_multipleModes)
318  {
319  // transform gradients into physical Fourier space
320  fields[i]->HomogeneousBwdTrans(grad[0], grad[0]);
321  fields[i]->HomogeneousBwdTrans(grad[1], grad[1]);
322  fields[i]->HomogeneousBwdTrans(grad[2], grad[2]);
323  }
324  }
325  break;
326  }
327 
328  // Calculate U_j du'_i/dx_j
329  Vmath::Vmul(nPointsTot,grad[0], 1, m_baseflow[0], 1, outarray[i], 1);
330  for( int j = 1; j < nDerivs; ++j)
331  {
332  Vmath::Vvtvp(nPointsTot,grad[j], 1,
333  m_baseflow[j], 1,
334  outarray[i], 1,
335  outarray[i], 1);
336  }
337 
338  // Add u'_j dU_i/dx_j
339  int lim = (m_halfMode || m_singleMode) ? 2 : ndim;
340  if (m_halfMode && i==2)
341  {
342  lim = 0;
343  }
344  for( int j = 0; j < lim; ++j)
345  {
346  Vmath::Vvtvp(nPointsTot,m_gradBase[i*nBaseDerivs + j], 1,
347  velocity[j], 1,
348  outarray[i], 1,
349  outarray[i], 1);
350  }
351 
352  if(m_multipleModes)
353  {
354  fields[i]->HomogeneousFwdTrans(outarray[i],outarray[i]);
355  }
356  Vmath::Neg(nPointsTot,outarray[i],1);
357  }
358 }
359 
361  const Array<OneD, Array<OneD, NekDouble> > &inarray,
363 {
364  if (m_session->GetSolverInfo("EqType") == "UnsteadyNavierStokes")
365  {
366  // The SFD method is only applied to the velocity variables in
367  // incompressible
368  ASSERTL1(inarray.num_elements() == (m_baseflow.num_elements() - 1),
369  "Number of base flow variables does not match what is "
370  "expected.");
371  }
372  else
373  {
374  ASSERTL1(inarray.num_elements() == (m_baseflow.num_elements()),
375  "Number of base flow variables does not match what is expected.");
376  }
377 
378  int npts = inarray[0].num_elements();
379 
380  for (int i = 0; i < inarray.num_elements(); ++i)
381  {
382  ASSERTL1(npts == m_baseflow[i].num_elements(),
383  "Size of base flow array does not match expected.");
384  Vmath::Vcopy(npts, inarray[i], 1, m_baseflow[i], 1);
385  UpdateGradBase(i, fields[i]);
386  }
387 }
388 
389 
390 /**
391  * Import field from infile and load into \a m_fields. This routine will
392  * also perform a \a BwdTrans to ensure data is in both the physical and
393  * coefficient storage.
394  * @param pInFile Filename to read.
395  * @param pFields Array of expansion lists
396  */
398  std::string pInfile,
400  int pSlice)
401 {
402  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
403  std::vector<std::vector<NekDouble> > FieldData;
404 
405  int nqtot = m_baseflow[0].num_elements();
406  Array<OneD, NekDouble> tmp_coeff(pFields[0]->GetNcoeffs(), 0.0);
407 
408  int numexp = pFields[0]->GetExpSize();
409  Array<OneD,int> ElementGIDs(numexp);
410 
411  // Define list of global element ids
412  for(int i = 0; i < numexp; ++i)
413  {
414  ElementGIDs[i] = pFields[0]->GetExp(i)->GetGeom()->GetGlobalID();
415  }
416 
417  //Get Homogeneous
419  m_session, pInfile);
420  fld->Import(pInfile, FieldDef, FieldData,
422  ElementGIDs);
423 
424  int nSessionVar = m_session->GetVariables().size();
425  int nSessionConvVar = nSessionVar - 1;
426  int nFileVar = FieldDef[0]->m_fields.size();
427  int nFileConvVar = nFileVar - 1; // Ignore pressure
428  if (m_halfMode)
429  {
430  ASSERTL0(nFileVar == 3, "For half mode, expect 2D2C base flow.");
431  nFileConvVar = 2;
432  }
433 
434  for(int j = 0; j < nFileConvVar; ++j)
435  {
436  for(int i = 0; i < FieldDef.size(); ++i)
437  {
438  bool flag = FieldDef[i]->m_fields[j] ==
439  m_session->GetVariable(j);
440 
441  ASSERTL0(flag, (std::string("Order of ") + pInfile
442  + std::string(" data and that defined in "
443  "the session file differs")).c_str());
444 
445  pFields[j]->ExtractDataToCoeffs(
446  FieldDef[i],
447  FieldData[i],
448  FieldDef[i]->m_fields[j],
449  tmp_coeff);
450  }
451 
452  if(m_singleMode || m_halfMode)
453  {
454  pFields[j]->GetPlane(0)->BwdTrans(tmp_coeff, m_baseflow[j]);
455 
456  if(m_singleMode)
457  {
458  //copy the bwd trans into the second plane for single
459  //Mode Analysis
460  int ncplane=(pFields[0]->GetNpoints())/m_npointsZ;
461  Vmath::Vcopy(ncplane,&m_baseflow[j][0],1,&m_baseflow[j][ncplane],1);
462  }
463  }
464  else // fully 3D base flow - put in physical space.
465  {
466  bool oldwavespace = pFields[j]->GetWaveSpace();
467  pFields[j]->SetWaveSpace(false);
468  pFields[j]->BwdTrans(tmp_coeff, m_baseflow[j]);
469  pFields[j]->SetWaveSpace(oldwavespace);
470  }
471  }
472 
473  // Zero unused fields (e.g. w in a 2D2C base flow).
474  for (int j = nFileConvVar; j < nSessionConvVar; ++j)
475  {
476  Vmath::Fill(nqtot, 0.0, m_baseflow[j], 1);
477  }
478 
479  // If time-periodic, put loaded data into the slice storage.
480  if(m_session->DefinesParameter("N_slices"))
481  {
482  for(int i = 0; i < nSessionConvVar; ++i)
483  {
484  Vmath::Vcopy(nqtot, &m_baseflow[i][0], 1, &m_interp[i][pSlice*nqtot], 1);
485  }
486  }
487 }
488 
489 
491  const NekDouble m_slices,
492  const Array<OneD, const NekDouble> &inarray,
493  Array<OneD, NekDouble> &outarray,
494  const NekDouble m_time,
495  const NekDouble m_period)
496 {
497  int npoints = m_baseflow[0].num_elements();
498  NekDouble BetaT = 2*M_PI*fmod (m_time, m_period) / m_period;
499  NekDouble phase;
500  Array<OneD, NekDouble> auxiliary(npoints);
501 
502  Vmath::Vcopy(npoints,&inarray[0],1,&outarray[0],1);
503  Vmath::Svtvp(npoints, cos(0.5*m_slices*BetaT),&inarray[npoints],1,&outarray[0],1,&outarray[0],1);
504 
505  for (int i = 2; i < m_slices; i += 2)
506  {
507  phase = (i>>1) * BetaT;
508 
509  Vmath::Svtvp(npoints, cos(phase),&inarray[i*npoints],1,&outarray[0],1,&outarray[0],1);
510  Vmath::Svtvp(npoints, -sin(phase), &inarray[(i+1)*npoints], 1, &outarray[0], 1,&outarray[0],1);
511  }
512 
513 }
514 
516  const int var,
517  const MultiRegions::ExpListSharedPtr &field)
518 {
519  int nBaseDerivs = (m_halfMode || m_singleMode) ? 2 : m_spacedim;
520  switch(m_spacedim)
521  {
522  case 1: // 1D
523  {
524  field->PhysDeriv(m_baseflow[var], m_gradBase[var*nBaseDerivs + 0]);
525  }
526  break;
527  case 2: //2D
528  {
529  field->PhysDeriv(m_baseflow[var], m_gradBase[var*nBaseDerivs + 0],
530  m_gradBase[var*nBaseDerivs + 1]);
531  }
532  break;
533  case 3:
534  {
535  if(m_halfMode) // can assume W = 0 and d/dz == 0
536  {
537  if( var < 2)
538  {
539  field->PhysDeriv(m_baseflow[var],
540  m_gradBase[var*nBaseDerivs + 0],
541  m_gradBase[var*nBaseDerivs + 1]);
542  }
543  }
544  else if(m_singleMode) // single mode where d/dz = 0
545  {
546  field->PhysDeriv(m_baseflow[var], m_gradBase[var*nBaseDerivs + 0],
547  m_gradBase[var*nBaseDerivs + 1]);
548  }
549  else
550  {
551  // Differentiate base flow in physical space
552  bool oldwavespace = field->GetWaveSpace();
553  field->SetWaveSpace(false);
554  field->PhysDeriv(m_baseflow[var], m_gradBase[var*nBaseDerivs + 0],
555  m_gradBase[var*nBaseDerivs + 1],
556  m_gradBase[var*nBaseDerivs + 2]);
557  field->SetWaveSpace(oldwavespace);
558  }
559  }
560  break;
561  }
562 }
563 
565 {
566  DNekMatSharedPtr loc_mat;
567  DNekBlkMatSharedPtr BlkMatrix;
568  int n_exp = 0;
569 
570  n_exp = m_baseflow[0].num_elements(); // will operatore on m_phys
571 
572  Array<OneD,unsigned int> nrows(n_exp);
573  Array<OneD,unsigned int> ncols(n_exp);
574 
575  nrows = Array<OneD, unsigned int>(n_exp,m_slices);
576  ncols = Array<OneD, unsigned int>(n_exp,m_slices);
577 
578  MatrixStorage blkmatStorage = eDIAGONAL;
579  BlkMatrix = MemoryManager<DNekBlkMat>
580  ::AllocateSharedPtr(nrows,ncols,blkmatStorage);
581 
582 
585  StdRegions::StdSegExp StdSeg(BK);
586 
588  StdSeg.DetShapeType(),
589  StdSeg);
590 
591  loc_mat = StdSeg.GetStdMatrix(matkey);
592 
593  // set up array of block matrices.
594  for(int i = 0; i < n_exp; ++i)
595  {
596  BlkMatrix->SetBlock(i,i,loc_mat);
597  }
598 
599  return BlkMatrix;
600 }
601 
602 //Discrete Fourier Transform for Floquet analysis
603 void LinearisedAdvection::DFT(const string file,
605  const NekDouble m_slices)
606 {
607  int ConvectedFields = m_baseflow.num_elements()-1;
608  int npoints = m_baseflow[0].num_elements();
609  m_interp = Array<OneD, Array<OneD, NekDouble> > (ConvectedFields);
610 
611  for (int i = 0; i < ConvectedFields; ++i)
612  {
613  m_interp[i] = Array<OneD,NekDouble>(npoints*m_slices, 0.0);
614  }
615 
616  // Import the slides into auxiliary vector
617  // The base flow should be stored in the form "filename_%d.ext"
618  // A subdirectory can also be included, such as "dir/filename_%d.ext"
619  size_t found = file.find("%d");
620  ASSERTL0(found != string::npos && file.find("%d", found+1) == string::npos,
621  "Since N_slices is specified, the filename provided for function "
622  "'BaseFlow' must include exactly one instance of the format "
623  "specifier '%d', to index the time-slices.");
624  char* buffer = new char[file.length() + 8];
625  for (int i = 0; i < m_slices; ++i)
626  {
627  sprintf(buffer, file.c_str(), i);
628  ImportFldBase(buffer,pFields,i);
629  }
630  delete[] buffer;
631 
632 
633  // Discrete Fourier Transform of the fields
634  for(int k=0; k< ConvectedFields;++k)
635  {
636 #ifdef NEKTAR_USING_FFTW
637 
638  //Discrete Fourier Transform using FFTW
639  Array<OneD, NekDouble> fft_in(npoints*m_slices);
640  Array<OneD, NekDouble> fft_out(npoints*m_slices);
641 
644 
645  //Shuffle the data
646  for(int j= 0; j < m_slices; ++j)
647  {
648  Vmath::Vcopy(npoints,&m_interp[k][j*npoints],1,&(fft_in[j]),m_slices);
649  }
650 
651  m_FFT = LibUtilities::GetNektarFFTFactory().CreateInstance("NekFFTW", m_slices);
652 
653  //FFT Transform
654  for(int i=0; i<npoints; i++)
655  {
656  m_FFT->FFTFwdTrans(m_tmpIN =fft_in + i*m_slices, m_tmpOUT =fft_out + i*m_slices);
657 
658  }
659 
660  //Reshuffle data
661  for(int s = 0; s < m_slices; ++s)
662  {
663  Vmath::Vcopy(npoints,&fft_out[s],m_slices,&m_interp[k][s*npoints],1);
664 
665  }
666 
667  Vmath::Zero(fft_in.num_elements(),&fft_in[0],1);
668  Vmath::Zero(fft_out.num_elements(),&fft_out[0],1);
669 #else
670  //Discrete Fourier Transform using MVM
671  DNekBlkMatSharedPtr blkmat;
673 
674  int nrows = blkmat->GetRows();
675  int ncols = blkmat->GetColumns();
676 
677  Array<OneD, NekDouble> sortedinarray(ncols);
678  Array<OneD, NekDouble> sortedoutarray(nrows);
679 
680  //Shuffle the data
681  for(int j= 0; j < m_slices; ++j)
682  {
683  Vmath::Vcopy(npoints,&m_interp[k][j*npoints],1,&(sortedinarray[j]),m_slices);
684  }
685 
686  // Create NekVectors from the given data arrays
687  NekVector<NekDouble> in (ncols,sortedinarray,eWrapper);
688  NekVector<NekDouble> out(nrows,sortedoutarray,eWrapper);
689 
690  // Perform matrix-vector multiply.
691  out = (*blkmat)*in;
692 
693  //Reshuffle data
694  for(int s = 0; s < m_slices; ++s)
695  {
696  Vmath::Vcopy(npoints,&sortedoutarray[s],m_slices,&m_interp[k][s*npoints],1);
697  }
698 
699  for(int r=0; r<sortedinarray.num_elements();++r)
700  {
701  sortedinarray[0]=0;
702  sortedoutarray[0]=0;
703  }
704 
705 #endif
706 
707  //scaling of the Fourier coefficients
708  NekDouble j=-1;
709  for (int i = 2; i < m_slices; i += 2)
710  {
711  Vmath::Smul(2*npoints,j,&m_interp[k][i*npoints],1,&m_interp[k][i*npoints],1);
712  j=-j;
713 
714  }
715 
716  }
717 
718 }
719 
720 } //end of namespace
721 
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:198
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.
virtual void v_SetBaseFlow(const Array< OneD, Array< OneD, NekDouble > > &inarray, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields)
Overrides the base flow used during linearised advection.
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:485
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:442
array buffer
Definition: GsLib.hpp:60
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
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:706
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:66
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:213
int m_npointsZ
number of points in Z direction (if homogeneous)
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
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:396
boost::shared_ptr< FieldIO > FieldIOSharedPtr
Definition: FieldIO.h:309
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)
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, const Array< OneD, Array< OneD, NekDouble > > &pFwd=NullNekDoubleArrayofArray, const Array< OneD, Array< OneD, NekDouble > > &pBwd=NullNekDoubleArrayofArray)
Advects a vector field.
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
static boost::shared_ptr< FieldIO > CreateForFile(const LibUtilities::SessionReaderSharedPtr session, const std::string &filename)
Construct a FieldIO object for a given input filename.
Definition: FieldIO.cpp:212
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 UpdateGradBase(const int var, const MultiRegions::ExpListSharedPtr &field)
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:373
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:228
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1061
static FieldMetaDataMap NullFieldMetaDataMap
Definition: FieldIO.h:55
bool m_useFFT
flag to determine if use or not the FFT for transformations
Describes the specification for a Basis.
Definition: Basis.h:50
Array< OneD, Array< OneD, NekDouble > > m_gradBase
virtual SOLVER_UTILS_EXPORT void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Initialises the advection object.
Definition: Advection.cpp:100
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:183
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