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

Last change on this file since 821 was 821, checked in by mhnguyen, 6 years ago

Making some improvements of transformation algorithm

+) Correct the way to enlisting transformations in an element (domain, axis)
+) Optimize generic transformation to make sure temporary grid to be created on demand
+) Update some mpi tag to prevent conflict
+) Correct some minor stuffs
+) Update documents

Test
+) On Curie
+) all test pass

File size: 34.9 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#include "mpi_tag.hpp"
21
22namespace xios {
23CGridTransformation::CGridTransformation(CGrid* destination, CGrid* source)
24: gridSource_(source), gridDestination_(destination), originalGridSource_(source),
25  algoTypes_(), nbAlgos_(0), currentGridIndexToOriginalGridIndex_()
26
27{
28  //Verify the compatibity between two grids
29  int numElement = gridDestination_->axis_domain_order.numElements();
30  if (numElement != gridSource_->axis_domain_order.numElements())
31    ERROR("CGridTransformation::CGridTransformation(CGrid* destination, CGrid* source)",
32       << "Two grids have different number of elements"
33       << "Number of elements of grid source " <<gridSource_->getId() << " is " << gridSource_->axis_domain_order.numElements()  << std::endl
34       << "Number of elements of grid destination " <<gridDestination_->getId() << " is " << numElement);
35
36  for (int i = 0; i < numElement; ++i)
37  {
38    if (gridDestination_->axis_domain_order(i) != gridSource_->axis_domain_order(i))
39      ERROR("CGridTransformation::CGridTransformation(CGrid* destination, CGrid* source)",
40         << "Transformed grid and its grid source have incompatible elements"
41         << "Grid source " <<gridSource_->getId() << std::endl
42         << "Grid destination " <<gridDestination_->getId());
43  }
44
45//  std::vector<CAxis*> axisSrcTmp = gridSource_->getAxis(), axisSrc;
46//  std::vector<CDomain*> domainSrcTmp = gridSource_->getDomains(), domainSrc;
47//  for (int idx = 0; idx < axisSrcTmp.size(); ++idx)
48//  {
49//    CAxis* axis = CAxis::createAxis();
50//    axis->axis_ref.setValue(axisSrcTmp[idx]->getId());
51//    axis->solveRefInheritance(true);
52//    axis->solveInheritanceTransformation();
53//    axis->checkAttributesOnClient();
54//    axisSrc.push_back(axis);
55//  }
56//
57//  for (int idx = 0; idx < domainSrcTmp.size(); ++idx)
58//  {
59//    CDomain* domain = CDomain::createDomain();
60//    domain->domain_ref.setValue(domainSrcTmp[idx]->getId());
61//    domain->solveRefInheritance(true);
62//    domain->solveInheritanceTransformation();
63//    domain->checkAttributesOnClient();
64//    domainSrc.push_back(domain);
65//  }
66//
67//  gridSource_ = CGrid::createGrid(domainSrc, axisSrc, gridDestination_->axis_domain_order);
68//  gridSource_->computeGridGlobalDimension(domainSrc, axisSrc, gridDestination_->axis_domain_order);
69
70  initializeMappingOfOriginalGridSource();
71//  initializeAlgorithms();
72}
73
74/*!
75  Initialize the mapping between the first grid source and the original one
76  In a series of transformation, for each step, there is a need to "create" a new grid that plays a role of "temporary" source.
77Because at the end of the series, we need to know about the index mapping between the final grid destination and original grid source,
78for each transformation, we need to make sure that the current "temporary source" maps its global index correctly to the original one.
79*/
80void CGridTransformation::initializeMappingOfOriginalGridSource()
81{
82  CContext* context = CContext::getCurrent();
83  CContextClient* client = context->client;
84
85  // Initialize algorithms
86  initializeAlgorithms();
87
88  ListAlgoType::const_iterator itb = listAlgos_.begin(),
89                               ite = listAlgos_.end(), it;
90
91  for (it = itb; it != ite; ++it)
92  {
93    ETranformationType transType = (it->second).first;
94    if (!isSpecialTransformation(transType)) ++nbAlgos_;
95  }
96
97  if (1<nbAlgos_)  // Only when there are more than 1 algorithm, will we create temporary grid source
98  {
99    std::vector<CAxis*> axisSrcTmp = gridSource_->getAxis(), axisSrc;
100    std::vector<CDomain*> domainSrcTmp = gridSource_->getDomains(), domainSrc;
101    for (int idx = 0; idx < axisSrcTmp.size(); ++idx)
102    {
103      CAxis* axis = CAxis::createAxis();
104      axis->axis_ref.setValue(axisSrcTmp[idx]->getId());
105      axis->solveRefInheritance(true);
106      axis->solveInheritanceTransformation();
107      axis->checkAttributesOnClient();
108      axisSrc.push_back(axis);
109    }
110
111    for (int idx = 0; idx < domainSrcTmp.size(); ++idx)
112    {
113      CDomain* domain = CDomain::createDomain();
114      domain->domain_ref.setValue(domainSrcTmp[idx]->getId());
115      domain->solveRefInheritance(true);
116      domain->solveInheritanceTransformation();
117      domain->checkAttributesOnClient();
118      domainSrc.push_back(domain);
119    }
120
121    gridSource_ = CGrid::createGrid(domainSrc, axisSrc, gridDestination_->axis_domain_order);
122    gridSource_->computeGridGlobalDimension(domainSrc, axisSrc, gridDestination_->axis_domain_order);
123  }
124}
125
126CGridTransformation::~CGridTransformation()
127{
128  std::list<CGenericAlgorithmTransformation*>::const_iterator itb = algoTransformation_.begin(), it,
129                                                              ite = algoTransformation_.end();
130  for (it = itb; it != ite; ++it) delete (*it);
131}
132
133/*!
134  Initialize the algorithms (transformations)
135*/
136void CGridTransformation::initializeAlgorithms()
137{
138  std::vector<int> axisPositionInGrid;
139  std::vector<int> domPositionInGrid;
140  std::vector<CAxis*> axisListDestP = gridDestination_->getAxis();
141  std::vector<CDomain*> domListDestP = gridDestination_->getDomains();
142
143  int idx = 0;
144  for (int i = 0; i < gridDestination_->axis_domain_order.numElements(); ++i)
145  {
146    if (false == (gridDestination_->axis_domain_order)(i))
147    {
148      axisPositionInGrid.push_back(idx);
149      ++idx;
150    }
151    else
152    {
153      ++idx;
154      domPositionInGrid.push_back(idx);
155      ++idx;
156    }
157  }
158
159  for (int i = 0; i < axisListDestP.size(); ++i)
160  {
161    elementPosition2AxisPositionInGrid_[axisPositionInGrid[i]] = i;
162  }
163
164  for (int i = 0; i < domListDestP.size(); ++i)
165  {
166    elementPosition2DomainPositionInGrid_[domPositionInGrid[i]] = i;
167  }
168
169  idx = 0;
170  for (int i = 0; i < gridDestination_->axis_domain_order.numElements(); ++i)
171  {
172    if (false == (gridDestination_->axis_domain_order)(i))
173    {
174      initializeAxisAlgorithms(idx);
175      ++idx;
176    }
177    else
178    {
179      ++idx;
180      initializeDomainAlgorithms(idx);
181      ++idx;
182    }
183  }
184}
185
186/*!
187  Initialize the algorithms corresponding to transformation info contained in each axis.
188If an axis has transformations, these transformations will be represented in form of vector of CTransformation pointers
189In general, each axis can have several transformations performed on itself. However, should they be done seperately or combinely (of course in order)?
190For now, one approach is to do these combinely but maybe this needs changing.
191\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)
192*/
193void CGridTransformation::initializeAxisAlgorithms(int axisPositionInGrid)
194{
195  std::vector<CAxis*> axisListDestP = gridDestination_->getAxis();
196  if (!axisListDestP.empty())
197  {
198    if (axisListDestP[elementPosition2AxisPositionInGrid_[axisPositionInGrid]]->hasTransformation())
199    {
200      CAxis::TransMapTypes trans = axisListDestP[elementPosition2AxisPositionInGrid_[axisPositionInGrid]]->getAllTransformations();
201      CAxis::TransMapTypes::const_iterator itb = trans.begin(), it,
202                                           ite = trans.end();
203      int transformationOrder = 0;
204      for (it = itb; it != ite; ++it)
205      {
206        listAlgos_.push_back(std::make_pair(axisPositionInGrid, std::make_pair(it->first, transformationOrder)));
207        algoTypes_.push_back(false);
208        ++transformationOrder;
209      }
210    }
211  }
212}
213
214/*!
215  Initialize the algorithms corresponding to transformation info contained in each domain.
216If a domain has transformations, they will be represented in form of vector of CTransformation pointers
217In general, each domain can have several transformations performed on itself.
218\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)
219*/
220void CGridTransformation::initializeDomainAlgorithms(int domPositionInGrid)
221{
222  std::vector<CDomain*> domListDestP = gridDestination_->getDomains();
223  if (!domListDestP.empty())
224  {
225    if (domListDestP[elementPosition2DomainPositionInGrid_[domPositionInGrid]]->hasTransformation())
226    {
227      CDomain::TransMapTypes trans = domListDestP[elementPosition2DomainPositionInGrid_[domPositionInGrid]]->getAllTransformations();
228      CDomain::TransMapTypes::const_iterator itb = trans.begin(), it,
229                                             ite = trans.end();
230      int transformationOrder = 0;
231      for (it = itb; it != ite; ++it)
232      {
233        listAlgos_.push_back(std::make_pair(domPositionInGrid, std::make_pair(it->first, transformationOrder)));
234        algoTypes_.push_back(true);
235        ++transformationOrder;
236      }
237    }
238  }
239
240}
241
242/*!
243  Select algorithm correspoding to its transformation type and its position in each element
244  \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)
245                                             and position of axis is 2
246  \param [in] transType transformation type, for now we have Zoom_axis, inverse_axis
247  \param [in] transformationOrder position of the transformation in an element (an element can have several transformation)
248  \param [in] isDomainAlgo flag to specify type of algorithm (for domain or axis)
249*/
250void CGridTransformation::selectAlgo(int elementPositionInGrid, ETranformationType transType, int transformationOrder, bool isDomainAlgo)
251{
252   if (isDomainAlgo) selectDomainAlgo(elementPositionInGrid, transType, transformationOrder);
253   else selectAxisAlgo(elementPositionInGrid, transType, transformationOrder);
254}
255
256/*!
257  Select algorithm of an axis correspoding to its transformation type and its position in each element
258  \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)
259                                             and position of axis is 2
260  \param [in] transType transformation type, for now we have Zoom_axis, inverse_axis
261  \param [in] transformationOrder position of the transformation in an element (an element can have several transformation)
262*/
263void CGridTransformation::selectAxisAlgo(int elementPositionInGrid, ETranformationType transType, int transformationOrder)
264{
265  std::vector<CAxis*> axisListDestP = gridDestination_->getAxis();
266  std::vector<CAxis*> axisListSrcP = gridSource_->getAxis();
267
268  int axisIndex =  elementPosition2AxisPositionInGrid_[elementPositionInGrid];
269  CAxis::TransMapTypes trans = axisListDestP[axisIndex]->getAllTransformations();
270  CAxis::TransMapTypes::const_iterator it = trans.begin();
271
272  for (int i = 0; i < transformationOrder; ++i, ++it) {}  // Find the correct transformation
273
274  CZoomAxis* zoomAxis = 0;
275  CInterpolateAxis* interpAxis = 0;
276  CGenericAlgorithmTransformation* algo = 0;
277  switch (transType)
278  {
279    case TRANS_INTERPOLATE_AXIS:
280      interpAxis = dynamic_cast<CInterpolateAxis*> (it->second);
281      algo = new CAxisAlgorithmInterpolate(axisListDestP[axisIndex], axisListSrcP[axisIndex], interpAxis);
282      break;
283    case TRANS_ZOOM_AXIS:
284      zoomAxis = dynamic_cast<CZoomAxis*> (it->second);
285      algo = new CAxisAlgorithmZoom(axisListDestP[axisIndex], axisListSrcP[axisIndex], zoomAxis);
286      break;
287    case TRANS_INVERSE_AXIS:
288      algo = new CAxisAlgorithmInverse(axisListDestP[axisIndex], axisListSrcP[axisIndex]);
289      break;
290    default:
291      break;
292  }
293  algoTransformation_.push_back(algo);
294
295}
296
297/*!
298  Select algorithm of a domain correspoding to its transformation type and its position in each element
299  \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)
300                                             and position of axis is 2
301  \param [in] transType transformation type, for now we have Zoom_axis, inverse_axis
302  \param [in] transformationOrder position of the transformation in an element (an element can have several transformation)
303*/
304void CGridTransformation::selectDomainAlgo(int elementPositionInGrid, ETranformationType transType, int transformationOrder)
305{
306  std::vector<CDomain*> domainListDestP = gridDestination_->getDomains();
307  std::vector<CDomain*> domainListSrcP = gridSource_->getDomains();
308
309  int domainIndex =  elementPosition2DomainPositionInGrid_[elementPositionInGrid];
310  CDomain::TransMapTypes trans = domainListDestP[domainIndex]->getAllTransformations();
311  CDomain::TransMapTypes::const_iterator it = trans.begin();
312
313  for (int i = 0; i < transformationOrder; ++i, ++it) {}  // Find the correct transformation
314
315  CZoomDomain* zoomDomain = 0;
316  CInterpolateDomain* interpFileDomain = 0;
317  CGenericAlgorithmTransformation* algo = 0;
318  switch (transType)
319  {
320    case TRANS_INTERPOLATE_DOMAIN:
321      interpFileDomain = dynamic_cast<CInterpolateDomain*> (it->second);
322      algo = new CDomainAlgorithmInterpolate(domainListDestP[domainIndex], domainListSrcP[domainIndex],interpFileDomain);
323      break;
324    case TRANS_ZOOM_DOMAIN:
325      zoomDomain = dynamic_cast<CZoomDomain*> (it->second);
326      algo = new CDomainAlgorithmZoom(domainListDestP[domainIndex], domainListSrcP[domainIndex], zoomDomain);
327      break;
328    default:
329      break;
330  }
331  algoTransformation_.push_back(algo);
332}
333
334/*!
335  Assign the current grid destination to the grid source in the new transformation.
336The current grid destination plays the role of grid source in next transformation (if any).
337Only element on which the transformation is performed is modified
338  \param [in] elementPositionInGrid position of element in grid
339  \param [in] transType transformation type
340*/
341void CGridTransformation::setUpGrid(int elementPositionInGrid, ETranformationType transType)
342{
343  std::vector<CAxis*> axisListDestP = gridDestination_->getAxis();
344  std::vector<CAxis*> axisListSrcP = gridSource_->getAxis();
345
346  std::vector<CDomain*> domListDestP = gridDestination_->getDomains();
347  std::vector<CDomain*> domListSrcP = gridSource_->getDomains();
348
349  int axisIndex, domainIndex;
350  switch (transType)
351  {
352    case TRANS_INTERPOLATE_DOMAIN:
353    case TRANS_ZOOM_DOMAIN:
354      domainIndex = elementPosition2DomainPositionInGrid_[elementPositionInGrid];
355      domListSrcP[domainIndex]->clearAllAttributes();
356      domListSrcP[domainIndex]->duplicateAttributes(domListDestP[domainIndex]);
357      break;
358
359    case TRANS_INTERPOLATE_AXIS:
360    case TRANS_ZOOM_AXIS:
361    case TRANS_INVERSE_AXIS:
362      axisIndex =  elementPosition2AxisPositionInGrid_[elementPositionInGrid];
363      axisListSrcP[axisIndex]->clearAllAttributes();
364      axisListSrcP[axisIndex]->duplicateAttributes(axisListDestP[axisIndex]);
365      break;
366    default:
367      break;
368  }
369  gridSource_->createMask();
370  gridSource_->computeGridGlobalDimension(domListSrcP, axisListSrcP, gridSource_->axis_domain_order);
371}
372
373/*!
374  Perform all transformations
375  For each transformation, there are some things to do:
376  -) Chose the correct algorithm by transformation type and position of element
377  -) Calculate the mapping of global index between the current grid source and grid destination
378  -) Calculate the mapping of global index between current grid DESTINATION and ORIGINAL grid SOURCE
379  -) Make current grid destination become grid source in the next transformation
380*/
381void CGridTransformation::computeAll()
382{
383  CContext* context = CContext::getCurrent();
384  CContextClient* client = context->client;
385
386  ListAlgoType::const_iterator itb = listAlgos_.begin(),
387                               ite = listAlgos_.end(), it;
388
389  CGenericAlgorithmTransformation* algo = 0;
390  int nbAgloTransformation = 0; // Only count for executed transformation. Generate domain is a special one, not executed in the list
391  for (it = itb; it != ite; ++it)
392  {
393    int elementPositionInGrid = it->first;
394    ETranformationType transType = (it->second).first;
395    int transformationOrder = (it->second).second;
396    std::map<size_t, std::vector<std::pair<size_t,double> > > globaIndexWeightFromDestToSource;
397
398    // First of all, select an algorithm
399    selectAlgo(elementPositionInGrid, transType, transformationOrder, algoTypes_[std::distance(itb, it)]);
400    algo = algoTransformation_.back();
401
402    if (0 != algo) // Only registered transformation can be executed
403    {
404      // Recalculate the distribution of grid destination
405      CDistributionClient distributionClientDest(client->clientRank, gridDestination_);
406      const std::vector<size_t>& globalIndexGridDestSendToServer = distributionClientDest.getGlobalDataIndexSendToServer();
407
408      // ComputeTransformation of global index of each element
409      std::vector<int> gridDestinationDimensionSize = gridDestination_->getGlobalDimension();
410      std::vector<int> gridSrcDimensionSize = gridSource_->getGlobalDimension();
411      int elementPosition = it->first;
412      algo->computeGlobalSourceIndex(elementPosition,
413                                     gridDestinationDimensionSize,
414                                     gridSrcDimensionSize,
415                                     globalIndexGridDestSendToServer,
416                                     globaIndexWeightFromDestToSource);
417
418      if (1 < nbAlgos_)
419      {
420        // Compute transformation of global indexes among grids
421        computeTransformationFromOriginalGridSource(globaIndexWeightFromDestToSource);
422
423        // Now grid destination becomes grid source in a new transformation
424        if (nbAgloTransformation != (nbAlgos_-1)) setUpGrid(elementPositionInGrid, transType);
425      }
426      else
427      {
428        currentGridIndexToOriginalGridIndex_.swap(globaIndexWeightFromDestToSource);
429      }
430
431      ++nbAgloTransformation;
432    }
433  }
434
435  if (0 != nbAgloTransformation)
436  {
437    updateFinalGridDestination();
438    computeFinalTransformationMapping();
439  }
440}
441
442
443/*!
444  After applying the algorithms, there are some informations on grid destination needing change, for now, there are:
445   +) mask
446*/
447void CGridTransformation::updateFinalGridDestination()
448{
449  CContext* context = CContext::getCurrent();
450  CContextClient* client = context->client;
451
452  //First of all, retrieve info of local mask of grid destination
453  CDistributionClient distributionClientDest(client->clientRank, gridDestination_);
454  const std::vector<int>& localMaskIndexOnClientDest = distributionClientDest.getLocalMaskIndexOnClient();
455  const std::vector<size_t>& globalIndexOnClientDest = distributionClientDest.getGlobalDataIndexSendToServer();
456
457  std::vector<size_t>::const_iterator itbArr, itArr, iteArr;
458  itbArr = globalIndexOnClientDest.begin();
459  iteArr = globalIndexOnClientDest.end();
460
461  GlobalIndexMap::const_iterator iteGlobalMap = currentGridIndexToOriginalGridIndex_.end();
462  const size_t sfmax = NumTraits<unsigned long>::sfmax();
463  int maskIndexNum = 0;
464  for (itArr = itbArr; itArr != iteArr; ++itArr)
465  {
466    if (iteGlobalMap != currentGridIndexToOriginalGridIndex_.find(*itArr))
467    {
468      const std::vector<std::pair<size_t,double> >& vecIndex = currentGridIndexToOriginalGridIndex_[*itArr];
469      for (int idx = 0; idx < vecIndex.size(); ++idx)
470      {
471        if (sfmax == vecIndex[idx].first)
472        {
473          ++maskIndexNum;
474          break;
475        }
476      }
477    }
478  }
479
480  CArray<int,1> maskIndexToModify(maskIndexNum);
481  maskIndexNum = 0;
482  for (itArr = itbArr; itArr != iteArr; ++itArr)
483  {
484    if (iteGlobalMap != currentGridIndexToOriginalGridIndex_.find(*itArr))
485    {
486      const std::vector<std::pair<size_t,double> >& vecIndex = currentGridIndexToOriginalGridIndex_[*itArr];
487      for (int idx = 0; idx < vecIndex.size(); ++idx)
488      {
489        if (sfmax == vecIndex[idx].first)
490        {
491          int localIdx = std::distance(itbArr, itArr);
492          maskIndexToModify(maskIndexNum) = localMaskIndexOnClientDest[localIdx];
493          ++maskIndexNum;
494          break;
495        }
496      }
497    }
498  }
499
500  gridDestination_->modifyMask(maskIndexToModify);
501}
502
503/*!
504  A transformation from a grid source to grid destination often passes several intermediate grids, which play a role of
505temporary grid source and/or grid destination. This function makes sure that global index of original grid source are mapped correctly to
506the final grid destination
507*/
508void CGridTransformation::computeTransformationFromOriginalGridSource(const std::map<size_t, std::vector<std::pair<size_t,double> > >& globaIndexMapFromDestToSource)
509{
510  CContext* context = CContext::getCurrent();
511  CContextClient* client = context->client;
512
513  if (currentGridIndexToOriginalGridIndex_.empty())
514  {
515    currentGridIndexToOriginalGridIndex_ = globaIndexMapFromDestToSource;
516    return;
517  }
518
519  CTransformationMapping transformationMap(gridDestination_, gridSource_);
520
521    // Then compute transformation mapping among clients
522  transformationMap.computeTransformationMapping(globaIndexMapFromDestToSource);
523
524  const std::map<int,std::vector<std::vector<std::pair<size_t,double> > > >& globalIndexToReceive = transformationMap.getGlobalIndexReceivedOnGridDestMapping();
525  const std::map<int,std::vector<size_t> >& globalIndexToSend = transformationMap.getGlobalIndexSendToGridDestMapping();
526
527 // Sending global index of original grid source
528  std::map<int,std::vector<size_t> >::const_iterator itbSend = globalIndexToSend.begin(), itSend,
529                                                     iteSend = globalIndexToSend.end();
530 int sendBuffSize = 0;
531 for (itSend = itbSend; itSend != iteSend; ++itSend) sendBuffSize += (itSend->second).size();
532 // We use the first element of each block to send number of element in this block
533 sendBuffSize += globalIndexToSend.size();
534
535
536 typedef unsigned long Scalar;
537 Scalar* sendBuff, *currentSendBuff;
538 if (0 != sendBuffSize) sendBuff = new Scalar [sendBuffSize];
539 for (StdSize idx = 0; idx < sendBuffSize; ++idx) sendBuff[idx] = NumTraits<Scalar>::sfmax();
540
541 std::map<int, MPI_Request> requestsCurrentGrid, requestsOriginalGrid, requestsWeightGrid;
542 GlobalIndexMap::const_iterator iteGlobalIndex = currentGridIndexToOriginalGridIndex_.end();
543
544  // Only send global index of original source corresponding to non-masked index
545  // Use first position of each block to specify the number of elemnt in this block
546 int globalIndexOriginalSrcSendBuffSize = 0;
547 int currentBuffPosition = 0;
548 for (itSend = itbSend; itSend != iteSend; ++itSend)
549 {
550   int destRank = itSend->first;
551   const std::vector<size_t>& globalIndexOfCurrentGridSourceToSend = itSend->second;
552   int countSize  = globalIndexOfCurrentGridSourceToSend.size();
553   size_t countBlock = 0;
554   for (int idx = 0; idx < countSize; ++idx)
555   {
556     size_t index = globalIndexOfCurrentGridSourceToSend[idx];
557     if (iteGlobalIndex != currentGridIndexToOriginalGridIndex_.find(index))// searchCurrentSrc.search(itbIndex, iteIndex, globalIndexOfCurrentGridSourceToSend[idx], itIndex))
558     {
559       globalIndexOriginalSrcSendBuffSize += currentGridIndexToOriginalGridIndex_[index].size() + 1; // 1 for number of elements in this block
560       sendBuff[idx+currentBuffPosition+1] = index;
561       countBlock += currentGridIndexToOriginalGridIndex_[index].size() + 1;
562     }
563   }
564   sendBuff[currentBuffPosition] = countBlock;
565   currentSendBuff = sendBuff + currentBuffPosition;
566   MPI_Isend(currentSendBuff, countSize +1, MPI_UNSIGNED_LONG, destRank, MPI_GRID_TRANSFORMATION_CURRENT_GRID_INDEX, client->intraComm, &requestsCurrentGrid[destRank]);
567   currentBuffPosition += countSize + 1;
568 }
569
570 Scalar* sendOriginalIndexBuff, *currentOriginalIndexSendBuff;
571 if (0 != globalIndexOriginalSrcSendBuffSize) sendOriginalIndexBuff = new Scalar [globalIndexOriginalSrcSendBuffSize];
572 double* sendOriginalWeightBuff, *currentOriginalWeightSendBuff;
573 if (0 != globalIndexOriginalSrcSendBuffSize) sendOriginalWeightBuff = new double [globalIndexOriginalSrcSendBuffSize];
574
575 currentBuffPosition = 0;
576 for (itSend = itbSend; itSend != iteSend; ++itSend)
577 {
578   int destRank = itSend->first;
579   const std::vector<size_t>& globalIndexOfCurrentGridSourceToSend = itSend->second;
580   int countSize = globalIndexOfCurrentGridSourceToSend.size();
581   int increaseStep = 0;
582   for (int idx = 0; idx < countSize; ++idx)
583   {
584     size_t index = globalIndexOfCurrentGridSourceToSend[idx];
585     if (iteGlobalIndex != currentGridIndexToOriginalGridIndex_.find(index))
586     {
587       size_t vectorSize = currentGridIndexToOriginalGridIndex_[index].size();
588       sendOriginalIndexBuff[currentBuffPosition+increaseStep]  = vectorSize;
589       sendOriginalWeightBuff[currentBuffPosition+increaseStep] = (double)vectorSize;
590       const std::vector<std::pair<size_t,double> >& indexWeightPair = currentGridIndexToOriginalGridIndex_[index];
591       for (size_t i = 0; i < vectorSize; ++i)
592       {
593         ++increaseStep;
594         sendOriginalIndexBuff[currentBuffPosition+increaseStep]  = indexWeightPair[i].first;
595         sendOriginalWeightBuff[currentBuffPosition+increaseStep] = indexWeightPair[i].second;
596       }
597       ++increaseStep;
598     }
599   }
600
601   currentOriginalIndexSendBuff = sendOriginalIndexBuff + currentBuffPosition;
602   currentOriginalWeightSendBuff = sendOriginalWeightBuff + currentBuffPosition;
603   if (0 != increaseStep)
604   {
605     MPI_Isend(currentOriginalIndexSendBuff, increaseStep, MPI_UNSIGNED_LONG, destRank, MPI_GRID_TRANSFORMATION_ORIGINAL_GRID_INDEX, client->intraComm, &requestsOriginalGrid[destRank]);
606     MPI_Isend(currentOriginalWeightSendBuff, increaseStep, MPI_DOUBLE, destRank, MPI_GRID_TRANSFORMATION_ORIGINAL_GRID_WEIGHT, client->intraComm, &requestsWeightGrid[destRank]);
607   }
608   currentBuffPosition += increaseStep;
609 }
610
611
612 // Receiving global index of grid source sending from current grid source
613 std::map<int,std::vector<std::vector<std::pair<size_t,double> > > >::const_iterator itbRecv = globalIndexToReceive.begin(), itRecv,
614                                                                                     iteRecv = globalIndexToReceive.end();
615 int recvBuffSize = 0;
616 for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv) recvBuffSize += (itRecv->second).size();
617 recvBuffSize += globalIndexToReceive.size();
618
619 Scalar* recvBuff, *currentRecvBuff;
620 if (0 != recvBuffSize) recvBuff = new Scalar [recvBuffSize];
621 for (StdSize idx = 0; idx < recvBuffSize; ++idx) recvBuff[idx] = NumTraits<Scalar>::sfmax();
622
623 std::map<int,int> countBlockMap;
624 int globalIndexOriginalSrcRecvBuffSize = 0;
625 int currentRecvBuffPosition = 0;
626 for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv)
627 {
628   MPI_Status status;
629   int srcRank = itRecv->first;
630   int countSize = (itRecv->second).size();
631   currentRecvBuff = recvBuff + currentRecvBuffPosition;
632   MPI_Recv(currentRecvBuff, countSize +1, MPI_UNSIGNED_LONG, srcRank, MPI_GRID_TRANSFORMATION_CURRENT_GRID_INDEX, client->intraComm, &status);
633   globalIndexOriginalSrcRecvBuffSize += *currentRecvBuff;
634   countBlockMap[srcRank] = *currentRecvBuff;
635   currentRecvBuffPosition += countSize +1;
636 }
637
638 Scalar* recvOriginalIndexBuff, *currentOriginalIndexRecvBuff;
639 if (0 != globalIndexOriginalSrcRecvBuffSize) recvOriginalIndexBuff = new Scalar [globalIndexOriginalSrcRecvBuffSize];
640 double* recvOriginalWeightBuff, *currentOriginalWeightRecvBuff;
641 if (0 != globalIndexOriginalSrcRecvBuffSize) recvOriginalWeightBuff = new double [globalIndexOriginalSrcRecvBuffSize];
642
643 int countBlock = 0;
644 currentRecvBuffPosition = 0;
645 currentBuffPosition = 0;
646 for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv)
647 {
648   MPI_Status statusIndex, statusWeight;
649   int srcRank = itRecv->first;
650   countBlock = countBlockMap[srcRank];
651   currentOriginalIndexRecvBuff = recvOriginalIndexBuff + currentBuffPosition;
652   currentOriginalWeightRecvBuff = recvOriginalWeightBuff + currentBuffPosition;
653   if (0 != countBlock)
654   {
655     MPI_Recv(currentOriginalIndexRecvBuff, countBlock, MPI_UNSIGNED_LONG, srcRank, MPI_GRID_TRANSFORMATION_ORIGINAL_GRID_INDEX, client->intraComm, &statusIndex);
656     MPI_Recv(currentOriginalWeightRecvBuff, countBlock, MPI_DOUBLE, srcRank, MPI_GRID_TRANSFORMATION_ORIGINAL_GRID_WEIGHT, client->intraComm, &statusWeight);
657   }
658   currentBuffPosition += countBlock;
659 }
660
661 // We process everything in here, even case of masked index
662 // The way to process masked index needs discussing
663 const size_t sfmax = NumTraits<unsigned long>::sfmax();
664 GlobalIndexMap currentToOriginalTmp;
665
666 currentRecvBuffPosition = 0;
667 currentRecvBuff = recvBuff;
668 currentOriginalIndexRecvBuff  = recvOriginalIndexBuff;
669 currentOriginalWeightRecvBuff = recvOriginalWeightBuff;
670 for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv)
671 {
672   int countBlockRank = countBlockMap[itRecv->first];
673
674   ++currentRecvBuff;  // it's very subtle here, pay attention
675   int countSize = (itRecv->second).size();
676   for (int idx = 0; idx < countSize; ++idx)
677   {
678      ++currentRecvBuff;
679     int ssize = (itRecv->second)[idx].size();
680     if (sfmax != *currentRecvBuff)
681     {
682       if (0 != countBlockRank)
683       {
684         countBlock = *(currentOriginalIndexRecvBuff+currentRecvBuffPosition);
685         for (int i = 0; i < ssize; ++i)
686         {
687           for (int j = 0; j < countBlock; ++j)
688           {
689             size_t globalOriginalIndex = *(currentOriginalIndexRecvBuff+currentRecvBuffPosition+j+1);
690             double weightGlobal = *(currentOriginalWeightRecvBuff+currentRecvBuffPosition+j+1) * (itRecv->second)[idx][i].second;
691             currentToOriginalTmp[(itRecv->second)[idx][i].first].push_back(make_pair(globalOriginalIndex,weightGlobal));
692           }
693         }
694         currentRecvBuffPosition += countBlock+1;
695       }
696     }
697//     else
698//     {
699//       for (int i = 0; i < ssize; ++i)
700//       {
701//         currentToOriginalTmp[(itRecv->second)[idx][i].first].push_back(make_pair(sfmax,1.0));
702//       }
703//     }
704   }
705 }
706
707 currentGridIndexToOriginalGridIndex_.swap(currentToOriginalTmp);
708
709 std::map<int, MPI_Request>::iterator itRequest;
710 for (itRequest = requestsCurrentGrid.begin(); itRequest != requestsCurrentGrid.end(); ++itRequest)
711   MPI_Wait(&itRequest->second, MPI_STATUS_IGNORE);
712 for (itRequest = requestsOriginalGrid.begin(); itRequest != requestsOriginalGrid.end(); ++itRequest)
713   MPI_Wait(&itRequest->second, MPI_STATUS_IGNORE);
714 for (itRequest = requestsWeightGrid.begin(); itRequest != requestsWeightGrid.end(); ++itRequest)
715   MPI_Wait(&itRequest->second, MPI_STATUS_IGNORE);
716
717 if (0 != sendBuffSize) delete [] sendBuff;
718 if (0 != recvBuffSize) delete [] recvBuff;
719 if (0 != globalIndexOriginalSrcSendBuffSize) delete [] sendOriginalIndexBuff;
720 if (0 != globalIndexOriginalSrcSendBuffSize) delete [] sendOriginalWeightBuff;
721 if (0 != globalIndexOriginalSrcRecvBuffSize) delete [] recvOriginalIndexBuff;
722 if (0 != globalIndexOriginalSrcRecvBuffSize) delete [] recvOriginalWeightBuff;
723}
724
725/*!
726  Compute transformation mapping between grid source and grid destination
727  The transformation between grid source and grid destination is represented in form of mapping between global index
728of two grids. Then local index mapping between data on each grid will be found out thanks to these global indexes
729*/
730void CGridTransformation::computeFinalTransformationMapping()
731{
732  CContext* context = CContext::getCurrent();
733  CContextClient* client = context->client;
734
735  CTransformationMapping transformationMap(gridDestination_, originalGridSource_);
736
737  transformationMap.computeTransformationMapping(currentGridIndexToOriginalGridIndex_);
738
739  const std::map<int,std::vector<std::vector<std::pair<size_t,double> > > >& globalIndexToReceive = transformationMap.getGlobalIndexReceivedOnGridDestMapping();
740  const std::map<int,std::vector<size_t> >& globalIndexToSend = transformationMap.getGlobalIndexSendToGridDestMapping();
741
742  CDistributionClient distributionClientDest(client->clientRank, gridDestination_);
743  CDistributionClient distributionClientSrc(client->clientRank, originalGridSource_);
744
745  const std::vector<size_t>& globalIndexOnClientDest = distributionClientDest.getGlobalDataIndexSendToServer();
746  const std::vector<size_t>& globalIndexOnClientSrc = distributionClientSrc.getGlobalDataIndexSendToServer();
747
748  std::vector<size_t>::const_iterator itbArr, itArr, iteArr;
749  std::vector<int>::const_iterator itIndex, itbIndex, iteIndex;
750  std::map<int,std::vector<std::vector<std::pair<size_t,double> > > >::const_iterator itbMapRecv, itMapRecv, iteMapRecv;
751
752  std::vector<int> permutIndex;
753  typedef XIOSBinarySearchWithIndex<size_t> BinarySearch;
754
755  // Find out local index on grid destination (received)
756  XIOSAlgorithms::fillInIndex(globalIndexOnClientDest.size(), permutIndex);
757  XIOSAlgorithms::sortWithIndex<size_t, CVectorStorage>(globalIndexOnClientDest, permutIndex);
758  itbIndex = permutIndex.begin();
759  iteIndex = permutIndex.end();
760  BinarySearch searchClientDest(globalIndexOnClientDest);
761  itbMapRecv = globalIndexToReceive.begin();
762  iteMapRecv = globalIndexToReceive.end();
763  for (itMapRecv = itbMapRecv; itMapRecv != iteMapRecv; ++itMapRecv)
764  {
765    int sourceRank = itMapRecv->first;
766    int numGlobalIndex = (itMapRecv->second).size();
767    for (int i = 0; i < numGlobalIndex; ++i)
768    {
769      int vecSize = ((itMapRecv->second)[i]).size();
770      std::vector<std::pair<int,double> > tmpVec;
771      for (int idx = 0; idx < vecSize; ++idx)
772      {
773        size_t globalIndex = (itMapRecv->second)[i][idx].first;
774        double weight = (itMapRecv->second)[i][idx].second;
775        if (searchClientDest.search(itbIndex, iteIndex, globalIndex, itIndex))
776        {
777          tmpVec.push_back(make_pair(*itIndex, weight));
778        }
779      }
780      localIndexToReceiveOnGridDest_[sourceRank].push_back(tmpVec);
781    }
782  }
783
784  // Find out local index on grid source (to send)
785  std::map<int,std::vector<size_t> >::const_iterator itbMap, itMap, iteMap;
786  XIOSAlgorithms::fillInIndex(globalIndexOnClientSrc.size(), permutIndex);
787  XIOSAlgorithms::sortWithIndex<size_t, CVectorStorage>(globalIndexOnClientSrc, permutIndex);
788  itbIndex = permutIndex.begin();
789  iteIndex = permutIndex.end();
790  BinarySearch searchClientSrc(globalIndexOnClientSrc);
791  itbMap = globalIndexToSend.begin();
792  iteMap = globalIndexToSend.end();
793  for (itMap = itbMap; itMap != iteMap; ++itMap)
794  {
795    int destRank = itMap->first;
796    int vecSize = itMap->second.size();
797    localIndexToSendFromGridSource_[destRank].resize(vecSize);
798    for (int idx = 0; idx < vecSize; ++idx)
799    {
800      if (searchClientSrc.search(itbIndex, iteIndex, itMap->second[idx], itIndex))
801      {
802        localIndexToSendFromGridSource_[destRank](idx) = *itIndex;
803      }
804    }
805  }
806}
807
808bool CGridTransformation::isSpecialTransformation(ETranformationType transType)
809{
810  bool res;
811  switch (transType)
812  {
813    case TRANS_GENERATE_RECTILINEAR_DOMAIN:
814     res = true;
815     break;
816    default:
817     res = false;
818     break;
819  }
820
821  return res;
822}
823
824/*!
825  Local index of data which need sending from the grid source
826  \return local index of data
827*/
828const std::map<int, CArray<int,1> >& CGridTransformation::getLocalIndexToSendFromGridSource() const
829{
830  return localIndexToSendFromGridSource_;
831}
832
833/*!
834  Local index of data which will be received on the grid destination
835  \return local index of data
836*/
837const std::map<int,std::vector<std::vector<std::pair<int,double> > > >& CGridTransformation::getLocalIndexToReceiveOnGridDest() const
838{
839  return localIndexToReceiveOnGridDest_;
840}
841
842}
Note: See TracBrowser for help on using the repository browser.