source: XIOS/trunk/src/transformation/grid_transformation.cpp @ 795

Last change on this file since 795 was 795, checked in by mhnguyen, 7 years ago

Fixng a minor bug

+) Remove function that shadows the one of ancestor class

Test
+) On Curie
+) test_client, test_complete are correct

File size: 28.6 KB
Line 
1/*!
2   \file grid_transformation.cpp
3   \author Ha NGUYEN
4   \since 14 May 2015
5   \date 02 Jul 2015
6
7   \brief Interface for all transformations.
8 */
9#include "grid_transformation.hpp"
10#include "axis_algorithm_inverse.hpp"
11#include "axis_algorithm_zoom.hpp"
12#include "axis_algorithm_interpolate.hpp"
13#include "domain_algorithm_zoom.hpp"
14#include "domain_algorithm_interpolate.hpp"
15#include "context.hpp"
16#include "context_client.hpp"
17#include "transformation_mapping.hpp"
18#include "axis_algorithm_transformation.hpp"
19#include "distribution_client.hpp"
20
21namespace xios {
22CGridTransformation::CGridTransformation(CGrid* destination, CGrid* source)
23: gridSource_(source), gridDestination_(destination), originalGridSource_(source),
24  globalIndexOfCurrentGridSource_(), globalIndexOfOriginalGridSource_(), weightOfGlobalIndexOfOriginalGridSource_(0), algoTypes_()
25{
26  //Verify the compatibity between two grids
27  int numElement = gridDestination_->axis_domain_order.numElements();
28  if (numElement != gridSource_->axis_domain_order.numElements())
29    ERROR("CGridTransformation::CGridTransformation(CGrid* destination, CGrid* source)",
30       << "Two grids have different number of elements"
31       << "Number of elements of grid source " <<gridSource_->getId() << " is " << gridSource_->axis_domain_order.numElements()  << std::endl
32       << "Number of elements of grid destination " <<gridDestination_->getId() << " is " << numElement);
33
34  for (int i = 0; i < numElement; ++i)
35  {
36    if (gridDestination_->axis_domain_order(i) != gridSource_->axis_domain_order(i))
37      ERROR("CGridTransformation::CGridTransformation(CGrid* destination, CGrid* source)",
38         << "Transformed grid and its grid source have incompatible elements"
39         << "Grid source " <<gridSource_->getId() << std::endl
40         << "Grid destination " <<gridDestination_->getId());
41  }
42
43  std::vector<CAxis*> axisSrcTmp = gridSource_->getAxis(), axisSrc;
44  std::vector<CDomain*> domainSrcTmp = gridSource_->getDomains(), domainSrc;
45  for (int idx = 0; idx < axisSrcTmp.size(); ++idx)
46  {
47    CAxis* axis = CAxis::createAxis();
48    axis->axis_ref.setValue(axisSrcTmp[idx]->getId());
49    axis->solveRefInheritance(true);
50    axis->solveInheritanceTransformation();
51    axis->checkAttributesOnClient();
52    axisSrc.push_back(axis);
53  }
54
55  for (int idx = 0; idx < domainSrcTmp.size(); ++idx)
56  {
57    CDomain* domain = CDomain::createDomain();
58    domain->domain_ref.setValue(domainSrcTmp[idx]->getId());
59    domain->solveRefInheritance(true);
60    domain->solveInheritanceTransformation();
61    domain->checkAttributesOnClient();
62    domainSrc.push_back(domain);
63  }
64
65  gridSource_ = CGrid::createGrid(domainSrc, axisSrc, gridDestination_->axis_domain_order);
66  gridSource_->computeGridGlobalDimension(domainSrc, axisSrc, gridDestination_->axis_domain_order);
67
68  initializeMappingOfOriginalGridSource();
69  initializeAlgorithms();
70}
71
72/*!
73  Initialize the mapping between the first grid source and the original one
74  In a series of transformation, for each step, there is a need to "create" a new grid that plays a role of "temporary" source.
75Because at the end of the series, we need to know about the index mapping between the final grid destination and original grid source,
76for each transformation, we need to make sure that the current "temporary source" maps its global index correctly to the original one.
77*/
78void CGridTransformation::initializeMappingOfOriginalGridSource()
79{
80  CContext* context = CContext::getCurrent();
81  CContextClient* client = context->client;
82
83  CDistributionClient distribution(client->clientRank, originalGridSource_);
84  const std::vector<size_t>& globalIndexGridSrcSendToServer = distribution.getGlobalDataIndexSendToServer();
85
86  weightOfGlobalIndexOfOriginalGridSource_.resize(globalIndexGridSrcSendToServer.size());
87  globalIndexOfCurrentGridSource_  = globalIndexGridSrcSendToServer;
88  globalIndexOfOriginalGridSource_ = globalIndexGridSrcSendToServer;
89  weightOfGlobalIndexOfOriginalGridSource_ = 1.0;
90}
91
92CGridTransformation::~CGridTransformation()
93{
94  std::list<CGenericAlgorithmTransformation*>::const_iterator itb = algoTransformation_.begin(), it,
95                                                              ite = algoTransformation_.end();
96  for (it = itb; it != ite; ++it) delete (*it);
97}
98
99/*!
100  Initialize the algorithms (transformations)
101*/
102void CGridTransformation::initializeAlgorithms()
103{
104  std::vector<int> axisPositionInGrid;
105  std::vector<int> domPositionInGrid;
106  std::vector<CAxis*> axisListDestP = gridDestination_->getAxis();
107  std::vector<CDomain*> domListDestP = gridDestination_->getDomains();
108
109  int idx = 0;
110  for (int i = 0; i < gridDestination_->axis_domain_order.numElements(); ++i)
111  {
112    if (false == (gridDestination_->axis_domain_order)(i))
113    {
114      axisPositionInGrid.push_back(idx);
115      ++idx;
116    }
117    else
118    {
119      ++idx;
120      domPositionInGrid.push_back(idx);
121      ++idx;
122    }
123  }
124
125  for (int i = 0; i < axisListDestP.size(); ++i)
126  {
127    elementPosition2AxisPositionInGrid_[axisPositionInGrid[i]] = i;
128  }
129
130  for (int i = 0; i < domListDestP.size(); ++i)
131  {
132    elementPosition2DomainPositionInGrid_[domPositionInGrid[i]] = i;
133  }
134
135  idx = 0;
136  for (int i = 0; i < gridDestination_->axis_domain_order.numElements(); ++i)
137  {
138    if (false == (gridDestination_->axis_domain_order)(i))
139    {
140      initializeAxisAlgorithms(idx);
141      ++idx;
142    }
143    else
144    {
145      ++idx;
146      initializeDomainAlgorithms(idx);
147      ++idx;
148    }
149  }
150}
151
152
153
154/*!
155  Initialize the algorithms corresponding to transformation info contained in each axis.
156If an axis has transformations, these transformations will be represented in form of vector of CTransformation pointers
157In general, each axis can have several transformations performed on itself. However, should they be done seperately or combinely (of course in order)?
158For now, one approach is to do these combinely but maybe this needs changing.
159\param [in] axisPositionInGrid position of an axis in grid. (for example: a grid with one domain and one axis, position of domain is 1, position of axis is 2)
160*/
161void CGridTransformation::initializeAxisAlgorithms(int axisPositionInGrid)
162{
163  std::vector<CAxis*> axisListDestP = gridDestination_->getAxis();
164  if (!axisListDestP.empty())
165  {
166    if (axisListDestP[elementPosition2AxisPositionInGrid_[axisPositionInGrid]]->hasTransformation())
167    {
168      CAxis::TransMapTypes trans = axisListDestP[elementPosition2AxisPositionInGrid_[axisPositionInGrid]]->getAllTransformations();
169      CAxis::TransMapTypes::const_iterator itb = trans.begin(), it,
170                                           ite = trans.end();
171      int transformationOrder = 0;
172      for (it = itb; it != ite; ++it)
173      {
174        listAlgos_.push_back(std::make_pair(axisPositionInGrid, std::make_pair(it->first, transformationOrder)));
175        algoTypes_.push_back(false);
176        ++transformationOrder;
177      }
178    }
179  }
180}
181
182/*!
183  Initialize the algorithms corresponding to transformation info contained in each domain.
184If a domain has transformations, they will be represented in form of vector of CTransformation pointers
185In general, each domain can have several transformations performed on itself.
186\param [in] domPositionInGrid position of a domain in grid. (for example: a grid with one domain and one axis, position of domain is 1, position of axis is 2)
187*/
188void CGridTransformation::initializeDomainAlgorithms(int domPositionInGrid)
189{
190  std::vector<CDomain*> domListDestP = gridDestination_->getDomains();
191  if (!domListDestP.empty())
192  {
193    if (domListDestP[elementPosition2DomainPositionInGrid_[domPositionInGrid]]->hasTransformation())
194    {
195      CDomain::TransMapTypes trans = domListDestP[elementPosition2DomainPositionInGrid_[domPositionInGrid]]->getAllTransformations();
196      CDomain::TransMapTypes::const_iterator itb = trans.begin(), it,
197                                             ite = trans.end();
198      int transformationOrder = 0;
199      for (it = itb; it != ite; ++it)
200      {
201        listAlgos_.push_back(std::make_pair(domPositionInGrid, std::make_pair(it->first, transformationOrder)));
202        algoTypes_.push_back(true);
203        ++transformationOrder;
204      }
205    }
206  }
207
208}
209
210/*!
211  Select algorithm correspoding to its transformation type and its position in each element
212  \param [in] elementPositionInGrid position of element in grid. e.g: a grid has 1 domain and 1 axis, then position of domain is 1 (because it contains 2 basic elements)
213                                             and position of axis is 2
214  \param [in] transType transformation type, for now we have Zoom_axis, inverse_axis
215  \param [in] transformationOrder position of the transformation in an element (an element can have several transformation)
216  \param [in] isDomainAlgo flag to specify type of algorithm (for domain or axis)
217*/
218void CGridTransformation::selectAlgo(int elementPositionInGrid, ETranformationType transType, int transformationOrder, bool isDomainAlgo)
219{
220   if (isDomainAlgo) selectDomainAlgo(elementPositionInGrid, transType, transformationOrder);
221   else selectAxisAlgo(elementPositionInGrid, transType, transformationOrder);
222}
223
224/*!
225  Select algorithm of an axis correspoding to its transformation type and its position in each element
226  \param [in] elementPositionInGrid position of element in grid. e.g: a grid has 1 domain and 1 axis, then position of domain is 1 (because it contains 2 basic elements)
227                                             and position of axis is 2
228  \param [in] transType transformation type, for now we have Zoom_axis, inverse_axis
229  \param [in] transformationOrder position of the transformation in an element (an element can have several transformation)
230*/
231void CGridTransformation::selectAxisAlgo(int elementPositionInGrid, ETranformationType transType, int transformationOrder)
232{
233  std::vector<CAxis*> axisListDestP = gridDestination_->getAxis();
234  std::vector<CAxis*> axisListSrcP = gridSource_->getAxis();
235
236  int axisIndex =  elementPosition2AxisPositionInGrid_[elementPositionInGrid];
237  CAxis::TransMapTypes trans = axisListDestP[axisIndex]->getAllTransformations();
238  CAxis::TransMapTypes::const_iterator it = trans.begin();
239
240  for (int i = 0; i < transformationOrder; ++i, ++it) {}  // Find the correct transformation
241
242  CZoomAxis* zoomAxis = 0;
243  CInterpolateAxis* interpAxis = 0;
244  CGenericAlgorithmTransformation* algo = 0;
245  switch (transType)
246  {
247    case TRANS_INTERPOLATE_AXIS:
248      interpAxis = dynamic_cast<CInterpolateAxis*> (it->second);
249      algo = new CAxisAlgorithmInterpolate(axisListDestP[axisIndex], axisListSrcP[axisIndex], interpAxis);
250      break;
251    case TRANS_ZOOM_AXIS:
252      zoomAxis = dynamic_cast<CZoomAxis*> (it->second);
253      algo = new CAxisAlgorithmZoom(axisListDestP[axisIndex], axisListSrcP[axisIndex], zoomAxis);
254      break;
255    case TRANS_INVERSE_AXIS:
256      algo = new CAxisAlgorithmInverse(axisListDestP[axisIndex], axisListSrcP[axisIndex]);
257      break;
258    default:
259      break;
260  }
261  algoTransformation_.push_back(algo);
262
263}
264
265/*!
266  Select algorithm of a domain correspoding to its transformation type and its position in each element
267  \param [in] elementPositionInGrid position of element in grid. e.g: a grid has 1 domain and 1 axis, then position of domain is 1 (because it contains 2 basic elements)
268                                             and position of axis is 2
269  \param [in] transType transformation type, for now we have Zoom_axis, inverse_axis
270  \param [in] transformationOrder position of the transformation in an element (an element can have several transformation)
271*/
272void CGridTransformation::selectDomainAlgo(int elementPositionInGrid, ETranformationType transType, int transformationOrder)
273{
274  std::vector<CDomain*> domainListDestP = gridDestination_->getDomains();
275  std::vector<CDomain*> domainListSrcP = gridSource_->getDomains();
276
277  int domainIndex =  elementPosition2DomainPositionInGrid_[elementPositionInGrid];
278  CDomain::TransMapTypes trans = domainListDestP[domainIndex]->getAllTransformations();
279  CDomain::TransMapTypes::const_iterator it = trans.begin();
280
281  for (int i = 0; i < transformationOrder; ++i, ++it) {}  // Find the correct transformation
282
283  CZoomDomain* zoomDomain = 0;
284  CInterpolateDomain* interpFileDomain = 0;
285  CGenericAlgorithmTransformation* algo = 0;
286  switch (transType)
287  {
288    case TRANS_INTERPOLATE_DOMAIN:
289      interpFileDomain = dynamic_cast<CInterpolateDomain*> (it->second);
290      algo = new CDomainAlgorithmInterpolate(domainListDestP[domainIndex], domainListSrcP[domainIndex],interpFileDomain);
291      break;
292    case TRANS_ZOOM_DOMAIN:
293      zoomDomain = dynamic_cast<CZoomDomain*> (it->second);
294      algo = new CDomainAlgorithmZoom(domainListDestP[domainIndex], domainListSrcP[domainIndex], zoomDomain);
295      break;
296    default:
297      break;
298  }
299  algoTransformation_.push_back(algo);
300}
301
302/*!
303  Assign the current grid destination to the grid source in the new transformation.
304The current grid destination plays the role of grid source in next transformation (if any).
305Only element on which the transformation is performed is modified
306  \param [in] elementPositionInGrid position of element in grid
307  \param [in] transType transformation type
308*/
309void CGridTransformation::setUpGrid(int elementPositionInGrid, ETranformationType transType)
310{
311  std::vector<CAxis*> axisListDestP = gridDestination_->getAxis();
312  std::vector<CAxis*> axisListSrcP = gridSource_->getAxis();
313
314  std::vector<CDomain*> domListDestP = gridDestination_->getDomains();
315  std::vector<CDomain*> domListSrcP = gridSource_->getDomains();
316
317  int axisIndex, domainIndex;
318  switch (transType)
319  {
320    case TRANS_INTERPOLATE_DOMAIN:
321    case TRANS_ZOOM_DOMAIN:
322      domainIndex = elementPosition2DomainPositionInGrid_[elementPositionInGrid];
323      domListSrcP[domainIndex]->clearAllAttributes();
324      domListSrcP[domainIndex]->duplicateAttributes(domListDestP[domainIndex]);
325      break;
326
327    case TRANS_INTERPOLATE_AXIS:
328    case TRANS_ZOOM_AXIS:
329    case TRANS_INVERSE_AXIS:
330      axisIndex =  elementPosition2AxisPositionInGrid_[elementPositionInGrid];
331      axisListSrcP[axisIndex]->clearAllAttributes();
332      axisListSrcP[axisIndex]->duplicateAttributes(axisListDestP[axisIndex]);
333      break;
334    default:
335      break;
336  }
337}
338
339/*!
340  Perform all transformations
341  For each transformation, there are some things to do:
342  -) Chose the correct algorithm by transformation type and position of element
343  -) Calculate the mapping of global index between the current grid source and grid destination
344  -) Calculate the mapping of global index between current grid DESTINATION and ORIGINAL grid SOURCE
345  -) Make current grid destination become grid source in the next transformation
346*/
347void CGridTransformation::computeAll()
348{
349  CContext* context = CContext::getCurrent();
350  CContextClient* client = context->client;
351
352  ListAlgoType::const_iterator itb = listAlgos_.begin(),
353                               ite = listAlgos_.end(), it;
354  CGenericAlgorithmTransformation* algo = 0;
355  int nbAgloTransformation = 0; // Only count for executed transformation. Generate domain is a special one, not executed in the list
356  for (it = itb; it != ite; ++it)
357  {
358    int elementPositionInGrid = it->first;
359    ETranformationType transType = (it->second).first;
360    int transformationOrder = (it->second).second;
361    std::map<size_t, std::vector<std::pair<size_t,double> > > globaIndexWeightFromDestToSource;
362
363    // First of all, select an algorithm
364    selectAlgo(elementPositionInGrid, transType, transformationOrder, algoTypes_[std::distance(itb, it)]);
365    algo = algoTransformation_.back();
366
367    if (0 != algo) // Only registered transformation can be executed
368    {
369      // Recalculate the distribution of grid destination
370      CDistributionClient distributionClientDest(client->clientRank, gridDestination_);
371      const std::vector<size_t>& globalIndexGridDestSendToServer = distributionClientDest.getGlobalDataIndexSendToServer();
372
373      // ComputeTransformation of global index of each element
374      std::vector<int> gridDestinationDimensionSize = gridDestination_->getGlobalDimension();
375      std::vector<int> gridSrcDimensionSize = gridSource_->getGlobalDimension();
376      int elementPosition = it->first;
377      algo->computeGlobalSourceIndex(elementPosition,
378                                     gridDestinationDimensionSize,
379                                     gridSrcDimensionSize,
380                                     globalIndexGridDestSendToServer,
381                                     globaIndexWeightFromDestToSource);
382
383      // Compute transformation of global indexes among grids
384      computeTransformationFromOriginalGridSource(globaIndexWeightFromDestToSource);
385
386      // Now grid destination becomes grid source in a new transformation
387      setUpGrid(elementPositionInGrid, transType);
388      ++nbAgloTransformation;
389    }
390  }
391
392  if (0 != nbAgloTransformation)
393  {
394    updateFinalGridDestination();
395    computeFinalTransformationMapping();
396  }
397}
398
399
400/*!
401  After applying the algorithms, there are some informations on grid destination needing change, for now, there are:
402   +) mask
403*/
404void CGridTransformation::updateFinalGridDestination()
405{
406  CContext* context = CContext::getCurrent();
407  CContextClient* client = context->client;
408
409  //First of all, retrieve info of local mask of grid destination
410  CDistributionClient distributionClientDest(client->clientRank, gridDestination_);
411  const std::vector<int>& localMaskIndexOnClientDest = distributionClientDest.getLocalMaskIndexOnClient();
412  const std::vector<size_t>& globalIndexOnClientDest = distributionClientDest.getGlobalDataIndexSendToServer();
413
414  std::vector<size_t>::const_iterator itbArr, itArr, iteArr;
415  itbArr = globalIndexOnClientDest.begin();
416  iteArr = globalIndexOnClientDest.end();
417
418  // Then find out which index became invalid (become masked after being applied the algorithms, or demande some masked points from grid source)
419  int num = globalIndexOfOriginalGridSource_.size();
420  const size_t sfmax = NumTraits<unsigned long>::sfmax();
421  int maskIndexNum = 0;
422  for (int idx = 0; idx < num; ++idx)
423  {
424    if (sfmax == globalIndexOfOriginalGridSource_[idx])
425    {
426      size_t maskedGlobalIndex = globalIndexOfCurrentGridSource_[idx];
427      itArr = std::find(itbArr, iteArr, maskedGlobalIndex);
428      if (iteArr != itArr) ++maskIndexNum;
429    }
430  }
431
432  CArray<int,1> maskIndexToModify(maskIndexNum);
433  maskIndexNum = 0;
434  for (int idx = 0; idx < num; ++idx)
435  {
436    if (sfmax == globalIndexOfOriginalGridSource_[idx])
437    {
438      size_t maskedGlobalIndex = globalIndexOfCurrentGridSource_[idx];
439      itArr = std::find(itbArr, iteArr, maskedGlobalIndex);
440      if (iteArr != itArr)
441      {
442        int localIdx = std::distance(itbArr, itArr);
443        maskIndexToModify(maskIndexNum) = localMaskIndexOnClientDest[localIdx];
444        ++maskIndexNum;
445      }
446    }
447  }
448
449  gridDestination_->modifyMask(maskIndexToModify);
450}
451
452/*!
453  A transformation from a grid source to grid destination often passes several intermediate grids, which play a role of
454temporary grid source and/or grid destination. This function makes sure that global index of original grid source are mapped correctly to
455the final grid destination
456*/
457void CGridTransformation::computeTransformationFromOriginalGridSource(const std::map<size_t, std::vector<std::pair<size_t,double> > >& globaIndexMapFromDestToSource)
458{
459  CContext* context = CContext::getCurrent();
460  CContextClient* client = context->client;
461
462  CTransformationMapping transformationMap(gridDestination_, gridSource_);
463
464    // Then compute transformation mapping among clients
465  transformationMap.computeTransformationMapping(globaIndexMapFromDestToSource);
466
467  const std::map<int,std::vector<std::vector<std::pair<size_t,double> > > >& globalIndexToReceive = transformationMap.getGlobalIndexReceivedOnGridDestMapping();
468  const std::map<int,std::vector<size_t> >& globalIndexToSend = transformationMap.getGlobalIndexSendToGridDestMapping();
469
470 // Sending global index of original grid source
471  std::map<int,std::vector<size_t> >::const_iterator itbSend = globalIndexToSend.begin(), itSend,
472                                                     iteSend = globalIndexToSend.end();
473  std::vector<size_t>::const_iterator itbArr = globalIndexOfCurrentGridSource_.begin(), itArr,
474                                      iteArr = globalIndexOfCurrentGridSource_.end();
475 int sendBuffSize = 0;
476 for (itSend = itbSend; itSend != iteSend; ++itSend) sendBuffSize += (itSend->second).size();
477
478 typedef unsigned long Scalar;
479 unsigned long* sendBuff, *currentSendBuff;
480 if (0 != sendBuffSize) sendBuff = new unsigned long [sendBuffSize];
481 for (StdSize idx = 0; idx < sendBuffSize; ++idx) sendBuff[idx] = NumTraits<Scalar>::sfmax();
482
483 std::map<int, MPI_Request> requests;
484
485 std::vector<int> permutIndex(globalIndexOfCurrentGridSource_.size());
486 typedef XIOSBinarySearchWithIndex<size_t> BinarySearch;
487 XIOSAlgorithms::fillInIndex(globalIndexOfCurrentGridSource_.size(), permutIndex);
488 XIOSAlgorithms::sortWithIndex<size_t, CVectorStorage>(globalIndexOfCurrentGridSource_, permutIndex);
489 BinarySearch searchCurrentSrc(globalIndexOfCurrentGridSource_);
490 std::vector<int>::iterator itbIndex = permutIndex.begin(), itIndex,
491                            iteIndex = permutIndex.end();
492
493  // Find out local index on grid destination (received)
494 int currentBuffPosition = 0;
495 for (itSend = itbSend; itSend != iteSend; ++itSend)
496 {
497   int destRank = itSend->first;
498   const std::vector<size_t>& globalIndexOfCurrentGridSourceToSend = itSend->second;
499   int countSize = globalIndexOfCurrentGridSourceToSend.size();
500   for (int idx = 0; idx < (countSize); ++idx)
501   {
502     if (searchCurrentSrc.search(itbIndex, iteIndex, globalIndexOfCurrentGridSourceToSend[idx], itIndex))
503     {
504       sendBuff[idx+currentBuffPosition] = globalIndexOfOriginalGridSource_[*itIndex];
505     }
506   }
507   currentSendBuff = sendBuff + currentBuffPosition;
508   MPI_Isend(currentSendBuff, countSize, MPI_UNSIGNED_LONG, destRank, 14, client->intraComm, &requests[destRank]);
509   currentBuffPosition += countSize;
510 }
511
512 // Receiving global index of grid source sending from current grid source
513 std::map<int,std::vector<std::vector<std::pair<size_t,double> > > >::const_iterator itbRecv = globalIndexToReceive.begin(), itRecv,
514                                                                                     iteRecv = globalIndexToReceive.end();
515 int recvBuffSize = 0;
516 for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv) recvBuffSize += (itRecv->second).size();
517
518 unsigned long* recvBuff, *currentRecvBuff;
519 if (0 != recvBuffSize) recvBuff = new unsigned long [recvBuffSize];
520 for (StdSize idx = 0; idx < recvBuffSize; ++idx) recvBuff[idx] = NumTraits<Scalar>::sfmax();
521
522 int currentRecvBuffPosition = 0;
523 for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv)
524 {
525   MPI_Status status;
526   int srcRank = itRecv->first;
527   int countSize = (itRecv->second).size();
528   currentRecvBuff = recvBuff + currentRecvBuffPosition;
529   MPI_Recv(currentRecvBuff, countSize, MPI_UNSIGNED_LONG, srcRank, 14, client->intraComm, &status);
530   currentRecvBuffPosition += countSize;
531 }
532
533 int nbCurrentGridSource = 0;
534 for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv)
535 {
536   int ssize = (itRecv->second).size();
537   for (int idx = 0; idx < ssize; ++idx)
538   {
539     nbCurrentGridSource += (itRecv->second)[idx].size();
540   }
541 }
542
543 if (globalIndexOfCurrentGridSource_.size() != nbCurrentGridSource)
544 {
545   globalIndexOfCurrentGridSource_.resize(nbCurrentGridSource);
546   globalIndexOfOriginalGridSource_.resize(nbCurrentGridSource);
547   weightOfGlobalIndexOfOriginalGridSource_.resize(nbCurrentGridSource);
548 }
549
550 int k = 0;
551 currentRecvBuff = recvBuff;
552 for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv)
553 {
554   int countSize = (itRecv->second).size();
555   for (int idx = 0; idx < countSize; ++idx, ++currentRecvBuff)
556   {
557     int ssize = (itRecv->second)[idx].size();
558     for (int i = 0; i < ssize; ++i)
559     {
560       globalIndexOfCurrentGridSource_[k] = ((itRecv->second)[idx][i]).first;
561       weightOfGlobalIndexOfOriginalGridSource_(k) = ((itRecv->second)[idx][i]).second;
562       globalIndexOfOriginalGridSource_[k] = *currentRecvBuff;
563       ++k;
564     }
565   }
566 }
567
568 std::map<int, MPI_Request>::iterator itRequest;
569 for (itRequest = requests.begin(); itRequest != requests.end(); ++itRequest)
570   MPI_Wait(&itRequest->second, MPI_STATUS_IGNORE);
571
572 if (0 != sendBuffSize) delete [] sendBuff;
573 if (0 != recvBuffSize) delete [] recvBuff;
574}
575
576/*!
577  Compute transformation mapping between grid source and grid destination
578  The transformation between grid source and grid destination is represented in form of mapping between global index
579of two grids. Then local index mapping between data on each grid will be found out thanks to these global indexes
580*/
581void CGridTransformation::computeFinalTransformationMapping()
582{
583  CContext* context = CContext::getCurrent();
584  CContextClient* client = context->client;
585
586  CTransformationMapping transformationMap(gridDestination_, originalGridSource_);
587
588  std::map<size_t, std::vector<std::pair<size_t,double> > > globaIndexWeightFromDestToSource;
589  int nb = globalIndexOfCurrentGridSource_.size();
590  const size_t sfmax = NumTraits<unsigned long>::sfmax();
591  for (int idx = 0; idx < nb; ++idx)
592  {
593    if (sfmax != globalIndexOfOriginalGridSource_[idx])
594      globaIndexWeightFromDestToSource[globalIndexOfCurrentGridSource_[idx]].push_back(make_pair(globalIndexOfOriginalGridSource_[idx], weightOfGlobalIndexOfOriginalGridSource_(idx))) ;
595  }
596
597  // Then compute transformation mapping among clients
598  transformationMap.computeTransformationMapping(globaIndexWeightFromDestToSource);
599
600  const std::map<int,std::vector<std::vector<std::pair<size_t,double> > > >& globalIndexToReceive = transformationMap.getGlobalIndexReceivedOnGridDestMapping();
601  const std::map<int,std::vector<size_t> >& globalIndexToSend = transformationMap.getGlobalIndexSendToGridDestMapping();
602
603  CDistributionClient distributionClientDest(client->clientRank, gridDestination_);
604  CDistributionClient distributionClientSrc(client->clientRank, originalGridSource_);
605
606  const std::vector<size_t>& globalIndexOnClientDest = distributionClientDest.getGlobalDataIndexSendToServer();
607  const std::vector<size_t>& globalIndexOnClientSrc = distributionClientSrc.getGlobalDataIndexSendToServer();
608
609  std::vector<size_t>::const_iterator itbArr, itArr, iteArr;
610  std::vector<int>::const_iterator itIndex, itbIndex, iteIndex;
611  std::map<int,std::vector<std::vector<std::pair<size_t,double> > > >::const_iterator itbMapRecv, itMapRecv, iteMapRecv;
612
613  std::vector<int> permutIndex;
614  typedef XIOSBinarySearchWithIndex<size_t> BinarySearch;
615
616  // Find out local index on grid destination (received)
617  XIOSAlgorithms::fillInIndex(globalIndexOnClientDest.size(), permutIndex);
618  XIOSAlgorithms::sortWithIndex<size_t, CVectorStorage>(globalIndexOnClientDest, permutIndex);
619  itbIndex = permutIndex.begin();
620  iteIndex = permutIndex.end();
621  BinarySearch searchClientDest(globalIndexOnClientDest);
622  itbMapRecv = globalIndexToReceive.begin();
623  iteMapRecv = globalIndexToReceive.end();
624  for (itMapRecv = itbMapRecv; itMapRecv != iteMapRecv; ++itMapRecv)
625  {
626    int sourceRank = itMapRecv->first;
627    int numGlobalIndex = (itMapRecv->second).size();
628    for (int i = 0; i < numGlobalIndex; ++i)
629    {
630      int vecSize = ((itMapRecv->second)[i]).size();
631      std::vector<std::pair<int,double> > tmpVec;
632      for (int idx = 0; idx < vecSize; ++idx)
633      {
634        size_t globalIndex = (itMapRecv->second)[i][idx].first;
635        double weight = (itMapRecv->second)[i][idx].second;
636        if (searchClientDest.search(itbIndex, iteIndex, globalIndex, itIndex))
637        {
638          tmpVec.push_back(make_pair(*itIndex, weight));
639        }
640      }
641      localIndexToReceiveOnGridDest_[sourceRank].push_back(tmpVec);
642    }
643  }
644
645  // Find out local index on grid source (to send)
646  std::map<int,std::vector<size_t> >::const_iterator itbMap, itMap, iteMap;
647  XIOSAlgorithms::fillInIndex(globalIndexOnClientSrc.size(), permutIndex);
648  XIOSAlgorithms::sortWithIndex<size_t, CVectorStorage>(globalIndexOnClientSrc, permutIndex);
649  itbIndex = permutIndex.begin();
650  iteIndex = permutIndex.end();
651  BinarySearch searchClientSrc(globalIndexOnClientSrc);
652  itbMap = globalIndexToSend.begin();
653  iteMap = globalIndexToSend.end();
654  for (itMap = itbMap; itMap != iteMap; ++itMap)
655  {
656    int destRank = itMap->first;
657    int vecSize = itMap->second.size();
658    localIndexToSendFromGridSource_[destRank].resize(vecSize);
659    for (int idx = 0; idx < vecSize; ++idx)
660    {
661      if (searchClientSrc.search(itbIndex, iteIndex, itMap->second[idx], itIndex))
662      {
663        localIndexToSendFromGridSource_[destRank](idx) = *itIndex;
664      }
665    }
666  }
667}
668
669/*!
670  Local index of data which need sending from the grid source
671  \return local index of data
672*/
673const std::map<int, CArray<int,1> >& CGridTransformation::getLocalIndexToSendFromGridSource() const
674{
675  return localIndexToSendFromGridSource_;
676}
677
678/*!
679  Local index of data which will be received on the grid destination
680  \return local index of data
681*/
682const std::map<int,std::vector<std::vector<std::pair<int,double> > > >& CGridTransformation::getLocalIndexToReceiveOnGridDest() const
683{
684  return localIndexToReceiveOnGridDest_;
685}
686
687}
Note: See TracBrowser for help on using the repository browser.