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 // 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: An ExpList which is homogeneous in 1-direction
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
37 #include <LibUtilities/Foundations/ManagerAccess.h> // for PointsManager, etc
38 #include <StdRegions/StdSegExp.h>
39 #include <StdRegions/StdPointExp.h>
40 #include <LocalRegions/Expansion.h>
42 
43 namespace Nektar
44 {
45  namespace MultiRegions
46  {
47  // Forward declaration for typedefs
49  ExpList(),
50  m_homogeneousBasis(LibUtilities::NullBasisSharedPtr),
51  m_lhom(1),
52  m_homogeneous1DBlockMat(MemoryManager<Homo1DBlockMatrixMap>::AllocateSharedPtr())
53  {
54  }
55 
57  &pSession,const LibUtilities::BasisKey &HomoBasis, const NekDouble lhom, const bool useFFT, const bool dealiasing):
58  ExpList(pSession),
59  m_useFFT(useFFT),
60  m_lhom(lhom),
61  m_homogeneous1DBlockMat(MemoryManager<Homo1DBlockMatrixMap>::AllocateSharedPtr()),
62  m_dealiasing(dealiasing)
63  {
64  ASSERTL2(HomoBasis != LibUtilities::NullBasisKey,"Homogeneous Basis is a null basis");
65 
67 
68  m_StripZcomm = m_session->DefinesSolverInfo("HomoStrip") ?
69  m_comm->GetColumnComm()->GetColumnComm() :
70  m_comm->GetColumnComm();
71 
74 
76  m_homogeneousBasis->GetNumPoints() /
77  m_StripZcomm->GetSize());
78 
79  if(m_useFFT)
80  {
82  "NekFFTW", m_homogeneousBasis->GetNumPoints());
83  }
84 
85  if(m_dealiasing)
86  {
87  if(m_useFFT)
88  {
89  NekDouble size = 1.5*m_homogeneousBasis->GetNumPoints();
90  m_padsize = int(size);
92  .CreateInstance("NekFFTW", m_padsize);
93  }
94  else
95  {
96  ASSERTL0(false, "Dealiasing available just in combination "
97  "with FFTW");
98  }
99  }
100  }
101 
102 
103  /**
104  * @param In ExpListHomogeneous1D object to copy.
105  */
107  ExpList(In,false),
108  m_transposition(In.m_transposition),
109  m_useFFT(In.m_useFFT),
110  m_FFT(In.m_FFT),
111  m_tmpIN(In.m_tmpIN),
112  m_tmpOUT(In.m_tmpOUT),
113  m_homogeneousBasis(In.m_homogeneousBasis),
114  m_lhom(In.m_lhom),
115  m_homogeneous1DBlockMat(In.m_homogeneous1DBlockMat),
116  m_dealiasing(In.m_dealiasing),
117  m_padsize(In.m_padsize)
118  {
119  m_planes = Array<OneD, ExpListSharedPtr>(In.m_planes.num_elements());
120  }
121 
122  /**
123  * Destructor
124  */
126  {
127  }
128 
130  Array<OneD, NekDouble> &outarray,
131  CoeffState coeffstate,
132  bool Shuff,
133  bool UnShuff)
134  {
135  // Forwards trans
136  Homogeneous1DTrans(inarray,outarray,true,coeffstate,Shuff,UnShuff);
137  }
138 
140  Array<OneD, NekDouble> &outarray,
141  CoeffState coeffstate,
142  bool Shuff,
143  bool UnShuff)
144  {
145  // Backwards trans
146  Homogeneous1DTrans(inarray,outarray,false,coeffstate,Shuff,UnShuff);
147  }
148 
149  /**
150  * Dealiasing routine
151  */
153  const Array<OneD, NekDouble> &inarray2,
154  Array<OneD, NekDouble> &outarray,
155  CoeffState coeffstate)
156  {
157  // inarray1 = first term of the product in full physical space
158  // inarray2 = second term of the product in full physical space
159  // dealiased product stored in outarray
160 
161  int num_dofs = inarray1.num_elements();
162 
163  int N = m_homogeneousBasis->GetNumPoints();
164 
165  Array<OneD, NekDouble> V1(num_dofs);
166  Array<OneD, NekDouble> V2(num_dofs);
167  Array<OneD, NekDouble> V1V2(num_dofs);
168 
169  HomogeneousFwdTrans(inarray1,V1,coeffstate);
170  HomogeneousFwdTrans(inarray2,V2,coeffstate);
171 
172  int num_points_per_plane = num_dofs/m_planes.num_elements();
173  int num_proc;
174  if(!m_session->DefinesSolverInfo("HomoStrip"))
175  {
176  num_proc = m_comm->GetColumnComm()->GetSize();
177  }
178  else
179  {
180  num_proc = m_StripZcomm->GetSize();
181  }
182  int num_dfts_per_proc = num_points_per_plane / num_proc
183  + (num_points_per_plane % num_proc > 0);
184 
185  Array<OneD, NekDouble> ShufV1(num_dfts_per_proc*N,0.0);
186  Array<OneD, NekDouble> ShufV2(num_dfts_per_proc*N,0.0);
187  Array<OneD, NekDouble> ShufV1V2(num_dfts_per_proc*N,0.0);
188 
189  Array<OneD, NekDouble> ShufV1_PAD_coef(m_padsize,0.0);
190  Array<OneD, NekDouble> ShufV2_PAD_coef(m_padsize,0.0);
191  Array<OneD, NekDouble> ShufV1_PAD_phys(m_padsize,0.0);
192  Array<OneD, NekDouble> ShufV2_PAD_phys(m_padsize,0.0);
193 
194  Array<OneD, NekDouble> ShufV1V2_PAD_coef(m_padsize,0.0);
195  Array<OneD, NekDouble> ShufV1V2_PAD_phys(m_padsize,0.0);
196 
197  m_transposition->Transpose(V1, ShufV1, false, LibUtilities::eXYtoZ);
198  m_transposition->Transpose(V2, ShufV2, false, LibUtilities::eXYtoZ);
199 
200  // Looping on the pencils
201  for(int i = 0 ; i < num_dfts_per_proc ; i++)
202  {
203  // Copying the i-th pencil pf lenght N into a bigger
204  // pencil of lenght 2N We are in Fourier space
205  Vmath::Vcopy(N, &(ShufV1[i*N]), 1, &(ShufV1_PAD_coef[0]), 1);
206  Vmath::Vcopy(N, &(ShufV2[i*N]), 1, &(ShufV2_PAD_coef[0]), 1);
207 
208  // Moving to physical space using the padded system
209  m_FFT_deal->FFTBwdTrans(ShufV1_PAD_coef, ShufV1_PAD_phys);
210  m_FFT_deal->FFTBwdTrans(ShufV2_PAD_coef, ShufV2_PAD_phys);
211 
212  // Perfroming the vectors multiplication in physical space on
213  // the padded system
214  Vmath::Vmul(m_padsize, ShufV1_PAD_phys, 1,
215  ShufV2_PAD_phys, 1,
216  ShufV1V2_PAD_phys, 1);
217 
218  // Moving back the result (V1*V2)_phys in Fourier space, padded
219  // system
220  m_FFT_deal->FFTFwdTrans(ShufV1V2_PAD_phys, ShufV1V2_PAD_coef);
221 
222  // Copying the first half of the padded pencil in the full
223  // vector (Fourier space)
224  Vmath::Vcopy(N, &(ShufV1V2_PAD_coef[0]), 1,
225  &(ShufV1V2[i*N]), 1);
226  }
227 
228  m_transposition->Transpose(ShufV1V2, V1V2, false,
230 
231  // Moving the results in physical space for the output
232  HomogeneousBwdTrans(V1V2, outarray, coeffstate);
233  }
234 
235  /**
236  * Forward transform
237  */
239  {
240  int cnt = 0, cnt1 = 0;
241  Array<OneD, NekDouble> tmparray;
242 
243  for(int n = 0; n < m_planes.num_elements(); ++n)
244  {
245  m_planes[n]->FwdTrans(inarray+cnt, tmparray = outarray + cnt1,
246  coeffstate);
247  cnt += m_planes[n]->GetTotPoints();
248 
249  cnt1 += m_planes[n]->GetNcoeffs(); // need to skip ncoeffs
250  }
251  if(!m_WaveSpace)
252  {
253  HomogeneousFwdTrans(outarray,outarray,coeffstate);
254  }
255  }
256 
257  /**
258  * Forward transform element by element
259  */
261  {
262  int cnt = 0, cnt1 = 0;
263  Array<OneD, NekDouble> tmparray;
264 
265  //spectral element FwdTrans plane by plane
266  for(int n = 0; n < m_planes.num_elements(); ++n)
267  {
268  m_planes[n]->FwdTrans_IterPerExp(inarray+cnt, tmparray = outarray + cnt1);
269  cnt += m_planes[n]->GetTotPoints();
270  cnt1 += m_planes[n]->GetNcoeffs();
271  }
272  if(!m_WaveSpace)
273  {
274  HomogeneousFwdTrans(outarray,outarray);
275  }
276  }
277 
278  /**
279  * Backward transform
280  */
282  {
283  int cnt = 0, cnt1 = 0;
284  Array<OneD, NekDouble> tmparray;
285 
286  for(int n = 0; n < m_planes.num_elements(); ++n)
287  {
288  m_planes[n]->BwdTrans(inarray+cnt, tmparray = outarray + cnt1,
289  coeffstate);
290  cnt += m_planes[n]->GetNcoeffs();
291  cnt1 += m_planes[n]->GetTotPoints();
292  }
293  if(!m_WaveSpace)
294  {
295  HomogeneousBwdTrans(outarray,outarray);
296  }
297  }
298 
299  /**
300  * Backward transform element by element
301  */
303  {
304  int cnt = 0, cnt1 = 0;
305  Array<OneD, NekDouble> tmparray;
306 
307  for(int n = 0; n < m_planes.num_elements(); ++n)
308  {
309  m_planes[n]->BwdTrans_IterPerExp(inarray+cnt, tmparray = outarray + cnt1);
310 
311  cnt += m_planes[n]->GetNcoeffs();
312  cnt1 += m_planes[n]->GetTotPoints();
313  }
314  if(!m_WaveSpace)
315  {
316  HomogeneousBwdTrans(outarray,outarray);
317  }
318  }
319 
320  /**
321  * Inner product
322  */
324  {
325  int cnt = 0, cnt1 = 0;
326  Array<OneD, NekDouble> tmparray;
327 
328  for(int n = 0; n < m_planes.num_elements(); ++n)
329  {
330  m_planes[n]->IProductWRTBase(inarray+cnt, tmparray = outarray + cnt1,coeffstate);
331 
332  cnt1 += m_planes[n]->GetNcoeffs();
333  cnt += m_planes[n]->GetTotPoints();
334  }
335  }
336 
337  /**
338  * Inner product element by element
339  */
341  {
342  int cnt = 0, cnt1 = 0;
343  Array<OneD, NekDouble> tmparray;
344 
345  for(int n = 0; n < m_planes.num_elements(); ++n)
346  {
347  m_planes[n]->IProductWRTBase_IterPerExp(inarray+cnt, tmparray = outarray + cnt1);
348 
349  cnt1 += m_planes[n]->GetNcoeffs();
350  cnt += m_planes[n]->GetTotPoints();
351  }
352  }
353 
354  /**
355  * Homogeneous transform Bwd/Fwd (MVM and FFT)
356  */
358  bool IsForwards,
359  CoeffState coeffstate,
360  bool Shuff,
361  bool UnShuff)
362  {
363  int num_dofs;
364 
365  if(IsForwards)
366  {
367  num_dofs = inarray.num_elements();
368  }
369  else
370  {
371  num_dofs = outarray.num_elements();
372  }
373 
374  if(m_useFFT)
375  {
376  int num_points_per_plane = num_dofs/m_planes.num_elements();
377  int num_dfts_per_proc;
378  if(!m_session->DefinesSolverInfo("HomoStrip"))
379  {
380  int nP = m_comm->GetColumnComm()->GetSize();
381  num_dfts_per_proc = num_points_per_plane / nP
382  + (num_points_per_plane % nP > 0);
383  }
384  else
385  {
386  int nP = m_StripZcomm->GetSize();
387  num_dfts_per_proc = num_points_per_plane / nP
388  + (num_points_per_plane % nP > 0);
389  }
390 
391  Array<OneD, NekDouble> fft_in (num_dfts_per_proc*m_homogeneousBasis->GetNumPoints(),0.0);
392  Array<OneD, NekDouble> fft_out(num_dfts_per_proc*m_homogeneousBasis->GetNumPoints(),0.0);
393 
394  if(Shuff)
395  {
396  m_transposition->Transpose(inarray,fft_in,false,LibUtilities::eXYtoZ);
397  }
398  else
399  {
400  Vmath::Vcopy(num_dfts_per_proc*m_homogeneousBasis->GetNumPoints(),
401  inarray,1,fft_in,1);
402  }
403 
404  if(IsForwards)
405  {
406  for(int i = 0 ; i < num_dfts_per_proc ; i++)
407  {
408  m_FFT->FFTFwdTrans(m_tmpIN = fft_in + i*m_homogeneousBasis->GetNumPoints(), m_tmpOUT = fft_out + i*m_homogeneousBasis->GetNumPoints());
409  }
410  }
411  else
412  {
413  for(int i = 0 ; i < num_dfts_per_proc ; i++)
414  {
415  m_FFT->FFTBwdTrans(m_tmpIN = fft_in + i*m_homogeneousBasis->GetNumPoints(), m_tmpOUT = fft_out + i*m_homogeneousBasis->GetNumPoints());
416  }
417  }
418 
419  if(UnShuff)
420  {
421  m_transposition->Transpose(fft_out,outarray,false,LibUtilities::eZtoXY);
422  }
423  else
424  {
425  Vmath::Vcopy(num_dfts_per_proc*m_homogeneousBasis->GetNumPoints(),
426  fft_out,1,outarray,1);
427  }
428  }
429  else
430  {
431  DNekBlkMatSharedPtr blkmat;
432 
433  if(num_dofs == m_npoints) //transform phys space
434  {
435  if(IsForwards)
436  {
438  }
439  else
440  {
442  }
443  }
444  else
445  {
446  if(IsForwards)
447  {
449  }
450  else
451  {
453  }
454  }
455 
456  int nrows = blkmat->GetRows();
457  int ncols = blkmat->GetColumns();
458 
459  Array<OneD, NekDouble> sortedinarray(ncols,0.0);
460  Array<OneD, NekDouble> sortedoutarray(nrows,0.0);
461 
462  if(Shuff)
463  {
464  m_transposition->Transpose(inarray,sortedinarray,!IsForwards,LibUtilities::eXYtoZ);
465  }
466  else
467  {
468  Vmath::Vcopy(ncols,inarray,1,sortedinarray,1);
469  }
470 
471  // Create NekVectors from the given data arrays
472  NekVector<NekDouble> in (ncols,sortedinarray,eWrapper);
473  NekVector<NekDouble> out(nrows,sortedoutarray,eWrapper);
474 
475  // Perform matrix-vector multiply.
476  out = (*blkmat)*in;
477 
478  if(UnShuff)
479  {
480  m_transposition->Transpose(sortedoutarray,outarray,IsForwards,LibUtilities::eZtoXY);
481  }
482  else
483  {
484  Vmath::Vcopy(nrows,sortedinarray,1,outarray,1);
485  }
486 
487  }
488  }
489 
491  {
492  Homo1DBlockMatrixMap::iterator matrixIter = m_homogeneous1DBlockMat->find(mattype);
493 
494  if(matrixIter == m_homogeneous1DBlockMat->end())
495  {
496  return ((*m_homogeneous1DBlockMat)[mattype] =
497  GenHomogeneous1DBlockMatrix(mattype,coeffstate));
498  }
499  else
500  {
501  return matrixIter->second;
502  }
503  }
504 
505 
507  {
508  DNekMatSharedPtr loc_mat;
509  DNekBlkMatSharedPtr BlkMatrix;
510  int n_exp = 0;
511  int num_trans_per_proc = 0;
512 
513  if((mattype == eForwardsCoeffSpace1D)
514  ||(mattype == eBackwardsCoeffSpace1D)) // will operate on m_coeffs
515  {
516  n_exp = m_planes[0]->GetNcoeffs();
517  }
518  else
519  {
520  n_exp = m_planes[0]->GetTotPoints(); // will operatore on m_phys
521  }
522 
523  num_trans_per_proc = n_exp/m_comm->GetColumnComm()->GetSize() + (n_exp%m_comm->GetColumnComm()->GetSize() > 0);
524 
525  Array<OneD,unsigned int> nrows(num_trans_per_proc);
526  Array<OneD,unsigned int> ncols(num_trans_per_proc);
527 
528  if((mattype == eForwardsCoeffSpace1D)||(mattype == eForwardsPhysSpace1D))
529  {
530  nrows = Array<OneD, unsigned int>(num_trans_per_proc,m_homogeneousBasis->GetNumModes());
531  ncols = Array<OneD, unsigned int>(num_trans_per_proc,m_homogeneousBasis->GetNumPoints());
532  }
533  else
534  {
535  nrows = Array<OneD, unsigned int>(num_trans_per_proc,m_homogeneousBasis->GetNumPoints());
536  ncols = Array<OneD, unsigned int>(num_trans_per_proc,m_homogeneousBasis->GetNumModes());
537  }
538 
539  MatrixStorage blkmatStorage = eDIAGONAL;
540  BlkMatrix = MemoryManager<DNekBlkMat>
541  ::AllocateSharedPtr(nrows,ncols,blkmatStorage);
542 
543  //Half Mode
545  {
546  StdRegions::StdPointExp StdPoint(m_homogeneousBasis->GetBasisKey());
547 
548  if((mattype == eForwardsCoeffSpace1D)||(mattype == eForwardsPhysSpace1D))
549  {
551  StdPoint.DetShapeType(),
552  StdPoint);
553 
554  loc_mat = StdPoint.GetStdMatrix(matkey);
555  }
556  else
557  {
559  StdPoint.DetShapeType(),
560  StdPoint);
561 
562  loc_mat = StdPoint.GetStdMatrix(matkey);
563  }
564  }
565  //other cases
566  else
567  {
568  StdRegions::StdSegExp StdSeg(m_homogeneousBasis->GetBasisKey());
569 
570  if((mattype == eForwardsCoeffSpace1D)||(mattype == eForwardsPhysSpace1D))
571  {
573  StdSeg.DetShapeType(),
574  StdSeg);
575 
576  loc_mat = StdSeg.GetStdMatrix(matkey);
577  }
578  else
579  {
581  StdSeg.DetShapeType(),
582  StdSeg);
583 
584  loc_mat = StdSeg.GetStdMatrix(matkey);
585  }
586 
587  }
588 
589  // set up array of block matrices.
590  for(int i = 0; i < num_trans_per_proc; ++i)
591  {
592  BlkMatrix->SetBlock(i,i,loc_mat);
593  }
594 
595  return BlkMatrix;
596  }
597 
598  std::vector<LibUtilities::FieldDefinitionsSharedPtr> ExpListHomogeneous1D::v_GetFieldDefinitions()
599  {
600  std::vector<LibUtilities::FieldDefinitionsSharedPtr> returnval;
601 
602  // Set up Homogeneous length details.
604 
605  std::vector<NekDouble> HomoLen;
606  HomoLen.push_back(m_lhom);
607 
608  std::vector<unsigned int> PlanesIDs;
609 
610  for(int i = 0; i < m_planes.num_elements(); i++)
611  {
612  PlanesIDs.push_back(m_transposition->GetPlaneID(i));
613  }
614 
615  int NumHomoStrip;
616  m_session->LoadParameter("Strip_Z",NumHomoStrip,1);
617 
618  m_planes[0]->GeneralGetFieldDefinitions(returnval, 1, NumHomoStrip, HomoBasis, HomoLen, PlanesIDs);
619 
620  return returnval;
621  }
622 
623  void ExpListHomogeneous1D::v_GetFieldDefinitions(std::vector<LibUtilities::FieldDefinitionsSharedPtr> &fielddef)
624  {
625  // Set up Homogeneous length details.
627 
628  std::vector<NekDouble> HomoLen;
629  HomoLen.push_back(m_lhom);
630 
631  std::vector<unsigned int> PlanesIDs;
632 
633  for(int i = 0; i < m_planes.num_elements(); i++)
634  {
635  PlanesIDs.push_back(m_transposition->GetPlaneID(i));
636  }
637 
638  int NumHomoStrip;
639  m_session->LoadParameter("Strip_Z",NumHomoStrip,1);
640 
641  // enforce NumHomoDir == 1 by direct call
642  m_planes[0]->GeneralGetFieldDefinitions(fielddef, 1, NumHomoStrip, HomoBasis,HomoLen,PlanesIDs);
643  }
644 
645 
646  /** This routine appends the data from the expansion list into
647  the output format where each element is given by looping
648  over its Fourier modes where as data in the expandion is
649  stored with all consecutive elements and then the Fourier
650  modes
651  */
653  {
654  int i,n;
655  int ncoeffs_per_plane = m_planes[0]->GetNcoeffs();
656 
657  // Determine mapping from element ids to location in
658  // expansion list
659  map<int, int> ElmtID_to_ExpID;
660  for(i = 0; i < m_planes[0]->GetExpSize(); ++i)
661  {
662  ElmtID_to_ExpID[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
663  }
664 
665  for(i = 0; i < fielddef->m_elementIDs.size(); ++i)
666  {
667  int eid = ElmtID_to_ExpID[fielddef->m_elementIDs[i]];
668  int datalen = (*m_exp)[eid]->GetNcoeffs();
669 
670  for(n = 0; n < m_planes.num_elements(); ++n)
671  {
672  fielddata.insert(fielddata.end(),&coeffs[m_coeff_offset[eid]+n*ncoeffs_per_plane],&coeffs[m_coeff_offset[eid]+n*ncoeffs_per_plane]+datalen);
673  }
674  }
675  }
676 
677  void ExpListHomogeneous1D::v_AppendFieldData(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector<NekDouble> &fielddata)
678  {
679  v_AppendFieldData(fielddef,fielddata,m_coeffs);
680  }
681 
682  //Extract the data in fielddata into the m_coeff list
685  std::vector<NekDouble> &fielddata,
686  std::string &field,
687  Array<OneD, NekDouble> &coeffs)
688  {
689  int i,n;
690  int offset = 0;
691  int nzmodes = 1;
692  int datalen = fielddata.size()/fielddef->m_fields.size();
693  std::vector<unsigned int> fieldDefHomoZids;
694 
695 
696  // Find data location according to field definition
697  for(i = 0; i < fielddef->m_fields.size(); ++i)
698  {
699  if(fielddef->m_fields[i] == field)
700  {
701  break;
702  }
703  offset += datalen;
704  }
705 
706  if(i == fielddef->m_fields.size())
707  {
708  cout << "Field "<< field<< "not found in data file. " << endl;
709  }
710  else
711  {
712 
713  int modes_offset = 0;
714  int planes_offset = 0;
715  Array<OneD, NekDouble> coeff_tmp;
717 
718 
719  // Build map of plane IDs lying on this processor.
720  std::map<int,int> homoZids;
721  for (i = 0; i < m_planes.num_elements(); ++i)
722  {
723  homoZids[m_transposition->GetPlaneID(i)] = i;
724  }
725 
726  if(fielddef->m_numHomogeneousDir)
727  {
728  for(i = 0; i < fielddef->m_basis.size(); ++i)
729  {
730  if(fielddef->m_basis[i] == m_homogeneousBasis->GetBasisType())
731  {
732  nzmodes = fielddef->m_homogeneousZIDs.size();
733  break;
734  }
735  }
736  ASSERTL1(i != fielddef->m_basis.size(),"Failed to determine number of Homogeneous modes");
737 
738  fieldDefHomoZids = fielddef->m_homogeneousZIDs;
739  }
740  else // input file is 2D and so set nzmodes to 1
741  {
742  nzmodes = 1;
743  fieldDefHomoZids.push_back(0);
744  }
745 
746  // Determine mapping from element ids to location in expansion list.
747  map<int, int> ElmtID_to_ExpID;
748  for(i = 0; i < m_planes[0]->GetExpSize(); ++i)
749  {
750  ElmtID_to_ExpID[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
751  }
752 
753 
754  // calculate number of modes in the current partition
755  int ncoeffs_per_plane = m_planes[0]->GetNcoeffs();
756 
757  for(i = 0; i < fielddef->m_elementIDs.size(); ++i)
758  {
759  if(fielddef->m_uniOrder == true) // reset modes_offset to zero
760  {
761  modes_offset = 0;
762  }
763 
764  int datalen = LibUtilities::GetNumberOfCoefficients(fielddef->m_shapeType,
765  fielddef->m_numModes,
766  modes_offset);
767 
768  it = ElmtID_to_ExpID.find(fielddef->m_elementIDs[i]);
769 
770  // ensure element is on this partition for parallel case.
771  if(it == ElmtID_to_ExpID.end())
772  {
773  // increase offset for correct FieldData access
774  offset += datalen*nzmodes;
775  continue;
776  }
777 
778  int eid = it->second;
779 
780 
781  for(n = 0; n < nzmodes; ++n, offset += datalen)
782  {
783 
784  it = homoZids.find(fieldDefHomoZids[n]);
785 
786  // Check to make sure this mode number lies in this field.
787  if (it == homoZids.end())
788  {
789  continue;
790  }
791 
792  planes_offset = it->second;
793  if(datalen == (*m_exp)[eid]->GetNcoeffs())
794  {
795  Vmath::Vcopy(datalen,&fielddata[offset],1,&coeffs[m_coeff_offset[eid]+planes_offset*ncoeffs_per_plane],1);
796  }
797  else // unpack data to new order
798  {
799  (*m_exp)[eid]->ExtractDataToCoeffs(&fielddata[offset], fielddef->m_numModes,modes_offset,&coeffs[m_coeff_offset[eid] + planes_offset*ncoeffs_per_plane]);
800  }
801  }
802  modes_offset += (*m_exp)[0]->GetNumBases();
803  }
804  }
805  }
806 
807  //Extract the data in fielddata into the m_coeff list
809  const boost::shared_ptr<ExpList> &fromExpList,const Array<OneD, const NekDouble> &fromCoeffs, Array<OneD, NekDouble> &toCoeffs)
810  {
811  int i;
812  int fromNcoeffs_per_plane = fromExpList->GetPlane(0)->GetNcoeffs();
813  Array<OneD, NekDouble> tocoeffs_tmp, fromcoeffs_tmp;
814 
815  for(i = 0; i < m_planes.num_elements(); ++i)
816  {
817  m_planes[i]->ExtractCoeffsToCoeffs(fromExpList->GetPlane(i),fromcoeffs_tmp = fromCoeffs + fromNcoeffs_per_plane*i, tocoeffs_tmp = toCoeffs + m_ncoeffs*i);
818  }
819  }
820 
821  void ExpListHomogeneous1D::v_WriteVtkPieceData(std::ostream &outfile, int expansion,
822  std::string var)
823  {
824  int i;
825  int nq = (*m_exp)[expansion]->GetTotPoints();
826  int npoints_per_plane = m_planes[0]->GetTotPoints();
827 
828  // printing the fields of that zone
829  outfile << " <DataArray type=\"Float32\" Name=\""
830  << var << "\">" << endl;
831  outfile << " ";
832  for (int n = 0; n < m_planes.num_elements(); ++n)
833  {
834  const Array<OneD, NekDouble> phys = m_phys + m_phys_offset[expansion] + n*npoints_per_plane;
835  for(i = 0; i < nq; ++i)
836  {
837  outfile << (fabs(phys[i]) < NekConstants::kNekZeroTol ? 0 : phys[i]) << " ";
838  }
839  }
840  outfile << endl;
841  outfile << " </DataArray>" << endl;
842  }
843 
845  {
846  int cnt,cnt1;
847  Array<OneD, NekDouble> tmparray;
848  cnt = m_planes[0]->GetTotPoints();
849  cnt1 = m_planes[0]->Get1DScaledTotPoints(scale);
850 
851  ASSERTL1(m_planes.num_elements()*cnt1 <= outarray.num_elements(),"size of outarray does not match internal estimage");
852 
853 
854  for(int i = 0; i < m_planes.num_elements(); i++)
855  {
856 
857  m_planes[i]->PhysInterp1DScaled(scale,inarray+i*cnt,
858  tmparray = outarray+i*cnt1);
859  }
860  }
861 
862 
864  {
865  int cnt,cnt1;
866  Array<OneD, NekDouble> tmparray;
867  cnt = m_planes[0]->Get1DScaledTotPoints(scale);
868  cnt1 = m_planes[0]->GetTotPoints();
869 
870  ASSERTL1(m_planes.num_elements()*cnt <= inarray.num_elements(),"size of outarray does not match internal estimage");
871 
872 
873  for(int i = 0; i < m_planes.num_elements(); i++)
874  {
875  m_planes[i]->PhysGalerkinProjection1DScaled(scale,inarray+i*cnt,
876  tmparray = outarray+i*cnt1);
877  }
878 
879  }
881  Array<OneD, NekDouble> &out_d0,
882  Array<OneD, NekDouble> &out_d1,
883  Array<OneD, NekDouble> &out_d2)
884  {
885  int nT_pts = inarray.num_elements(); //number of total points = n. of Fourier points * n. of points per plane (nT_pts)
886  int nP_pts = nT_pts/m_planes.num_elements(); //number of points per plane = n of Fourier transform required (nP_pts)
887 
888  Array<OneD, NekDouble> temparray(nT_pts);
889  Array<OneD, NekDouble> outarray(nT_pts);
893 
894  for(int i = 0; i < m_planes.num_elements(); i++)
895  {
896  m_planes[i]->PhysDeriv(inarray + i*nP_pts ,tmp2 = out_d0 + i*nP_pts , tmp3 = out_d1 + i*nP_pts );
897  }
898 
899  if(out_d2 != NullNekDouble1DArray)
900  {
903  {
904  if(m_WaveSpace)
905  {
906  temparray = inarray;
907  }
908  else
909  {
910  HomogeneousFwdTrans(inarray,temparray);
911  }
912 
913  NekDouble sign = -1.0;
914  NekDouble beta;
915 
916  //Half Mode
918  {
919  beta = sign*2*M_PI*(m_transposition->GetK(0))/m_lhom;
920 
921  Vmath::Smul(nP_pts,beta,temparray,1,outarray,1);
922  }
923  else if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeIm)
924  {
925  beta = -sign*2*M_PI*(m_transposition->GetK(0))/m_lhom;
926 
927  Vmath::Smul(nP_pts,beta,temparray,1,outarray,1);
928  }
929 
930  //Fully complex
931  else
932  {
933  for(int i = 0; i < m_planes.num_elements(); i++)
934  {
935  beta = -sign*2*M_PI*(m_transposition->GetK(i))/m_lhom;
936 
937  Vmath::Smul(nP_pts,beta,tmp1 = temparray + i*nP_pts,1,tmp2 = outarray + (i-int(sign))*nP_pts,1);
938 
939  sign = -1.0*sign;
940  }
941  }
942 
943  if(m_WaveSpace)
944  {
945  out_d2 = outarray;
946  }
947  else
948  {
949  HomogeneousBwdTrans(outarray,out_d2);
950  }
951  }
952  else
953  {
954  if(!m_session->DefinesSolverInfo("HomoStrip"))
955  {
956  ASSERTL0(m_comm->GetColumnComm()->GetSize() == 1,
957  "Parallelisation in the homogeneous direction "
958  "implemented just for Fourier basis");
959  }
960  else
961  {
962  ASSERTL0(m_StripZcomm->GetSize() == 1,
963  "Parallelisation in the homogeneous direction "
964  "implemented just for Fourier basis");
965  }
966 
967  if(m_WaveSpace)
968  {
969  ASSERTL0(false, "Semi-phyisical time-stepping not "
970  "implemented yet for non-Fourier "
971  "basis");
972  }
973  else
974  {
975  StdRegions::StdSegExp StdSeg(m_homogeneousBasis->GetBasisKey());
976 
977  m_transposition->Transpose(inarray,temparray,false,LibUtilities::eXYtoZ);
978 
979  for(int i = 0; i < nP_pts; i++)
980  {
981  StdSeg.PhysDeriv(temparray + i*m_planes.num_elements(), tmp2 = outarray + i*m_planes.num_elements());
982  }
983 
984  m_transposition->Transpose(outarray,out_d2,false,LibUtilities::eZtoXY);
985 
986  Vmath::Smul(nT_pts,2.0/m_lhom,out_d2,1,out_d2,1);
987  }
988  }
989  }
990  }
991 
994 
995  {
996  int nT_pts = inarray.num_elements(); //number of total points = n. of Fourier points * n. of points per plane (nT_pts)
997  int nP_pts = nT_pts/m_planes.num_elements(); //number of points per plane = n of Fourier transform required (nP_pts)
998 
999  int dir= (int)edir;
1000 
1001  Array<OneD, NekDouble> temparray(nT_pts);
1002  Array<OneD, NekDouble> outarray(nT_pts);
1005 
1006  if (dir < 2)
1007  {
1008  for(int i=0; i < m_planes.num_elements(); i++)
1009  {
1010  m_planes[i]->PhysDeriv(edir, inarray + i*nP_pts ,tmp2 = out_d + i*nP_pts);
1011  }
1012  }
1013  else
1014  {
1017  {
1018  if(m_WaveSpace)
1019  {
1020  temparray = inarray;
1021  }
1022  else
1023  {
1024  HomogeneousFwdTrans(inarray,temparray);
1025  }
1026 
1027  NekDouble sign = -1.0;
1028  NekDouble beta;
1029 
1030  //HalfMode
1032  {
1033  beta = 2*sign*M_PI*(m_transposition->GetK(0))/m_lhom;
1034 
1035  Vmath::Smul(nP_pts,beta,temparray,1,outarray,1);
1036  }
1037  else if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeIm)
1038  {
1039  beta = -2*sign*M_PI*(m_transposition->GetK(0))/m_lhom;
1040 
1041  Vmath::Smul(nP_pts,beta,temparray,1,outarray,1);
1042  }
1043  //Fully complex
1044  else
1045  {
1046  for(int i = 0; i < m_planes.num_elements(); i++)
1047  {
1048  beta = -sign*2*M_PI*(m_transposition->GetK(i))/m_lhom;
1049 
1050  Vmath::Smul(nP_pts,beta,tmp1 = temparray + i*nP_pts,1,tmp2 = outarray + (i-int(sign))*nP_pts,1);
1051 
1052  sign = -1.0*sign;
1053  }
1054  }
1055  if(m_WaveSpace)
1056  {
1057  out_d = outarray;
1058  }
1059  else
1060  {
1061  HomogeneousBwdTrans(outarray,out_d);
1062  }
1063  }
1064  else
1065  {
1066  if(!m_session->DefinesSolverInfo("HomoStrip"))
1067  {
1068  ASSERTL0(m_comm->GetColumnComm()->GetSize() == 1,
1069  "Parallelisation in the homogeneous direction "
1070  "implemented just for Fourier basis");
1071  }
1072  else
1073  {
1074  ASSERTL0(m_StripZcomm->GetSize() == 1,
1075  "Parallelisation in the homogeneous direction "
1076  "implemented just for Fourier basis");
1077  }
1078 
1079  if(m_WaveSpace)
1080  {
1081  ASSERTL0(false,"Semi-phyisical time-stepping not implemented yet for non-Fourier basis");
1082  }
1083  else
1084  {
1085  StdRegions::StdSegExp StdSeg(m_homogeneousBasis->GetBasisKey());
1086 
1087  m_transposition->Transpose(inarray,temparray,false,LibUtilities::eXYtoZ);
1088 
1089  for(int i = 0; i < nP_pts; i++)
1090  {
1091  StdSeg.PhysDeriv(temparray + i*m_planes.num_elements(), tmp2 = outarray + i*m_planes.num_elements());
1092  }
1093 
1094  m_transposition->Transpose(outarray,out_d,false,LibUtilities::eZtoXY);
1095 
1096  Vmath::Smul(nT_pts,2.0/m_lhom,out_d,1,out_d,1);
1097  }
1098  }
1099  }
1100  }
1101 
1103  Array<OneD, NekDouble> &out_d0,
1104  Array<OneD, NekDouble> &out_d1,
1105  Array<OneD, NekDouble> &out_d2)
1106 
1107  {
1108  v_PhysDeriv(inarray,out_d0,out_d1,out_d2);
1109  }
1110 
1112  const Array<OneD, const NekDouble> &inarray,
1113  Array<OneD, NekDouble> &out_d)
1114  {
1115  v_PhysDeriv(edir,inarray,out_d);
1116  }
1117 
1119  {
1120  return m_transposition;
1121  }
1122 
1124  {
1125  return m_lhom;
1126  }
1127 
1129  {
1130  return m_transposition->GetPlanesIDs();
1131  }
1132  } //end of namespace
1133 } //end of namespace
1134 
1135 
1136 /**
1137 * $Log: v $
1138 *
1139 **/
1140 
Abstraction of a two-dimensional multi-elemental expansion which is merely a collection of local expa...
void HomogeneousBwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate=eLocal, bool Shuff=true, bool UnShuff=true)
virtual void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate)
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
DNekBlkMatSharedPtr GetHomogeneous1DBlockMatrix(Homogeneous1DMatType mattype, CoeffState coeffstate=eLocal) const
virtual void v_PhysGalerkinProjection1DScaled(const NekDouble scale, const Array< OneD, NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
boost::shared_ptr< Transposition > TranspositionSharedPtr
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 Array< OneD, NekDouble > NullNekDouble1DArray
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:22
LibUtilities::TranspositionSharedPtr m_transposition
virtual void v_WriteVtkPieceData(std::ostream &outfile, int expansion, std::string var)
virtual void v_PhysInterp1DScaled(const NekDouble scale, const Array< OneD, NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
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_ExtractCoeffsToCoeffs(const boost::shared_ptr< ExpList > &fromExpList, const Array< OneD, const NekDouble > &fromCoeffs, Array< OneD, NekDouble > &toCoeffs)
boost::shared_ptr< FieldDefinitions > FieldDefinitionsSharedPtr
Definition: FieldIO.h:131
LibUtilities::NektarFFTSharedPtr m_FFT_deal
static BasisSharedPtr NullBasisSharedPtr
Definition: Basis.h:358
NekDouble m_lhom
Width of homogeneous direction.
Array< OneD, NekDouble > m_phys
The global expansion evaluated at the quadrature points.
Definition: ExpList.h:926
Fourier Expansion .
Definition: BasisType.h:52
virtual void v_HomogeneousFwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate=eLocal, bool Shuff=true, bool UnShuff=true)
Array< OneD, NekDouble > m_coeffs
Concatenation of all local expansion coefficients.
Definition: ExpList.h:909
boost::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:70
NektarFFTFactory & GetNektarFFTFactory()
Definition: NektarFFT.cpp:69
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:50
virtual LibUtilities::TranspositionSharedPtr v_GetTransposition(void)
int GetNumberOfCoefficients(ShapeType shape, std::vector< unsigned int > &modes, int offset)
Definition: ShapeType.hpp:312
virtual void v_AppendFieldData(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector< NekDouble > &fielddata)
BasisManagerT & BasisManager(void)
Class representing a segment element in reference space.
Definition: StdSegExp.h:54
virtual void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate)
Array< OneD, int > m_coeff_offset
Offset of elemental data into the array m_coeffs.
Definition: ExpList.h:958
LibUtilities::BasisSharedPtr m_homogeneousBasis
Definition of the total number of degrees of freedom and quadrature points. Sets up the storage for m...
Base class for all multi-elemental spectral/hp expansions.
Definition: ExpList.h:101
static const NekDouble kNekZeroTol
Fourier Modified expansions with just the real part of the first mode .
Definition: BasisType.h:59
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:199
boost::shared_ptr< LocalRegions::ExpansionVector > m_exp
The list of local expansions.
Definition: ExpList.h:947
Array< OneD, int > m_phys_offset
Offset of elemental data into the array m_phys.
Definition: ExpList.h:961
virtual void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate)
int m_ncoeffs
The total number of local degrees of freedom. m_ncoeffs .
Definition: ExpList.h:887
virtual void v_BwdTrans_IterPerExp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
LibUtilities::SessionReaderSharedPtr m_session
Session.
Definition: ExpList.h:880
Array< OneD, ExpListSharedPtr > m_planes
double NekDouble
void HomogeneousFwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate=eLocal, bool Shuff=true, bool UnShuff=true)
virtual void v_FwdTrans_IterPerExp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void v_HomogeneousBwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate=eLocal, bool Shuff=true, bool UnShuff=true)
Fourier Modified expansions with just the imaginary part of the first mode .
Definition: BasisType.h:60
virtual Array< OneD, const unsigned int > v_GetZIDs(void)
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.
boost::shared_ptr< DNekBlkMat > DNekBlkMatSharedPtr
Definition: NekTypeDefs.hpp:72
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
virtual void v_DealiasedProd(const Array< OneD, NekDouble > &inarray1, const Array< OneD, NekDouble > &inarray2, Array< OneD, NekDouble > &outarray, CoeffState coeffstate=eLocal)
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition: ExpList.h:877
void Homogeneous1DTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool IsForwards, CoeffState coeffstate=eLocal, bool Shuff=true, bool UnShuff=true)
virtual std::vector< LibUtilities::FieldDefinitionsSharedPtr > v_GetFieldDefinitions(void)
DNekBlkMatSharedPtr GenHomogeneous1DBlockMatrix(Homogeneous1DMatType mattype, CoeffState coeffstate=eLocal) const
map< Homogeneous1DMatType, DNekBlkMatSharedPtr > Homo1DBlockMatrixMap
A map between homo matrix keys and their associated block matrices.
virtual void v_IProductWRTBase_IterPerExp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Fourier ModifiedExpansion with just the first mode .
Definition: BasisType.h:58
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed t...
Definition: ErrorUtil.hpp:213
static const BasisKey NullBasisKey(eNoBasisType, 0, NullPointsKey)
Defines a null basis with no type or points.
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:191
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1038
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)
Describes the specification for a Basis.
Definition: Basis.h:50
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:169
LibUtilities::NektarFFTSharedPtr m_FFT