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

Last change on this file since 842 was 842, checked in by mhnguyen, 6 years ago

Small improvements on transformation mapping

+) Remove complex structure of transformation mapping by a simpler one

Test
+) On Curie
+) All tests pass

File size: 20.4 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      if (1 < nbAlgos_)
416      {
417        // Now grid destination becomes grid source in a new transformation
418        if (nbAgloTransformation != (nbAlgos_-1)) setUpGrid(elementPositionInGrid, transType, nbAgloTransformation);
419      }
420      ++nbAgloTransformation;
421    }
422  }
423}
424
425/*!
426  Compute exchange index between grid source and grid destination
427  \param [in] globalIndexWeightFromDestToSource global index mapping between grid destination and grid source
428*/
429void CGridTransformation::computeTransformationMapping(const DestinationIndexMap& globalIndexWeightFromDestToSource)
430{
431  CContext* context = CContext::getCurrent();
432  CContextClient* client = context->client;
433
434  CTransformationMapping transformationMap(gridDestination_, gridSource_);
435
436  transformationMap.computeTransformationMapping(globalIndexWeightFromDestToSource);
437
438  const CTransformationMapping::ReceivedIndexMap& globalIndexToReceive = transformationMap.getGlobalIndexReceivedOnGridDestMapping();
439  CTransformationMapping::ReceivedIndexMap::const_iterator itbMapRecv, itMapRecv, iteMapRecv;
440  itbMapRecv = globalIndexToReceive.begin();
441  iteMapRecv = globalIndexToReceive.end();
442  nbLocalIndexOnGridDest_.push_back(globalIndexWeightFromDestToSource.size());
443  localIndexToReceiveOnGridDest_.push_back(RecvIndexGridDestinationMap());
444  RecvIndexGridDestinationMap& recvTmp = localIndexToReceiveOnGridDest_.back();
445  for (itMapRecv = itbMapRecv; itMapRecv != iteMapRecv; ++itMapRecv)
446  {
447    int sourceRank = itMapRecv->first;
448    int numGlobalIndex = (itMapRecv->second).size();
449    recvTmp[sourceRank].resize(numGlobalIndex);
450    for (int i = 0; i < numGlobalIndex; ++i)
451    {
452      recvTmp[sourceRank][i] = make_pair((itMapRecv->second)[i].localIndex,(itMapRecv->second)[i].weight);
453    }
454  }
455
456  // Find out local index on grid source (to send)
457  const CTransformationMapping::SentIndexMap& globalIndexToSend = transformationMap.getGlobalIndexSendToGridDestMapping();
458  CTransformationMapping::SentIndexMap::const_iterator itbMap, itMap, iteMap;
459  itbMap = globalIndexToSend.begin();
460  iteMap = globalIndexToSend.end();
461  localIndexToSendFromGridSource_.push_back(SendingIndexGridSourceMap());
462  SendingIndexGridSourceMap& tmpSend = localIndexToSendFromGridSource_.back();
463  for (itMap = itbMap; itMap != iteMap; ++itMap)
464  {
465    int destRank = itMap->first;
466    int vecSize = itMap->second.size();
467    tmpSend[destRank].resize(vecSize);
468    for (int idx = 0; idx < vecSize; ++idx)
469    {
470      tmpSend[destRank](idx) = itMap->second[idx].first;
471    }
472  }
473}
474
475bool CGridTransformation::isSpecialTransformation(ETranformationType transType)
476{
477  bool res;
478  switch (transType)
479  {
480    case TRANS_GENERATE_RECTILINEAR_DOMAIN:
481     res = true;
482     break;
483    default:
484     res = false;
485     break;
486  }
487
488  return res;
489}
490
491/*!
492  Local index of data which need sending from the grid source
493  \return local index of data
494*/
495const std::list<CGridTransformation::SendingIndexGridSourceMap>& CGridTransformation::getLocalIndexToSendFromGridSource() const
496{
497  return localIndexToSendFromGridSource_;
498}
499
500/*!
501  Local index of data which will be received on the grid destination
502  \return local index of data
503*/
504const std::list<CGridTransformation::RecvIndexGridDestinationMap>& CGridTransformation::getLocalIndexToReceiveOnGridDest() const
505{
506  return localIndexToReceiveOnGridDest_;
507}
508
509const std::list<size_t>& CGridTransformation::getNbLocalIndexToReceiveOnGridDest() const
510{
511  return nbLocalIndexOnGridDest_;
512}
513
514}
Note: See TracBrowser for help on using the repository browser.