Nektar++
LinearisedAdvection.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File LinearisedAdvection.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: Evaluation of the linearised advective term
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
36 #include <StdRegions/StdSegExp.h>
37 
38 using namespace std;
39 
40 namespace Nektar
41 {
42 string LinearisedAdvection::className =
44  "Linearised", LinearisedAdvection::create,
45  "Linearised Non-Conservative");
46 
47 /**
48  * Constructor. Creates ...
49  *
50  * \param
51  * \param
52  */
53 
54 LinearisedAdvection::LinearisedAdvection() : Advection()
55 {
56 }
57 
61 {
62  Advection::v_InitObject(pSession, pFields);
63 
64  m_session = pSession;
67  m_session, pFields[0]->GetGraph());
68  m_spacedim = pFields[0]->GetGraph()->GetSpaceDimension();
69  m_expdim = pFields[0]->GetGraph()->GetMeshDimension();
70 
71  // Setting parameters for homogeneous problems
72  m_HomoDirec = 0;
73  m_useFFT = false;
75  m_singleMode = false;
76  m_halfMode = false;
77  m_multipleModes = false;
78 
79  if (m_session->DefinesSolverInfo("HOMOGENEOUS"))
80  {
81  std::string HomoStr = m_session->GetSolverInfo("HOMOGENEOUS");
82  m_spacedim = 3;
83 
84  if ((HomoStr == "HOMOGENEOUS1D") || (HomoStr == "Homogeneous1D") ||
85  (HomoStr == "1D") || (HomoStr == "Homo1D"))
86  {
88  m_LhomZ = m_session->GetParameter("LZ");
89  m_HomoDirec = 1;
90 
91  ASSERTL0(m_session->DefinesSolverInfo("ModeType"),
92  "Need to specify ModeType as HalfMode,SingleMode or "
93  "MultipleModes");
94 
95  m_session->MatchSolverInfo("ModeType", "SingleMode", m_singleMode,
96  false);
97  m_session->MatchSolverInfo("ModeType", "HalfMode", m_halfMode,
98  false);
99  m_session->MatchSolverInfo("ModeType", "MultipleModes",
100  m_multipleModes, false);
101 
102  if (m_singleMode)
103  {
104  m_npointsZ = 2;
105  }
106  else if (m_halfMode)
107  {
108  m_npointsZ = 1;
109  }
110  else if (m_multipleModes)
111  {
112  m_npointsZ = m_session->GetParameter("HomModesZ");
113  }
114  }
115 
116  if ((HomoStr == "HOMOGENEOUS2D") || (HomoStr == "Homogeneous2D") ||
117  (HomoStr == "2D") || (HomoStr == "Homo2D"))
118  {
120  m_session->LoadParameter("HomModesY", m_npointsY);
121  m_session->LoadParameter("LY", m_LhomY);
122  m_session->LoadParameter("HomModesZ", m_npointsZ);
123  m_session->LoadParameter("LZ", m_LhomZ);
124  m_HomoDirec = 2;
125  }
126 
127  if ((HomoStr == "HOMOGENEOUS3D") || (HomoStr == "Homogeneous3D") ||
128  (HomoStr == "3D") || (HomoStr == "Homo3D"))
129  {
131  m_session->LoadParameter("HomModesX", m_npointsX);
132  m_session->LoadParameter("LX", m_LhomX);
133  m_session->LoadParameter("HomModesY", m_npointsY);
134  m_session->LoadParameter("LY", m_LhomY);
135  m_session->LoadParameter("HomModesZ", m_npointsZ);
136  m_session->LoadParameter("LZ", m_LhomZ);
137  m_HomoDirec = 3;
138  }
139 
140  if (m_session->DefinesSolverInfo("USEFFT"))
141  {
142  m_useFFT = true;
143  }
144  }
145  else
146  {
147  m_npointsZ = 1; // set to default value so can use to identify 2d or 3D
148  // (homogeneous) expansions
149  }
150 
151  int nvar = m_session->GetVariables().size();
153  for (int i = 0; i < nvar; ++i)
154  {
155  m_baseflow[i] = Array<OneD, NekDouble>(pFields[i]->GetTotPoints(), 0.0);
156  }
157 
158  int nBaseDerivs = (m_halfMode || m_singleMode) ? 2 : m_spacedim;
159  m_gradBase = Array<OneD, Array<OneD, NekDouble>>(nvar * nBaseDerivs);
160  for (int i = 0; i < nvar; ++i)
161  {
162  for (int j = 0; j < nBaseDerivs; ++j)
163  {
164  m_gradBase[i * nBaseDerivs + j] =
165  Array<OneD, NekDouble>(pFields[i]->GetTotPoints(), 0.0);
166  }
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  // Periodic base flows
174  if (m_session->DefinesParameter("N_slices"))
175  {
176  m_session->LoadParameter("N_slices", m_slices);
177  if (m_slices > 1)
178  {
179  ASSERTL0(m_session->GetFunctionType("BaseFlow", 0) ==
181  "Base flow should be a sequence of files.");
182  m_session->LoadParameter("BaseFlow_interporder", m_interporder, 0);
184  "BaseFlow_interporder should be smaller than or equal to "
185  "N_slices.");
187  m_session->LoadParameter("N_start", m_start, 0);
188  m_session->LoadParameter("N_skip", m_skip, 1);
189  DFT(file, pFields, m_slices);
190  }
191  else
192  {
193  ASSERTL0(false, "Number of slices must be a positive number "
194  "greater than 1");
195  }
196  }
197  // Steady base-flow
198  else
199  {
200  m_slices = 1;
201 
202  // BaseFlow from file
203  if (m_session->GetFunctionType("BaseFlow", m_session->GetVariable(0)) ==
205  {
206  ImportFldBase(file, pFields, 1);
207  }
208  // analytic base flow
209  else
210  {
211  int nq = pFields[0]->GetNpoints();
212  Array<OneD, NekDouble> x0(nq);
213  Array<OneD, NekDouble> x1(nq);
214  Array<OneD, NekDouble> x2(nq);
215 
216  // get the coordinates (assuming all fields have the same
217  // discretisation)
218  pFields[0]->GetCoords(x0, x1, x2);
219 
220  for (unsigned int i = 0; i < pFields.size(); i++)
221  {
223  m_session->GetFunction("BaseFlow", i);
224 
225  ifunc->Evaluate(x0, x1, x2, m_baseflow[i]);
226  }
227  }
228  }
229 
230  for (int i = 0; i < nvar; ++i)
231  {
232  UpdateGradBase(i, pFields[i]);
233  }
234 
235  if (m_session->DefinesParameter("period"))
236  {
237  m_period = m_session->GetParameter("period");
238  }
239  else
240  {
241  m_period =
242  (m_session->GetParameter("TimeStep") * m_slices) / (m_slices - 1.);
243  }
244  if (m_session->GetComm()->GetRank() == 0)
245  {
246  cout << "baseflow info : interpolation order " << m_interporder
247  << ", period " << m_period << ", periodicity ";
248  if (m_isperiodic)
249  {
250  cout << "yes\n";
251  }
252  else
253  {
254  cout << "no\n";
255  }
256  cout << "baseflow info : files from " << m_start << " to "
257  << (m_start + (m_slices - 1) * m_skip) << " (skip " << m_skip
258  << ") with " << (m_slices - (m_interporder > 1))
259  << " time intervals" << endl;
260  }
261 }
262 
264 {
265 }
266 
267 // Advection function
268 
270  const int nConvectiveFields,
272  const Array<OneD, Array<OneD, NekDouble>> &advVel,
273  const Array<OneD, Array<OneD, NekDouble>> &inarray,
274  Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble &time,
275  const Array<OneD, Array<OneD, NekDouble>> &pFwd,
276  const Array<OneD, Array<OneD, NekDouble>> &pBwd)
277 {
278  ASSERTL1(nConvectiveFields == inarray.size(),
279  "Number of convective fields and Inarray are not compatible");
280 
281  int nPointsTot = fields[0]->GetNpoints();
282  int ndim = advVel.size();
283  int nBaseDerivs = (m_halfMode || m_singleMode) ? 2 : m_spacedim;
284  int nDerivs = (m_halfMode) ? 2 : m_spacedim;
285 
286  Array<OneD, Array<OneD, NekDouble>> velocity(ndim);
287  for (int i = 0; i < ndim; ++i)
288  {
289  if (fields[i]->GetWaveSpace() && !m_singleMode && !m_halfMode)
290  {
291  velocity[i] = Array<OneD, NekDouble>(nPointsTot, 0.0);
292  fields[i]->HomogeneousBwdTrans(advVel[i], velocity[i]);
293  }
294  else
295  {
296  velocity[i] = advVel[i];
297  }
298  }
299 
300  Array<OneD, Array<OneD, NekDouble>> grad(nDerivs);
301  for (int i = 0; i < nDerivs; ++i)
302  {
303  grad[i] = Array<OneD, NekDouble>(nPointsTot);
304  }
305 
306  // Evaluation of the base flow for periodic cases
307  if (m_slices > 1)
308  {
309  for (int i = 0; i < ndim; ++i)
310  {
312  UpdateGradBase(i, fields[i]);
313  }
314  }
315 
316  // Evaluate the linearised advection term
317  for (int i = 0; i < nConvectiveFields; ++i)
318  {
319  // Calculate gradient
320  switch (nDerivs)
321  {
322  case 1:
323  {
324  fields[i]->PhysDeriv(inarray[i], grad[0]);
325  }
326  break;
327  case 2:
328  {
329  fields[i]->PhysDeriv(inarray[i], grad[0], grad[1]);
330  }
331  break;
332  case 3:
333  {
334  fields[i]->PhysDeriv(inarray[i], grad[0], grad[1], grad[2]);
335  if (m_multipleModes)
336  {
337  // transform gradients into physical Fourier space
338  fields[i]->HomogeneousBwdTrans(grad[0], grad[0]);
339  fields[i]->HomogeneousBwdTrans(grad[1], grad[1]);
340  fields[i]->HomogeneousBwdTrans(grad[2], grad[2]);
341  }
342  }
343  break;
344  }
345 
346  // Calculate U_j du'_i/dx_j
347  Vmath::Vmul(nPointsTot, grad[0], 1, m_baseflow[0], 1, outarray[i], 1);
348  for (int j = 1; j < nDerivs; ++j)
349  {
350  Vmath::Vvtvp(nPointsTot, grad[j], 1, m_baseflow[j], 1, outarray[i],
351  1, outarray[i], 1);
352  }
353 
354  // Add u'_j dU_i/dx_j
355  int lim = (m_halfMode || m_singleMode) ? 2 : ndim;
356  if (m_halfMode && i == 2)
357  {
358  lim = 0;
359  }
360  for (int j = 0; j < lim; ++j)
361  {
362  Vmath::Vvtvp(nPointsTot, m_gradBase[i * nBaseDerivs + j], 1,
363  velocity[j], 1, outarray[i], 1, outarray[i], 1);
364  }
365 
366  if (m_multipleModes)
367  {
368  fields[i]->HomogeneousFwdTrans(outarray[i], outarray[i]);
369  }
370  Vmath::Neg(nPointsTot, outarray[i], 1);
371  }
372 }
373 
375  const Array<OneD, Array<OneD, NekDouble>> &inarray,
377 {
378  if (m_session->GetSolverInfo("EqType") == "UnsteadyNavierStokes")
379  {
380  // The SFD method is only applied to the velocity variables in
381  // incompressible
382  ASSERTL1(inarray.size() == (m_baseflow.size() - 1),
383  "Number of base flow variables does not match what is "
384  "expected.");
385  }
386  else
387  {
388  ASSERTL1(
389  inarray.size() == (m_baseflow.size()),
390  "Number of base flow variables does not match what is expected.");
391  }
392 
393  int npts = inarray[0].size();
394 
395  for (int i = 0; i < inarray.size(); ++i)
396  {
397  ASSERTL1(npts == m_baseflow[i].size(),
398  "Size of base flow array does not match expected.");
399  Vmath::Vcopy(npts, inarray[i], 1, m_baseflow[i], 1);
400  UpdateGradBase(i, fields[i]);
401  }
402 }
403 
404 /**
405  * Import field from infile and load into \a m_fields. This routine will
406  * also perform a \a BwdTrans to ensure data is in both the physical and
407  * coefficient storage.
408  * @param pInFile Filename to read.
409  * @param pFields Array of expansion lists
410  */
412  std::string pInfile, Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
413  int pSlice)
414 {
415  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
416  std::vector<std::vector<NekDouble>> FieldData;
417 
418  int nqtot = m_baseflow[0].size();
419  Array<OneD, NekDouble> tmp_coeff(pFields[0]->GetNcoeffs(), 0.0);
420 
421  int numexp = pFields[0]->GetExpSize();
422  Array<OneD, int> ElementGIDs(numexp);
423 
424  // Define list of global element ids
425  for (int i = 0; i < numexp; ++i)
426  {
427  ElementGIDs[i] = pFields[0]->GetExp(i)->GetGeom()->GetGlobalID();
428  }
429 
430  // Get Homogeneous
433  fld->Import(pInfile, FieldDef, FieldData,
435 
436  int nSessionVar = m_session->GetVariables().size();
437  int nSessionConvVar = nSessionVar - 1;
438  int nFileVar = FieldDef[0]->m_fields.size();
439  int nFileConvVar = nFileVar - 1; // Ignore pressure
440  if (m_halfMode)
441  {
442  ASSERTL0(nFileVar == 3, "For half mode, expect 2D2C base flow.");
443  nFileConvVar = 2;
444  }
445 
446  for (int j = 0; j < nFileConvVar; ++j)
447  {
448  for (int i = 0; i < FieldDef.size(); ++i)
449  {
450  bool flag = FieldDef[i]->m_fields[j] == m_session->GetVariable(j);
451 
452  ASSERTL0(flag, (std::string("Order of ") + pInfile +
453  std::string(" data and that defined in "
454  "the session file differs"))
455  .c_str());
456 
457  pFields[j]->ExtractDataToCoeffs(
458  FieldDef[i], FieldData[i], FieldDef[i]->m_fields[j], tmp_coeff);
459  }
460 
461  if (m_singleMode || m_halfMode)
462  {
463  pFields[j]->GetPlane(0)->BwdTrans(tmp_coeff, m_baseflow[j]);
464 
465  if (m_singleMode)
466  {
467  // copy the bwd trans into the second plane for single
468  // Mode Analysis
469  int ncplane = (pFields[0]->GetNpoints()) / m_npointsZ;
470  Vmath::Vcopy(ncplane, &m_baseflow[j][0], 1,
471  &m_baseflow[j][ncplane], 1);
472  }
473  }
474  else // fully 3D base flow - put in physical space.
475  {
476  bool oldwavespace = pFields[j]->GetWaveSpace();
477  pFields[j]->SetWaveSpace(false);
478  pFields[j]->BwdTrans(tmp_coeff, m_baseflow[j]);
479  pFields[j]->SetWaveSpace(oldwavespace);
480  }
481  }
482 
483  // Zero unused fields (e.g. w in a 2D2C base flow).
484  for (int j = nFileConvVar; j < nSessionConvVar; ++j)
485  {
486  Vmath::Fill(nqtot, 0.0, m_baseflow[j], 1);
487  }
488 
489  // If time-periodic, put loaded data into the slice storage.
490  if (m_slices > 1)
491  {
492  for (int i = 0; i < nSessionConvVar; ++i)
493  {
494  Vmath::Vcopy(nqtot, &m_baseflow[i][0], 1,
495  &m_interp[i][pSlice * nqtot], 1);
496  }
497  }
498 }
499 
501  const NekDouble m_slices, const Array<OneD, const NekDouble> &inarray,
502  Array<OneD, NekDouble> &outarray, const NekDouble m_time,
503  const NekDouble m_period)
504 {
505  int npoints = m_baseflow[0].size();
506  if (m_isperiodic)
507  {
508  NekDouble BetaT = 2 * M_PI * fmod(m_time, m_period) / m_period;
509  NekDouble phase;
510  Array<OneD, NekDouble> auxiliary(npoints);
511 
512  Vmath::Vcopy(npoints, &inarray[0], 1, &outarray[0], 1);
513  Vmath::Svtvp(npoints, cos(0.5 * m_slices * BetaT), &inarray[npoints], 1,
514  &outarray[0], 1, &outarray[0], 1);
515 
516  for (int i = 2; i < m_slices; i += 2)
517  {
518  phase = (i >> 1) * BetaT;
519 
520  Vmath::Svtvp(npoints, cos(phase), &inarray[i * npoints], 1,
521  &outarray[0], 1, &outarray[0], 1);
522  Vmath::Svtvp(npoints, -sin(phase), &inarray[(i + 1) * npoints], 1,
523  &outarray[0], 1, &outarray[0], 1);
524  }
525  }
526  else
527  {
528  NekDouble x = m_time;
529  x = x / m_period * (m_slices - 1);
530  int ix = x;
531  if (ix < 0)
532  {
533  ix = 0;
534  }
535  if (ix > m_slices - 2)
536  {
537  ix = m_slices - 2;
538  }
539  int padleft = m_interporder / 2 - 1;
540  if (padleft > ix)
541  {
542  padleft = ix;
543  }
544  int padright = m_interporder - 1 - padleft;
545  if (padright > m_slices - 1 - ix)
546  {
547  padright = m_slices - 1 - ix;
548  }
549  padleft = m_interporder - 1 - padright;
551  for (int i = 0; i < m_interporder; ++i)
552  {
553  for (int j = 0; j < m_interporder; ++j)
554  {
555  if (i != j)
556  {
557  coeff[i] *= (x - ix + padleft - (NekDouble)j) /
558  ((NekDouble)i - (NekDouble)j);
559  }
560  }
561  }
562  Vmath::Zero(npoints, &outarray[0], 1);
563  for (int i = ix - padleft; i < ix + padright + 1; ++i)
564  {
565  Vmath::Svtvp(npoints, coeff[i - ix + padleft],
566  &inarray[i * npoints], 1, &outarray[0], 1,
567  &outarray[0], 1);
568  }
569  }
570 }
571 
573  const int var, const MultiRegions::ExpListSharedPtr &field)
574 {
575  int nBaseDerivs = (m_halfMode || m_singleMode) ? 2 : m_spacedim;
576  switch (m_spacedim)
577  {
578  case 1: // 1D
579  {
580  field->PhysDeriv(m_baseflow[var],
581  m_gradBase[var * nBaseDerivs + 0]);
582  }
583  break;
584  case 2: // 2D
585  {
586  field->PhysDeriv(m_baseflow[var], m_gradBase[var * nBaseDerivs + 0],
587  m_gradBase[var * nBaseDerivs + 1]);
588  }
589  break;
590  case 3:
591  {
592  if (m_halfMode) // can assume W = 0 and d/dz == 0
593  {
594  if (var < 2)
595  {
596  field->PhysDeriv(m_baseflow[var],
597  m_gradBase[var * nBaseDerivs + 0],
598  m_gradBase[var * nBaseDerivs + 1]);
599  }
600  }
601  else if (m_singleMode) // single mode where d/dz = 0
602  {
603  field->PhysDeriv(m_baseflow[var],
604  m_gradBase[var * nBaseDerivs + 0],
605  m_gradBase[var * nBaseDerivs + 1]);
606  }
607  else
608  {
609  // Differentiate base flow in physical space
610  bool oldwavespace = field->GetWaveSpace();
611  field->SetWaveSpace(false);
612  field->PhysDeriv(m_baseflow[var],
613  m_gradBase[var * nBaseDerivs + 0],
614  m_gradBase[var * nBaseDerivs + 1],
615  m_gradBase[var * nBaseDerivs + 2]);
616  field->SetWaveSpace(oldwavespace);
617  }
618  }
619  break;
620  }
621 }
622 
624  FloquetMatType mattype, bool UseContCoeffs) const
625 {
626  DNekMatSharedPtr loc_mat;
627  DNekBlkMatSharedPtr BlkMatrix;
628  int n_exp = 0;
629 
630  n_exp = m_baseflow[0].size(); // will operatore on m_phys
631 
632  Array<OneD, unsigned int> nrows(n_exp);
633  Array<OneD, unsigned int> ncols(n_exp);
634 
635  nrows = Array<OneD, unsigned int>(n_exp, m_slices);
636  ncols = Array<OneD, unsigned int>(n_exp, m_slices);
637 
638  MatrixStorage blkmatStorage = eDIAGONAL;
639  BlkMatrix = MemoryManager<DNekBlkMat>::AllocateSharedPtr(nrows, ncols,
640  blkmatStorage);
641 
645  StdRegions::StdSegExp StdSeg(BK);
646 
648  StdSeg.DetShapeType(), StdSeg);
649 
650  loc_mat = StdSeg.GetStdMatrix(matkey);
651 
652  // set up array of block matrices.
653  for (int i = 0; i < n_exp; ++i)
654  {
655  BlkMatrix->SetBlock(i, i, loc_mat);
656  }
657 
658  return BlkMatrix;
659 }
660 
661 // Discrete Fourier Transform for Floquet analysis
663  const string file, Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
664  const NekDouble m_slices)
665 {
666  int ConvectedFields = m_baseflow.size() - 1;
667  int npoints = m_baseflow[0].size();
668  m_interp = Array<OneD, Array<OneD, NekDouble>>(ConvectedFields);
669 
670  for (int i = 0; i < ConvectedFields; ++i)
671  {
672  m_interp[i] = Array<OneD, NekDouble>(npoints * m_slices, 0.0);
673  }
674 
675  // Import the slides into auxiliary vector
676  // The base flow should be stored in the form "filename_%d.ext"
677  // A subdirectory can also be included, such as "dir/filename_%d.ext"
678  size_t found = file.find("%d");
679  ASSERTL0(found != string::npos &&
680  file.find("%d", found + 1) == string::npos,
681  "Since N_slices is specified, the filename provided for function "
682  "'BaseFlow' must include exactly one instance of the format "
683  "specifier '%d', to index the time-slices.");
684  char *buffer = new char[file.length() + 8];
685  int nstart = m_start;
686  for (int i = nstart; i < nstart + m_slices * m_skip; i += m_skip)
687  {
688  sprintf(buffer, file.c_str(), i);
689  ImportFldBase(buffer, pFields, (i - nstart) / m_skip);
690  if (m_session->GetComm()->GetRank() == 0)
691  {
692  cout << "read base flow file " << buffer << endl;
693  }
694  }
695  delete[] buffer;
696  if (!m_isperiodic)
697  {
698  return;
699  }
700 
701  // Discrete Fourier Transform of the fields
702  for (int k = 0; k < ConvectedFields; ++k)
703  {
704 #ifdef NEKTAR_USING_FFTW
705 
706  // Discrete Fourier Transform using FFTW
707  Array<OneD, NekDouble> fft_in(npoints * m_slices);
708  Array<OneD, NekDouble> fft_out(npoints * m_slices);
709 
712 
713  // Shuffle the data
714  for (int j = 0; j < m_slices; ++j)
715  {
716  Vmath::Vcopy(npoints, &m_interp[k][j * npoints], 1, &(fft_in[j]),
717  m_slices);
718  }
719 
721  m_slices);
722 
723  // FFT Transform
724  for (int i = 0; i < npoints; i++)
725  {
726  m_FFT->FFTFwdTrans(m_tmpIN = fft_in + i * m_slices,
727  m_tmpOUT = fft_out + i * m_slices);
728  }
729 
730  // Reshuffle data
731  for (int s = 0; s < m_slices; ++s)
732  {
733  Vmath::Vcopy(npoints, &fft_out[s], m_slices,
734  &m_interp[k][s * npoints], 1);
735  }
736 
737  Vmath::Zero(fft_in.size(), &fft_in[0], 1);
738  Vmath::Zero(fft_out.size(), &fft_out[0], 1);
739 #else
740  // Discrete Fourier Transform using MVM
741  DNekBlkMatSharedPtr blkmat;
743 
744  int nrows = blkmat->GetRows();
745  int ncols = blkmat->GetColumns();
746 
747  Array<OneD, NekDouble> sortedinarray(ncols);
748  Array<OneD, NekDouble> sortedoutarray(nrows);
749 
750  // Shuffle the data
751  for (int j = 0; j < m_slices; ++j)
752  {
753  Vmath::Vcopy(npoints, &m_interp[k][j * npoints], 1,
754  &(sortedinarray[j]), m_slices);
755  }
756 
757  // Create NekVectors from the given data arrays
758  NekVector<NekDouble> in(ncols, sortedinarray, eWrapper);
759  NekVector<NekDouble> out(nrows, sortedoutarray, eWrapper);
760 
761  // Perform matrix-vector multiply.
762  out = (*blkmat) * in;
763 
764  // Reshuffle data
765  for (int s = 0; s < m_slices; ++s)
766  {
767  Vmath::Vcopy(npoints, &sortedoutarray[s], m_slices,
768  &m_interp[k][s * npoints], 1);
769  }
770 
771  for (int r = 0; r < sortedinarray.size(); ++r)
772  {
773  sortedinarray[0] = 0;
774  sortedoutarray[0] = 0;
775  }
776 
777 #endif
778  }
779 }
780 
781 } // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:249
Describes the specification for a Basis.
Definition: Basis.h:50
static std::shared_ptr< FieldIO > CreateForFile(const LibUtilities::SessionReaderSharedPtr session, const std::string &filename)
Construct a FieldIO object for a given input filename.
Definition: FieldIO.cpp:227
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:198
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:144
Defines a specification for a set of points.
Definition: Points.h:59
SpatialDomains::BoundaryConditionsSharedPtr m_boundaryConditions
int m_npointsZ
number of points in Z direction (if homogeneous)
NekDouble m_LhomZ
physical length in Z direction (if homogeneous)
enum HomogeneousType m_HomogeneousType
NekDouble m_period
period length
int m_npointsX
number of points in X direction (if homogeneous)
int m_npointsY
number of points in Y direction (if homogeneous)
int m_HomoDirec
number of homogenous directions
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.
bool m_useFFT
flag to determine if use or not the FFT for transformations
LibUtilities::SessionReaderSharedPtr m_session
Array< OneD, NekDouble > m_tmpIN
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)
DNekBlkMatSharedPtr GetFloquetBlockMatrix(FloquetMatType mattype, bool UseContCoeffs=false) const
void DFT(const std::string file, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const NekDouble m_slices)
bool m_singleMode
flag to determine if use single mode or not
Array< OneD, Array< OneD, NekDouble > > m_baseflow
Storage for base flow.
Array< OneD, NekDouble > m_tmpOUT
Array< OneD, Array< OneD, NekDouble > > m_gradBase
void UpdateBase(const NekDouble m_slices, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const NekDouble m_time, const NekDouble m_period)
bool m_multipleModes
flag to determine if use multiple mode or not
LibUtilities::NektarFFTSharedPtr m_FFT
auxiliary variables
bool m_halfMode
flag to determine if use half mode or not
NekDouble m_LhomY
physical length in Y direction (if homogeneous)
void ImportFldBase(std::string pInfile, Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, int slice)
Import Base flow.
NekDouble m_LhomX
physical length in X direction (if homogeneous)
Array< OneD, Array< OneD, NekDouble > > m_interp
interpolation vector
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.
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
An abstract base class encapsulating the concept of advection of a vector field.
Definition: Advection.h:72
virtual SOLVER_UTILS_EXPORT void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
Initialises the advection object.
Definition: Advection.cpp:353
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:611
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:375
Class representing a segment element in reference space.
Definition: StdSegExp.h:54
array buffer
Definition: GsLib.hpp:83
NektarFFTFactory & GetNektarFFTFactory()
Definition: NektarFFT.cpp:69
std::shared_ptr< FieldIO > FieldIOSharedPtr
Definition: FieldIO.h:301
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< Equation > EquationSharedPtr
Definition: Equation.h:130
static FieldMetaDataMap NullFieldMetaDataMap
Definition: FieldIO.h:53
@ eFourierEvenlySpaced
1D Evenly-spaced points using Fourier Fit
Definition: PointsType.h:76
@ eFourier
Fourier Expansion .
Definition: BasisType.h:57
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Definition: Advection.cpp:47
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
std::shared_ptr< DNekBlkMat > DNekBlkMatSharedPtr
Definition: NekTypeDefs.hpp:77
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:75
double NekDouble
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:209
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:622
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:518
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:574
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:492
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:45
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255