Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Transposition.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File Transposition.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: Data manager for homogeneous transpositions
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
37 
38 #include <LibUtilities/BasicUtils/ErrorUtil.hpp> // for ASSERTL0, etc
39 #include <LibUtilities/BasicUtils/SharedArray.hpp> // for Array
40 #include <LibUtilities/BasicUtils/Vmath.hpp> // for Vcopy
41 #include <LibUtilities/Foundations/Basis.h> // for BasisKey
43 
44 namespace Nektar
45 {
46 namespace LibUtilities
47 {
48 /**
49  * Constructor for 1D transform.
50  */
54 {
55  m_hcomm = hcomm1;
57 
62 
63  m_num_homogeneous_points[0] = HomoBasis0.GetNumPoints();
64  m_num_homogeneous_coeffs[0] = HomoBasis0.GetNumModes();
65  m_num_processes[0] = m_hcomm->GetSize();
66  m_num_points_per_proc[0] = m_num_homogeneous_points[0] / m_num_processes[0];
67  m_rank_id = m_hcomm->GetRank();
68 
69  //================================================================
70  // TODO: Need to be generalised for 1D, 2D and 3D
73 
74  for (int i = 0; i < m_num_points_per_proc[0]; i++)
75  {
76  m_planes_IDs[i] = m_rank_id * m_num_points_per_proc[0] + i;
77  }
78 
79  int global_rank_id = hcomm0->GetColumnComm()->GetRank();
80  int NumStrips = hcomm0->GetColumnComm()->GetSize() / m_hcomm->GetSize();
81  m_strip_ID = 0;
82 
83  if (NumStrips > 1)
84  {
85  m_strip_ID = (NumStrips > global_rank_id)
86  ? global_rank_id
87  : (global_rank_id - NumStrips);
88  }
89 
90  if (HomoBasis0.GetBasisType() == LibUtilities::eFourier)
91  {
92  for (int i = 0; i < m_num_points_per_proc[0]; i++)
93  {
94  m_K[i] = m_planes_IDs[i] / 2;
95  }
96  }
97 
99  {
100  m_K[0] = 1;
101  m_K[1] = 1;
102  }
103 
104  if (HomoBasis0.GetBasisType() == LibUtilities::eFourierHalfModeRe ||
106  {
107  m_K[0] = 1;
108  }
109  //================================================================
110 }
111 
112 /**
113  * Constructor for 2D transform.
114  */
116  const LibUtilities::BasisKey &HomoBasis1,
118 {
119  m_hcomm = hcomm;
121 
126 
127  m_num_homogeneous_points[0] = HomoBasis0.GetNumPoints();
128  m_num_homogeneous_coeffs[0] = HomoBasis0.GetNumModes();
129  m_num_homogeneous_points[1] = HomoBasis1.GetNumPoints();
130  m_num_homogeneous_coeffs[1] = HomoBasis1.GetNumModes();
131 
132  m_num_processes[0] = m_hcomm->GetRowComm()->GetSize();
133  m_num_processes[1] = m_hcomm->GetColumnComm()->GetSize();
134 
135  m_num_points_per_proc[0] = m_num_homogeneous_points[0] / m_num_processes[0];
136  m_num_points_per_proc[1] = m_num_homogeneous_points[1] / m_num_processes[1];
137 
138  //================================================================
139  // TODO: Need set up for 2D lines IDs and Ks if Fourier
140  //================================================================
141 }
142 
143 /**
144  * Constructor for 3D transform.
145  */
147  const LibUtilities::BasisKey &HomoBasis1,
148  const LibUtilities::BasisKey &HomoBasis2,
150 {
151  m_hcomm = hcomm;
153 
158 
159  //================================================================
160  // TODO: Need set up for 3D
161  ASSERTL0(false, "Transposition is not set up for 3D.");
162  //================================================================
163 }
164 
165 /**
166  * Destructor
167  */
169 {
170 }
171 
172 //====================================================================
173 // TODO: Need to generalise the following methods for 1D, 2D and 3D
174 unsigned int Transposition::GetK(int i)
175 {
176  return m_K[i];
177 }
178 
180 {
181  return m_K;
182 }
183 
184 unsigned int Transposition::GetPlaneID(int i)
185 {
186  return m_planes_IDs[i];
187 }
188 
190 {
191  return m_planes_IDs;
192 }
193 
194 unsigned int Transposition::GetStripID(void)
195 {
196  return m_strip_ID;
197 }
198 
199 /**
200  * Main method: General transposition, the dir parameters define if
201  * 1D,2D,3D and which transposition is required at the same time
202  */
204  Array<OneD, NekDouble> &outarray, bool UseNumMode,
205  TranspositionDir dir)
206 {
207  switch (dir)
208  {
209  case eXYtoZ:
210  {
211  TransposeXYtoZ(inarray, outarray, UseNumMode);
212  }
213  break;
214  case eZtoXY:
215  {
216  TransposeZtoXY(inarray, outarray, UseNumMode);
217  }
218  break;
219  case eXtoYZ:
220  {
221  TransposeXtoYZ(inarray, outarray, UseNumMode);
222  }
223  break;
224  case eYZtoX:
225  {
226  TransposeYZtoX(inarray, outarray, UseNumMode);
227  }
228  break;
229  case eYZtoZY:
230  {
231  TransposeYZtoZY(inarray, outarray, UseNumMode);
232  }
233  break;
234  case eZYtoYZ:
235  {
236  TransposeZYtoYZ(inarray, outarray, UseNumMode);
237  }
238  break;
239  case eXtoY:
240  {
241  ASSERTL0(false, "Transposition not implemented yet.");
242  }
243  break;
244  case eYtoZ:
245  {
246  ASSERTL0(false, "Transposition not implemented yet.");
247  }
248  break;
249  case eZtoX:
250  {
251  ASSERTL0(false, "Transposition not implemented yet.");
252  }
253  break;
254  default:
255  {
256  ASSERTL0(false, "Transposition type does not exist.");
257  }
258  }
259 }
260 
261 /**
262  * Homogeneous 1D transposition from SEM to Homogeneous ordering.
263  */
265  Array<OneD, NekDouble> &outarray,
266  bool UseNumMode)
267 {
268  if (m_num_processes[0] > 1)
269  {
270  // Paramerers set up
271  int i, packed_len;
272  int copy_len = 0;
273  int index = 0;
274  int cnt = 0;
275 
276  int num_dofs = inarray.num_elements();
277  int num_points_per_plane = num_dofs / m_num_points_per_proc[0];
278  int num_pencil_per_proc =
279  (num_points_per_plane / m_num_processes[0]) +
280  (num_points_per_plane % m_num_processes[0] > 0);
281 
283  m_OffsetMap = Array<OneD, int>(m_num_processes[0], 0);
284 
285  for (i = 0; i < m_num_processes[0]; i++)
286  {
287  m_SizeMap[i] = num_pencil_per_proc * m_num_points_per_proc[0];
288  m_OffsetMap[i] = i * num_pencil_per_proc * m_num_points_per_proc[0];
289  }
290 
291  Array<OneD, NekDouble> tmp_outarray(
292  num_pencil_per_proc * m_num_homogeneous_points[0], 0.0);
293 
294  if (UseNumMode)
295  {
296  packed_len = m_num_homogeneous_coeffs[0];
297  }
298  else
299  {
300  packed_len = m_num_homogeneous_points[0];
301  }
302 
303  // Start Transposition
304  while (index < num_points_per_plane)
305  {
306  copy_len = num_pencil_per_proc < (num_points_per_plane - index)
307  ? num_pencil_per_proc
308  : (num_points_per_plane - index);
309 
310  for (i = 0; i < m_num_points_per_proc[0]; i++)
311  {
312  Vmath::Vcopy(copy_len,
313  &(inarray[index + (i * num_points_per_plane)]), 1,
314  &(outarray[cnt]), 1);
315 
316  cnt += num_pencil_per_proc;
317  }
318 
319  index += copy_len;
320  }
321 
322  m_hcomm->AlltoAllv(outarray, m_SizeMap, m_OffsetMap, tmp_outarray,
323  m_SizeMap, m_OffsetMap);
324 
325  for (i = 0; i < packed_len; ++i)
326  {
327  Vmath::Vcopy(num_pencil_per_proc,
328  &(tmp_outarray[i * num_pencil_per_proc]), 1,
329  &(outarray[i]), packed_len);
330  }
331  // End Transposition
332  }
333 
334  // Serial case implementation (more efficient then MPI 1 processor
335  // implemenation)
336  else
337  {
338  int i, pts_per_plane;
339  int n = inarray.num_elements();
340  int packed_len;
341 
342  pts_per_plane = n / m_num_points_per_proc[0];
343 
344  if (UseNumMode)
345  {
346  packed_len = m_num_homogeneous_coeffs[0];
347  }
348  else
349  {
350  packed_len = m_num_homogeneous_points[0];
351  }
352 
353  ASSERTL1(&inarray[0] != &outarray[0],
354  "Inarray and outarray cannot be the same");
355 
356  for (i = 0; i < packed_len; ++i)
357  {
358  Vmath::Vcopy(pts_per_plane, &(inarray[i * pts_per_plane]), 1,
359  &(outarray[i]), packed_len);
360  }
361  }
362 }
363 
364 /**
365  * Homogeneous 1D transposition from Homogeneous to SEM ordering.
366  */
368  Array<OneD, NekDouble> &outarray,
369  bool UseNumMode)
370 {
371  if (m_num_processes[0] > 1)
372  {
373  // Paramerers set up
374  int i, packed_len;
375  int copy_len = 0;
376  int index = 0;
377  int cnt = 0;
378 
379  int num_dofs = outarray.num_elements();
380  int num_points_per_plane = num_dofs / m_num_points_per_proc[0];
381  int num_pencil_per_proc =
382  (num_points_per_plane / m_num_processes[0]) +
383  (num_points_per_plane % m_num_processes[0] > 0);
384 
386  m_OffsetMap = Array<OneD, int>(m_num_processes[0], 0);
387 
388  for (i = 0; i < m_num_processes[0]; i++)
389  {
390  m_SizeMap[i] = num_pencil_per_proc * m_num_points_per_proc[0];
391  m_OffsetMap[i] = i * num_pencil_per_proc * m_num_points_per_proc[0];
392  }
393 
394  Array<OneD, NekDouble> tmp_inarray(
395  num_pencil_per_proc * m_num_homogeneous_points[0], 0.0);
396  Array<OneD, NekDouble> tmp_outarray(
397  num_pencil_per_proc * m_num_homogeneous_points[0], 0.0);
398 
399  if (UseNumMode)
400  {
401  packed_len = m_num_homogeneous_coeffs[0];
402  }
403  else
404  {
405  packed_len = m_num_homogeneous_points[0];
406  }
407 
408  // Start Transposition
409  for (i = 0; i < packed_len; ++i)
410  {
411  Vmath::Vcopy(num_pencil_per_proc, &(inarray[i]), packed_len,
412  &(tmp_inarray[i * num_pencil_per_proc]), 1);
413  }
414 
415  m_hcomm->AlltoAllv(tmp_inarray, m_SizeMap, m_OffsetMap, tmp_outarray,
416  m_SizeMap, m_OffsetMap);
417 
418  while (index < num_points_per_plane)
419  {
420  copy_len = num_pencil_per_proc < (num_points_per_plane - index)
421  ? num_pencil_per_proc
422  : (num_points_per_plane - index);
423 
424  for (i = 0; i < m_num_points_per_proc[0]; i++)
425  {
426  Vmath::Vcopy(copy_len, &(tmp_outarray[cnt]), 1,
427  &(outarray[index + (i * num_points_per_plane)]),
428  1);
429 
430  cnt += num_pencil_per_proc;
431  }
432 
433  index += copy_len;
434  }
435  // End Transposition
436  }
437 
438  // Serial case implementation (more efficient then MPI 1 processor
439  // implemenation)
440  else
441  {
442  int i, pts_per_plane;
443  int n = inarray.num_elements();
444  int packed_len;
445 
446  // use length of inarray to determine data storage type
447  // (i.e.modal or physical).
448  pts_per_plane = n / m_num_points_per_proc[0];
449 
450  if (UseNumMode)
451  {
452  packed_len = m_num_homogeneous_coeffs[0];
453  }
454  else
455  {
456  packed_len = m_num_homogeneous_points[0];
457  }
458 
459  ASSERTL1(&inarray[0] != &outarray[0],
460  "Inarray and outarray cannot be the same");
461 
462  for (i = 0; i < packed_len; ++i)
463  {
464  Vmath::Vcopy(pts_per_plane, &(inarray[i]), packed_len,
465  &(outarray[i * pts_per_plane]), 1);
466  }
467  }
468 }
469 
470 /**
471  * Homogeneous 2D transposition from SEM to Homogeneous(YZ) ordering.
472  */
474  Array<OneD, NekDouble> &outarray,
475  bool UseNumMode)
476 {
477  if (m_num_processes[0] > 1 || m_num_processes[1] > 1)
478  {
479  ASSERTL0(false, "Parallel transposition not implemented yet for "
480  "3D-Homo-2D approach.");
481  }
482  else
483  {
484  int i, pts_per_line;
485  int n = inarray.num_elements();
486  int packed_len;
487 
488  pts_per_line =
490 
491  if (UseNumMode)
492  {
493  packed_len =
495  }
496  else
497  {
498  packed_len =
500  }
501 
502  ASSERTL1(&inarray[0] != &outarray[0],
503  "Inarray and outarray cannot be the same");
504 
505  for (i = 0; i < packed_len; ++i)
506  {
507  Vmath::Vcopy(pts_per_line, &(inarray[i * pts_per_line]), 1,
508  &(outarray[i]), packed_len);
509  }
510  }
511 }
512 
513 /**
514  * Homogeneous 2D transposition from Homogeneous (YZ) ordering to SEM.
515  */
517  Array<OneD, NekDouble> &outarray,
518  bool UseNumMode)
519 {
520  if (m_num_processes[0] > 1 || m_num_processes[1] > 1)
521  {
522  ASSERTL0(false, "Parallel transposition not implemented yet for "
523  "3D-Homo-2D approach.");
524  }
525  else
526  {
527  int i, pts_per_line;
528  int n = inarray.num_elements();
529  int packed_len;
530 
531  pts_per_line =
533 
534  if (UseNumMode)
535  {
536  packed_len =
538  }
539  else
540  {
541  packed_len =
543  }
544 
545  ASSERTL1(&inarray[0] != &outarray[0],
546  "Inarray and outarray cannot be the same");
547 
548  for (i = 0; i < packed_len; ++i)
549  {
550  Vmath::Vcopy(pts_per_line, &(inarray[i]), packed_len,
551  &(outarray[i * pts_per_line]), 1);
552  }
553  }
554 }
555 
556 /**
557  * Homogeneous 2D transposition from Y ordering to Z.
558  */
560  Array<OneD, NekDouble> &outarray,
561  bool UseNumMode)
562 {
563  if (m_num_processes[0] > 1 || m_num_processes[1] > 1)
564  {
565  ASSERTL0(false, "Parallel transposition not implemented yet for "
566  "3D-Homo-2D approach.");
567  }
568  else
569  {
571  int s = inarray.num_elements();
572 
573  int pts_per_line = s / n;
574 
575  int packed_len = pts_per_line * m_num_homogeneous_points[1];
576 
577  for (int i = 0; i < m_num_homogeneous_points[0]; ++i)
578  {
579  Vmath::Vcopy(packed_len, &(inarray[i]), m_num_homogeneous_points[0],
580  &(outarray[i * packed_len]), 1);
581  }
582  }
583 }
584 
585 /**
586  * Homogeneous 2D transposition from Z ordering to Y.
587  */
589  Array<OneD, NekDouble> &outarray,
590  bool UseNumMode)
591 {
592  if (m_num_processes[0] > 1 || m_num_processes[1] > 1)
593  {
594  ASSERTL0(false, "Parallel transposition not implemented yet for "
595  "3D-Homo-2D approach.");
596  }
597  else
598  {
600  int s = inarray.num_elements();
601 
602  int pts_per_line = s / n;
603 
604  int packed_len = pts_per_line * m_num_homogeneous_points[1];
605 
606  for (int i = 0; i < packed_len; ++i)
607  {
608  Vmath::Vcopy(m_num_homogeneous_points[0], &(inarray[i]), packed_len,
609  &(outarray[i * m_num_homogeneous_points[0]]), 1);
610  }
611  }
612 }
613 
614 // TODO: Impelement 2D and 3D transposition routines
615 }
616 }
Array< OneD, int > m_num_homogeneous_points
Total homogeneous points per direction.
void TransposeZYtoYZ(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool UseNumMode=false)
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
Array< OneD, unsigned int > GetPlanesIDs(void)
Array< OneD, int > m_num_points_per_proc
Number of homogeneous points on each processor per direction.
void TransposeXYtoZ(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool UseNumMode=false)
BasisType GetBasisType() const
Return type of expansion basis.
Definition: Basis.h:139
void TransposeYZtoX(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool UseNumMode=false)
Transposition(const LibUtilities::BasisKey &HomoBasis0, LibUtilities::CommSharedPtr hcomm0, LibUtilities::CommSharedPtr hcomm1)
void TransposeZtoXY(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool UseNumMode=false)
void TransposeXtoYZ(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool UseNumMode=false)
Fourier Expansion .
Definition: BasisType.h:52
Array< OneD, unsigned int > m_planes_IDs
IDs of the planes on the processes.
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:55
void TransposeYZtoZY(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool UseNumMode=false)
Fourier Modified expansions with just the real part of the first mode .
Definition: BasisType.h:59
int GetNumPoints() const
Return points order at which basis is defined.
Definition: Basis.h:128
Array< OneD, int > m_SizeMap
MPI_Alltoallv map containing size of send/recv buffer.
unsigned int m_strip_ID
IDs of the strips on the processes.
Fourier Modified expansions with just the imaginary part of the first mode .
Definition: BasisType.h:60
Array< OneD, unsigned int > m_K
Fourier wave numbers associated with the planes.
Fourier ModifiedExpansion with just the first mode .
Definition: BasisType.h:58
void Transpose(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool UseNumMode=false, TranspositionDir dir=eNoTrans)
Array< OneD, int > m_OffsetMap
MPI_Alltoallv offset map of send/recv buffer in global vector.
Array< OneD, unsigned int > GetKs(void)
Array< OneD, int > m_num_homogeneous_coeffs
Total number of homogeneous coefficients.
int GetNumModes() const
Returns the order of the basis.
Definition: Basis.h:84
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1061
Describes the specification for a Basis.
Definition: Basis.h:50