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

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

Weight computation of dynamic transformation is done only one for each time stamp

+) Each weight computation of dynamic transformation attached to timestamp
+) Remove some redundant codes

Test
+) On Curie
+) All tests pass

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