source: XIOS/trunk/src/filter/grid_transformation.cpp @ 622

Last change on this file since 622 was 622, checked in by mhnguyen, 9 years ago

Final testing transfomation algorithm: inverse axis (local commit)

+) Make some minor change to make sure one element (axis or domain) be able to have several similar transformation

Test
+) On Curie
+) test_new_feature: test passed with correct data written

File size: 19.9 KB
Line 
1#include "grid_transformation.hpp"
2#include "axis_inverse.hpp"
3#include "axis_zoom.hpp"
4#include "context.hpp"
5#include "context_client.hpp"
6#include "transformation_mapping.hpp"
7
8#include "axis_algorithm_transformation.hpp"
9
10namespace xios {
11CGridTransformation::CGridTransformation(CGrid* destination, CGrid* source)
12: gridSource_(source), gridDestination_(destination), originalGridSource_(source),
13  globalIndexOfCurrentGridSource_(0), globalIndexOfOriginalGridSource_(0)
14{
15  //Verify the compatibity between two grids
16  int numElement = gridDestination_->axis_domain_order.numElements();
17  if (numElement != gridSource_->axis_domain_order.numElements())
18    ERROR("CGridTransformation::CGridTransformation(CGrid* destination, CGrid* source)",
19       << "Two grids have different number of elements"
20       << "Number of elements of grid source " <<gridSource_->getId() << " is " << gridSource_->axis_domain_order.numElements()  << std::endl
21       << "Number of elements of grid destination " <<gridDestination_->getId() << " is " << numElement);
22
23  for (int i = 0; i < numElement; ++i)
24  {
25    if (gridDestination_->axis_domain_order(i) != gridSource_->axis_domain_order(i))
26      ERROR("CGridTransformation::CGridTransformation(CGrid* destination, CGrid* source)",
27         << "Transformed grid and its grid source have incompatible elements"
28         << "Grid source " <<gridSource_->getId() << std::endl
29         << "Grid destination " <<gridDestination_->getId());
30  }
31
32  std::vector<CAxis*> axisSrcTmp = gridSource_->getAxis(), axisSrc;
33  std::vector<CDomain*> domainSrcTmp = gridSource_->getDomains(), domainSrc;
34  for (int idx = 0; idx < axisSrcTmp.size(); ++idx)
35  {
36    CAxis* axis = CAxis::createAxis();
37    axisSrcTmp[idx]->duplicateAttributes(axis);
38    axisSrc.push_back(axis);
39  }
40
41  for (int idx = 0; idx < domainSrcTmp.size(); ++idx)
42  {
43    CDomain* domain = CDomain::createDomain();
44    domainSrcTmp[idx]->duplicateAttributes(domain);
45    domainSrc.push_back(domain);
46  }
47
48  gridSource_ = CGrid::createGrid(domainSrc, axisSrc, gridDestination_->axis_domain_order);
49  gridSourceDimensionSize_ = gridSource_->getGlobalDimension();
50  gridDestinationDimensionSize_ = gridDestination_->getGlobalDimension();
51
52  initializeMappingOfOriginalGridSource();
53  initializeAlgorithms();
54}
55
56void CGridTransformation::initializeMappingOfOriginalGridSource()
57{
58  CContext* context = CContext::getCurrent();
59  CContextClient* client=context->client;
60
61  CDistributionClient distribution(client->clientRank, originalGridSource_);
62  const CArray<size_t,1>& globalIndexGridDestSendToServer = distribution.getGlobalDataIndexSendToServer();
63
64  globalIndexOfCurrentGridSource_   = new CArray<size_t,1>(globalIndexGridDestSendToServer.numElements());
65  globalIndexOfOriginalGridSource_  = new CArray<size_t,1>(globalIndexGridDestSendToServer.numElements());
66  *globalIndexOfCurrentGridSource_  = globalIndexGridDestSendToServer;
67  *globalIndexOfOriginalGridSource_ = globalIndexGridDestSendToServer;
68}
69
70CGridTransformation::~CGridTransformation()
71{
72  std::list<CGenericAlgorithmTransformation*>::const_iterator itb = algoTransformation_.begin(), it,
73                                                              ite = algoTransformation_.end();
74  for (it = itb; it != ite; ++it) delete (*it);
75
76  std::map<int, std::vector<CArray<int,1>* > >::const_iterator itMapRecv, iteMapRecv;
77  itMapRecv = localIndexToReceiveOnGridDest_.begin();
78  iteMapRecv = localIndexToReceiveOnGridDest_.end();
79  for (; itMapRecv != iteMapRecv; ++itMapRecv)
80  {
81    int numVec = (itMapRecv->second).size();
82    for (int idx = 0; idx < numVec; ++idx) delete (itMapRecv->second)[idx];
83  }
84
85  std::map<int, CArray<int,1>* >::const_iterator itMap, iteMap;
86  itMap = localIndexToSendFromGridSource_.begin();
87  iteMap = localIndexToSendFromGridSource_.end();
88  for (; itMap != iteMap; ++itMap) delete (itMap->second);
89
90  if (0 != globalIndexOfCurrentGridSource_) delete globalIndexOfCurrentGridSource_;
91  if (0 != globalIndexOfOriginalGridSource_) delete globalIndexOfOriginalGridSource_;
92}
93
94void CGridTransformation::initializeAlgorithms()
95{
96  initializeAxisAlgorithms();
97  initializeDomainAlgorithms();
98}
99
100/*!
101  Initialize the algorithms corresponding to transformation info contained in each axis.
102If an axis has transformations, these transformations will be represented in form of vector of CTransformation pointers
103In general, each axis can have several transformations performed on itself. However, should they be done seperately or combinely (of course in order)?
104For now, one approach is to do these combinely but maybe this needs changing.
105*/
106void CGridTransformation::initializeAxisAlgorithms()
107{
108  std::vector<int> axisPositionInGrid;
109  std::vector<CAxis*> axisListDestP = gridDestination_->getAxis();
110  std::vector<CAxis*> axisListSrcP = gridSource_->getAxis();
111  if (!axisListDestP.empty())
112  {
113    int idx = 0;
114    for (int i = 0; i < gridDestination_->axis_domain_order.numElements(); ++i)
115    {
116      if (false == (gridDestination_->axis_domain_order)(i))
117      {
118        axisPositionInGrid.push_back(idx);
119        ++idx;
120      }
121      else idx += 2;
122    }
123
124    for (int i = 0; i < axisListDestP.size(); ++i)
125    {
126      elementPosition2AxisPositionInGrid_[axisPositionInGrid[i]] = i;
127      if (axisListDestP[i]->hasTransformation())
128      {
129        CAxis::TransMapTypes trans = axisListDestP[i]->getAllTransformations();
130        CAxis::TransMapTypes::const_iterator itb = trans.begin(), it,
131                                             ite = trans.end();
132        int transformationOrder = 0;
133        for (it = itb; it != ite; ++it)
134        {
135          listAlgos_.push_back(std::make_pair(axisPositionInGrid[i], std::make_pair(it->first, transformationOrder)));
136          ++transformationOrder;
137        }
138      }
139    }
140  }
141}
142
143void CGridTransformation::initializeDomainAlgorithms()
144{
145
146}
147
148void CGridTransformation::selectAlgo(int elementPositionInGrid, ETranformationType transType, int transformationOrder)
149{
150   selectAxisAlgo(elementPositionInGrid, transType, transformationOrder);
151}
152
153void CGridTransformation::selectAxisAlgo(int elementPositionInGrid, ETranformationType transType, int transformationOrder)
154{
155  std::vector<CAxis*> axisListDestP = gridDestination_->getAxis();
156  std::vector<CAxis*> axisListSrcP = gridSource_->getAxis();
157
158  int axisIndex =  elementPosition2AxisPositionInGrid_[elementPositionInGrid];
159  CAxis::TransMapTypes trans = axisListDestP[axisIndex]->getAllTransformations();
160  CAxis::TransMapTypes::const_iterator it = trans.begin();
161
162  for (int i = 0; i < transformationOrder; ++i, ++it) {}  // Find the correct transformation
163
164  CZoomAxis* zoomAxis = 0;
165  CGenericAlgorithmTransformation* algo = 0;
166  switch (transType)
167  {
168    case TRANS_ZOOM_AXIS:
169      zoomAxis = dynamic_cast<CZoomAxis*> (it->second);
170      algo = new CAxisZoom(axisListDestP[axisIndex], axisListSrcP[axisIndex], zoomAxis);
171      break;
172    case TRANS_INVERSE_AXIS:
173      algo = new CAxisInverse(axisListDestP[axisIndex], axisListSrcP[axisIndex]);
174      break;
175    default:
176      break;
177  }
178  algoTransformation_.push_back(algo);
179
180}
181
182void CGridTransformation::selectDomainAlgo(int elementPositionInGrid, ETranformationType transType, int transformationOrder)
183{
184
185}
186
187void CGridTransformation::setUpGrid(int elementPositionInGrid, ETranformationType transType)
188{
189  std::vector<CAxis*> axisListDestP = gridDestination_->getAxis();
190  std::vector<CAxis*> axisListSrcP = gridSource_->getAxis();
191
192  int axisIndex;
193  switch (transType)
194  {
195    case TRANS_ZOOM_AXIS:
196    case TRANS_INVERSE_AXIS:
197      axisIndex =  elementPosition2AxisPositionInGrid_[elementPositionInGrid];
198      axisListDestP[axisIndex]->duplicateAttributes(axisListSrcP[axisIndex]);
199      break;
200    default:
201      break;
202  }
203}
204
205void CGridTransformation::computeAll()
206{
207  CContext* context = CContext::getCurrent();
208  CContextClient* client=context->client;
209
210  ListAlgoType::const_iterator itb = listAlgos_.begin(),
211                               ite = listAlgos_.end(), it;
212  CGenericAlgorithmTransformation* algo = 0;
213  for (it = itb; it != ite; ++it)
214  {
215    std::map<size_t, std::set<size_t> > globaIndexMapFromDestToSource;
216    int elementPositionInGrid = it->first;
217    ETranformationType transType = (it->second).first;
218    int transformationOrder = (it->second).second;
219
220    // First of all, select an algorithm
221    selectAlgo(elementPositionInGrid, transType, transformationOrder);
222    algo = algoTransformation_.back();
223
224    // Recalculate the distribution of grid destination
225    CDistributionClient distributionClientDest(client->clientRank, gridDestination_);
226    const CArray<size_t,1>& globalIndexGridDestSendToServer = distributionClientDest.getGlobalDataIndexSendToServer();
227
228    // ComputeTransformation of global index of each element
229    std::vector<int> gridDestinationDimensionSize = gridDestination_->getGlobalDimension();
230    int elementPosition = it->first;
231    algo->computeGlobalSourceIndex(elementPosition,
232                                   gridDestinationDimensionSize,
233                                   globalIndexGridDestSendToServer,
234                                   globaIndexMapFromDestToSource);
235
236    // Compute transformation of global indexes among grids
237    computeTransformationFromOriginalGridSource(globaIndexMapFromDestToSource);
238
239    // Now grid destination becomes grid source in a new transformation
240    setUpGrid(elementPositionInGrid, transType);
241  }
242
243 std::cout << "global index destination 0 " << *globalIndexOfCurrentGridSource_ << std::endl;
244 std::cout << "global index destination 1 " << *globalIndexOfOriginalGridSource_ << std::endl;
245  computeFinalTransformationMapping();
246}
247
248
249/*!
250  Compute index mapping representing transformation between two grids
251  Each domain and each axis can contain some information of transformation, these information are then converted into
252form of global index mapping reprensenting transformation between two grids.
253*/
254void CGridTransformation::computeTransformation()
255{
256//  const CArray<size_t,1>& globalIndexGridDestSendToServer = gridDestination_->getDistributionClient()->getGlobalDataIndexSendToServer();
257//  std::list<std::pair<int,CGenericAlgorithmTransformation*> >::const_iterator itbMap, itMap, iteMap;
258//  itbMap = algoTransformation_.begin();
259//  iteMap = algoTransformation_.end();
260//
261//  std::vector<CGenericAlgorithmTransformation*>::const_iterator itbVec, itVec, iteVec;
262//
263//  for (itMap = itbMap; itMap != iteMap; ++itMap)
264//  {
265//    int elementPosition = itMap->first;
266//    itbVec = (itMap->second).begin();
267//    iteVec = (itMap->second).end();
268//    for (itVec = itbVec; itVec != iteVec; ++itVec)
269//    {
270//      (*itVec)->computeGlobalSourceIndex(elementPosition,
271//                                         gridDestinationDimensionSize_,
272//                                         globalIndexGridDestSendToServer,
273//                                         globaIndexMapFromDestToSource_);
274//    }
275//  }
276}
277
278/*!
279  A transformation from a grid source to grid destination often passes several intermediate grids, which play a role of
280temporary grid source and/or grid destination. This function makes sure that global index of original grid source are mapped correctly to
281the final grid destination
282*/
283void CGridTransformation::computeTransformationFromOriginalGridSource(const std::map<size_t, std::set<size_t> >& globaIndexMapFromDestToSource)
284{
285  CContext* context = CContext::getCurrent();
286  CContextClient* client=context->client;
287
288  CTransformationMapping transformationMap(gridDestination_, gridSource_);
289
290    // Then compute transformation mapping among clients
291  transformationMap.computeTransformationMapping(globaIndexMapFromDestToSource);
292
293  const std::map<int,std::vector<std::vector<size_t> > >& globalIndexToReceive = transformationMap.getGlobalIndexReceivedOnGridDestMapping();
294  const std::map<int,std::vector<size_t> >& globalIndexToSend = transformationMap.getGlobalIndexSendToGridDestMapping();
295
296 // Sending global index of original grid source
297  std::map<int,std::vector<size_t> >::const_iterator itbSend = globalIndexToSend.begin(), itSend,
298                                                     iteSend = globalIndexToSend.end();
299  CArray<size_t,1>::const_iterator itbArr = globalIndexOfCurrentGridSource_->begin(), itArr,
300                                   iteArr = globalIndexOfCurrentGridSource_->end();
301 int sendBuffSize = 0;
302 for (itSend = itbSend; itSend != iteSend; ++itSend) sendBuffSize += (itSend->second).size();
303
304 unsigned long* sendBuff, *currentSendBuff;
305 if (0 != sendBuffSize) sendBuff = new unsigned long [sendBuffSize];
306 int currentBuffPosition = 0;
307 for (itSend = itbSend; itSend != iteSend; ++itSend)
308 {
309   int destRank = itSend->first;
310   const std::vector<size_t>& globalIndexOfCurrentGridSourceToSend = itSend->second;
311   int countSize = globalIndexOfCurrentGridSourceToSend.size();
312   for (int idx = 0; idx < (countSize); ++idx)
313   {
314     itArr = std::find(itbArr, iteArr, globalIndexOfCurrentGridSourceToSend[idx]);
315     if (iteArr != itArr)
316     {
317       int index = std::distance(itbArr, itArr);
318       sendBuff[idx+currentBuffPosition] = (*globalIndexOfOriginalGridSource_)(index);
319     }
320   }
321   currentSendBuff = sendBuff + currentBuffPosition;
322   MPI_Send(currentSendBuff, countSize, MPI_UNSIGNED_LONG, destRank, 14, client->intraComm);
323   currentBuffPosition += countSize;
324 }
325
326 // Receiving global index of grid source sending from current grid source
327 std::map<int,std::vector<std::vector<size_t> > >::const_iterator itbRecv = globalIndexToReceive.begin(), itRecv,
328                                                                  iteRecv = globalIndexToReceive.end();
329 int recvBuffSize = 0;
330 for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv) recvBuffSize += (itRecv->second).size();
331
332 unsigned long* recvBuff, *currentRecvBuff;
333 if (0 != recvBuffSize) recvBuff = new unsigned long [recvBuffSize];
334 int currentRecvBuffPosition = 0;
335 for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv)
336 {
337   MPI_Status status;
338   int srcRank = itRecv->first;
339   int countSize = (itRecv->second).size();
340   currentRecvBuff = recvBuff + currentRecvBuffPosition;
341   MPI_Recv(currentRecvBuff, countSize, MPI_UNSIGNED_LONG, srcRank, 14, client->intraComm, &status);
342   currentRecvBuffPosition += countSize;
343 }
344
345 int nbCurrentGridSource = 0;
346 for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv)
347 {
348   int ssize = (itRecv->second).size();
349   for (int idx = 0; idx < ssize; ++idx)
350   {
351     nbCurrentGridSource += (itRecv->second)[idx].size();
352   }
353 }
354
355 globalIndexOfCurrentGridSource_->resize(nbCurrentGridSource);
356 globalIndexOfOriginalGridSource_->resize(nbCurrentGridSource);
357 int k = 0;
358 currentRecvBuff = recvBuff;
359 for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv)
360 {
361   int countSize = (itRecv->second).size();
362   for (int idx = 0; idx < countSize; ++idx, ++currentRecvBuff)
363   {
364     int ssize = (itRecv->second)[idx].size();
365     for (int i = 0; i < ssize; ++i)
366     {
367       (*globalIndexOfCurrentGridSource_)(k) = (itRecv->second)[idx][i];
368       (*globalIndexOfOriginalGridSource_)(k) = *currentRecvBuff;
369       ++k;
370     }
371   }
372 }
373
374 if (0 != sendBuffSize) delete [] sendBuff;
375 if (0 != recvBuffSize) delete [] recvBuff;
376}
377
378/*!
379  Compute transformation mapping between grid source and grid destination
380  The transformation between grid source and grid destination is represented in form of mapping between global index
381of two grids. Then local index mapping between data on each grid will be found out thanks to these global indexes
382*/
383void CGridTransformation::computeFinalTransformationMapping()
384{
385  CContext* context = CContext::getCurrent();
386  CContextClient* client=context->client;
387
388  CTransformationMapping transformationMap(gridDestination_, originalGridSource_);
389
390  std::map<size_t, std::set<size_t> > globaIndexMapFromDestToSource;
391
392  int nb = globalIndexOfCurrentGridSource_->numElements();
393  for (int idx = 0; idx < nb; ++idx)
394  {
395    globaIndexMapFromDestToSource[(*globalIndexOfCurrentGridSource_)(idx)].insert((*globalIndexOfOriginalGridSource_)(idx));
396  }
397
398  // Then compute transformation mapping among clients
399  transformationMap.computeTransformationMapping(globaIndexMapFromDestToSource);
400
401  const std::map<int,std::vector<std::vector<size_t> > >& globalIndexToReceive = transformationMap.getGlobalIndexReceivedOnGridDestMapping();
402  const std::map<int,std::vector<size_t> >& globalIndexToSend = transformationMap.getGlobalIndexSendToGridDestMapping();
403
404  CDistributionClient distributionClientDest(client->clientRank, gridDestination_);
405  CDistributionClient distributionClientSrc(client->clientRank, gridSource_);
406
407  CArray<int, 1> localIndexOnClientDest = distributionClientDest.getLocalDataIndexOnClient(); //gridDestination_->getDistributionClient()->getLocalDataIndexOnClient();
408  CArray<size_t,1> globalIndexOnClientDest = distributionClientDest.getGlobalDataIndexSendToServer(); //gridDestination_->getDistributionClient()->getGlobalDataIndexSendToServer();
409
410 std::cout << "dest: local index " << localIndexOnClientDest << std::endl;
411 std::cout << "dest: global index " << globalIndexOnClientDest << std::endl;
412  CArray<int, 1> localIndexOnClientSrc = distributionClientSrc.getLocalDataIndexOnClient(); //gridSource_->getDistributionClient()->getLocalDataIndexOnClient();
413  CArray<size_t,1> globalIndexOnClientSrc = distributionClientSrc.getGlobalDataIndexSendToServer(); //gridSource_->getDistributionClient()->getGlobalDataIndexSendToServer();
414 std::cout << "src: local index " << localIndexOnClientSrc << std::endl;
415 std::cout << "src: global index " << globalIndexOnClientSrc << std::endl;
416  std::vector<size_t>::const_iterator itbVec, itVec, iteVec;
417  CArray<size_t, 1>::iterator itbArr, itArr, iteArr;
418
419  std::map<int,std::vector<std::vector<size_t> > >::const_iterator itbMapRecv, itMapRecv, iteMapRecv;
420
421  // Find out local index on grid destination (received)
422  itbMapRecv = globalIndexToReceive.begin();
423  iteMapRecv = globalIndexToReceive.end();
424  itbArr = globalIndexOnClientDest.begin();
425  iteArr = globalIndexOnClientDest.end();
426  for (itMapRecv = itbMapRecv; itMapRecv != iteMapRecv; ++itMapRecv)
427  {
428    int sourceRank = itMapRecv->first;
429    int numGlobalIndex = (itMapRecv->second).size();
430    for (int i = 0; i < numGlobalIndex; ++i)
431    {
432      int vecSize = ((itMapRecv->second)[i]).size();
433      CArray<int,1>* ptr = new CArray<int,1>(vecSize);
434      localIndexToReceiveOnGridDest_[sourceRank].push_back(ptr);
435      for (int idx = 0; idx < vecSize; ++idx)
436      {
437        itArr = std::find(itbArr, iteArr, (itMapRecv->second)[i][idx]);
438        if (iteArr != itArr)
439        {
440          int localIdx = std::distance(itbArr, itArr);
441//          (*localIndexToReceiveOnGridDest_[sourceRank][i])(idx) = localIndexOnClientDest(localIdx); // Local index of un-extracted data (only domain)
442          (*localIndexToReceiveOnGridDest_[sourceRank][i])(idx) = (localIdx); // Local index of extracted data
443        }
444      }
445    }
446  }
447
448  std::map<int,std::vector<size_t> >::const_iterator itbMap, itMap, iteMap;
449  // Find out local index on grid source (to send)
450  itbMap = globalIndexToSend.begin();
451  iteMap = globalIndexToSend.end();
452  itbArr = globalIndexOnClientSrc.begin();
453  iteArr = globalIndexOnClientSrc.end();
454  for (itMap = itbMap; itMap != iteMap; ++itMap)
455  {
456    CArray<int,1>* ptr = new CArray<int,1>((itMap->second).size());
457    localIndexToSendFromGridSource_[itMap->first] = ptr;
458    int destRank = itMap->first;
459    int vecSize = (itMap->second).size();
460    for (int idx = 0; idx < vecSize; ++idx)
461    {
462      itArr = std::find(itbArr, iteArr, (itMap->second)[idx]);
463      if (iteArr != itArr)
464      {
465        int localIdx = std::distance(itbArr, itArr);
466//        (*localIndexToSendFromGridSource_[destRank])(idx) = localIndexOnClientSrc(localIdx);
467        (*localIndexToSendFromGridSource_[destRank])(idx) = (localIdx);
468      }
469    }
470  }
471}
472
473/*!
474  Local index of data which need sending from the grid source
475  \return local index of data
476*/
477std::map<int, CArray<int,1>* > CGridTransformation::getLocalIndexToSendFromGridSource()
478{
479  return localIndexToSendFromGridSource_;
480}
481
482/*!
483  Local index of data which will be received on the grid destination
484  \return local index of data
485*/
486std::map<int, std::vector<CArray<int,1>* > > CGridTransformation::getLocalIndexToReceiveOnGridDest()
487{
488  return localIndexToReceiveOnGridDest_;
489}
490
491}
Note: See TracBrowser for help on using the repository browser.