Nektar++
ExpListHomogeneous1D.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File ExpListHomogeneous1D.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: An ExpList which is homogeneous in 1-direction
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #include <LibUtilities/Foundations/ManagerAccess.h> // for PointsManager, etc
36 #include <LocalRegions/Expansion.h>
39 #include <StdRegions/StdPointExp.h>
40 #include <StdRegions/StdSegExp.h>
41 
42 using namespace std;
43 
44 namespace Nektar
45 {
46 namespace MultiRegions
47 {
48 // Forward declaration for typedefs
49 ExpListHomogeneous1D::ExpListHomogeneous1D(const ExpansionType type)
50  : ExpList(type), m_homogeneousBasis(LibUtilities::NullBasisSharedPtr),
51  m_lhom(1), m_homogeneous1DBlockMat(
52  MemoryManager<Homo1DBlockMatrixMap>::AllocateSharedPtr())
53 {
54 }
55 
57  const ExpansionType type,
59  const LibUtilities::BasisKey &HomoBasis, const NekDouble lhom,
60  const bool useFFT, const bool dealiasing)
61  : ExpList(type), m_useFFT(useFFT), m_lhom(lhom),
62  m_homogeneous1DBlockMat(
63  MemoryManager<Homo1DBlockMatrixMap>::AllocateSharedPtr()),
64  m_dealiasing(dealiasing)
65 {
66  m_session = pSession;
67  m_comm = pSession->GetComm();
68 
70  "Homogeneous Basis is a null basis");
71 
73 
74  m_StripZcomm = m_session->DefinesSolverInfo("HomoStrip")
75  ? m_comm->GetColumnComm()->GetColumnComm()
76  : m_comm->GetColumnComm();
77 
78  ASSERTL0(m_homogeneousBasis->GetNumPoints() % m_StripZcomm->GetSize() == 0,
79  "HomModesZ should be a multiple of npz.");
80 
81  if ((m_homogeneousBasis->GetBasisType() !=
83  (m_homogeneousBasis->GetBasisType() !=
85  {
86  ASSERTL0(
87  (m_homogeneousBasis->GetNumPoints() / m_StripZcomm->GetSize()) %
88  2 ==
89  0,
90  "HomModesZ/npz should be an even integer.");
91  }
92 
95  HomoBasis, m_comm, m_StripZcomm);
96 
98  m_homogeneousBasis->GetNumPoints() / m_StripZcomm->GetSize());
99 
100  if (m_useFFT)
101  {
103  "NekFFTW", m_homogeneousBasis->GetNumPoints());
104  }
105 
106  if (m_dealiasing)
107  {
108  if (m_useFFT)
109  {
110  int size = m_homogeneousBasis->GetNumPoints() +
111  m_homogeneousBasis->GetNumPoints() / 2;
112  m_padsize = size + (size % 2);
114  "NekFFTW", m_padsize);
115  }
116  else
117  {
118  ASSERTL0(false, "Dealiasing available just in combination "
119  "with FFTW");
120  }
121  }
122 }
123 
124 /**
125  * @param In ExpListHomogeneous1D object to copy.
126  */
128  : ExpList(In, false), m_transposition(In.m_transposition),
129  m_StripZcomm(In.m_StripZcomm), m_useFFT(In.m_useFFT), m_FFT(In.m_FFT),
130  m_tmpIN(In.m_tmpIN), m_tmpOUT(In.m_tmpOUT),
131  m_homogeneousBasis(In.m_homogeneousBasis), m_lhom(In.m_lhom),
132  m_homogeneous1DBlockMat(In.m_homogeneous1DBlockMat),
133  m_dealiasing(In.m_dealiasing), m_padsize(In.m_padsize)
134 {
136 }
137 
139  const ExpListHomogeneous1D &In, const std::vector<unsigned int> &eIDs,
140  const Collections::ImplementationType ImpType)
141  : ExpList(In, eIDs, false, ImpType), m_transposition(In.m_transposition),
142  m_useFFT(In.m_useFFT), m_FFT(In.m_FFT), m_tmpIN(In.m_tmpIN),
143  m_tmpOUT(In.m_tmpOUT), m_homogeneousBasis(In.m_homogeneousBasis),
144  m_lhom(In.m_lhom),
145  m_homogeneous1DBlockMat(
146  MemoryManager<Homo1DBlockMatrixMap>::AllocateSharedPtr()),
147  m_dealiasing(In.m_dealiasing), m_padsize(In.m_padsize)
148 {
150 }
151 
152 /**
153  * Destructor
154  */
156 {
157 }
158 
160  const Array<OneD, const NekDouble> &inarray,
161  Array<OneD, NekDouble> &outarray, bool Shuff, bool UnShuff)
162 {
163  // Forwards trans
164  Homogeneous1DTrans(inarray, outarray, true, Shuff, UnShuff);
165 }
166 
168  const Array<OneD, const NekDouble> &inarray,
169  Array<OneD, NekDouble> &outarray, bool Shuff, bool UnShuff)
170 {
171  // Backwards trans
172  Homogeneous1DTrans(inarray, outarray, false, Shuff, UnShuff);
173 }
174 
175 /**
176  * Dealiasing routine
177  *
178  * @param inarray1 First term of the product
179  * @param inarray2 Second term of the product
180  * @param outarray Dealiased product
181  */
183  const Array<OneD, NekDouble> &inarray1,
184  const Array<OneD, NekDouble> &inarray2, Array<OneD, NekDouble> &outarray)
185 {
186  int num_dofs = inarray1.size();
187  int N = m_homogeneousBasis->GetNumPoints();
188 
189  Array<OneD, NekDouble> V1(num_dofs);
190  Array<OneD, NekDouble> V2(num_dofs);
191  Array<OneD, NekDouble> V1V2(num_dofs);
192 
193  if (m_WaveSpace)
194  {
195  V1 = inarray1;
196  V2 = inarray2;
197  }
198  else
199  {
200  HomogeneousFwdTrans(inarray1, V1);
201  HomogeneousFwdTrans(inarray2, V2);
202  }
203 
204  int num_points_per_plane = num_dofs / m_planes.size();
205  int num_proc;
206  if (!m_session->DefinesSolverInfo("HomoStrip"))
207  {
208  num_proc = m_comm->GetColumnComm()->GetSize();
209  }
210  else
211  {
212  num_proc = m_StripZcomm->GetSize();
213  }
214  int num_dfts_per_proc =
215  num_points_per_plane / num_proc + (num_points_per_plane % num_proc > 0);
216 
217  Array<OneD, NekDouble> ShufV1(num_dfts_per_proc * N, 0.0);
218  Array<OneD, NekDouble> ShufV2(num_dfts_per_proc * N, 0.0);
219  Array<OneD, NekDouble> ShufV1V2(num_dfts_per_proc * N, 0.0);
220 
221  Array<OneD, NekDouble> ShufV1_PAD_coef(m_padsize, 0.0);
222  Array<OneD, NekDouble> ShufV2_PAD_coef(m_padsize, 0.0);
223  Array<OneD, NekDouble> ShufV1_PAD_phys(m_padsize, 0.0);
224  Array<OneD, NekDouble> ShufV2_PAD_phys(m_padsize, 0.0);
225 
226  Array<OneD, NekDouble> ShufV1V2_PAD_coef(m_padsize, 0.0);
227  Array<OneD, NekDouble> ShufV1V2_PAD_phys(m_padsize, 0.0);
228 
229  m_transposition->Transpose(V1, ShufV1, false, LibUtilities::eXYtoZ);
230  m_transposition->Transpose(V2, ShufV2, false, LibUtilities::eXYtoZ);
231 
232  // Looping on the pencils
233  for (int i = 0; i < num_dfts_per_proc; i++)
234  {
235  // Copying the i-th pencil pf lenght N into a bigger
236  // pencil of lenght 2N We are in Fourier space
237  Vmath::Vcopy(N, &(ShufV1[i * N]), 1, &(ShufV1_PAD_coef[0]), 1);
238  Vmath::Vcopy(N, &(ShufV2[i * N]), 1, &(ShufV2_PAD_coef[0]), 1);
239 
240  // Moving to physical space using the padded system
241  m_FFT_deal->FFTBwdTrans(ShufV1_PAD_coef, ShufV1_PAD_phys);
242  m_FFT_deal->FFTBwdTrans(ShufV2_PAD_coef, ShufV2_PAD_phys);
243 
244  // Perfroming the vectors multiplication in physical space on
245  // the padded system
246  Vmath::Vmul(m_padsize, ShufV1_PAD_phys, 1, ShufV2_PAD_phys, 1,
247  ShufV1V2_PAD_phys, 1);
248 
249  // Moving back the result (V1*V2)_phys in Fourier space, padded
250  // system
251  m_FFT_deal->FFTFwdTrans(ShufV1V2_PAD_phys, ShufV1V2_PAD_coef);
252 
253  // Copying the first half of the padded pencil in the full
254  // vector (Fourier space)
255  Vmath::Vcopy(N, &(ShufV1V2_PAD_coef[0]), 1, &(ShufV1V2[i * N]), 1);
256  }
257 
258  // Moving the results to the output
259  if (m_WaveSpace)
260  {
261  m_transposition->Transpose(ShufV1V2, outarray, false,
263  }
264  else
265  {
266  m_transposition->Transpose(ShufV1V2, V1V2, false, LibUtilities::eZtoXY);
267  HomogeneousBwdTrans(V1V2, outarray);
268  }
269 }
270 
271 /**
272  * Dealiasing routine for dot product
273  *
274  * @param inarray1 First term of the product with dimension ndim
275  * (e.g. U)
276  * @param inarray2 Second term of the product with dimension ndim*nvec
277  * (e.g. grad U)
278  * @param outarray Dealiased product with dimension nvec
279  */
281  const Array<OneD, Array<OneD, NekDouble>> &inarray1,
282  const Array<OneD, Array<OneD, NekDouble>> &inarray2,
283  Array<OneD, Array<OneD, NekDouble>> &outarray)
284 {
285  int ndim = inarray1.size();
286  ASSERTL1(inarray2.size() % ndim == 0,
287  "Wrong dimensions for DealiasedDotProd.");
288  int nvec = inarray2.size() / ndim;
289 
290  int num_dofs = inarray1[0].size();
291  int N = m_homogeneousBasis->GetNumPoints();
292 
293  int num_points_per_plane = num_dofs / m_planes.size();
294  int num_proc;
295  if (!m_session->DefinesSolverInfo("HomoStrip"))
296  {
297  num_proc = m_comm->GetColumnComm()->GetSize();
298  }
299  else
300  {
301  num_proc = m_StripZcomm->GetSize();
302  }
303  int num_dfts_per_proc =
304  num_points_per_plane / num_proc + (num_points_per_plane % num_proc > 0);
305 
306  // Get inputs in Fourier space
308  Array<OneD, Array<OneD, NekDouble>> V2(ndim * nvec);
309  if (m_WaveSpace)
310  {
311  for (int i = 0; i < ndim; i++)
312  {
313  V1[i] = inarray1[i];
314  }
315  for (int i = 0; i < ndim * nvec; i++)
316  {
317  V2[i] = inarray2[i];
318  }
319  }
320  else
321  {
322  for (int i = 0; i < ndim; i++)
323  {
324  V1[i] = Array<OneD, NekDouble>(num_dofs);
325  HomogeneousFwdTrans(inarray1[i], V1[i]);
326  }
327  for (int i = 0; i < ndim * nvec; i++)
328  {
329  V2[i] = Array<OneD, NekDouble>(num_dofs);
330  HomogeneousFwdTrans(inarray2[i], V2[i]);
331  }
332  }
333 
334  // Allocate variables for ffts
336  Array<OneD, NekDouble> ShufV1_PAD_coef(m_padsize, 0.0);
337  Array<OneD, Array<OneD, NekDouble>> ShufV1_PAD_phys(ndim);
338  for (int i = 0; i < ndim; i++)
339  {
340  ShufV1[i] = Array<OneD, NekDouble>(num_dfts_per_proc * N, 0.0);
341  ShufV1_PAD_phys[i] = Array<OneD, NekDouble>(m_padsize, 0.0);
342  }
343 
344  Array<OneD, Array<OneD, NekDouble>> ShufV2(ndim * nvec);
345  Array<OneD, NekDouble> ShufV2_PAD_coef(m_padsize, 0.0);
346  Array<OneD, Array<OneD, NekDouble>> ShufV2_PAD_phys(ndim * nvec);
347  for (int i = 0; i < ndim * nvec; i++)
348  {
349  ShufV2[i] = Array<OneD, NekDouble>(num_dfts_per_proc * N, 0.0);
350  ShufV2_PAD_phys[i] = Array<OneD, NekDouble>(m_padsize, 0.0);
351  }
352 
353  Array<OneD, Array<OneD, NekDouble>> ShufV1V2(nvec);
354  Array<OneD, NekDouble> ShufV1V2_PAD_coef(m_padsize, 0.0);
355  Array<OneD, NekDouble> ShufV1V2_PAD_phys(m_padsize, 0.0);
356  for (int i = 0; i < nvec; i++)
357  {
358  ShufV1V2[i] = Array<OneD, NekDouble>(num_dfts_per_proc * N, 0.0);
359  }
360 
361  // Transpose inputs
362  for (int i = 0; i < ndim; i++)
363  {
364  m_transposition->Transpose(V1[i], ShufV1[i], false,
366  }
367  for (int i = 0; i < ndim * nvec; i++)
368  {
369  m_transposition->Transpose(V2[i], ShufV2[i], false,
371  }
372 
373  // Looping on the pencils
374  for (int i = 0; i < num_dfts_per_proc; i++)
375  {
376  for (int j = 0; j < ndim; j++)
377  {
378  // Copying the i-th pencil pf lenght N into a bigger
379  // pencil of lenght 1.5N We are in Fourier space
380  Vmath::Vcopy(N, &(ShufV1[j][i * N]), 1, &(ShufV1_PAD_coef[0]), 1);
381  // Moving to physical space using the padded system
382  m_FFT_deal->FFTBwdTrans(ShufV1_PAD_coef, ShufV1_PAD_phys[j]);
383  }
384  for (int j = 0; j < ndim * nvec; j++)
385  {
386  Vmath::Vcopy(N, &(ShufV2[j][i * N]), 1, &(ShufV2_PAD_coef[0]), 1);
387  m_FFT_deal->FFTBwdTrans(ShufV2_PAD_coef, ShufV2_PAD_phys[j]);
388  }
389 
390  // Performing the vectors multiplication in physical space on
391  // the padded system
392  for (int j = 0; j < nvec; j++)
393  {
394  Vmath::Zero(m_padsize, ShufV1V2_PAD_phys, 1);
395  for (int k = 0; k < ndim; k++)
396  {
397  Vmath::Vvtvp(m_padsize, ShufV1_PAD_phys[k], 1,
398  ShufV2_PAD_phys[j * ndim + k], 1,
399  ShufV1V2_PAD_phys, 1, ShufV1V2_PAD_phys, 1);
400  }
401  // Moving back the result (V1*V2)_phys in Fourier space,
402  // padded system
403  m_FFT_deal->FFTFwdTrans(ShufV1V2_PAD_phys, ShufV1V2_PAD_coef);
404  // Copying the first half of the padded pencil in the full
405  // vector (Fourier space)
406  Vmath::Vcopy(N, &(ShufV1V2_PAD_coef[0]), 1, &(ShufV1V2[j][i * N]),
407  1);
408  }
409  }
410 
411  // Moving the results to the output
412  if (m_WaveSpace)
413  {
414  for (int j = 0; j < nvec; j++)
415  {
416  m_transposition->Transpose(ShufV1V2[j], outarray[j], false,
418  }
419  }
420  else
421  {
422  Array<OneD, NekDouble> V1V2(num_dofs);
423  for (int j = 0; j < nvec; j++)
424  {
425  m_transposition->Transpose(ShufV1V2[j], V1V2, false,
427  HomogeneousBwdTrans(V1V2, outarray[j]);
428  }
429  }
430 }
431 
432 /**
433  * Forward transform
434  */
436  const Array<OneD, const NekDouble> &inarray,
437  Array<OneD, NekDouble> &outarray)
438 {
439  int cnt = 0, cnt1 = 0;
440  Array<OneD, NekDouble> tmparray;
441 
442  for (int n = 0; n < m_planes.size(); ++n)
443  {
444  m_planes[n]->FwdTrans(inarray + cnt, tmparray = outarray + cnt1);
445  cnt += m_planes[n]->GetTotPoints();
446 
447  cnt1 += m_planes[n]->GetNcoeffs(); // need to skip ncoeffs
448  }
449  if (!m_WaveSpace)
450  {
451  HomogeneousFwdTrans(outarray, outarray);
452  }
453 }
454 
455 /**
456  * Forward transform element by element
457  */
459  const Array<OneD, const NekDouble> &inarray,
460  Array<OneD, NekDouble> &outarray)
461 {
462  int cnt = 0, cnt1 = 0;
463  Array<OneD, NekDouble> tmparray;
464 
465  // spectral element FwdTrans plane by plane
466  for (int n = 0; n < m_planes.size(); ++n)
467  {
468  m_planes[n]->FwdTransLocalElmt(inarray + cnt,
469  tmparray = outarray + cnt1);
470  cnt += m_planes[n]->GetTotPoints();
471  cnt1 += m_planes[n]->GetNcoeffs();
472  }
473  if (!m_WaveSpace)
474  {
475  HomogeneousFwdTrans(outarray, outarray);
476  }
477 }
478 
479 /**
480  * Forward transform element by element with boundaries constrained
481  */
483  const Array<OneD, const NekDouble> &inarray,
484  Array<OneD, NekDouble> &outarray)
485 {
486  int cnt = 0, cnt1 = 0;
487  Array<OneD, NekDouble> tmparray;
488 
489  // spectral element FwdTrans plane by plane
490  for (int n = 0; n < m_planes.size(); ++n)
491  {
492  m_planes[n]->FwdTransBndConstrained(inarray + cnt,
493  tmparray = outarray + cnt1);
494  cnt += m_planes[n]->GetTotPoints();
495  cnt1 += m_planes[n]->GetNcoeffs();
496  }
497  if (!m_WaveSpace)
498  {
499  HomogeneousFwdTrans(outarray, outarray);
500  }
501 }
502 
503 /**
504  * Backward transform
505  */
507  const Array<OneD, const NekDouble> &inarray,
508  Array<OneD, NekDouble> &outarray)
509 {
510  int cnt = 0, cnt1 = 0;
511  Array<OneD, NekDouble> tmparray;
512 
513  for (int n = 0; n < m_planes.size(); ++n)
514  {
515  m_planes[n]->BwdTrans(inarray + cnt, tmparray = outarray + cnt1);
516  cnt += m_planes[n]->GetNcoeffs();
517  cnt1 += m_planes[n]->GetTotPoints();
518  }
519  if (!m_WaveSpace)
520  {
521  HomogeneousBwdTrans(outarray, outarray);
522  }
523 }
524 
525 /**
526  * Inner product
527  */
529  const Array<OneD, const NekDouble> &inarray,
530  Array<OneD, NekDouble> &outarray)
531 {
532  int cnt = 0, cnt1 = 0;
533  Array<OneD, NekDouble> tmparray, tmpIn;
534 
535  if (m_WaveSpace)
536  {
537  tmpIn = inarray;
538  }
539  else
540  {
541  tmpIn = Array<OneD, NekDouble>(inarray.size(), 0.0);
542  HomogeneousFwdTrans(inarray, tmpIn);
543  }
544 
545  for (int n = 0; n < m_planes.size(); ++n)
546  {
547  m_planes[n]->IProductWRTBase(tmpIn + cnt, tmparray = outarray + cnt1);
548 
549  cnt1 += m_planes[n]->GetNcoeffs();
550  cnt += m_planes[n]->GetTotPoints();
551  }
552 }
553 
554 /**
555  * Homogeneous transform Bwd/Fwd (MVM and FFT)
556  */
558  const Array<OneD, const NekDouble> &inarray,
559  Array<OneD, NekDouble> &outarray, bool IsForwards, bool Shuff, bool UnShuff)
560 {
561  int num_dofs;
562 
563  if (IsForwards)
564  {
565  num_dofs = inarray.size();
566  }
567  else
568  {
569  num_dofs = outarray.size();
570  }
571 
572  if (m_useFFT)
573  {
574  int num_points_per_plane = num_dofs / m_planes.size();
575  int num_dfts_per_proc;
576  if (!m_session->DefinesSolverInfo("HomoStrip"))
577  {
578  int nP = m_comm->GetColumnComm()->GetSize();
579  num_dfts_per_proc =
580  num_points_per_plane / nP + (num_points_per_plane % nP > 0);
581  }
582  else
583  {
584  int nP = m_StripZcomm->GetSize();
585  num_dfts_per_proc =
586  num_points_per_plane / nP + (num_points_per_plane % nP > 0);
587  }
588 
589  Array<OneD, NekDouble> fft_in(
590  num_dfts_per_proc * m_homogeneousBasis->GetNumPoints(), 0.0);
591  Array<OneD, NekDouble> fft_out(
592  num_dfts_per_proc * m_homogeneousBasis->GetNumPoints(), 0.0);
593 
594  if (Shuff)
595  {
596  m_transposition->Transpose(inarray, fft_in, false,
598  }
599  else
600  {
601  Vmath::Vcopy(num_dfts_per_proc * m_homogeneousBasis->GetNumPoints(),
602  inarray, 1, fft_in, 1);
603  }
604 
605  if (IsForwards)
606  {
607  for (int i = 0; i < num_dfts_per_proc; i++)
608  {
609  m_FFT->FFTFwdTrans(
610  m_tmpIN = fft_in + i * m_homogeneousBasis->GetNumPoints(),
611  m_tmpOUT =
612  fft_out + i * m_homogeneousBasis->GetNumPoints());
613  }
614  }
615  else
616  {
617  for (int i = 0; i < num_dfts_per_proc; i++)
618  {
619  m_FFT->FFTBwdTrans(
620  m_tmpIN = fft_in + i * m_homogeneousBasis->GetNumPoints(),
621  m_tmpOUT =
622  fft_out + i * m_homogeneousBasis->GetNumPoints());
623  }
624  }
625 
626  if (UnShuff)
627  {
628  m_transposition->Transpose(fft_out, outarray, false,
630  }
631  else
632  {
633  Vmath::Vcopy(num_dfts_per_proc * m_homogeneousBasis->GetNumPoints(),
634  fft_out, 1, outarray, 1);
635  }
636  }
637  else
638  {
639  DNekBlkMatSharedPtr blkmat;
640 
641  if (num_dofs == m_npoints) // transform phys space
642  {
643  if (IsForwards)
644  {
646  }
647  else
648  {
650  }
651  }
652  else
653  {
654  if (IsForwards)
655  {
657  }
658  else
659  {
661  }
662  }
663 
664  int nrows = blkmat->GetRows();
665  int ncols = blkmat->GetColumns();
666 
667  Array<OneD, NekDouble> sortedinarray(ncols, 0.0);
668  Array<OneD, NekDouble> sortedoutarray(nrows, 0.0);
669 
670  if (Shuff)
671  {
672  m_transposition->Transpose(inarray, sortedinarray, !IsForwards,
674  }
675  else
676  {
677  Vmath::Vcopy(ncols, inarray, 1, sortedinarray, 1);
678  }
679 
680  // Create NekVectors from the given data arrays
681  NekVector<NekDouble> in(ncols, sortedinarray, eWrapper);
682  NekVector<NekDouble> out(nrows, sortedoutarray, eWrapper);
683 
684  // Perform matrix-vector multiply.
685  out = (*blkmat) * in;
686 
687  if (UnShuff)
688  {
689  m_transposition->Transpose(sortedoutarray, outarray, IsForwards,
691  }
692  else
693  {
694  Vmath::Vcopy(nrows, sortedoutarray, 1, outarray, 1);
695  }
696  }
697 }
698 
700  Homogeneous1DMatType mattype) const
701 {
702  auto matrixIter = m_homogeneous1DBlockMat->find(mattype);
703 
704  if (matrixIter == m_homogeneous1DBlockMat->end())
705  {
706  return ((*m_homogeneous1DBlockMat)[mattype] =
707  GenHomogeneous1DBlockMatrix(mattype));
708  }
709  else
710  {
711  return matrixIter->second;
712  }
713 }
714 
716  Homogeneous1DMatType mattype) const
717 {
718  DNekMatSharedPtr loc_mat;
719  DNekBlkMatSharedPtr BlkMatrix;
720  int n_exp = 0;
721  int num_trans_per_proc = 0;
722 
723  if ((mattype == eForwardsCoeffSpace1D) ||
724  (mattype == eBackwardsCoeffSpace1D)) // will operate on m_coeffs
725  {
726  n_exp = m_planes[0]->GetNcoeffs();
727  }
728  else
729  {
730  n_exp = m_planes[0]->GetTotPoints(); // will operatore on m_phys
731  }
732 
733  num_trans_per_proc = n_exp / m_comm->GetColumnComm()->GetSize() +
734  (n_exp % m_comm->GetColumnComm()->GetSize() > 0);
735 
736  Array<OneD, unsigned int> nrows(num_trans_per_proc);
737  Array<OneD, unsigned int> ncols(num_trans_per_proc);
738 
739  if ((mattype == eForwardsCoeffSpace1D) || (mattype == eForwardsPhysSpace1D))
740  {
741  nrows = Array<OneD, unsigned int>(num_trans_per_proc,
742  m_homogeneousBasis->GetNumModes());
743  ncols = Array<OneD, unsigned int>(num_trans_per_proc,
744  m_homogeneousBasis->GetNumPoints());
745  }
746  else
747  {
748  nrows = Array<OneD, unsigned int>(num_trans_per_proc,
749  m_homogeneousBasis->GetNumPoints());
750  ncols = Array<OneD, unsigned int>(num_trans_per_proc,
751  m_homogeneousBasis->GetNumModes());
752  }
753 
754  MatrixStorage blkmatStorage = eDIAGONAL;
755  BlkMatrix = MemoryManager<DNekBlkMat>::AllocateSharedPtr(nrows, ncols,
756  blkmatStorage);
757 
758  // Half Mode
759  if (m_homogeneousBasis->GetBasisType() ==
762  {
763  StdRegions::StdPointExp StdPoint(m_homogeneousBasis->GetBasisKey());
764 
765  if ((mattype == eForwardsCoeffSpace1D) ||
766  (mattype == eForwardsPhysSpace1D))
767  {
769  StdPoint.DetShapeType(), StdPoint);
770 
771  loc_mat = StdPoint.GetStdMatrix(matkey);
772  }
773  else
774  {
776  StdPoint.DetShapeType(), StdPoint);
777 
778  loc_mat = StdPoint.GetStdMatrix(matkey);
779  }
780  }
781  // other cases
782  else
783  {
784  StdRegions::StdSegExp StdSeg(m_homogeneousBasis->GetBasisKey());
785 
786  if ((mattype == eForwardsCoeffSpace1D) ||
787  (mattype == eForwardsPhysSpace1D))
788  {
790  StdSeg.DetShapeType(), StdSeg);
791 
792  loc_mat = StdSeg.GetStdMatrix(matkey);
793  }
794  else
795  {
797  StdSeg.DetShapeType(), StdSeg);
798 
799  loc_mat = StdSeg.GetStdMatrix(matkey);
800  }
801  }
802 
803  // set up array of block matrices.
804  for (int i = 0; i < num_trans_per_proc; ++i)
805  {
806  BlkMatrix->SetBlock(i, i, loc_mat);
807  }
808 
809  return BlkMatrix;
810 }
811 
812 std::vector<LibUtilities::FieldDefinitionsSharedPtr> ExpListHomogeneous1D::
814 {
815  std::vector<LibUtilities::FieldDefinitionsSharedPtr> returnval;
816 
817  // Set up Homogeneous length details.
819 
820  std::vector<NekDouble> HomoLen;
821  HomoLen.push_back(m_lhom);
822 
823  std::vector<unsigned int> StripsIDs;
824 
825  bool strips;
826  m_session->MatchSolverInfo("HomoStrip", "True", strips, false);
827  if (strips)
828  {
829  StripsIDs.push_back(m_transposition->GetStripID());
830  }
831 
832  std::vector<unsigned int> PlanesIDs;
833  int IDoffset = 0;
834 
835  // introduce a 2 plane offset for single mode case so can
836  // be post-processed or used in MultiMode expansion.
838  {
839  IDoffset = 2;
840  }
841 
842  for (int i = 0; i < m_planes.size(); i++)
843  {
844  PlanesIDs.push_back(m_transposition->GetPlaneID(i) + IDoffset);
845  }
846 
847  m_planes[0]->GeneralGetFieldDefinitions(returnval, 1, HomoBasis, HomoLen,
848  strips, StripsIDs, PlanesIDs);
849  return returnval;
850 }
851 
853  std::vector<LibUtilities::FieldDefinitionsSharedPtr> &fielddef)
854 {
855  // Set up Homogeneous length details.
857 
858  std::vector<NekDouble> HomoLen;
859  HomoLen.push_back(m_lhom);
860 
861  std::vector<unsigned int> StripsIDs;
862 
863  bool strips;
864  m_session->MatchSolverInfo("HomoStrip", "True", strips, false);
865  if (strips)
866  {
867  StripsIDs.push_back(m_transposition->GetStripID());
868  }
869 
870  std::vector<unsigned int> PlanesIDs;
871  int IDoffset = 0;
872 
874  {
875  IDoffset = 2;
876  }
877 
878  for (int i = 0; i < m_planes.size(); i++)
879  {
880  PlanesIDs.push_back(m_transposition->GetPlaneID(i) + IDoffset);
881  }
882 
883  // enforce NumHomoDir == 1 by direct call
884  m_planes[0]->GeneralGetFieldDefinitions(fielddef, 1, HomoBasis, HomoLen,
885  strips, StripsIDs, PlanesIDs);
886 }
887 
888 /** This routine appends the data from the expansion list into
889  the output format where each element is given by looping
890  over its Fourier modes where as data in the expandion is
891  stored with all consecutive elements and then the Fourier
892  modes
893  */
896  std::vector<NekDouble> &fielddata, Array<OneD, NekDouble> &coeffs)
897 {
898  int i, n;
899  int ncoeffs_per_plane = m_planes[0]->GetNcoeffs();
900 
901  // Determine mapping from element ids to location in
902  // expansion list
903  if (m_elmtToExpId.size() == 0)
904  {
905  for (i = 0; i < m_planes[0]->GetExpSize(); ++i)
906  {
907  m_elmtToExpId[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
908  }
909  }
910 
911  for (i = 0; i < fielddef->m_elementIDs.size(); ++i)
912  {
913  int eid = m_elmtToExpId[fielddef->m_elementIDs[i]];
914  int datalen = (*m_exp)[eid]->GetNcoeffs();
915 
916  for (n = 0; n < m_planes.size(); ++n)
917  {
918  fielddata.insert(
919  fielddata.end(),
920  &coeffs[m_coeff_offset[eid] + n * ncoeffs_per_plane],
921  &coeffs[m_coeff_offset[eid] + n * ncoeffs_per_plane] + datalen);
922  }
923  }
924 }
925 
928  std::vector<NekDouble> &fielddata)
929 {
930  v_AppendFieldData(fielddef, fielddata, m_coeffs);
931 }
932 
933 // Extract the data in fielddata into the m_coeff list
936  std::vector<NekDouble> &fielddata, std::string &field,
937  Array<OneD, NekDouble> &coeffs)
938 {
939  int i, n;
940  int offset = 0;
941  int nzmodes = 1;
942  int datalen = fielddata.size() / fielddef->m_fields.size();
943  std::vector<unsigned int> fieldDefHomoZids;
944 
945  // Find data location according to field definition
946  for (i = 0; i < fielddef->m_fields.size(); ++i)
947  {
948  if (fielddef->m_fields[i] == field)
949  {
950  break;
951  }
952  offset += datalen;
953  }
954 
955  if (i == fielddef->m_fields.size())
956  {
957  cout << "Field " << field << "not found in data file. " << endl;
958  }
959  else
960  {
961 
962  int modes_offset = 0;
963  int planes_offset = 0;
964  Array<OneD, NekDouble> coeff_tmp;
965 
966  // Build map of plane IDs lying on this processor and determine
967  // mapping from element ids to location in expansion list.
968  if (m_zIdToPlane.size() == 0)
969  {
970  for (i = 0; i < m_planes.size(); ++i)
971  {
972  m_zIdToPlane[m_transposition->GetPlaneID(i)] = i;
973  }
974 
975  for (i = 0; i < m_planes[0]->GetExpSize(); ++i)
976  {
977  m_elmtToExpId[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
978  }
979  }
980 
981  if (fielddef->m_numHomogeneousDir)
982  {
983  nzmodes = fielddef->m_homogeneousZIDs.size();
984  fieldDefHomoZids = fielddef->m_homogeneousZIDs;
985  }
986  else // input file is 2D and so set nzmodes to 1
987  {
988  nzmodes = 1;
989  fieldDefHomoZids.push_back(0);
990  }
991 
992  // calculate number of modes in the current partition
993  int ncoeffs_per_plane = m_planes[0]->GetNcoeffs();
994 
995  for (i = 0; i < fielddef->m_elementIDs.size(); ++i)
996  {
997  if (fielddef->m_uniOrder == true) // reset modes_offset to zero
998  {
999  modes_offset = 0;
1000  }
1001 
1003  fielddef->m_shapeType, fielddef->m_numModes, modes_offset);
1004 
1005  auto it = m_elmtToExpId.find(fielddef->m_elementIDs[i]);
1006 
1007  // ensure element is on this partition for parallel case.
1008  if (it == m_elmtToExpId.end())
1009  {
1010  // increase offset for correct FieldData access
1011  offset += datalen * nzmodes;
1012  modes_offset +=
1013  (*m_exp)[0]->GetNumBases() + fielddef->m_numHomogeneousDir;
1014  continue;
1015  }
1016 
1017  int eid = it->second;
1018  bool sameBasis = true;
1019  for (int j = 0; j < fielddef->m_basis.size() - 1; ++j)
1020  {
1021  if (fielddef->m_basis[j] != (*m_exp)[eid]->GetBasisType(j))
1022  {
1023  sameBasis = false;
1024  break;
1025  }
1026  }
1027 
1028  for (n = 0; n < nzmodes; ++n, offset += datalen)
1029  {
1030 
1031  it = m_zIdToPlane.find(fieldDefHomoZids[n]);
1032 
1033  // Check to make sure this mode number lies in this field.
1034  if (it == m_zIdToPlane.end())
1035  {
1036  continue;
1037  }
1038 
1039  planes_offset = it->second;
1040  if (datalen == (*m_exp)[eid]->GetNcoeffs() && sameBasis)
1041  {
1042  Vmath::Vcopy(datalen, &fielddata[offset], 1,
1043  &coeffs[m_coeff_offset[eid] +
1044  planes_offset * ncoeffs_per_plane],
1045  1);
1046  }
1047  else // unpack data to new order
1048  {
1049  (*m_exp)[eid]->ExtractDataToCoeffs(
1050  &fielddata[offset], fielddef->m_numModes, modes_offset,
1051  &coeffs[m_coeff_offset[eid] +
1052  planes_offset * ncoeffs_per_plane],
1053  fielddef->m_basis);
1054  }
1055  }
1056  modes_offset +=
1057  (*m_exp)[0]->GetNumBases() + fielddef->m_numHomogeneousDir;
1058  }
1059  }
1060 }
1061 
1062 // Extract the data in fielddata into the m_coeff list
1064  const std::shared_ptr<ExpList> &fromExpList,
1065  const Array<OneD, const NekDouble> &fromCoeffs,
1066  Array<OneD, NekDouble> &toCoeffs)
1067 {
1068  int i;
1069  int fromNcoeffs_per_plane = fromExpList->GetPlane(0)->GetNcoeffs();
1070  int toNcoeffs_per_plane = m_planes[0]->GetNcoeffs();
1071  Array<OneD, NekDouble> tocoeffs_tmp, fromcoeffs_tmp;
1072 
1073  for (i = 0; i < m_planes.size(); ++i)
1074  {
1075  m_planes[i]->ExtractCoeffsToCoeffs(
1076  fromExpList->GetPlane(i),
1077  fromcoeffs_tmp = fromCoeffs + fromNcoeffs_per_plane * i,
1078  tocoeffs_tmp = toCoeffs + toNcoeffs_per_plane * i);
1079  }
1080 }
1081 
1083  int expansion, std::string var)
1084 {
1085  // If there is only one plane (e.g. HalfMode), we write a 2D plane.
1086  if (m_planes.size() == 1)
1087  {
1088  m_planes[0]->WriteVtkPieceData(outfile, expansion, var);
1089  return;
1090  }
1091 
1092  int i;
1093  int nq = (*m_exp)[expansion]->GetTotPoints();
1094  int npoints_per_plane = m_planes[0]->GetTotPoints();
1095 
1096  // If we are using Fourier points, output extra plane to fill domain
1097  int outputExtraPlane = 0;
1098  Array<OneD, NekDouble> extraPlane;
1099  if (m_homogeneousBasis->GetBasisType() == LibUtilities::eFourier &&
1100  m_homogeneousBasis->GetPointsType() ==
1102  {
1103  outputExtraPlane = 1;
1104  // Get extra plane data
1105  if (m_StripZcomm->GetSize() == 1)
1106  {
1107  extraPlane = m_phys + m_phys_offset[expansion];
1108  }
1109  else
1110  {
1111  // Determine to and from rank for communication
1112  int size = m_StripZcomm->GetSize();
1113  int rank = m_StripZcomm->GetRank();
1114  int fromRank = (rank + 1) % size;
1115  int toRank = (rank == 0) ? size - 1 : rank - 1;
1116  // Communicate using SendRecv
1117  extraPlane = Array<OneD, NekDouble>(nq);
1118  Array<OneD, NekDouble> send(nq, m_phys + m_phys_offset[expansion]);
1119  m_StripZcomm->SendRecv(toRank, send, fromRank, extraPlane);
1120  }
1121  }
1122 
1123  // printing the fields of that zone
1124  outfile << " <DataArray type=\"Float64\" Name=\"" << var << "\">"
1125  << endl;
1126  outfile << " ";
1127  for (int n = 0; n < m_planes.size(); ++n)
1128  {
1129  const Array<OneD, NekDouble> phys =
1130  m_phys + m_phys_offset[expansion] + n * npoints_per_plane;
1131  for (i = 0; i < nq; ++i)
1132  {
1133  outfile << (fabs(phys[i]) < NekConstants::kNekZeroTol ? 0 : phys[i])
1134  << " ";
1135  }
1136  }
1137  if (outputExtraPlane)
1138  {
1139  for (i = 0; i < nq; ++i)
1140  {
1141  outfile << (fabs(extraPlane[i]) < NekConstants::kNekZeroTol
1142  ? 0
1143  : extraPlane[i])
1144  << " ";
1145  }
1146  }
1147  outfile << endl;
1148  outfile << " </DataArray>" << endl;
1149 }
1150 
1152  const NekDouble scale, const Array<OneD, NekDouble> &inarray,
1153  Array<OneD, NekDouble> &outarray)
1154 {
1155  int cnt, cnt1;
1156  Array<OneD, NekDouble> tmparray;
1157  cnt = m_planes[0]->GetTotPoints();
1158  cnt1 = m_planes[0]->Get1DScaledTotPoints(scale);
1159 
1160  ASSERTL1(m_planes.size() * cnt1 <= outarray.size(),
1161  "size of outarray does not match internal estimage");
1162 
1163  for (int i = 0; i < m_planes.size(); i++)
1164  {
1165 
1166  m_planes[i]->PhysInterp1DScaled(scale, inarray + i * cnt,
1167  tmparray = outarray + i * cnt1);
1168  }
1169 }
1170 
1172  const NekDouble scale, const Array<OneD, NekDouble> &inarray,
1173  Array<OneD, NekDouble> &outarray)
1174 {
1175  int cnt, cnt1;
1176  Array<OneD, NekDouble> tmparray;
1177  cnt = m_planes[0]->Get1DScaledTotPoints(scale);
1178  cnt1 = m_planes[0]->GetTotPoints();
1179  ASSERTL1(m_planes.size() * cnt <= inarray.size(),
1180  "size of outarray does not match internal estimage");
1181 
1182  for (int i = 0; i < m_planes.size(); i++)
1183  {
1184  m_planes[i]->PhysGalerkinProjection1DScaled(
1185  scale, inarray + i * cnt, tmparray = outarray + i * cnt1);
1186  }
1187 }
1189  const Array<OneD, const NekDouble> &inarray, Array<OneD, NekDouble> &out_d0,
1191 {
1192  int nT_pts = inarray.size(); // number of total points = n. of Fourier
1193  // points * n. of points per plane (nT_pts)
1194  int nP_pts =
1195  nT_pts / m_planes.size(); // number of points per plane = n of
1196  // Fourier transform required (nP_pts)
1197 
1198  Array<OneD, NekDouble> temparray(nT_pts);
1199  Array<OneD, NekDouble> outarray(nT_pts);
1203 
1204  for (int i = 0; i < m_planes.size(); i++)
1205  {
1206  m_planes[i]->PhysDeriv(inarray + i * nP_pts, tmp2 = out_d0 + i * nP_pts,
1207  tmp3 = out_d1 + i * nP_pts);
1208  }
1209 
1210  if (out_d2 != NullNekDouble1DArray)
1211  {
1212  if (m_homogeneousBasis->GetBasisType() == LibUtilities::eFourier ||
1213  m_homogeneousBasis->GetBasisType() ==
1215  m_homogeneousBasis->GetBasisType() ==
1217  m_homogeneousBasis->GetBasisType() ==
1219  {
1220  if (m_WaveSpace)
1221  {
1222  temparray = inarray;
1223  }
1224  else
1225  {
1226  HomogeneousFwdTrans(inarray, temparray);
1227  }
1228 
1229  NekDouble sign = -1.0;
1230  NekDouble beta;
1231 
1232  // Half Mode
1233  if (m_homogeneousBasis->GetBasisType() ==
1235  {
1236  beta = sign * 2 * M_PI * (m_transposition->GetK(0)) / m_lhom;
1237 
1238  Vmath::Smul(nP_pts, beta, temparray, 1, outarray, 1);
1239  }
1240  else if (m_homogeneousBasis->GetBasisType() ==
1242  {
1243  beta = -sign * 2 * M_PI * (m_transposition->GetK(0)) / m_lhom;
1244 
1245  Vmath::Smul(nP_pts, beta, temparray, 1, outarray, 1);
1246  }
1247 
1248  // Fully complex
1249  else
1250  {
1251  for (int i = 0; i < m_planes.size(); i++)
1252  {
1253  beta =
1254  -sign * 2 * M_PI * (m_transposition->GetK(i)) / m_lhom;
1255 
1256  Vmath::Smul(nP_pts, beta, tmp1 = temparray + i * nP_pts, 1,
1257  tmp2 = outarray + (i - int(sign)) * nP_pts, 1);
1258 
1259  sign = -1.0 * sign;
1260  }
1261  }
1262 
1263  if (m_WaveSpace)
1264  {
1265  out_d2 = outarray;
1266  }
1267  else
1268  {
1269  HomogeneousBwdTrans(outarray, out_d2);
1270  }
1271  }
1272  else
1273  {
1274  if (!m_session->DefinesSolverInfo("HomoStrip"))
1275  {
1276  ASSERTL0(m_comm->GetColumnComm()->GetSize() == 1,
1277  "Parallelisation in the homogeneous direction "
1278  "implemented just for Fourier basis");
1279  }
1280  else
1281  {
1282  ASSERTL0(m_StripZcomm->GetSize() == 1,
1283  "Parallelisation in the homogeneous direction "
1284  "implemented just for Fourier basis");
1285  }
1286 
1287  if (m_WaveSpace)
1288  {
1289  ASSERTL0(false, "Semi-phyisical time-stepping not "
1290  "implemented yet for non-Fourier "
1291  "basis");
1292  }
1293  else
1294  {
1295  StdRegions::StdSegExp StdSeg(m_homogeneousBasis->GetBasisKey());
1296 
1297  m_transposition->Transpose(inarray, temparray, false,
1299 
1300  for (int i = 0; i < nP_pts; i++)
1301  {
1302  StdSeg.PhysDeriv(temparray + i * m_planes.size(),
1303  tmp2 = outarray + i * m_planes.size());
1304  }
1305 
1306  m_transposition->Transpose(outarray, out_d2, false,
1308 
1309  Vmath::Smul(nT_pts, 2.0 / m_lhom, out_d2, 1, out_d2, 1);
1310  }
1311  }
1312  }
1313 }
1314 
1316  Direction edir, const Array<OneD, const NekDouble> &inarray,
1317  Array<OneD, NekDouble> &out_d)
1318 
1319 {
1320  int nT_pts = inarray.size(); // number of total points = n. of Fourier
1321  // points * n. of points per plane (nT_pts)
1322  int nP_pts =
1323  nT_pts / m_planes.size(); // number of points per plane = n of
1324  // Fourier transform required (nP_pts)
1325 
1326  int dir = (int)edir;
1327 
1328  Array<OneD, NekDouble> temparray(nT_pts);
1329  Array<OneD, NekDouble> outarray(nT_pts);
1332 
1333  if (dir < 2)
1334  {
1335  for (int i = 0; i < m_planes.size(); i++)
1336  {
1337  m_planes[i]->PhysDeriv(edir, inarray + i * nP_pts,
1338  tmp2 = out_d + i * nP_pts);
1339  }
1340  }
1341  else
1342  {
1343  if (m_homogeneousBasis->GetBasisType() == LibUtilities::eFourier ||
1344  m_homogeneousBasis->GetBasisType() ==
1346  m_homogeneousBasis->GetBasisType() ==
1348  m_homogeneousBasis->GetBasisType() ==
1350  {
1351  if (m_WaveSpace)
1352  {
1353  temparray = inarray;
1354  }
1355  else
1356  {
1357  HomogeneousFwdTrans(inarray, temparray);
1358  }
1359 
1360  NekDouble sign = -1.0;
1361  NekDouble beta;
1362 
1363  // HalfMode
1364  if (m_homogeneousBasis->GetBasisType() ==
1366  {
1367  beta = 2 * sign * M_PI * (m_transposition->GetK(0)) / m_lhom;
1368 
1369  Vmath::Smul(nP_pts, beta, temparray, 1, outarray, 1);
1370  }
1371  else if (m_homogeneousBasis->GetBasisType() ==
1373  {
1374  beta = -2 * sign * M_PI * (m_transposition->GetK(0)) / m_lhom;
1375 
1376  Vmath::Smul(nP_pts, beta, temparray, 1, outarray, 1);
1377  }
1378  // Fully complex
1379  else
1380  {
1381  for (int i = 0; i < m_planes.size(); i++)
1382  {
1383  beta =
1384  -sign * 2 * M_PI * (m_transposition->GetK(i)) / m_lhom;
1385 
1386  Vmath::Smul(nP_pts, beta, tmp1 = temparray + i * nP_pts, 1,
1387  tmp2 = outarray + (i - int(sign)) * nP_pts, 1);
1388 
1389  sign = -1.0 * sign;
1390  }
1391  }
1392  if (m_WaveSpace)
1393  {
1394  out_d = outarray;
1395  }
1396  else
1397  {
1398  HomogeneousBwdTrans(outarray, out_d);
1399  }
1400  }
1401  else
1402  {
1403  if (!m_session->DefinesSolverInfo("HomoStrip"))
1404  {
1405  ASSERTL0(m_comm->GetColumnComm()->GetSize() == 1,
1406  "Parallelisation in the homogeneous direction "
1407  "implemented just for Fourier basis");
1408  }
1409  else
1410  {
1411  ASSERTL0(m_StripZcomm->GetSize() == 1,
1412  "Parallelisation in the homogeneous direction "
1413  "implemented just for Fourier basis");
1414  }
1415 
1416  if (m_WaveSpace)
1417  {
1418  ASSERTL0(false, "Semi-phyisical time-stepping not implemented "
1419  "yet for non-Fourier basis");
1420  }
1421  else
1422  {
1423  StdRegions::StdSegExp StdSeg(m_homogeneousBasis->GetBasisKey());
1424 
1425  m_transposition->Transpose(inarray, temparray, false,
1427 
1428  for (int i = 0; i < nP_pts; i++)
1429  {
1430  StdSeg.PhysDeriv(temparray + i * m_planes.size(),
1431  tmp2 = outarray + i * m_planes.size());
1432  }
1433 
1434  m_transposition->Transpose(outarray, out_d, false,
1436 
1437  Vmath::Smul(nT_pts, 2.0 / m_lhom, out_d, 1, out_d, 1);
1438  }
1439  }
1440  }
1441 }
1442 
1444  const Array<OneD, const NekDouble> &inarray, Array<OneD, NekDouble> &out_d0,
1446 
1447 {
1448  v_PhysDeriv(inarray, out_d0, out_d1, out_d2);
1449 }
1450 
1452  Direction edir, const Array<OneD, const NekDouble> &inarray,
1453  Array<OneD, NekDouble> &out_d)
1454 {
1455  v_PhysDeriv(edir, inarray, out_d);
1456 }
1457 
1459  void)
1460 {
1461  return m_transposition;
1462 }
1463 
1465 {
1466  return m_lhom;
1467 }
1468 
1470 {
1471  m_lhom = lhom;
1472 }
1473 
1475 {
1476  return m_transposition->GetPlanesIDs();
1477 }
1478 
1480  const Array<OneD, const NekDouble> &inarray)
1481 {
1482  NekDouble val = 0.0;
1483  int i = 0;
1484 
1485  for (i = 0; i < (*m_exp).size(); ++i)
1486  {
1487  val += (*m_exp)[i]->Integral(inarray + m_phys_offset[i]);
1488  }
1489  val *= m_lhom / m_homogeneousBasis->GetNumModes();
1490 
1491  m_comm->AllReduce(val, LibUtilities::ReduceSum);
1492 
1493  return val;
1494 }
1495 } // namespace MultiRegions
1496 } // 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
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed to...
Definition: ErrorUtil.hpp:272
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:15
Describes the specification for a Basis.
Definition: Basis.h:50
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:144
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Abstraction of a two-dimensional multi-elemental expansion which is merely a collection of local expa...
void PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2)
virtual void v_AppendFieldData(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector< NekDouble > &fielddata)
virtual void v_PhysInterp1DScaled(const NekDouble scale, const Array< OneD, NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
NekDouble m_lhom
Width of homogeneous direction.
LibUtilities::TranspositionSharedPtr m_transposition
virtual void v_FwdTransBndConstrained(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void v_ExtractCoeffsToCoeffs(const std::shared_ptr< ExpList > &fromExpList, const Array< OneD, const NekDouble > &fromCoeffs, Array< OneD, NekDouble > &toCoeffs)
virtual void v_DealiasedProd(const Array< OneD, NekDouble > &inarray1, const Array< OneD, NekDouble > &inarray2, Array< OneD, NekDouble > &outarray)
void HomogeneousBwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool Shuff=true, bool UnShuff=true)
virtual void v_FwdTransLocalElmt(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Array< OneD, ExpListSharedPtr > m_planes
std::unordered_map< int, int > m_zIdToPlane
virtual void v_DealiasedDotProd(const Array< OneD, Array< OneD, NekDouble >> &inarray1, const Array< OneD, Array< OneD, NekDouble >> &inarray2, Array< OneD, Array< OneD, NekDouble >> &outarray)
NekDouble v_Integral(const Array< OneD, const NekDouble > &inarray)
LibUtilities::BasisSharedPtr m_homogeneousBasis
Definition of the total number of degrees of freedom and quadrature points. Sets up the storage for m...
void HomogeneousFwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool Shuff=true, bool UnShuff=true)
LibUtilities::NektarFFTSharedPtr m_FFT_deal
virtual void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual Array< OneD, const unsigned int > v_GetZIDs(void)
void Homogeneous1DTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool IsForwards, bool Shuff=true, bool UnShuff=true)
virtual void v_SetHomoLen(const NekDouble lhom)
LibUtilities::NektarFFTSharedPtr m_FFT
virtual void v_PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2)
virtual std::vector< LibUtilities::FieldDefinitionsSharedPtr > v_GetFieldDefinitions(void)
DNekBlkMatSharedPtr GenHomogeneous1DBlockMatrix(Homogeneous1DMatType mattype) const
virtual void v_WriteVtkPieceData(std::ostream &outfile, int expansion, std::string var)
virtual void v_HomogeneousFwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool Shuff=true, bool UnShuff=true)
virtual void v_ExtractDataToCoeffs(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector< NekDouble > &fielddata, std::string &field, Array< OneD, NekDouble > &coeffs)
Extract data from raw field data into expansion list.
virtual void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
DNekBlkMatSharedPtr GetHomogeneous1DBlockMatrix(Homogeneous1DMatType mattype) const
virtual void v_PhysGalerkinProjection1DScaled(const NekDouble scale, const Array< OneD, NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void v_HomogeneousBwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool Shuff=true, bool UnShuff=true)
virtual LibUtilities::TranspositionSharedPtr v_GetTransposition(void)
ExpListHomogeneous1D(const ExpansionType type)
Default constructor.
Base class for all multi-elemental spectral/hp expansions.
Definition: ExpList.h:101
Array< OneD, NekDouble > m_coeffs
Concatenation of all local expansion coefficients.
Definition: ExpList.h:1158
int GetNcoeffs(void) const
Returns the total number of local degrees of freedom .
Definition: ExpList.h:1641
Array< OneD, int > m_coeff_offset
Offset of elemental data into the array m_coeffs.
Definition: ExpList.h:1210
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition: ExpList.h:1126
std::shared_ptr< LocalRegions::ExpansionVector > m_exp
The list of local expansions.
Definition: ExpList.h:1196
std::unordered_map< int, int > m_elmtToExpId
Mapping from geometry ID of element to index inside m_exp.
Definition: ExpList.h:1227
LibUtilities::SessionReaderSharedPtr m_session
Session.
Definition: ExpList.h:1129
Array< OneD, int > m_phys_offset
Offset of elemental data into the array m_phys.
Definition: ExpList.h:1213
Array< OneD, NekDouble > m_phys
The global expansion evaluated at the quadrature points.
Definition: ExpList.h:1175
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
void PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1=NullNekDouble1DArray, Array< OneD, NekDouble > &out_d2=NullNekDouble1DArray)
Definition: StdExpansion.h:850
Class representing a segment element in reference space.
Definition: StdSegExp.h:54
static BasisSharedPtr NullBasisSharedPtr
Definition: Basis.h:369
BasisManagerT & BasisManager(void)
NektarFFTFactory & GetNektarFFTFactory()
Definition: NektarFFT.cpp:69
static const BasisKey NullBasisKey(eNoBasisType, 0, NullPointsKey)
Defines a null basis with no type or points.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< FieldDefinitions > FieldDefinitionsSharedPtr
Definition: FieldIO.h:177
int GetNumberOfCoefficients(ShapeType shape, std::vector< unsigned int > &modes, int offset=0)
Definition: ShapeType.hpp:308
std::shared_ptr< Transposition > TranspositionSharedPtr
@ beta
Gauss Radau pinned at x=-1,.
Definition: PointsType.h:61
@ eFourierEvenlySpaced
1D Evenly-spaced points using Fourier Fit
Definition: PointsType.h:76
@ eFourierSingleMode
Fourier ModifiedExpansion with just the first mode .
Definition: BasisType.h:66
@ eFourierHalfModeIm
Fourier Modified expansions with just the imaginary part of the first mode .
Definition: BasisType.h:70
@ eFourierHalfModeRe
Fourier Modified expansions with just the real part of the first mode .
Definition: BasisType.h:68
@ eFourier
Fourier Expansion .
Definition: BasisType.h:57
std::map< Homogeneous1DMatType, DNekBlkMatSharedPtr > Homo1DBlockMatrixMap
A map between homo matrix keys and their associated block matrices.
static const NekDouble kNekZeroTol
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
std::shared_ptr< DNekBlkMat > DNekBlkMatSharedPtr
Definition: NekTypeDefs.hpp:77
static Array< OneD, NekDouble > NullNekDouble1DArray
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 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 Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:248
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:492
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255