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

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

Cleaning up some redundant coeds and making some improvements

+) Remove some XIOS Search to make code run faster
+) Remove some commented codes

Test
+) On Curie
+) All tests pass

File size: 33.3 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 CDistributionClient::GlobalLocalDataMap& globalLocalIndexGridDestSendToServer = distributionClientDest.getGlobalLocalDataSendToServer();
395
396      // ComputeTransformation of global index of each element
397      std::vector<int> gridDestinationDimensionSize = gridDestination_->getGlobalDimension();
398      std::vector<int> gridSrcDimensionSize = gridSource_->getGlobalDimension();
399      int elementPosition = it->first;
400      algo->computeGlobalSourceIndex(elementPosition,
401                                     gridDestinationDimensionSize,
402                                     gridSrcDimensionSize,
403                                     globalLocalIndexGridDestSendToServer,
404                                     globaIndexWeightFromDestToSource);
405
406      if (1 < nbAlgos_)
407      {
408        // Compute transformation of global indexes among grids
409        computeTransformationFromOriginalGridSource(globaIndexWeightFromDestToSource);
410
411        // Now grid destination becomes grid source in a new transformation
412        if (nbAgloTransformation != (nbAlgos_-1)) setUpGrid(elementPositionInGrid, transType, nbAgloTransformation);
413      }
414      else
415      {
416        currentGridIndexToOriginalGridIndex_.swap(globaIndexWeightFromDestToSource);
417      }
418
419      ++nbAgloTransformation;
420    }
421  }
422
423  if (0 != nbAgloTransformation)
424  {
425    updateFinalGridDestination();
426    computeFinalTransformationMapping();
427  }
428}
429
430
431/*!
432  After applying the algorithms, there are some informations on grid destination needing change, for now, there are:
433   +) mask
434*/
435void CGridTransformation::updateFinalGridDestination()
436{
437  CContext* context = CContext::getCurrent();
438  CContextClient* client = context->client;
439
440  //First of all, retrieve info of local mask of grid destination
441  CDistributionClient distributionClientDest(client->clientRank, gridDestination_);
442  const std::vector<int>& localMaskIndexOnClientDest = distributionClientDest.getLocalMaskIndexOnClient();
443  const CDistributionClient::GlobalLocalDataMap& globalIndexOnClientDest = distributionClientDest.getGlobalLocalDataSendToServer();
444
445  CDistributionClient::GlobalLocalDataMap::const_iterator itbArr, itArr, iteArr;
446  itbArr = globalIndexOnClientDest.begin();
447  iteArr = globalIndexOnClientDest.end();
448
449  DestinationIndexMap::const_iterator iteGlobalMap = currentGridIndexToOriginalGridIndex_.end();
450  const size_t sfmax = NumTraits<unsigned long>::sfmax();
451  int maskIndexNum = 0;
452  for (itArr = itbArr; itArr != iteArr; ++itArr)
453  {
454    if (iteGlobalMap != currentGridIndexToOriginalGridIndex_.find(itArr->first))
455    {
456      const std::vector<std::pair<int, std::pair<size_t,double> > >& vecIndex = currentGridIndexToOriginalGridIndex_[itArr->first];
457      for (int idx = 0; idx < vecIndex.size(); ++idx)
458      {
459        if (sfmax == (vecIndex[idx].second).first)
460        {
461          ++maskIndexNum;
462          break;
463        }
464      }
465    }
466  }
467
468  CArray<int,1> maskIndexToModify(maskIndexNum);
469  maskIndexNum = 0;
470  for (itArr = itbArr; itArr != iteArr; ++itArr)
471  {
472    if (iteGlobalMap != currentGridIndexToOriginalGridIndex_.find(itArr->first))
473    {
474      const std::vector<std::pair<int, std::pair<size_t,double> > >& vecIndex = currentGridIndexToOriginalGridIndex_[itArr->first];
475      for (int idx = 0; idx < vecIndex.size(); ++idx)
476      {
477        if (sfmax == (vecIndex[idx].second).first)
478        {
479          int localIdx = std::distance(itbArr, itArr);
480          maskIndexToModify(maskIndexNum) = localMaskIndexOnClientDest[localIdx];
481          ++maskIndexNum;
482          break;
483        }
484      }
485    }
486  }
487
488  gridDestination_->modifyMask(maskIndexToModify);
489}
490
491/*!
492  A transformation from a grid source to grid destination often passes several intermediate grids, which play a role of
493temporary grid source and/or grid destination. This function makes sure that global index of original grid source are mapped correctly to
494the final grid destination
495*/
496void CGridTransformation::computeTransformationFromOriginalGridSource(const DestinationIndexMap& globaIndexMapFromDestToSource)
497{
498  CContext* context = CContext::getCurrent();
499  CContextClient* client = context->client;
500
501  if (currentGridIndexToOriginalGridIndex_.empty())
502  {
503    currentGridIndexToOriginalGridIndex_ = globaIndexMapFromDestToSource;
504    return;
505  }
506
507  CTransformationMapping transformationMap(gridDestination_, gridSource_);
508
509    // Then compute transformation mapping among clients
510  transformationMap.computeTransformationMapping(globaIndexMapFromDestToSource);
511
512  const CTransformationMapping::ReceivedIndexMap& globalIndexToReceive = transformationMap.getGlobalIndexReceivedOnGridDestMapping();
513  const CTransformationMapping::SentIndexMap& globalIndexToSend = transformationMap.getGlobalIndexSendToGridDestMapping();
514
515 // Sending global index of original grid source
516  CTransformationMapping::SentIndexMap::const_iterator itbSend = globalIndexToSend.begin(), itSend,
517                                                       iteSend = globalIndexToSend.end();
518 int sendBuffSize = 0;
519 for (itSend = itbSend; itSend != iteSend; ++itSend) sendBuffSize += (itSend->second).size();
520 // We use the first element of each block to send number of element in this block
521 sendBuffSize += globalIndexToSend.size();
522
523
524 typedef unsigned long Scalar;
525 Scalar* sendBuff, *currentSendBuff;
526 if (0 != sendBuffSize) sendBuff = new Scalar [sendBuffSize];
527 for (StdSize idx = 0; idx < sendBuffSize; ++idx) sendBuff[idx] = NumTraits<Scalar>::sfmax();
528
529 std::map<int, MPI_Request> requestsCurrentGrid, requestsOriginalGridGlobalIndex, requestsOriginalGridLocalIndex, requestsWeightGrid;
530 DestinationIndexMap::const_iterator iteGlobalIndex = currentGridIndexToOriginalGridIndex_.end();
531
532  // Only send global index of original source corresponding to non-masked index
533  // Use first position of each block to specify the number of elemnt in this block
534 int globalIndexOriginalSrcSendBuffSize = 0;
535 int currentBuffPosition = 0;
536 for (itSend = itbSend; itSend != iteSend; ++itSend)
537 {
538   int destRank = itSend->first;
539   const std::vector<std::pair<int, size_t> >& globalIndexOfCurrentGridSourceToSend = itSend->second;
540   int countSize  = globalIndexOfCurrentGridSourceToSend.size();
541   size_t countBlock = 0;
542   for (int idx = 0; idx < countSize; ++idx)
543   {
544     size_t index = globalIndexOfCurrentGridSourceToSend[idx].second;
545     if (iteGlobalIndex != currentGridIndexToOriginalGridIndex_.find(index))
546     {
547       globalIndexOriginalSrcSendBuffSize += currentGridIndexToOriginalGridIndex_[index].size() + 1; // 1 for number of elements in this block
548       sendBuff[idx+currentBuffPosition+1] = index;
549       countBlock += currentGridIndexToOriginalGridIndex_[index].size() + 1;
550     }
551   }
552   sendBuff[currentBuffPosition] = countBlock;
553   currentSendBuff = sendBuff + currentBuffPosition;
554   MPI_Isend(currentSendBuff, countSize +1, MPI_UNSIGNED_LONG, destRank, MPI_GRID_TRANSFORMATION_CURRENT_GRID_INDEX, client->intraComm, &requestsCurrentGrid[destRank]);
555   currentBuffPosition += countSize + 1;
556 }
557
558 Scalar* sendOriginalGlobalIndexBuff, *currentOriginalGlobalIndexSendBuff;
559 if (0 != globalIndexOriginalSrcSendBuffSize) sendOriginalGlobalIndexBuff = new Scalar [globalIndexOriginalSrcSendBuffSize];
560 double* sendOriginalWeightBuff, *currentOriginalWeightSendBuff;
561 if (0 != globalIndexOriginalSrcSendBuffSize) sendOriginalWeightBuff = new double [globalIndexOriginalSrcSendBuffSize];
562
563 currentBuffPosition = 0;
564 for (itSend = itbSend; itSend != iteSend; ++itSend)
565 {
566   int destRank = itSend->first;
567   const std::vector<std::pair<int, size_t> >& globalIndexOfCurrentGridSourceToSend = itSend->second;
568   int countSize = globalIndexOfCurrentGridSourceToSend.size();
569   int increaseStep = 0;
570   for (int idx = 0; idx < countSize; ++idx)
571   {
572     size_t index = globalIndexOfCurrentGridSourceToSend[idx].second;
573     if (iteGlobalIndex != currentGridIndexToOriginalGridIndex_.find(index))
574     {
575       size_t vectorSize = currentGridIndexToOriginalGridIndex_[index].size();
576       sendOriginalGlobalIndexBuff[currentBuffPosition+increaseStep]  = vectorSize;
577       sendOriginalWeightBuff[currentBuffPosition+increaseStep] = (double)vectorSize;
578       const std::vector<std::pair<int, std::pair<size_t,double> > >& indexWeightPair = currentGridIndexToOriginalGridIndex_[index];
579       for (size_t i = 0; i < vectorSize; ++i)
580       {
581         ++increaseStep;
582         sendOriginalGlobalIndexBuff[currentBuffPosition+increaseStep]  = (indexWeightPair[i].second).first;
583         sendOriginalWeightBuff[currentBuffPosition+increaseStep] = (indexWeightPair[i].second).second;
584       }
585       ++increaseStep;
586     }
587   }
588
589   currentOriginalGlobalIndexSendBuff = sendOriginalGlobalIndexBuff + currentBuffPosition;
590   currentOriginalWeightSendBuff = sendOriginalWeightBuff + currentBuffPosition;
591   if (0 != increaseStep)
592   {
593     MPI_Isend(currentOriginalGlobalIndexSendBuff, increaseStep, MPI_UNSIGNED_LONG, destRank,
594               MPI_GRID_TRANSFORMATION_ORIGINAL_GRID_GLOBAL_INDEX, client->intraComm, &requestsOriginalGridGlobalIndex[destRank]);
595     MPI_Isend(currentOriginalWeightSendBuff, increaseStep, MPI_DOUBLE, destRank,
596               MPI_GRID_TRANSFORMATION_ORIGINAL_GRID_WEIGHT, client->intraComm, &requestsWeightGrid[destRank]);
597   }
598   currentBuffPosition += increaseStep;
599 }
600
601
602 // Receiving global index of grid source sending from current grid source
603 CTransformationMapping::ReceivedIndexMap::const_iterator itbRecv = globalIndexToReceive.begin(), itRecv,
604                                                          iteRecv = globalIndexToReceive.end();
605 int recvBuffSize = 0;
606 for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv) recvBuffSize += (itRecv->second).size();
607 recvBuffSize += globalIndexToReceive.size();
608
609 Scalar* recvBuff, *currentRecvBuff;
610 if (0 != recvBuffSize) recvBuff = new Scalar [recvBuffSize];
611 for (StdSize idx = 0; idx < recvBuffSize; ++idx) recvBuff[idx] = NumTraits<Scalar>::sfmax();
612
613 std::map<int,int> countBlockMap;
614 int globalIndexOriginalSrcRecvBuffSize = 0;
615 int currentRecvBuffPosition = 0;
616 for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv)
617 {
618   MPI_Status status;
619   int srcRank = itRecv->first;
620   int countSize = (itRecv->second).size();
621   currentRecvBuff = recvBuff + currentRecvBuffPosition;
622   MPI_Recv(currentRecvBuff, countSize +1, MPI_UNSIGNED_LONG, srcRank, MPI_GRID_TRANSFORMATION_CURRENT_GRID_INDEX, client->intraComm, &status);
623   globalIndexOriginalSrcRecvBuffSize += *currentRecvBuff;
624   countBlockMap[srcRank] = *currentRecvBuff;
625   currentRecvBuffPosition += countSize +1;
626 }
627
628 Scalar* recvOriginalGlobalIndexBuff, *currentOriginalGlobalIndexRecvBuff;
629 if (0 != globalIndexOriginalSrcRecvBuffSize) recvOriginalGlobalIndexBuff = new Scalar [globalIndexOriginalSrcRecvBuffSize];
630 double* recvOriginalWeightBuff, *currentOriginalWeightRecvBuff;
631 if (0 != globalIndexOriginalSrcRecvBuffSize) recvOriginalWeightBuff = new double [globalIndexOriginalSrcRecvBuffSize];
632
633 int countBlock = 0;
634 currentRecvBuffPosition = 0;
635 currentBuffPosition = 0;
636 for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv)
637 {
638   MPI_Status statusGlobalIndex, statusLocalIndex, statusWeight;
639   int srcRank = itRecv->first;
640   countBlock = countBlockMap[srcRank];
641   currentOriginalGlobalIndexRecvBuff = recvOriginalGlobalIndexBuff + currentBuffPosition;
642   currentOriginalWeightRecvBuff = recvOriginalWeightBuff + currentBuffPosition;
643   if (0 != countBlock)
644   {
645     MPI_Recv(currentOriginalGlobalIndexRecvBuff, countBlock, MPI_UNSIGNED_LONG, srcRank, MPI_GRID_TRANSFORMATION_ORIGINAL_GRID_GLOBAL_INDEX, client->intraComm, &statusGlobalIndex);
646     MPI_Recv(currentOriginalWeightRecvBuff, countBlock, MPI_DOUBLE, srcRank, MPI_GRID_TRANSFORMATION_ORIGINAL_GRID_WEIGHT, client->intraComm, &statusWeight);
647   }
648   currentBuffPosition += countBlock;
649 }
650
651 // We process everything in here, even case of masked index
652 // The way to process masked index needs discussing
653 const size_t sfmax = NumTraits<unsigned long>::sfmax();
654 DestinationIndexMap currentToOriginalTmp;
655
656 currentRecvBuffPosition = 0;
657 currentRecvBuff = recvBuff;
658 currentOriginalGlobalIndexRecvBuff  = recvOriginalGlobalIndexBuff;
659 currentOriginalWeightRecvBuff = recvOriginalWeightBuff;
660 for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv)
661 {
662   int countBlockRank = countBlockMap[itRecv->first];
663
664   ++currentRecvBuff;  // it's very subtle here, pay attention
665   int countSize = (itRecv->second).size();
666   for (int idx = 0; idx < countSize; ++idx)
667   {
668      ++currentRecvBuff;
669     int ssize = (itRecv->second)[idx].size();
670     if (sfmax != *currentRecvBuff)
671     {
672       if (0 != countBlockRank)
673       {
674         countBlock = *(currentOriginalGlobalIndexRecvBuff+currentRecvBuffPosition);
675         for (int i = 0; i < ssize; ++i)
676         {
677           for (int j = 0; j < countBlock; ++j)
678           {
679             size_t globalOriginalIndex = *(currentOriginalGlobalIndexRecvBuff+currentRecvBuffPosition+j+1);
680             int currentGridLocalIndex = (itRecv->second)[idx][i].first;
681             double weightGlobal = *(currentOriginalWeightRecvBuff+currentRecvBuffPosition+j+1) * (itRecv->second)[idx][i].second.second;
682             currentToOriginalTmp[(itRecv->second)[idx][i].second.first].push_back(make_pair(currentGridLocalIndex,make_pair(globalOriginalIndex,weightGlobal)));
683           }
684         }
685         currentRecvBuffPosition += countBlock+1;
686       }
687     }
688//     else
689//     {
690//       for (int i = 0; i < ssize; ++i)
691//       {
692//         currentToOriginalTmp[(itRecv->second)[idx][i].first].push_back(make_pair(sfmax,1.0));
693//       }
694//     }
695   }
696 }
697
698 currentGridIndexToOriginalGridIndex_.swap(currentToOriginalTmp);
699
700 std::map<int, MPI_Request>::iterator itRequest;
701 for (itRequest = requestsCurrentGrid.begin(); itRequest != requestsCurrentGrid.end(); ++itRequest)
702   MPI_Wait(&itRequest->second, MPI_STATUS_IGNORE);
703 for (itRequest = requestsOriginalGridGlobalIndex.begin(); itRequest != requestsOriginalGridGlobalIndex.end(); ++itRequest)
704   MPI_Wait(&itRequest->second, MPI_STATUS_IGNORE);
705 for (itRequest = requestsWeightGrid.begin(); itRequest != requestsWeightGrid.end(); ++itRequest)
706   MPI_Wait(&itRequest->second, MPI_STATUS_IGNORE);
707
708 if (0 != sendBuffSize) delete [] sendBuff;
709 if (0 != recvBuffSize) delete [] recvBuff;
710 if (0 != globalIndexOriginalSrcSendBuffSize) delete [] sendOriginalGlobalIndexBuff;
711 if (0 != globalIndexOriginalSrcSendBuffSize) delete [] sendOriginalWeightBuff;
712 if (0 != globalIndexOriginalSrcRecvBuffSize) delete [] recvOriginalGlobalIndexBuff;
713 if (0 != globalIndexOriginalSrcRecvBuffSize) delete [] recvOriginalWeightBuff;
714}
715
716/*!
717  Compute transformation mapping between grid source and grid destination
718  The transformation between grid source and grid destination is represented in form of mapping between global index
719of two grids. Then local index mapping between data on each grid will be found out thanks to these global indexes
720*/
721void CGridTransformation::computeFinalTransformationMapping()
722{
723  CContext* context = CContext::getCurrent();
724  CContextClient* client = context->client;
725
726  CTransformationMapping transformationMap(gridDestination_, originalGridSource_);
727
728  transformationMap.computeTransformationMapping(currentGridIndexToOriginalGridIndex_);
729
730  const CTransformationMapping::ReceivedIndexMap& globalIndexToReceive = transformationMap.getGlobalIndexReceivedOnGridDestMapping();
731  CTransformationMapping::ReceivedIndexMap::const_iterator itbMapRecv, itMapRecv, iteMapRecv;
732  itbMapRecv = globalIndexToReceive.begin();
733  iteMapRecv = globalIndexToReceive.end();
734  for (itMapRecv = itbMapRecv; itMapRecv != iteMapRecv; ++itMapRecv)
735  {
736    int sourceRank = itMapRecv->first;
737    int numGlobalIndex = (itMapRecv->second).size();
738    localIndexToReceiveOnGridDest_[sourceRank].resize(numGlobalIndex);
739    for (int i = 0; i < numGlobalIndex; ++i)
740    {
741      int vecSize = ((itMapRecv->second)[i]).size();
742      for (int idx = 0; idx < vecSize; ++idx)
743      {
744        const std::pair<int, std::pair<size_t,double> >& tmpPair = (itMapRecv->second)[i][idx];
745        localIndexToReceiveOnGridDest_[sourceRank][i].push_back(make_pair(tmpPair.first, tmpPair.second.second));
746      }
747    }
748  }
749
750  // Find out local index on grid source (to send)
751  const CTransformationMapping::SentIndexMap& globalIndexToSend = transformationMap.getGlobalIndexSendToGridDestMapping();
752  CTransformationMapping::SentIndexMap::const_iterator itbMap, itMap, iteMap;
753  itbMap = globalIndexToSend.begin();
754  iteMap = globalIndexToSend.end();
755  for (itMap = itbMap; itMap != iteMap; ++itMap)
756  {
757    int destRank = itMap->first;
758    int vecSize = itMap->second.size();
759    localIndexToSendFromGridSource_[destRank].resize(vecSize);
760    for (int idx = 0; idx < vecSize; ++idx)
761    {
762      localIndexToSendFromGridSource_[destRank](idx) = itMap->second[idx].first;
763    }
764  }
765}
766
767bool CGridTransformation::isSpecialTransformation(ETranformationType transType)
768{
769  bool res;
770  switch (transType)
771  {
772    case TRANS_GENERATE_RECTILINEAR_DOMAIN:
773     res = true;
774     break;
775    default:
776     res = false;
777     break;
778  }
779
780  return res;
781}
782
783/*!
784  Local index of data which need sending from the grid source
785  \return local index of data
786*/
787const std::map<int, CArray<int,1> >& CGridTransformation::getLocalIndexToSendFromGridSource() const
788{
789  return localIndexToSendFromGridSource_;
790}
791
792/*!
793  Local index of data which will be received on the grid destination
794  \return local index of data
795*/
796const std::map<int,std::vector<std::vector<std::pair<int,double> > > >& CGridTransformation::getLocalIndexToReceiveOnGridDest() const
797{
798  return localIndexToReceiveOnGridDest_;
799}
800
801}
Note: See TracBrowser for help on using the repository browser.