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

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

Chaning the way to process transformation to improve the performance.
Instead of exchanging global index and weights on full GRID, each process only
sends and receives the global index and weights on each ELEMENT, which can reduce
the message size of DHT.

+) Domain and axis now have their own exchange function to transfer global index and weight
+) Generic transformation now plays the role of "synthesizer" for all elements
+) Grid transformation now plays the role of transformation mapping, e.x: exchange final global index and weight
among processes.

Test
+) On Curie
+) Pass on all basic tests
+) Dynamic interpolation on axis hasn't been tested (and it seems to need more change to make it rework)

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