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

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

Adding rectilinear and curvilinear domain for expand_domain transformation
-) Rectilinear/curvilinear is expanded not only locally but also globally, its global size ni_glo, nj_glo become ni_glo+2 and nj_glo+2
-) Two attributes i_periodic, j_periodic are only used for rectilinear/curvilinear to process priodic condition

+) Do some minor modification

Test
+) Add test_connectivity_expand
+) On Curie
+) Work (but need more real tests)

  • Property svn:executable set to *
File size: 28.7 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 "grid_transformation_factory_impl.hpp"
11#include "algo_types.hpp"
12#include "context.hpp"
13#include "context_client.hpp"
14#include "distribution_client.hpp"
15#include "mpi_tag.hpp"
16#include "grid.hpp"
17#include <boost/unordered_map.hpp>
18
19namespace xios {
20CGridTransformation::CGridTransformation(CGrid* destination, CGrid* source)
21 : CGridTransformationSelector(destination, source),
22  tmpGridDestination_(destination), originalGridSource_(source),
23  tempGridSrcs_(), tempGridDests_(),
24  dynamicalTransformation_(false), timeStamp_()
25{
26}
27
28CGridTransformation::~CGridTransformation()
29{
30}
31
32/*!
33  Select algorithm of a scalar corresponding to its transformation type and its position in each element
34  \param [in] elementPositionInGrid position of element in grid. e.g: a grid has 1 domain and 1 axis, then position of domain is 0 and position of axis is 1
35  \param [in] transType transformation type, for now we have
36  \param [in] transformationOrder position of the transformation in an element (an element can have several transformation)
37*/
38void CGridTransformation::selectScalarAlgo(int elementPositionInGrid, ETranformationType transType, int transformationOrder)
39{
40  std::vector<CScalar*> scaListDestP = gridDestination_->getScalars();
41  int scalarDstIndex =  elementPositionInGridDst2ScalarPosition_[elementPositionInGrid];
42  CScalar::TransMapTypes trans = scaListDestP[scalarDstIndex]->getAllTransformations();
43  CScalar::TransMapTypes::const_iterator it = trans.begin();
44
45  for (int i = 0; i < transformationOrder; ++i, ++it) {}  // Find the correct transformation
46  CGenericAlgorithmTransformation* algo = 0;
47  algo = CGridTransformationFactory<CScalar>::createTransformation(transType,
48                                                                  gridDestination_,
49                                                                  gridSource_,
50                                                                  it->second,
51                                                                  elementPositionInGrid,
52                                                                  elementPositionInGridSrc2ScalarPosition_,
53                                                                  elementPositionInGridSrc2AxisPosition_,
54                                                                  elementPositionInGridSrc2DomainPosition_,
55                                                                  elementPositionInGridDst2ScalarPosition_,
56                                                                  elementPositionInGridDst2AxisPosition_,
57                                                                  elementPositionInGridDst2DomainPosition_);
58  algoTransformation_.push_back(algo);
59}
60
61/*!
62  Select algorithm of an axis corresponding to its transformation type and its position in each element
63  \param [in] elementPositionInGrid position of element in grid. e.g: a grid has 1 domain and 1 axis, then position of domain is 0 and position of axis is 1
64  \param [in] transType transformation type, for now we have zoom_axis, inverse_axis, interpolate_axis
65  \param [in] transformationOrder position of the transformation in an element (an element can have several transformation)
66*/
67void CGridTransformation::selectAxisAlgo(int elementPositionInGrid, ETranformationType transType, int transformationOrder)
68{
69  std::vector<CAxis*> axisListDestP = gridDestination_->getAxis();
70  int axisDstIndex =  elementPositionInGridDst2AxisPosition_[elementPositionInGrid];
71  CAxis::TransMapTypes trans = axisListDestP[axisDstIndex]->getAllTransformations();
72  CAxis::TransMapTypes::const_iterator it = trans.begin();
73  for (int i = 0; i < transformationOrder; ++i, ++it) {}  // Find the correct transformation
74
75  CGenericAlgorithmTransformation* algo = 0;
76  algo = CGridTransformationFactory<CAxis>::createTransformation(transType,
77                                                                 gridDestination_,
78                                                                 gridSource_,
79                                                                 it->second,
80                                                                 elementPositionInGrid,
81                                                                 elementPositionInGridSrc2ScalarPosition_,
82                                                                 elementPositionInGridSrc2AxisPosition_,
83                                                                 elementPositionInGridSrc2DomainPosition_,
84                                                                 elementPositionInGridDst2ScalarPosition_,
85                                                                 elementPositionInGridDst2AxisPosition_,
86                                                                 elementPositionInGridDst2DomainPosition_);
87  algoTransformation_.push_back(algo);
88}
89
90/*!
91  Select algorithm of a domain corresponding to its transformation type and its position in each element
92  \param [in] elementPositionInGrid position of element in grid. e.g: a grid has 1 domain and 1 axis, then position of domain is 0 and position of axis is 1
93  \param [in] transType transformation type, for now we have zoom_domain, interpolate_domain
94  \param [in] transformationOrder position of the transformation in an element (an element can have several transformation)
95*/
96void CGridTransformation::selectDomainAlgo(int elementPositionInGrid, ETranformationType transType, int transformationOrder)
97{
98  std::vector<CDomain*> domainListDestP = gridDestination_->getDomains();
99  int domainIndex =  elementPositionInGridDst2DomainPosition_[elementPositionInGrid];
100  CDomain::TransMapTypes trans = domainListDestP[domainIndex]->getAllTransformations();
101  CDomain::TransMapTypes::const_iterator it = trans.begin();
102  for (int i = 0; i < transformationOrder; ++i, ++it) {}  // Find the correct transformation 
103
104  CGenericAlgorithmTransformation* algo = 0;
105  algo = CGridTransformationFactory<CDomain>::createTransformation(transType,
106                                                                   gridDestination_,
107                                                                   gridSource_,
108                                                                   it->second,
109                                                                   elementPositionInGrid,
110                                                                   elementPositionInGridSrc2ScalarPosition_,
111                                                                   elementPositionInGridSrc2AxisPosition_,
112                                                                   elementPositionInGridSrc2DomainPosition_,
113                                                                   elementPositionInGridDst2ScalarPosition_,
114                                                                   elementPositionInGridDst2AxisPosition_,
115                                                                   elementPositionInGridDst2DomainPosition_);
116  algoTransformation_.push_back(algo);
117}
118
119/*!
120  Find position of element in a grid as well as its type (domain, axis, scalar) and position in its own element list
121  \return element position: map<int,<int,int> > corresponds to <element position in grid, <element type, element position in element list> >
122*/
123std::map<int,std::pair<int,int> > CGridTransformation::getElementPosition(CGrid* grid)
124{
125  std::vector<CScalar*> scalarListP = grid->getScalars(); 
126  std::vector<CAxis*> axisListP = grid->getAxis();
127  std::vector<CDomain*> domListP = grid->getDomains(); 
128  CArray<int,1> axisDomainOrder = grid->axis_domain_order;
129  int scalarIndex = 0, axisIndex = 0, domainIndex = 0;
130  int nbElement = axisDomainOrder.numElements(), elementDim;
131  std::map<int,std::pair<int,int> > elementPosition;
132  for (int idx = 0; idx < nbElement; ++idx)
133  {
134    elementDim = axisDomainOrder(idx);
135    switch (elementDim)
136    {
137      case 2:
138        elementPosition[idx] = std::make_pair(elementDim, domainIndex);
139        ++domainIndex;
140        break;
141      case 1:
142        elementPosition[idx] = std::make_pair(elementDim, axisIndex);
143        ++axisIndex;
144        break;
145      case 0:
146        elementPosition[idx] = std::make_pair(elementDim, scalarIndex);
147        ++scalarIndex;
148        break;
149      default:
150        break;       
151    }
152  }
153
154  return elementPosition; 
155}
156
157/*!
158  If there are more than one transformation, a new temporary grid will be created and it will play the role of grid destination.
159This new created one keeps a pointer to the real transformed element of grid destination and generate new copies of other elements from grid source.
160  \param [in] elementPositionInGrid position of element in grid
161  \param [in] transType transformation type
162*/
163void CGridTransformation::setUpGridDestination(int elementPositionInGrid, ETranformationType transType)
164{
165  if (isSpecialTransformation(transType)) return;
166
167  if (!tempGridDests_.empty() && (getNbAlgo() == tempGridDests_.size()))
168  {
169    tempGridDests_.resize(0);
170  }
171
172  if (1 == getNbAlgo()) 
173  {
174    tmpGridDestination_ = gridDestination_;
175    return;
176  }
177
178  std::vector<CScalar*> scalarListDestP = gridDestination_->getScalars();
179  std::vector<CScalar*> scalarListSrcP = gridSource_->getScalars(), scalarDst;
180
181  std::vector<CAxis*> axisListDestP = gridDestination_->getAxis();
182  std::vector<CAxis*> axisListSrcP = gridSource_->getAxis(), axisDst;
183
184  std::vector<CDomain*> domListDestP = gridDestination_->getDomains();
185  std::vector<CDomain*> domListSrcP = gridSource_->getDomains(), domainDst;
186
187  CArray<int,1> axisDomainOrderSrc = gridSource_->axis_domain_order;
188  CArray<int,1> axisDomainOrderDst = gridDestination_->axis_domain_order;
189
190  std::map<int,std::pair<int,int> > elementPositionSrc = getElementPosition(gridSource_);
191  std::map<int,std::pair<int,int> > elementPositionDst = getElementPosition(gridDestination_);
192
193  CArray<int,1> elementOrder(axisDomainOrderDst.numElements());
194  for (int idx = 0; idx < axisDomainOrderDst.numElements(); ++idx)
195  {
196    if (elementPositionInGrid == idx)
197    {
198      int dimElementDst = elementPositionDst[idx].first;
199      int elementIndex  = elementPositionDst[idx].second;
200      switch (dimElementDst) 
201      {
202        case 2:
203          domainDst.push_back(domListDestP[elementIndex]);
204          break;
205        case 1:
206          axisDst.push_back(axisListDestP[elementIndex]);
207          break;
208        case 0:
209          scalarDst.push_back(scalarListDestP[elementIndex]);
210          break;
211        default:
212          break;
213      }
214      elementOrder(idx) = dimElementDst;
215    }
216    else
217    {     
218      int dimElementSrc = elementPositionSrc[idx].first;
219      int elementIndex  = elementPositionSrc[idx].second;
220      switch (dimElementSrc)
221      {
222        case 2:
223          domainDst.push_back(domListSrcP[elementIndex]);
224          break;
225        case 1:
226          axisDst.push_back(axisListSrcP[elementIndex]);
227          break;
228        case 0:
229          scalarDst.push_back(scalarListSrcP[elementIndex]);
230          break;
231        default:
232          break;
233      }
234      elementOrder(idx) = dimElementSrc;
235    }
236  }
237
238  tmpGridDestination_ = CGrid::createGrid(domainDst, axisDst, scalarDst, elementOrder); 
239  tempGridDests_.push_back(tmpGridDestination_);
240}
241
242/*!
243  Assign the current grid destination to the grid source in the new transformation.
244The current grid destination plays the role of grid source in next transformation (if any).
245Only element on which the transformation is performed is modified
246  \param [in] elementPositionInGrid position of element in grid
247  \param [in] transType transformation type
248*/
249void CGridTransformation::setUpGridSource(int elementPositionInGrid)
250{
251  if (!tempGridSrcs_.empty() && (getNbAlgo()-1) == tempGridSrcs_.size())
252  {
253    tempGridSrcs_.resize(0);
254  }
255
256  std::vector<CScalar*> scalarListDestP = tmpGridDestination_->getScalars();
257  std::vector<CScalar*> scalarListSrcP = gridSource_->getScalars(), scalarSrc;
258
259  std::vector<CAxis*> axisListDestP = tmpGridDestination_->getAxis();
260  std::vector<CAxis*> axisListSrcP = gridSource_->getAxis(), axisSrc;
261
262  std::vector<CDomain*> domListDestP = tmpGridDestination_->getDomains();
263  std::vector<CDomain*> domListSrcP = gridSource_->getDomains(), domainSrc;
264
265  CArray<int,1> axisDomainOrderSrc = gridSource_->axis_domain_order;
266  CArray<int,1> axisDomainOrderDst = tmpGridDestination_->axis_domain_order;
267
268  std::map<int,std::pair<int,int> > elementPositionSrc = getElementPosition(gridSource_);
269  std::map<int,std::pair<int,int> > elementPositionDst = getElementPosition(tmpGridDestination_);
270
271  for (int idx = 0; idx < axisDomainOrderDst.numElements(); ++idx)
272  {   
273    if (elementPositionInGrid == idx)
274    {
275      int dimElementDst = elementPositionDst[idx].first;
276      int elementIndex  = elementPositionDst[idx].second;
277      if (2 == dimElementDst)
278      {
279        CDomain* domain = CDomain::createDomain();
280        domain->domain_ref.setValue(domListDestP[elementIndex]->getId());
281        domain->solveRefInheritance(true);
282        domain->checkAttributesOnClient();
283        domainSrc.push_back(domain);
284      }
285      else if (1 == dimElementDst)
286      {
287        CAxis* axis = CAxis::createAxis();
288        axis->axis_ref.setValue(axisListDestP[elementIndex]->getId());
289        axis->solveRefInheritance(true);
290        axis->checkAttributesOnClient();
291        axisSrc.push_back(axis); 
292      }
293      else
294      {
295        CScalar* scalar = CScalar::createScalar();
296        scalar->scalar_ref.setValue(scalarListDestP[elementIndex]->getId());
297        scalar->solveRefInheritance(true);
298        scalar->checkAttributesOnClient();
299        scalarSrc.push_back(scalar);
300      }
301    }
302    else
303    {     
304      int dimElementDst = elementPositionDst[idx].first;
305      int elementIndex  = elementPositionDst[idx].second;
306      switch (dimElementDst)
307      {
308        case 2:
309          domainSrc.push_back(domListDestP[elementIndex]);
310          break;
311        case 1:
312          axisSrc.push_back(axisListDestP[elementIndex]);
313          break;
314        case 0:
315          scalarSrc.push_back(scalarListDestP[elementIndex]);
316          break;
317        default:
318          break;
319      }
320    }
321  }
322
323  gridSource_ = CGrid::createGrid(domainSrc, axisSrc, scalarSrc, tmpGridDestination_->axis_domain_order); 
324
325  tempGridSrcs_.push_back(gridSource_);
326}
327
328/*!
329  Perform all transformations
330  For each transformation, there are some things to do:
331  -) Chose the correct algorithm by transformation type and position of element
332  -) Calculate the mapping of global index between the current grid source and grid destination
333  -) Calculate the mapping of global index between current grid DESTINATION and grid SOURCE
334  -) Make current grid destination become grid source in the next transformation
335*/
336void CGridTransformation::computeAll(const std::vector<CArray<double,1>* >& dataAuxInputs, Time timeStamp)
337{
338  if (nbNormalAlgos_ < 1) return;
339  if (!auxInputs_.empty() && !dynamicalTransformation_) { dynamicalTransformation_ = true; return; }
340  if (dynamicalTransformation_)
341  {
342    if (timeStamp_.insert(timeStamp).second)   //Reset map
343    {
344      std::list<SendingIndexGridSourceMap>().swap(localIndexToSendFromGridSource_);
345      std::list<RecvIndexGridDestinationMap>().swap(localIndexToReceiveOnGridDest_);
346      std::list<size_t>().swap(nbLocalIndexOnGridDest_);
347      std::list<std::vector<bool> >().swap(localMaskOnGridDest_);
348    }
349    else
350      return;
351  }
352
353  CContext* context = CContext::getCurrent();
354  CContextClient* client = context->client;
355
356  ListAlgoType::const_iterator itb = listAlgos_.begin(),
357                               ite = listAlgos_.end(), it;
358
359  CGenericAlgorithmTransformation* algo = 0;
360  int nbAgloTransformation = 0; // Only count for executed transformation. Generate domain is a special one, not executed in the list
361  for (it = itb; it != ite; ++it)
362  {
363    int elementPositionInGrid = it->first;
364    ETranformationType transType = (it->second).first;
365    int transformationOrder = (it->second).second.first;
366    int algoType = ((it->second).second.second); //algoTypes_[std::distance(itb, it)];
367    SourceDestinationIndexMap globaIndexWeightFromSrcToDst;
368   
369
370    // Create a temporary grid destination which contains transformed element of grid destination and
371    // non-transformed elements to grid source
372    setUpGridDestination(elementPositionInGrid, transType);
373
374    // First of all, select an algorithm
375    if (!dynamicalTransformation_ || (algoTransformation_.size() < listAlgos_.size()))
376    {
377      selectAlgo(elementPositionInGrid, transType, transformationOrder, algoType);
378      algo = algoTransformation_.back();
379    }
380    else
381      algo = algoTransformation_[std::distance(itb, it)];
382
383    if ((0 != algo) &&
384        ((CGenericAlgorithmTransformation::ELEMENT_NO_MODIFICATION_WITH_DATA == algo->type()) ||
385        (CGenericAlgorithmTransformation::ELEMENT_MODIFICATION_WITH_DATA == algo->type()))) // Only registered transformation can be executed
386    {
387      algo->computeIndexSourceMapping(dataAuxInputs);
388
389      // ComputeTransformation of global index of each element
390      int elementPosition = it->first;
391      algo->computeGlobalSourceIndex(elementPosition,
392                                     gridSource_,
393                                     tmpGridDestination_,
394                                     globaIndexWeightFromSrcToDst);
395
396      // Compute transformation of global indexes among grids
397      computeTransformationMapping(globaIndexWeightFromSrcToDst);
398
399      if (1 < nbNormalAlgos_)
400      {
401        // Now grid destination becomes grid source in a new transformation
402        if (nbAgloTransformation != (nbNormalAlgos_-1)) setUpGridSource(elementPositionInGrid);
403      }
404      ++nbAgloTransformation;
405    }
406  }
407}
408
409/*!
410  Compute exchange index between grid source and grid destination
411  \param [in] globalIndexWeightFromDestToSource global index mapping between grid destination and grid source
412*/
413void CGridTransformation::computeTransformationMapping(const SourceDestinationIndexMap& globaIndexWeightFromSrcToDst)
414{
415  CContext* context = CContext::getCurrent();
416  CContextClient* client = context->client;
417  int nbClient = client->clientSize;
418  int clientRank = client->clientRank;
419
420  // Recalculate the distribution of grid destination
421  CDistributionClient distributionClientDest(client->clientRank, tmpGridDestination_);
422  CDistributionClient::GlobalLocalDataMap& globalLocalIndexGridDestSendToServer = distributionClientDest.getGlobalLocalDataSendToServer();
423
424  // Update number of local index on each transformation
425  size_t nbLocalIndex = globalLocalIndexGridDestSendToServer.size();
426  nbLocalIndexOnGridDest_.push_back(nbLocalIndex);
427  localMaskOnGridDest_.push_back(std::vector<bool>());
428  std::vector<bool>& tmpMask = localMaskOnGridDest_.back();
429  tmpMask.resize(nbLocalIndex,false);
430
431  // Find out number of index sent from grid source and number of index received on grid destination
432  SourceDestinationIndexMap::const_iterator itbIndex = globaIndexWeightFromSrcToDst.begin(),
433                                            iteIndex = globaIndexWeightFromSrcToDst.end(), itIndex;
434  typedef boost::unordered_map<size_t, std::vector<std::pair<size_t,double> > > SendIndexMap;
435  std::map<int,int> sendRankSizeMap,recvRankSizeMap;
436  int connectedClient = globaIndexWeightFromSrcToDst.size();
437  int* recvCount=new int[nbClient];
438  int* displ=new int[nbClient];
439  int* sendRankBuff=new int[connectedClient];
440  int* sendSizeBuff=new int[connectedClient];
441  int n = 0;
442  for (itIndex = itbIndex; itIndex != iteIndex; ++itIndex, ++n)
443  {
444    sendRankBuff[n] = itIndex->first;
445    const SendIndexMap& sendIndexMap = itIndex->second;
446    SendIndexMap::const_iterator itbSend = sendIndexMap.begin(), iteSend = sendIndexMap.end(), itSend;
447    int sendSize = 0;
448    for (itSend = itbSend; itSend != iteSend; ++itSend)
449    {
450      sendSize += itSend->second.size();
451    }
452    sendSizeBuff[n] = sendSize;
453    sendRankSizeMap[itIndex->first] = sendSize;
454  }
455  MPI_Allgather(&connectedClient,1,MPI_INT,recvCount,1,MPI_INT,client->intraComm);
456
457  displ[0]=0 ;
458  for(int n=1;n<nbClient;n++) displ[n]=displ[n-1]+recvCount[n-1];
459  int recvSize=displ[nbClient-1]+recvCount[nbClient-1];
460  int* recvRankBuff=new int[recvSize];
461  int* recvSizeBuff=new int[recvSize];
462  MPI_Allgatherv(sendRankBuff,connectedClient,MPI_INT,recvRankBuff,recvCount,displ,MPI_INT,client->intraComm);
463  MPI_Allgatherv(sendSizeBuff,connectedClient,MPI_INT,recvSizeBuff,recvCount,displ,MPI_INT,client->intraComm);
464  for (int i = 0; i < nbClient; ++i)
465  {
466    int currentPos = displ[i];
467    for (int j = 0; j < recvCount[i]; ++j)
468      if (recvRankBuff[currentPos+j] == clientRank)
469      {
470        recvRankSizeMap[i] = recvSizeBuff[currentPos+j];
471      }
472  }
473
474  // Sending global index of grid source to corresponding process as well as the corresponding mask
475  std::vector<MPI_Request> requests;
476  std::vector<MPI_Status> status;
477  boost::unordered_map<int, unsigned char* > recvMaskDst;
478  boost::unordered_map<int, unsigned long* > recvGlobalIndexSrc;
479  for (std::map<int,int>::const_iterator itRecv = recvRankSizeMap.begin(); itRecv != recvRankSizeMap.end(); ++itRecv)
480  {
481    int recvRank = itRecv->first;
482    int recvSize = itRecv->second;
483    recvMaskDst[recvRank] = new unsigned char [recvSize];
484    recvGlobalIndexSrc[recvRank] = new unsigned long [recvSize];
485
486    requests.push_back(MPI_Request());
487    MPI_Irecv(recvGlobalIndexSrc[recvRank], recvSize, MPI_UNSIGNED_LONG, recvRank, 46, client->intraComm, &requests.back());
488    requests.push_back(MPI_Request());
489    MPI_Irecv(recvMaskDst[recvRank], recvSize, MPI_UNSIGNED_CHAR, recvRank, 47, client->intraComm, &requests.back());
490  }
491
492  boost::unordered_map<int, CArray<size_t,1> > globalIndexDst;
493  boost::unordered_map<int, CArray<double,1> > weightDst;
494  boost::unordered_map<int, unsigned char* > sendMaskDst;
495  boost::unordered_map<int, unsigned long* > sendGlobalIndexSrc;
496  for (itIndex = itbIndex; itIndex != iteIndex; ++itIndex)
497  {
498    int sendRank = itIndex->first;
499    int sendSize = sendRankSizeMap[sendRank];
500    const SendIndexMap& sendIndexMap = itIndex->second;
501    SendIndexMap::const_iterator itbSend = sendIndexMap.begin(), iteSend = sendIndexMap.end(), itSend;
502    globalIndexDst[sendRank].resize(sendSize);
503    weightDst[sendRank].resize(sendSize);
504    sendMaskDst[sendRank] = new unsigned char [sendSize];
505    sendGlobalIndexSrc[sendRank] = new unsigned long [sendSize];
506    int countIndex = 0;
507    for (itSend = itbSend; itSend != iteSend; ++itSend)
508    {
509      const std::vector<std::pair<size_t,double> >& dstWeight = itSend->second;
510      for (int idx = 0; idx < dstWeight.size(); ++idx)
511      {
512        globalIndexDst[sendRank](countIndex) = dstWeight[idx].first;
513        weightDst[sendRank](countIndex) = dstWeight[idx].second;
514        if (0 < globalLocalIndexGridDestSendToServer.count(dstWeight[idx].first))
515          sendMaskDst[sendRank][countIndex] = 1;
516        else
517          sendMaskDst[sendRank][countIndex] = 0;
518        sendGlobalIndexSrc[sendRank][countIndex] = itSend->first;
519        ++countIndex;
520      }
521    }
522
523    // Send global index source and mask
524    requests.push_back(MPI_Request());
525    MPI_Isend(sendGlobalIndexSrc[sendRank], sendSize, MPI_UNSIGNED_LONG, sendRank, 46, client->intraComm, &requests.back());
526    requests.push_back(MPI_Request());
527    MPI_Isend(sendMaskDst[sendRank], sendSize, MPI_UNSIGNED_CHAR, sendRank, 47, client->intraComm, &requests.back());
528  }
529
530  status.resize(requests.size());
531  MPI_Waitall(requests.size(), &requests[0], &status[0]);
532
533  // 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
534  std::vector<MPI_Request>().swap(requests);
535  std::vector<MPI_Status>().swap(status);
536  // Okie, on destination side, we will wait for information of masked index of source
537  for (std::map<int,int>::const_iterator itSend = sendRankSizeMap.begin(); itSend != sendRankSizeMap.end(); ++itSend)
538  {
539    int recvRank = itSend->first;
540    int recvSize = itSend->second;
541
542    requests.push_back(MPI_Request());
543    MPI_Irecv(sendMaskDst[recvRank], recvSize, MPI_UNSIGNED_CHAR, recvRank, 48, client->intraComm, &requests.back());
544  }
545
546  // Ok, now we fill in local index of grid source (we even count for masked index)
547  CDistributionClient distributionClientSrc(client->clientRank, gridSource_);
548  CDistributionClient::GlobalLocalDataMap& globalLocalIndexGridSrcSendToServer = distributionClientSrc.getGlobalLocalDataSendToServer();
549  localIndexToSendFromGridSource_.push_back(SendingIndexGridSourceMap());
550  SendingIndexGridSourceMap& tmpSend = localIndexToSendFromGridSource_.back();
551  for (std::map<int,int>::const_iterator itRecv = recvRankSizeMap.begin(); itRecv != recvRankSizeMap.end(); ++itRecv)
552  {
553    int recvRank = itRecv->first;
554    int recvSize = itRecv->second;
555    unsigned char* recvMask = recvMaskDst[recvRank];
556    unsigned long* recvIndexSrc = recvGlobalIndexSrc[recvRank];
557    int realSendSize = 0;
558    for (int idx = 0; idx < recvSize; ++idx)
559    {
560      if (0 != (*(recvMask+idx))) // OKie, now we have a demand from non-masked index destination
561        if (0 < globalLocalIndexGridSrcSendToServer.count(*(recvIndexSrc+idx))) // check whether index source is masked
562         ++realSendSize;
563        else // inform the destination that this index is masked
564         *(recvMask+idx) = 0;
565    }
566
567    tmpSend[recvRank].resize(realSendSize);
568    realSendSize = 0;
569    for (int idx = 0; idx < recvSize; ++idx)
570    {
571      if (0 != (*(recvMask+idx))) // OKie, now we have a demand from non-masked index destination
572      {
573        tmpSend[recvRank](realSendSize) = globalLocalIndexGridSrcSendToServer[*(recvIndexSrc+idx)];
574         ++realSendSize;
575      }
576    }
577
578    // Okie, now inform the destination which source index are masked
579    requests.push_back(MPI_Request());
580    MPI_Isend(recvMaskDst[recvRank], recvSize, MPI_UNSIGNED_CHAR, recvRank, 48, client->intraComm, &requests.back());
581  }
582  status.resize(requests.size());
583  MPI_Waitall(requests.size(), &requests[0], &status[0]);
584
585  // Cool, now we can fill in local index of grid destination (counted for masked index)
586  localIndexToReceiveOnGridDest_.push_back(RecvIndexGridDestinationMap());
587  RecvIndexGridDestinationMap& recvTmp = localIndexToReceiveOnGridDest_.back();
588  for (std::map<int,int>::const_iterator itSend = sendRankSizeMap.begin(); itSend != sendRankSizeMap.end(); ++itSend)
589  {
590    int recvRank = itSend->first;
591    int recvSize = itSend->second;
592    unsigned char* recvMask = sendMaskDst[recvRank];
593
594    CArray<size_t,1>& recvIndexDst = globalIndexDst[recvRank];
595    CArray<double,1>& recvWeightDst = weightDst[recvRank];
596    int realRecvSize = 0;
597    for (int idx = 0; idx < recvSize; ++idx)
598    {
599      if (0 != *(recvMask+idx)) // OKie, now we have a non-masked index destination
600         ++realRecvSize;
601    }
602
603    int localIndexDst;
604    recvTmp[recvRank].resize(realRecvSize);
605    realRecvSize = 0;
606    for (int idx = 0; idx < recvSize; ++idx)
607    {
608      if (0 != *(recvMask+idx)) // OKie, now we have a demand from non-masked index destination
609      {
610        recvTmp[recvRank][realRecvSize].first = globalLocalIndexGridDestSendToServer[recvIndexDst(idx)];
611        recvTmp[recvRank][realRecvSize].second = recvWeightDst(idx);
612        tmpMask[globalLocalIndexGridDestSendToServer[recvIndexDst(idx)]] = true;
613         ++realRecvSize;
614      }
615    }
616  }
617
618  delete [] recvCount;
619  delete [] displ;
620  delete [] sendRankBuff;
621  delete [] recvRankBuff;
622  delete [] sendSizeBuff;
623  delete [] recvSizeBuff;
624
625  boost::unordered_map<int, unsigned char* >::const_iterator itChar;
626  for (itChar = sendMaskDst.begin(); itChar != sendMaskDst.end(); ++itChar)
627    delete [] itChar->second;
628  for (itChar = recvMaskDst.begin(); itChar != recvMaskDst.end(); ++itChar)
629    delete [] itChar->second;
630  boost::unordered_map<int, unsigned long* >::const_iterator itLong;
631  for (itLong = sendGlobalIndexSrc.begin(); itLong != sendGlobalIndexSrc.end(); ++itLong)
632    delete [] itLong->second;
633  for (itLong = recvGlobalIndexSrc.begin(); itLong != recvGlobalIndexSrc.end(); ++itLong)
634    delete [] itLong->second;
635
636}
637
638/*!
639  Local index of data which need sending from the grid source
640  \return local index of data
641*/
642const std::list<CGridTransformation::SendingIndexGridSourceMap>& CGridTransformation::getLocalIndexToSendFromGridSource() const
643{
644  return localIndexToSendFromGridSource_;
645}
646
647/*!
648  Local index of data which will be received on the grid destination
649  \return local index of data
650*/
651const std::list<CGridTransformation::RecvIndexGridDestinationMap>& CGridTransformation::getLocalIndexToReceiveOnGridDest() const
652{
653  return localIndexToReceiveOnGridDest_;
654}
655
656/*!
657 Number of index will be received on the grid destination
658  \return number of index of data
659*/
660const std::list<size_t>& CGridTransformation::getNbLocalIndexToReceiveOnGridDest() const
661{
662  return nbLocalIndexOnGridDest_;
663}
664
665/*!
666  Local mask of data which will be received on the grid destination
667  \return local mask of data
668*/
669const std::list<std::vector<bool> >& CGridTransformation::getLocalMaskIndexOnGridDest() const
670{
671  return localMaskOnGridDest_;
672}
673
674}
Note: See TracBrowser for help on using the repository browser.