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

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

Fixing a bug in transformation in case of masking grid

+) Number of local index in grid destination does not depend on number of masked index in grid source

Test
+) On Curie
+) Simple test with NEMO
+) Up to 20 cores

File size: 20.5 KB
Line 
1/*!
2   \file grid_transformation.cpp
3   \author Ha NGUYEN
4   \since 14 May 2015
5   \date 02 Jul 2015
6
7   \brief Interface for all transformations.
8 */
9#include "grid_transformation.hpp"
10#include "axis_algorithm_inverse.hpp"
11#include "axis_algorithm_zoom.hpp"
12#include "axis_algorithm_interpolate.hpp"
13#include "domain_algorithm_zoom.hpp"
14#include "domain_algorithm_interpolate.hpp"
15#include "context.hpp"
16#include "context_client.hpp"
17#include "transformation_mapping.hpp"
18#include "axis_algorithm_transformation.hpp"
19#include "distribution_client.hpp"
20#include "mpi_tag.hpp"
21#include <boost/unordered_map.hpp>
22
23namespace xios {
24CGridTransformation::CGridTransformation(CGrid* destination, CGrid* source)
25: gridSource_(source), gridDestination_(destination), originalGridSource_(source),
26  algoTypes_(), nbAlgos_(0), currentGridIndexToOriginalGridIndex_(), tempGrids_(),
27  auxInputs_(), dynamicalTransformation_(false), timeStamp_()
28
29{
30  //Verify the compatibity between two grids
31  int numElement = gridDestination_->axis_domain_order.numElements();
32  if (numElement != gridSource_->axis_domain_order.numElements())
33    ERROR("CGridTransformation::CGridTransformation(CGrid* destination, CGrid* source)",
34       << "Two grids have different number of elements"
35       << "Number of elements of grid source " <<gridSource_->getId() << " is " << gridSource_->axis_domain_order.numElements()  << std::endl
36       << "Number of elements of grid destination " <<gridDestination_->getId() << " is " << numElement);
37
38  for (int i = 0; i < numElement; ++i)
39  {
40    if (gridDestination_->axis_domain_order(i) != gridSource_->axis_domain_order(i))
41      ERROR("CGridTransformation::CGridTransformation(CGrid* destination, CGrid* source)",
42         << "Transformed grid and its grid source have incompatible elements"
43         << "Grid source " <<gridSource_->getId() << std::endl
44         << "Grid destination " <<gridDestination_->getId());
45  }
46
47  initializeMappingOfOriginalGridSource();
48}
49
50/*!
51  Initialize the mapping between the first grid source and the original one
52  In a series of transformation, for each step, there is a need to "create" a new grid that plays a role of "temporary" source.
53Because at the end of the series, we need to know about the index mapping between the final grid destination and original grid source,
54for each transformation, we need to make sure that the current "temporary source" maps its global index correctly to the original one.
55*/
56void CGridTransformation::initializeMappingOfOriginalGridSource()
57{
58  CContext* context = CContext::getCurrent();
59  CContextClient* client = context->client;
60
61  // Initialize algorithms
62  initializeAlgorithms();
63
64  ListAlgoType::const_iterator itb = listAlgos_.begin(),
65                               ite = listAlgos_.end(), it;
66
67  for (it = itb; it != ite; ++it)
68  {
69    ETranformationType transType = (it->second).first;
70    if (!isSpecialTransformation(transType)) ++nbAlgos_;
71  }
72}
73
74CGridTransformation::~CGridTransformation()
75{
76  std::vector<CGenericAlgorithmTransformation*>::const_iterator itb = algoTransformation_.begin(), it,
77                                                                ite = algoTransformation_.end();
78  for (it = itb; it != ite; ++it) delete (*it);
79}
80
81/*!
82  Initialize the algorithms (transformations)
83*/
84void CGridTransformation::initializeAlgorithms()
85{
86  std::vector<int> axisPositionInGrid;
87  std::vector<int> domPositionInGrid;
88  std::vector<CAxis*> axisListDestP = gridDestination_->getAxis();
89  std::vector<CDomain*> domListDestP = gridDestination_->getDomains();
90
91  int idx = 0;
92  for (int i = 0; i < gridDestination_->axis_domain_order.numElements(); ++i)
93  {
94    if (false == (gridDestination_->axis_domain_order)(i))
95    {
96      axisPositionInGrid.push_back(idx);
97      ++idx;
98    }
99    else
100    {
101      ++idx;
102      domPositionInGrid.push_back(idx);
103      ++idx;
104    }
105  }
106
107  for (int i = 0; i < axisListDestP.size(); ++i)
108  {
109    elementPosition2AxisPositionInGrid_[axisPositionInGrid[i]] = i;
110  }
111
112  for (int i = 0; i < domListDestP.size(); ++i)
113  {
114    elementPosition2DomainPositionInGrid_[domPositionInGrid[i]] = i;
115  }
116
117  idx = 0;
118  for (int i = 0; i < gridDestination_->axis_domain_order.numElements(); ++i)
119  {
120    if (false == (gridDestination_->axis_domain_order)(i))
121    {
122      initializeAxisAlgorithms(idx);
123      ++idx;
124    }
125    else
126    {
127      ++idx;
128      initializeDomainAlgorithms(idx);
129      ++idx;
130    }
131  }
132}
133
134/*!
135  Initialize the algorithms corresponding to transformation info contained in each axis.
136If an axis has transformations, these transformations will be represented in form of vector of CTransformation pointers
137In general, each axis can have several transformations performed on itself. However, should they be done seperately or combinely (of course in order)?
138For now, one approach is to do these combinely but maybe this needs changing.
139\param [in] axisPositionInGrid position of an axis in grid. (for example: a grid with one domain and one axis, position of domain is 1, position of axis is 2)
140*/
141void CGridTransformation::initializeAxisAlgorithms(int axisPositionInGrid)
142{
143  std::vector<CAxis*> axisListDestP = gridDestination_->getAxis();
144  if (!axisListDestP.empty())
145  {
146    if (axisListDestP[elementPosition2AxisPositionInGrid_[axisPositionInGrid]]->hasTransformation())
147    {
148      CAxis::TransMapTypes trans = axisListDestP[elementPosition2AxisPositionInGrid_[axisPositionInGrid]]->getAllTransformations();
149      CAxis::TransMapTypes::const_iterator itb = trans.begin(), it,
150                                           ite = trans.end();
151      int transformationOrder = 0;
152      for (it = itb; it != ite; ++it)
153      {
154        listAlgos_.push_back(std::make_pair(axisPositionInGrid, std::make_pair(it->first, transformationOrder)));
155        algoTypes_.push_back(false);
156        ++transformationOrder;
157        std::vector<StdString> auxInput = (it->second)->checkAuxInputs();
158        for (int idx = 0; idx < auxInput.size(); ++idx) auxInputs_.push_back(auxInput[idx]);
159      }
160    }
161  }
162}
163
164/*!
165  Initialize the algorithms corresponding to transformation info contained in each domain.
166If a domain has transformations, they will be represented in form of vector of CTransformation pointers
167In general, each domain can have several transformations performed on itself.
168\param [in] domPositionInGrid position of a domain in grid. (for example: a grid with one domain and one axis, position of domain is 1, position of axis is 2)
169*/
170void CGridTransformation::initializeDomainAlgorithms(int domPositionInGrid)
171{
172  std::vector<CDomain*> domListDestP = gridDestination_->getDomains();
173  if (!domListDestP.empty())
174  {
175    if (domListDestP[elementPosition2DomainPositionInGrid_[domPositionInGrid]]->hasTransformation())
176    {
177      CDomain::TransMapTypes trans = domListDestP[elementPosition2DomainPositionInGrid_[domPositionInGrid]]->getAllTransformations();
178      CDomain::TransMapTypes::const_iterator itb = trans.begin(), it,
179                                             ite = trans.end();
180      int transformationOrder = 0;
181      for (it = itb; it != ite; ++it)
182      {
183        listAlgos_.push_back(std::make_pair(domPositionInGrid, std::make_pair(it->first, transformationOrder)));
184        algoTypes_.push_back(true);
185        ++transformationOrder;
186        std::vector<StdString> auxInput = (it->second)->checkAuxInputs();
187        for (int idx = 0; idx < auxInput.size(); ++idx) auxInputs_.push_back(auxInput[idx]);
188      }
189    }
190  }
191
192}
193
194/*!
195  Select algorithm correspoding to its transformation type and its position in each element
196  \param [in] elementPositionInGrid position of element in grid. e.g: a grid has 1 domain and 1 axis, then position of domain is 1 (because it contains 2 basic elements)
197                                             and position of axis is 2
198  \param [in] transType transformation type, for now we have Zoom_axis, inverse_axis
199  \param [in] transformationOrder position of the transformation in an element (an element can have several transformation)
200  \param [in] isDomainAlgo flag to specify type of algorithm (for domain or axis)
201*/
202void CGridTransformation::selectAlgo(int elementPositionInGrid, ETranformationType transType, int transformationOrder, bool isDomainAlgo)
203{
204   if (isDomainAlgo) selectDomainAlgo(elementPositionInGrid, transType, transformationOrder);
205   else selectAxisAlgo(elementPositionInGrid, transType, transformationOrder);
206}
207
208/*!
209  Select algorithm of an axis correspoding to its transformation type and its position in each element
210  \param [in] elementPositionInGrid position of element in grid. e.g: a grid has 1 domain and 1 axis, then position of domain is 1 (because it contains 2 basic elements)
211                                             and position of axis is 2
212  \param [in] transType transformation type, for now we have Zoom_axis, inverse_axis
213  \param [in] transformationOrder position of the transformation in an element (an element can have several transformation)
214*/
215void CGridTransformation::selectAxisAlgo(int elementPositionInGrid, ETranformationType transType, int transformationOrder)
216{
217  std::vector<CAxis*> axisListDestP = gridDestination_->getAxis();
218  std::vector<CAxis*> axisListSrcP = gridSource_->getAxis();
219
220  int axisIndex =  elementPosition2AxisPositionInGrid_[elementPositionInGrid];
221  CAxis::TransMapTypes trans = axisListDestP[axisIndex]->getAllTransformations();
222  CAxis::TransMapTypes::const_iterator it = trans.begin();
223
224  for (int i = 0; i < transformationOrder; ++i, ++it) {}  // Find the correct transformation
225
226  CZoomAxis* zoomAxis = 0;
227  CInterpolateAxis* interpAxis = 0;
228  CGenericAlgorithmTransformation* algo = 0;
229  switch (transType)
230  {
231    case TRANS_INTERPOLATE_AXIS:
232      interpAxis = dynamic_cast<CInterpolateAxis*> (it->second);
233      algo = new CAxisAlgorithmInterpolate(axisListDestP[axisIndex], axisListSrcP[axisIndex], interpAxis);
234      break;
235    case TRANS_ZOOM_AXIS:
236      zoomAxis = dynamic_cast<CZoomAxis*> (it->second);
237      algo = new CAxisAlgorithmZoom(axisListDestP[axisIndex], axisListSrcP[axisIndex], zoomAxis);
238      break;
239    case TRANS_INVERSE_AXIS:
240      algo = new CAxisAlgorithmInverse(axisListDestP[axisIndex], axisListSrcP[axisIndex]);
241      break;
242    default:
243      break;
244  }
245  algoTransformation_.push_back(algo);
246
247}
248
249/*!
250  Select algorithm of a domain correspoding to its transformation type and its position in each element
251  \param [in] elementPositionInGrid position of element in grid. e.g: a grid has 1 domain and 1 axis, then position of domain is 1 (because it contains 2 basic elements)
252                                             and position of axis is 2
253  \param [in] transType transformation type, for now we have Zoom_axis, inverse_axis
254  \param [in] transformationOrder position of the transformation in an element (an element can have several transformation)
255*/
256void CGridTransformation::selectDomainAlgo(int elementPositionInGrid, ETranformationType transType, int transformationOrder)
257{
258  std::vector<CDomain*> domainListDestP = gridDestination_->getDomains();
259  std::vector<CDomain*> domainListSrcP = gridSource_->getDomains();
260
261  int domainIndex =  elementPosition2DomainPositionInGrid_[elementPositionInGrid];
262  CDomain::TransMapTypes trans = domainListDestP[domainIndex]->getAllTransformations();
263  CDomain::TransMapTypes::const_iterator it = trans.begin();
264
265  for (int i = 0; i < transformationOrder; ++i, ++it) {}  // Find the correct transformation
266
267  CZoomDomain* zoomDomain = 0;
268  CInterpolateDomain* interpFileDomain = 0;
269  CGenericAlgorithmTransformation* algo = 0;
270  switch (transType)
271  {
272    case TRANS_INTERPOLATE_DOMAIN:
273      interpFileDomain = dynamic_cast<CInterpolateDomain*> (it->second);
274      algo = new CDomainAlgorithmInterpolate(domainListDestP[domainIndex], domainListSrcP[domainIndex],interpFileDomain);
275      break;
276    case TRANS_ZOOM_DOMAIN:
277      zoomDomain = dynamic_cast<CZoomDomain*> (it->second);
278      algo = new CDomainAlgorithmZoom(domainListDestP[domainIndex], domainListSrcP[domainIndex], zoomDomain);
279      break;
280    default:
281      break;
282  }
283  algoTransformation_.push_back(algo);
284}
285
286/*!
287  Assign the current grid destination to the grid source in the new transformation.
288The current grid destination plays the role of grid source in next transformation (if any).
289Only element on which the transformation is performed is modified
290  \param [in] elementPositionInGrid position of element in grid
291  \param [in] transType transformation type
292*/
293void CGridTransformation::setUpGrid(int elementPositionInGrid, ETranformationType transType, int nbTransformation)
294{
295  if (!tempGrids_.empty() && (getNbAlgo()-1) == tempGrids_.size())
296  {
297    gridSource_ = tempGrids_[nbTransformation];
298    return;
299  }
300
301  std::vector<CAxis*> axisListDestP = gridDestination_->getAxis();
302  std::vector<CAxis*> axisListSrcP = gridSource_->getAxis(), axisSrc;
303
304  std::vector<CDomain*> domListDestP = gridDestination_->getDomains();
305  std::vector<CDomain*> domListSrcP = gridSource_->getDomains(), domainSrc;
306
307  int axisIndex = -1, domainIndex = -1;
308  switch (transType)
309  {
310    case TRANS_INTERPOLATE_DOMAIN:
311    case TRANS_ZOOM_DOMAIN:
312      domainIndex = elementPosition2DomainPositionInGrid_[elementPositionInGrid];
313      break;
314
315    case TRANS_INTERPOLATE_AXIS:
316    case TRANS_ZOOM_AXIS:
317    case TRANS_INVERSE_AXIS:
318      axisIndex =  elementPosition2AxisPositionInGrid_[elementPositionInGrid];
319      break;
320    default:
321      break;
322  }
323
324  for (int idx = 0; idx < axisListSrcP.size(); ++idx)
325  {
326    CAxis* axis = CAxis::createAxis();
327    if (axisIndex != idx) axis->axis_ref.setValue(axisListSrcP[idx]->getId());
328    else axis->axis_ref.setValue(axisListDestP[idx]->getId());
329    axis->solveRefInheritance(true);
330    axis->checkAttributesOnClient();
331    axisSrc.push_back(axis);
332  }
333
334  for (int idx = 0; idx < domListSrcP.size(); ++idx)
335  {
336    CDomain* domain = CDomain::createDomain();
337    if (domainIndex != idx) domain->domain_ref.setValue(domListSrcP[idx]->getId());
338    else domain->domain_ref.setValue(domListDestP[idx]->getId());
339    domain->solveRefInheritance(true);
340    domain->checkAttributesOnClient();
341    domainSrc.push_back(domain);
342  }
343
344  gridSource_ = CGrid::createGrid(domainSrc, axisSrc, gridDestination_->axis_domain_order);
345  gridSource_->computeGridGlobalDimension(domainSrc, axisSrc, gridDestination_->axis_domain_order);
346
347  tempGrids_.push_back(gridSource_);
348}
349
350/*!
351  Perform all transformations
352  For each transformation, there are some things to do:
353  -) Chose the correct algorithm by transformation type and position of element
354  -) Calculate the mapping of global index between the current grid source and grid destination
355  -) Calculate the mapping of global index between current grid DESTINATION and ORIGINAL grid SOURCE
356  -) Make current grid destination become grid source in the next transformation
357*/
358void CGridTransformation::computeAll(const std::vector<CArray<double,1>* >& dataAuxInputs, Time timeStamp)
359{
360  if (nbAlgos_ < 1) return;
361  if (!auxInputs_.empty() && !dynamicalTransformation_) { dynamicalTransformation_ = true; return; }
362  if (dynamicalTransformation_)
363  {
364    if (timeStamp_.insert(timeStamp).second)
365      DestinationIndexMap().swap(currentGridIndexToOriginalGridIndex_);  // Reset map
366    else
367      return;
368  }
369
370  CContext* context = CContext::getCurrent();
371  CContextClient* client = context->client;
372
373  ListAlgoType::const_iterator itb = listAlgos_.begin(),
374                               ite = listAlgos_.end(), it;
375
376  CGenericAlgorithmTransformation* algo = 0;
377  int nbAgloTransformation = 0; // Only count for executed transformation. Generate domain is a special one, not executed in the list
378  for (it = itb; it != ite; ++it)
379  {
380    int elementPositionInGrid = it->first;
381    ETranformationType transType = (it->second).first;
382    int transformationOrder = (it->second).second;
383    DestinationIndexMap globaIndexWeightFromDestToSource;
384
385    // First of all, select an algorithm
386    if (!dynamicalTransformation_ || (algoTransformation_.size() < listAlgos_.size()))
387    {
388      selectAlgo(elementPositionInGrid, transType, transformationOrder, algoTypes_[std::distance(itb, it)]);
389      algo = algoTransformation_.back();
390    }
391    else
392      algo = algoTransformation_[std::distance(itb, it)];
393
394    if (0 != algo) // Only registered transformation can be executed
395    {
396      algo->computeIndexSourceMapping(dataAuxInputs);
397
398      // Recalculate the distribution of grid destination
399      CDistributionClient distributionClientDest(client->clientRank, gridDestination_);
400      const CDistributionClient::GlobalLocalDataMap& globalLocalIndexGridDestSendToServer = distributionClientDest.getGlobalLocalDataSendToServer();
401
402      // ComputeTransformation of global index of each element
403      std::vector<int> gridDestinationDimensionSize = gridDestination_->getGlobalDimension();
404      std::vector<int> gridSrcDimensionSize = gridSource_->getGlobalDimension();
405      int elementPosition = it->first;
406      algo->computeGlobalSourceIndex(elementPosition,
407                                     gridDestinationDimensionSize,
408                                     gridSrcDimensionSize,
409                                     globalLocalIndexGridDestSendToServer,
410                                     globaIndexWeightFromDestToSource);
411
412      // Compute transformation of global indexes among grids
413      computeTransformationMapping(globaIndexWeightFromDestToSource);
414
415      // Update number of local index on each transformation
416      nbLocalIndexOnGridDest_.push_back(globalLocalIndexGridDestSendToServer.size());
417
418      if (1 < nbAlgos_)
419      {
420        // Now grid destination becomes grid source in a new transformation
421        if (nbAgloTransformation != (nbAlgos_-1)) setUpGrid(elementPositionInGrid, transType, nbAgloTransformation);
422      }
423      ++nbAgloTransformation;
424    }
425  }
426}
427
428/*!
429  Compute exchange index between grid source and grid destination
430  \param [in] globalIndexWeightFromDestToSource global index mapping between grid destination and grid source
431*/
432void CGridTransformation::computeTransformationMapping(const DestinationIndexMap& globalIndexWeightFromDestToSource)
433{
434  CContext* context = CContext::getCurrent();
435  CContextClient* client = context->client;
436
437  CTransformationMapping transformationMap(gridDestination_, gridSource_);
438
439  transformationMap.computeTransformationMapping(globalIndexWeightFromDestToSource);
440
441  const CTransformationMapping::ReceivedIndexMap& globalIndexToReceive = transformationMap.getGlobalIndexReceivedOnGridDestMapping();
442  CTransformationMapping::ReceivedIndexMap::const_iterator itbMapRecv, itMapRecv, iteMapRecv;
443  itbMapRecv = globalIndexToReceive.begin();
444  iteMapRecv = globalIndexToReceive.end();
445  localIndexToReceiveOnGridDest_.push_back(RecvIndexGridDestinationMap());
446  RecvIndexGridDestinationMap& recvTmp = localIndexToReceiveOnGridDest_.back();
447  for (itMapRecv = itbMapRecv; itMapRecv != iteMapRecv; ++itMapRecv)
448  {
449    int sourceRank = itMapRecv->first;
450    int numGlobalIndex = (itMapRecv->second).size();
451    recvTmp[sourceRank].resize(numGlobalIndex);
452    for (int i = 0; i < numGlobalIndex; ++i)
453    {
454      recvTmp[sourceRank][i] = make_pair((itMapRecv->second)[i].localIndex,(itMapRecv->second)[i].weight);
455    }
456  }
457
458  // Find out local index on grid source (to send)
459  const CTransformationMapping::SentIndexMap& globalIndexToSend = transformationMap.getGlobalIndexSendToGridDestMapping();
460  CTransformationMapping::SentIndexMap::const_iterator itbMap, itMap, iteMap;
461  itbMap = globalIndexToSend.begin();
462  iteMap = globalIndexToSend.end();
463  localIndexToSendFromGridSource_.push_back(SendingIndexGridSourceMap());
464  SendingIndexGridSourceMap& tmpSend = localIndexToSendFromGridSource_.back();
465  for (itMap = itbMap; itMap != iteMap; ++itMap)
466  {
467    int destRank = itMap->first;
468    int vecSize = itMap->second.size();
469    tmpSend[destRank].resize(vecSize);
470    for (int idx = 0; idx < vecSize; ++idx)
471    {
472      tmpSend[destRank](idx) = itMap->second[idx].first;
473    }
474  }
475}
476
477bool CGridTransformation::isSpecialTransformation(ETranformationType transType)
478{
479  bool res;
480  switch (transType)
481  {
482    case TRANS_GENERATE_RECTILINEAR_DOMAIN:
483     res = true;
484     break;
485    default:
486     res = false;
487     break;
488  }
489
490  return res;
491}
492
493/*!
494  Local index of data which need sending from the grid source
495  \return local index of data
496*/
497const std::list<CGridTransformation::SendingIndexGridSourceMap>& CGridTransformation::getLocalIndexToSendFromGridSource() const
498{
499  return localIndexToSendFromGridSource_;
500}
501
502/*!
503  Local index of data which will be received on the grid destination
504  \return local index of data
505*/
506const std::list<CGridTransformation::RecvIndexGridDestinationMap>& CGridTransformation::getLocalIndexToReceiveOnGridDest() const
507{
508  return localIndexToReceiveOnGridDest_;
509}
510
511const std::list<size_t>& CGridTransformation::getNbLocalIndexToReceiveOnGridDest() const
512{
513  return nbLocalIndexOnGridDest_;
514}
515
516}
Note: See TracBrowser for help on using the repository browser.