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

Last change on this file since 829 was 829, checked in by mhnguyen, 5 years ago

Refactoring transformation code

+) On exchanging information during transformation, not only global index are sent but also local index
+) Correct a bug in distributed hash table (dht)
+) Add new type for dht
+) Clean up some redundant codes

Test
+) On Curie
+) Every test passes
+) Code runs faster in some cases (up to 30%)

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