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

Last change on this file since 824 was 824, checked in by mhnguyen, 8 years ago

Correcting the behavior of field_ref

+) If a field refers to another one via field_ref and there is no transformation between
grid of this field and one of refered field, this field will inherit attributes from field_ref.

Test
+) On Curie
+) All tests pass

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