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

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

Implementing dynamic interpolation on axis

+) Change grid transformation to make it more flexible
+) Make some small improvements

Test
+) On Curie
+) All test pass

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