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

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

Various clean up

+) Remove some redundant codes

Test
+) On Curie
+) tests pass

  • Property svn:executable set to *
File size: 27.6 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 "axis_algorithm_transformation.hpp"
18#include "distribution_client.hpp"
19#include "mpi_tag.hpp"
20#include "grid.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), 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  initializeTransformations();
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::initializeTransformations()
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(i);
97    }
98    else
99    {
100      domPositionInGrid.push_back(i);
101    }
102  }
103
104  for (int i = 0; i < axisListDestP.size(); ++i)
105  {
106    elementPosition2AxisPositionInGrid_[axisPositionInGrid[i]] = i;
107  }
108
109  for (int i = 0; i < domListDestP.size(); ++i)
110  {
111    elementPosition2DomainPositionInGrid_[domPositionInGrid[i]] = i;
112  }
113
114  idx = 0;
115  for (int i = 0; i < gridDestination_->axis_domain_order.numElements(); ++i)
116  {
117    if (false == (gridDestination_->axis_domain_order)(i))
118    {
119      initializeAxisAlgorithms(i);
120    }
121    else
122    {
123      initializeDomainAlgorithms(i);
124    }
125  }
126}
127
128/*!
129  Initialize the algorithms corresponding to transformation info contained in each axis.
130If an axis has transformations, these transformations will be represented in form of vector of CTransformation pointers
131In general, each axis can have several transformations performed on itself. However, should they be done seperately or combinely (of course in order)?
132For now, one approach is to do these combinely but maybe this needs changing.
133\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)
134*/
135void CGridTransformation::initializeAxisAlgorithms(int axisPositionInGrid)
136{
137  std::vector<CAxis*> axisListDestP = gridDestination_->getAxis();
138  if (!axisListDestP.empty())
139  {
140    if (axisListDestP[elementPosition2AxisPositionInGrid_[axisPositionInGrid]]->hasTransformation())
141    {
142      CAxis::TransMapTypes trans = axisListDestP[elementPosition2AxisPositionInGrid_[axisPositionInGrid]]->getAllTransformations();
143      CAxis::TransMapTypes::const_iterator itb = trans.begin(), it,
144                                           ite = trans.end();
145      int transformationOrder = 0;
146      for (it = itb; it != ite; ++it)
147      {
148        listAlgos_.push_back(std::make_pair(axisPositionInGrid, std::make_pair(it->first, transformationOrder)));
149        algoTypes_.push_back(false);
150        ++transformationOrder;
151        std::vector<StdString> auxInput = (it->second)->checkAuxInputs();
152        for (int idx = 0; idx < auxInput.size(); ++idx) auxInputs_.push_back(auxInput[idx]);
153      }
154    }
155  }
156}
157
158/*!
159  Initialize the algorithms corresponding to transformation info contained in each domain.
160If a domain has transformations, they will be represented in form of vector of CTransformation pointers
161In general, each domain can have several transformations performed on itself.
162\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)
163*/
164void CGridTransformation::initializeDomainAlgorithms(int domPositionInGrid)
165{
166  std::vector<CDomain*> domListDestP = gridDestination_->getDomains();
167  if (!domListDestP.empty())
168  {
169    if (domListDestP[elementPosition2DomainPositionInGrid_[domPositionInGrid]]->hasTransformation())
170    {
171      CDomain::TransMapTypes trans = domListDestP[elementPosition2DomainPositionInGrid_[domPositionInGrid]]->getAllTransformations();
172      CDomain::TransMapTypes::const_iterator itb = trans.begin(), it,
173                                             ite = trans.end();
174      int transformationOrder = 0;
175      for (it = itb; it != ite; ++it)
176      {
177        listAlgos_.push_back(std::make_pair(domPositionInGrid, std::make_pair(it->first, transformationOrder)));
178        algoTypes_.push_back(true);
179        ++transformationOrder;
180        std::vector<StdString> auxInput = (it->second)->checkAuxInputs();
181        for (int idx = 0; idx < auxInput.size(); ++idx) auxInputs_.push_back(auxInput[idx]);
182      }
183    }
184  }
185
186}
187
188/*!
189  Select algorithm correspoding to its transformation type and its position in each element
190  \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)
191                                             and position of axis is 2
192  \param [in] transType transformation type, for now we have Zoom_axis, inverse_axis
193  \param [in] transformationOrder position of the transformation in an element (an element can have several transformation)
194  \param [in] isDomainAlgo flag to specify type of algorithm (for domain or axis)
195*/
196void CGridTransformation::selectAlgo(int elementPositionInGrid, ETranformationType transType, int transformationOrder, bool isDomainAlgo)
197{
198   if (isDomainAlgo) selectDomainAlgo(elementPositionInGrid, transType, transformationOrder);
199   else selectAxisAlgo(elementPositionInGrid, transType, transformationOrder);
200}
201
202/*!
203  Select algorithm of an axis correspoding to its transformation type and its position in each element
204  \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)
205                                             and position of axis is 2
206  \param [in] transType transformation type, for now we have Zoom_axis, inverse_axis
207  \param [in] transformationOrder position of the transformation in an element (an element can have several transformation)
208*/
209void CGridTransformation::selectAxisAlgo(int elementPositionInGrid, ETranformationType transType, int transformationOrder)
210{
211  std::vector<CAxis*> axisListDestP = gridDestination_->getAxis();
212  std::vector<CAxis*> axisListSrcP = gridSource_->getAxis();
213
214  int axisIndex =  elementPosition2AxisPositionInGrid_[elementPositionInGrid];
215  CAxis::TransMapTypes trans = axisListDestP[axisIndex]->getAllTransformations();
216  CAxis::TransMapTypes::const_iterator it = trans.begin();
217
218  for (int i = 0; i < transformationOrder; ++i, ++it) {}  // Find the correct transformation
219
220  CZoomAxis* zoomAxis = 0;
221  CInterpolateAxis* interpAxis = 0;
222  CGenericAlgorithmTransformation* algo = 0;
223  switch (transType)
224  {
225    case TRANS_INTERPOLATE_AXIS:
226      interpAxis = dynamic_cast<CInterpolateAxis*> (it->second);
227      algo = new CAxisAlgorithmInterpolate(axisListDestP[axisIndex], axisListSrcP[axisIndex], interpAxis);
228      break;
229    case TRANS_ZOOM_AXIS:
230      zoomAxis = dynamic_cast<CZoomAxis*> (it->second);
231      algo = new CAxisAlgorithmZoom(axisListDestP[axisIndex], axisListSrcP[axisIndex], zoomAxis);
232      break;
233    case TRANS_INVERSE_AXIS:
234      algo = new CAxisAlgorithmInverse(axisListDestP[axisIndex], axisListSrcP[axisIndex]);
235      break;
236    default:
237      break;
238  }
239  algoTransformation_.push_back(algo);
240
241}
242
243/*!
244  Select algorithm of a domain correspoding to its transformation type and its position in each element
245  \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)
246                                             and position of axis is 2
247  \param [in] transType transformation type, for now we have Zoom_axis, inverse_axis
248  \param [in] transformationOrder position of the transformation in an element (an element can have several transformation)
249*/
250void CGridTransformation::selectDomainAlgo(int elementPositionInGrid, ETranformationType transType, int transformationOrder)
251{
252  std::vector<CDomain*> domainListDestP = gridDestination_->getDomains();
253  std::vector<CDomain*> domainListSrcP = gridSource_->getDomains();
254
255  int domainIndex =  elementPosition2DomainPositionInGrid_[elementPositionInGrid];
256  CDomain::TransMapTypes trans = domainListDestP[domainIndex]->getAllTransformations();
257  CDomain::TransMapTypes::const_iterator it = trans.begin();
258
259  for (int i = 0; i < transformationOrder; ++i, ++it) {}  // Find the correct transformation
260
261  CZoomDomain* zoomDomain = 0;
262  CInterpolateDomain* interpFileDomain = 0;
263  CGenericAlgorithmTransformation* algo = 0;
264  switch (transType)
265  {
266    case TRANS_INTERPOLATE_DOMAIN:
267      interpFileDomain = dynamic_cast<CInterpolateDomain*> (it->second);
268      algo = new CDomainAlgorithmInterpolate(domainListDestP[domainIndex], domainListSrcP[domainIndex],interpFileDomain);
269      break;
270    case TRANS_ZOOM_DOMAIN:
271      zoomDomain = dynamic_cast<CZoomDomain*> (it->second);
272      algo = new CDomainAlgorithmZoom(domainListDestP[domainIndex], domainListSrcP[domainIndex], zoomDomain);
273      break;
274    default:
275      break;
276  }
277  algoTransformation_.push_back(algo);
278}
279
280/*!
281  Assign the current grid destination to the grid source in the new transformation.
282The current grid destination plays the role of grid source in next transformation (if any).
283Only element on which the transformation is performed is modified
284  \param [in] elementPositionInGrid position of element in grid
285  \param [in] transType transformation type
286*/
287void CGridTransformation::setUpGrid(int elementPositionInGrid, ETranformationType transType, int nbTransformation)
288{
289  if (!tempGrids_.empty() && (getNbAlgo()-1) == tempGrids_.size())
290  {
291    gridSource_ = tempGrids_[nbTransformation];
292    return;
293  }
294
295  std::vector<CAxis*> axisListDestP = gridDestination_->getAxis();
296  std::vector<CAxis*> axisListSrcP = gridSource_->getAxis(), axisSrc;
297
298  std::vector<CDomain*> domListDestP = gridDestination_->getDomains();
299  std::vector<CDomain*> domListSrcP = gridSource_->getDomains(), domainSrc;
300
301  int axisIndex = -1, domainIndex = -1;
302  switch (transType)
303  {
304    case TRANS_INTERPOLATE_DOMAIN:
305    case TRANS_ZOOM_DOMAIN:
306      domainIndex = elementPosition2DomainPositionInGrid_[elementPositionInGrid];
307      break;
308
309    case TRANS_INTERPOLATE_AXIS:
310    case TRANS_ZOOM_AXIS:
311    case TRANS_INVERSE_AXIS:
312      axisIndex =  elementPosition2AxisPositionInGrid_[elementPositionInGrid];
313      break;
314    default:
315      break;
316  }
317
318  for (int idx = 0; idx < axisListSrcP.size(); ++idx)
319  {
320    CAxis* axis = CAxis::createAxis();
321    if (axisIndex != idx) axis->axis_ref.setValue(axisListSrcP[idx]->getId());
322    else axis->axis_ref.setValue(axisListDestP[idx]->getId());
323    axis->solveRefInheritance(true);
324    axis->checkAttributesOnClient();
325    axisSrc.push_back(axis);
326  }
327
328  for (int idx = 0; idx < domListSrcP.size(); ++idx)
329  {
330    CDomain* domain = CDomain::createDomain();
331    if (domainIndex != idx) domain->domain_ref.setValue(domListSrcP[idx]->getId());
332    else domain->domain_ref.setValue(domListDestP[idx]->getId());
333    domain->solveRefInheritance(true);
334    domain->checkAttributesOnClient();
335    domainSrc.push_back(domain);
336  }
337
338  gridSource_ = CGrid::createGrid(domainSrc, axisSrc, gridDestination_->axis_domain_order);
339  gridSource_->computeGridGlobalDimension(domainSrc, axisSrc, gridDestination_->axis_domain_order);
340
341  tempGrids_.push_back(gridSource_);
342}
343
344/*!
345  Perform all transformations
346  For each transformation, there are some things to do:
347  -) Chose the correct algorithm by transformation type and position of element
348  -) Calculate the mapping of global index between the current grid source and grid destination
349  -) Calculate the mapping of global index between current grid DESTINATION and grid SOURCE
350  -) Make current grid destination become grid source in the next transformation
351*/
352void CGridTransformation::computeAll(const std::vector<CArray<double,1>* >& dataAuxInputs, Time timeStamp)
353{
354  if (nbAlgos_ < 1) return;
355  if (!auxInputs_.empty() && !dynamicalTransformation_) { dynamicalTransformation_ = true; return; }
356  if (dynamicalTransformation_)
357  {
358    if (timeStamp_.insert(timeStamp).second)   //Reset map
359    {
360      std::list<SendingIndexGridSourceMap>().swap(localIndexToSendFromGridSource_);
361      std::list<RecvIndexGridDestinationMap>().swap(localIndexToReceiveOnGridDest_);
362      std::list<size_t>().swap(nbLocalIndexOnGridDest_);
363    }
364    else
365      return;
366  }
367
368  CContext* context = CContext::getCurrent();
369  CContextClient* client = context->client;
370
371  ListAlgoType::const_iterator itb = listAlgos_.begin(),
372                               ite = listAlgos_.end(), it;
373
374  CGenericAlgorithmTransformation* algo = 0;
375  int nbAgloTransformation = 0; // Only count for executed transformation. Generate domain is a special one, not executed in the list
376  for (it = itb; it != ite; ++it)
377  {
378    int elementPositionInGrid = it->first;
379    ETranformationType transType = (it->second).first;
380    int transformationOrder = (it->second).second;
381    SourceDestinationIndexMap globaIndexWeightFromSrcToDst;
382
383    // First of all, select an algorithm
384    if (!dynamicalTransformation_ || (algoTransformation_.size() < listAlgos_.size()))
385    {
386      selectAlgo(elementPositionInGrid, transType, transformationOrder, algoTypes_[std::distance(itb, it)]);
387      algo = algoTransformation_.back();
388    }
389    else
390      algo = algoTransformation_[std::distance(itb, it)];
391
392    if (0 != algo) // Only registered transformation can be executed
393    {
394      algo->computeIndexSourceMapping(dataAuxInputs);
395
396      // ComputeTransformation of global index of each element
397      int elementPosition = it->first;
398      algo->computeGlobalSourceIndex(elementPosition,
399                                     gridSource_,
400                                     gridDestination_,
401                                     globaIndexWeightFromSrcToDst);
402
403      // Compute transformation of global indexes among grids
404      computeTransformationMapping(globaIndexWeightFromSrcToDst);
405
406      if (1 < nbAlgos_)
407      {
408        // Now grid destination becomes grid source in a new transformation
409        if (nbAgloTransformation != (nbAlgos_-1)) setUpGrid(elementPositionInGrid, transType, nbAgloTransformation);
410      }
411      ++nbAgloTransformation;
412    }
413  }
414}
415
416/*!
417  Compute exchange index between grid source and grid destination
418  \param [in] globalIndexWeightFromDestToSource global index mapping between grid destination and grid source
419*/
420void CGridTransformation::computeTransformationMapping(const SourceDestinationIndexMap& globaIndexWeightFromSrcToDst)
421{
422  CContext* context = CContext::getCurrent();
423  CContextClient* client = context->client;
424  int nbClient = client->clientSize;
425  int clientRank = client->clientRank;
426
427  // Recalculate the distribution of grid destination
428  CDistributionClient distributionClientDest(client->clientRank, gridDestination_);
429  CDistributionClient::GlobalLocalDataMap& globalLocalIndexGridDestSendToServer = distributionClientDest.getGlobalLocalDataSendToServer();
430  // Update number of local index on each transformation
431  nbLocalIndexOnGridDest_.push_back(globalLocalIndexGridDestSendToServer.size());
432
433  // Find out number of index sent from grid source and number of index received on grid destination
434  SourceDestinationIndexMap::const_iterator itbIndex = globaIndexWeightFromSrcToDst.begin(),
435                                            iteIndex = globaIndexWeightFromSrcToDst.end(), itIndex;
436  typedef boost::unordered_map<size_t, std::vector<std::pair<size_t,double> > > SendIndexMap;
437  std::map<int,int> sendRankSizeMap,recvRankSizeMap;
438  int connectedClient = globaIndexWeightFromSrcToDst.size();
439  int* recvCount=new int[nbClient];
440  int* displ=new int[nbClient];
441  int* sendRankBuff=new int[connectedClient];
442  int* sendSizeBuff=new int[connectedClient];
443  int n = 0;
444  for (itIndex = itbIndex; itIndex != iteIndex; ++itIndex, ++n)
445  {
446    sendRankBuff[n] = itIndex->first;
447    const SendIndexMap& sendIndexMap = itIndex->second;
448    SendIndexMap::const_iterator itbSend = sendIndexMap.begin(), iteSend = sendIndexMap.end(), itSend;
449    int sendSize = 0;
450    for (itSend = itbSend; itSend != iteSend; ++itSend)
451    {
452      sendSize += itSend->second.size();
453    }
454    sendSizeBuff[n] = sendSize;
455    sendRankSizeMap[itIndex->first] = sendSize;
456  }
457  MPI_Allgather(&connectedClient,1,MPI_INT,recvCount,1,MPI_INT,client->intraComm);
458
459  displ[0]=0 ;
460  for(int n=1;n<nbClient;n++) displ[n]=displ[n-1]+recvCount[n-1];
461  int recvSize=displ[nbClient-1]+recvCount[nbClient-1];
462  int* recvRankBuff=new int[recvSize];
463  int* recvSizeBuff=new int[recvSize];
464  MPI_Allgatherv(sendRankBuff,connectedClient,MPI_INT,recvRankBuff,recvCount,displ,MPI_INT,client->intraComm);
465  MPI_Allgatherv(sendSizeBuff,connectedClient,MPI_INT,recvSizeBuff,recvCount,displ,MPI_INT,client->intraComm);
466  for (int i = 0; i < nbClient; ++i)
467  {
468    int currentPos = displ[i];
469    for (int j = 0; j < recvCount[i]; ++j)
470      if (recvRankBuff[currentPos+j] == clientRank)
471      {
472        recvRankSizeMap[i] = recvSizeBuff[currentPos+j];
473      }
474  }
475
476  // Sending global index of grid source to corresponding process as well as the corresponding mask
477  std::vector<MPI_Request> requests;
478  std::vector<MPI_Status> status;
479  boost::unordered_map<int, unsigned char* > recvMaskDst;
480  boost::unordered_map<int, unsigned long* > recvGlobalIndexSrc;
481  for (std::map<int,int>::const_iterator itRecv = recvRankSizeMap.begin(); itRecv != recvRankSizeMap.end(); ++itRecv)
482  {
483    int recvRank = itRecv->first;
484    int recvSize = itRecv->second;
485    recvMaskDst[recvRank] = new unsigned char [recvSize];
486    recvGlobalIndexSrc[recvRank] = new unsigned long [recvSize];
487
488    requests.push_back(MPI_Request());
489    MPI_Irecv(recvGlobalIndexSrc[recvRank], recvSize, MPI_UNSIGNED_LONG, recvRank, 46, client->intraComm, &requests.back());
490    requests.push_back(MPI_Request());
491    MPI_Irecv(recvMaskDst[recvRank], recvSize, MPI_UNSIGNED_CHAR, recvRank, 47, client->intraComm, &requests.back());
492  }
493
494  boost::unordered_map<int, CArray<size_t,1> > globalIndexDst;
495  boost::unordered_map<int, CArray<double,1> > weightDst;
496  boost::unordered_map<int, unsigned char* > sendMaskDst;
497  boost::unordered_map<int, unsigned long* > sendGlobalIndexSrc;
498  for (itIndex = itbIndex; itIndex != iteIndex; ++itIndex)
499  {
500    int sendRank = itIndex->first;
501    int sendSize = sendRankSizeMap[sendRank];
502    const SendIndexMap& sendIndexMap = itIndex->second;
503    SendIndexMap::const_iterator itbSend = sendIndexMap.begin(), iteSend = sendIndexMap.end(), itSend;
504    globalIndexDst[sendRank].resize(sendSize);
505    weightDst[sendRank].resize(sendSize);
506    sendMaskDst[sendRank] = new unsigned char [sendSize];
507    sendGlobalIndexSrc[sendRank] = new unsigned long [sendSize];
508    int countIndex = 0;
509    for (itSend = itbSend; itSend != iteSend; ++itSend)
510    {
511      const std::vector<std::pair<size_t,double> >& dstWeight = itSend->second;
512      for (int idx = 0; idx < dstWeight.size(); ++idx)
513      {
514        globalIndexDst[sendRank](countIndex) = dstWeight[idx].first;
515        weightDst[sendRank](countIndex) = dstWeight[idx].second;
516        if (0 < globalLocalIndexGridDestSendToServer.count(dstWeight[idx].first))
517          sendMaskDst[sendRank][countIndex] = 1;
518        else
519          sendMaskDst[sendRank][countIndex] = 0;
520        sendGlobalIndexSrc[sendRank][countIndex] = itSend->first;
521        ++countIndex;
522      }
523    }
524
525    // Send global index source and mask
526    requests.push_back(MPI_Request());
527    MPI_Isend(sendGlobalIndexSrc[sendRank], sendSize, MPI_UNSIGNED_LONG, sendRank, 46, client->intraComm, &requests.back());
528    requests.push_back(MPI_Request());
529    MPI_Isend(sendMaskDst[sendRank], sendSize, MPI_UNSIGNED_CHAR, sendRank, 47, client->intraComm, &requests.back());
530  }
531
532  status.resize(requests.size());
533  MPI_Waitall(requests.size(), &requests[0], &status[0]);
534
535  // Okie, now use the mask to identify which index source we need to send, then also signal the destination which masked index we will return
536  std::vector<MPI_Request>().swap(requests);
537  std::vector<MPI_Status>().swap(status);
538  // Okie, on destination side, we will wait for information of masked index of source
539  for (std::map<int,int>::const_iterator itSend = sendRankSizeMap.begin(); itSend != sendRankSizeMap.end(); ++itSend)
540  {
541    int recvRank = itSend->first;
542    int recvSize = itSend->second;
543
544    requests.push_back(MPI_Request());
545    MPI_Irecv(sendMaskDst[recvRank], recvSize, MPI_UNSIGNED_CHAR, recvRank, 48, client->intraComm, &requests.back());
546  }
547
548  // Ok, now we fill in local index of grid source (we even count for masked index)
549  CDistributionClient distributionClientSrc(client->clientRank, gridSource_);
550  CDistributionClient::GlobalLocalDataMap& globalLocalIndexGridSrcSendToServer = distributionClientSrc.getGlobalLocalDataSendToServer();
551  localIndexToSendFromGridSource_.push_back(SendingIndexGridSourceMap());
552  SendingIndexGridSourceMap& tmpSend = localIndexToSendFromGridSource_.back();
553  for (std::map<int,int>::const_iterator itRecv = recvRankSizeMap.begin(); itRecv != recvRankSizeMap.end(); ++itRecv)
554  {
555    int recvRank = itRecv->first;
556    int recvSize = itRecv->second;
557    unsigned char* recvMask = recvMaskDst[recvRank];
558    unsigned long* recvIndexSrc = recvGlobalIndexSrc[recvRank];
559    int realSendSize = 0;
560    for (int idx = 0; idx < recvSize; ++idx)
561    {
562      if (0 != (*(recvMask+idx))) // OKie, now we have a demand from non-masked index destination
563        if (0 < globalLocalIndexGridSrcSendToServer.count(*(recvIndexSrc+idx))) // check whether index source is masked
564         ++realSendSize;
565        else // inform the destination that this index is masked
566         *(recvMask+idx) = 0;
567    }
568
569    tmpSend[recvRank].resize(realSendSize);
570    realSendSize = 0;
571    for (int idx = 0; idx < recvSize; ++idx)
572    {
573      if (0 != (*(recvMask+idx))) // OKie, now we have a demand from non-masked index destination
574      {
575        tmpSend[recvRank](realSendSize) = globalLocalIndexGridSrcSendToServer[*(recvIndexSrc+idx)];
576         ++realSendSize;
577      }
578    }
579
580    // Okie, now inform the destination which source index are masked
581    requests.push_back(MPI_Request());
582    MPI_Isend(recvMaskDst[recvRank], recvSize, MPI_UNSIGNED_CHAR, recvRank, 48, client->intraComm, &requests.back());
583  }
584  status.resize(requests.size());
585  MPI_Waitall(requests.size(), &requests[0], &status[0]);
586
587  // Cool, now we can fill in local index of grid destination (counted for masked index)
588  localIndexToReceiveOnGridDest_.push_back(RecvIndexGridDestinationMap());
589  RecvIndexGridDestinationMap& recvTmp = localIndexToReceiveOnGridDest_.back();
590  for (std::map<int,int>::const_iterator itSend = sendRankSizeMap.begin(); itSend != sendRankSizeMap.end(); ++itSend)
591  {
592    int recvRank = itSend->first;
593    int recvSize = itSend->second;
594    unsigned char* recvMask = sendMaskDst[recvRank];
595
596    CArray<size_t,1>& recvIndexDst = globalIndexDst[recvRank];
597    CArray<double,1>& recvWeightDst = weightDst[recvRank];
598    int realRecvSize = 0;
599    for (int idx = 0; idx < recvSize; ++idx)
600    {
601      if (0 != *(recvMask+idx)) // OKie, now we have a non-masked index destination
602         ++realRecvSize;
603    }
604
605    int localIndexDst;
606    recvTmp[recvRank].resize(realRecvSize);
607    realRecvSize = 0;
608    for (int idx = 0; idx < recvSize; ++idx)
609    {
610      if (0 != *(recvMask+idx)) // OKie, now we have a demand from non-masked index destination
611      {
612        recvTmp[recvRank][realRecvSize].first = globalLocalIndexGridDestSendToServer[recvIndexDst(idx)];
613        recvTmp[recvRank][realRecvSize].second = recvWeightDst(idx);
614         ++realRecvSize;
615      }
616    }
617  }
618
619  delete [] recvCount;
620  delete [] displ;
621  delete [] sendRankBuff;
622  delete [] recvRankBuff;
623  delete [] sendSizeBuff;
624  delete [] recvSizeBuff;
625
626  boost::unordered_map<int, unsigned char* >::const_iterator itChar;
627  for (itChar = sendMaskDst.begin(); itChar != sendMaskDst.end(); ++itChar)
628    delete [] itChar->second;
629  for (itChar = recvMaskDst.begin(); itChar != recvMaskDst.end(); ++itChar)
630    delete [] itChar->second;
631  boost::unordered_map<int, unsigned long* >::const_iterator itLong;
632  for (itLong = sendGlobalIndexSrc.begin(); itLong != sendGlobalIndexSrc.end(); ++itLong)
633    delete [] itLong->second;
634  for (itLong = recvGlobalIndexSrc.begin(); itLong != recvGlobalIndexSrc.end(); ++itLong)
635    delete [] itLong->second;
636
637}
638
639bool CGridTransformation::isSpecialTransformation(ETranformationType transType)
640{
641  bool res;
642  switch (transType)
643  {
644    case TRANS_GENERATE_RECTILINEAR_DOMAIN:
645     res = true;
646     break;
647    default:
648     res = false;
649     break;
650  }
651
652  return res;
653}
654
655/*!
656  Local index of data which need sending from the grid source
657  \return local index of data
658*/
659const std::list<CGridTransformation::SendingIndexGridSourceMap>& CGridTransformation::getLocalIndexToSendFromGridSource() const
660{
661  return localIndexToSendFromGridSource_;
662}
663
664/*!
665  Local index of data which will be received on the grid destination
666  \return local index of data
667*/
668const std::list<CGridTransformation::RecvIndexGridDestinationMap>& CGridTransformation::getLocalIndexToReceiveOnGridDest() const
669{
670  return localIndexToReceiveOnGridDest_;
671}
672
673const std::list<size_t>& CGridTransformation::getNbLocalIndexToReceiveOnGridDest() const
674{
675  return nbLocalIndexOnGridDest_;
676}
677
678}
Note: See TracBrowser for help on using the repository browser.