source: XIOS/trunk/src/transformation/generic_algorithm_transformation.cpp @ 887

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

Adding a new type of element into grid: Scalar

+) Add a new node Scalar for xml
+) Make some change on writing scalar value
+) Reorganize some codes
+) Remove some redundant codes

Test
+) On Curie
+) All tests pass

File size: 21.3 KB
Line 
1/*!
2   \file generic_algorithm_transformation.hpp
3   \author Ha NGUYEN
4   \since 14 May 2015
5   \date 21 Mars 2016
6
7   \brief Interface for all transformation algorithms.
8 */
9#include "generic_algorithm_transformation.hpp"
10#include "context.hpp"
11#include "context_client.hpp"
12#include "client_client_dht_template.hpp"
13
14namespace xios {
15
16CGenericAlgorithmTransformation::CGenericAlgorithmTransformation()
17 : transformationMapping_(), transformationWeight_(), transformationPosition_(), idAuxInputs_()
18{
19}
20
21/*!
22  This function computes the global indexes of grid source, which the grid destination is in demand.
23  \param[in] elementPositionInGrid position of an element in a grid .E.g: if grid is composed of domain and axis (in order),
24                then position of axis in grid is 1 and domain is positioned at 0.
25  \param[in] gridSrc Grid source
26  \param[in] gridDst Grid destination
27  \param[in\out] globaIndexWeightFromSrcToDst mapping of each global index source and weight to index destination
28*/
29void CGenericAlgorithmTransformation::computeGlobalSourceIndex(int elementPositionInGrid,
30                                                               CGrid* gridSrc,
31                                                               CGrid* gridDst,
32                                                               SourceDestinationIndexMap& globaIndexWeightFromSrcToDst)
33 {
34  CContext* context = CContext::getCurrent();
35  CContextClient* client = context->client;
36  int nbClient = client->clientSize;
37
38  typedef boost::unordered_map<int, std::vector<std::pair<int,double> > > SrcToDstMap;
39
40  size_t indexSrcSize = 0;
41  for (size_t idxTrans = 0; idxTrans < transformationMapping_.size(); ++idxTrans)
42  {
43    TransformationIndexMap::const_iterator itbTransMap = transformationMapping_[idxTrans].begin(), itTransMap,
44                                           iteTransMap = transformationMapping_[idxTrans].end();
45    TransformationWeightMap::const_iterator itbTransWeight = transformationWeight_[idxTrans].begin(), itTransWeight;
46
47    itTransWeight = itbTransWeight;
48    for (itTransMap = itbTransMap; itTransMap != iteTransMap; ++itTransMap, ++itTransWeight)
49    {
50       indexSrcSize += (itTransMap->second).size();
51    }
52  }
53
54  bool isTransPosEmpty = transformationPosition_.empty();
55  CArray<size_t,1> transPos;
56  if (!isTransPosEmpty) transPos.resize(transformationMapping_.size());
57  CArray<size_t,1> indexSrc(indexSrcSize);
58  indexSrcSize = 0;
59  for (size_t idxTrans = 0; idxTrans < transformationMapping_.size(); ++idxTrans)
60  {
61    TransformationIndexMap::const_iterator itbTransMap = transformationMapping_[idxTrans].begin(), itTransMap,
62                                           iteTransMap = transformationMapping_[idxTrans].end();
63    TransformationWeightMap::const_iterator itbTransWeight = transformationWeight_[idxTrans].begin(), itTransWeight;
64
65    // Build mapping between global source element index and global destination element index.
66    itTransWeight = itbTransWeight;
67    for (itTransMap = itbTransMap; itTransMap != iteTransMap; ++itTransMap, ++itTransWeight)
68    {
69      const std::vector<int>& srcIndex = itTransMap->second;
70      for (int idx = 0; idx < srcIndex.size(); ++idx)
71      {
72        indexSrc(indexSrcSize) = srcIndex[idx];
73        ++indexSrcSize;
74      }
75    }
76
77    if (!isTransPosEmpty)
78    {
79      TransformationPositionMap::const_iterator itPosMap = transformationPosition_[idxTrans].begin();
80      transPos(idxTrans) = itPosMap->second[0];
81    }
82  }
83
84  std::vector<CAxis*> axisListDestP = gridDst->getAxis();
85  std::vector<CDomain*> domainListDestP = gridDst->getDomains();
86  CArray<int,1> axisDomainDstOrder = gridDst->axis_domain_order;
87  std::vector<CAxis*> axisListSrcP = gridSrc->getAxis();
88  std::vector<CDomain*> domainListSrcP = gridSrc->getDomains();
89
90  // Find out global index source of transformed element on corresponding process.
91  std::vector<boost::unordered_map<int,std::vector<size_t> > > globalElementIndexOnProc(axisDomainDstOrder.numElements());
92  CClientClientDHTInt::Index2VectorInfoTypeMap globalIndexOfTransformedElementOnProc;
93  int axisIndex = 0, domainIndex = 0;
94  for (int idx = 0; idx < axisDomainDstOrder.numElements(); ++idx)
95  {
96    if (idx == elementPositionInGrid)
97      computeExchangeGlobalIndex(indexSrc, globalIndexOfTransformedElementOnProc); //globalElementIndexOnProc[idx]);
98    if (2 == axisDomainDstOrder(idx)) // It's domain
99    {
100      if (idx != elementPositionInGrid)
101        computeExchangeDomainIndex(domainListDestP[domainIndex],
102                                   domainListSrcP[domainIndex],
103                                   transPos,
104                                   globalElementIndexOnProc[idx]);
105      ++domainIndex;
106
107    }
108    else if (1 == axisDomainDstOrder(idx))//it's an axis
109    {
110      if (idx != elementPositionInGrid)
111        computeExchangeAxisIndex(axisListDestP[axisIndex],
112                                 axisListSrcP[axisIndex],
113                                 transPos,
114                                 globalElementIndexOnProc[idx]);
115      ++axisIndex;
116
117    }
118  }
119
120  if (!isTransPosEmpty)
121  {
122    for (int idx = 0; idx < globalElementIndexOnProc.size(); ++idx)
123    {
124      if (idx != elementPositionInGrid)
125      {
126        boost::unordered_map<int,std::vector<size_t> >::iterator itb = globalElementIndexOnProc[idx].begin(), it,
127                                                                 ite = globalElementIndexOnProc[idx].end();
128        for (it = itb; it != ite; ++it) it->second.resize(1);
129      }
130    }
131  }
132
133  for (size_t idxTrans = 0; idxTrans < transformationMapping_.size(); ++idxTrans)
134  {
135    TransformationIndexMap::const_iterator itbTransMap = transformationMapping_[idxTrans].begin(), itTransMap,
136                                           iteTransMap = transformationMapping_[idxTrans].end();
137    TransformationWeightMap::const_iterator itbTransWeight = transformationWeight_[idxTrans].begin(), itTransWeight;
138    SrcToDstMap src2DstMap;
139    src2DstMap.rehash(std::ceil(transformationMapping_[idxTrans].size()/src2DstMap.max_load_factor()));
140
141    // Build mapping between global source element index and global destination element index.
142    boost::unordered_map<int,std::vector<size_t> >().swap(globalElementIndexOnProc[elementPositionInGrid]);
143    boost::unordered_map<int,int> tmpCounter;
144    itTransWeight = itbTransWeight;
145    for (itTransMap = itbTransMap; itTransMap != iteTransMap; ++itTransMap, ++itTransWeight)
146    {
147      const std::vector<int>& srcIndex = itTransMap->second;
148      const std::vector<double>& weight = itTransWeight->second;
149      for (int idx = 0; idx < srcIndex.size(); ++idx)
150      {
151        src2DstMap[srcIndex[idx]].push_back(make_pair(itTransMap->first, weight[idx]));
152        if (1 == globalIndexOfTransformedElementOnProc.count(srcIndex[idx]) && (0 == tmpCounter.count(srcIndex[idx])))
153        {
154          tmpCounter[srcIndex[idx]] = 1;
155          std::vector<int>& srcProc = globalIndexOfTransformedElementOnProc[srcIndex[idx]];
156          for (int j = 0; j < srcProc.size(); ++j)
157            globalElementIndexOnProc[elementPositionInGrid][srcProc[j]].push_back(srcIndex[idx]);
158        }
159      }
160    }
161
162    if (!isTransPosEmpty)
163    {
164      for (int idx = 0; idx < globalElementIndexOnProc.size(); ++idx)
165      {
166        if (idx != elementPositionInGrid)
167        {
168          boost::unordered_map<int,std::vector<size_t> >::iterator itb = globalElementIndexOnProc[idx].begin(), it,
169                                                                   ite = globalElementIndexOnProc[idx].end();
170          for (it = itb; it != ite; ++it) it->second[0] = transPos(idxTrans);
171        }
172      }
173    }
174
175    std::vector<std::vector<bool> > elementOnProc(axisDomainDstOrder.numElements(), std::vector<bool>(nbClient, false));
176    boost::unordered_map<int,std::vector<size_t> >::const_iterator it, itb, ite;
177    for (int idx = 0; idx < globalElementIndexOnProc.size(); ++idx)
178    {
179      itb = globalElementIndexOnProc[idx].begin();
180      ite = globalElementIndexOnProc[idx].end();
181      for (it = itb; it != ite; ++it) elementOnProc[idx][it->first] = true;
182    }
183
184    // Determine procs which contain global source index
185    std::vector<bool> intersectedProc(nbClient, true);
186    for (int idx = 0; idx < axisDomainDstOrder.numElements(); ++idx)
187    {
188      std::transform(elementOnProc[idx].begin(), elementOnProc[idx].end(),
189                     intersectedProc.begin(), intersectedProc.begin(),
190                     std::logical_and<bool>());
191    }
192
193    std::vector<int> srcRank;
194    for (int idx = 0; idx < nbClient; ++idx)
195    {
196      if (intersectedProc[idx]) srcRank.push_back(idx);
197    }
198
199    // Ok, now compute global index of grid source and ones of grid destination
200    computeGlobalGridIndexMapping(elementPositionInGrid,
201                                  srcRank,
202                                  src2DstMap,
203                                  gridSrc,
204                                  gridDst,
205                                  globalElementIndexOnProc,
206                                  globaIndexWeightFromSrcToDst);
207  }
208 }
209
210/*!
211  Compute mapping of global index of grid source and grid destination
212  \param [in] elementPositionInGrid position of element in grid. E.x: grid composed of domain and axis, domain has position 0 and axis 1.
213  \param [in] srcRank rank of client from which we demand global index of element source
214  \param [in] src2DstMap mapping of global index of element source and global index of element destination
215  \param[in] gridSrc Grid source
216  \param[in] gridDst Grid destination
217  \param[in] globalElementIndexOnProc Global index of element source on different client rank
218  \param[out] globaIndexWeightFromSrcToDst Mapping of global index of grid source and grid destination
219*/
220void CGenericAlgorithmTransformation::computeGlobalGridIndexMapping(int elementPositionInGrid,
221                                                                   const std::vector<int>& srcRank,
222                                                                   boost::unordered_map<int, std::vector<std::pair<int,double> > >& src2DstMap,
223                                                                   CGrid* gridSrc,
224                                                                   CGrid* gridDst,
225                                                                   std::vector<boost::unordered_map<int,std::vector<size_t> > >& globalElementIndexOnProc,
226                                                                   SourceDestinationIndexMap& globaIndexWeightFromSrcToDst)
227{
228  std::vector<CDomain*> domainListSrcP = gridSrc->getDomains();
229  std::vector<CAxis*> axisListSrcP = gridSrc->getAxis();
230  std::vector<CScalar*> scalarListSrcP = gridSrc->getScalars();
231  CArray<int,1> axisDomainSrcOrder = gridSrc->axis_domain_order;
232
233  size_t nbElement = axisDomainSrcOrder.numElements();
234  std::vector<size_t> nGlobSrc(nbElement);
235  size_t globalSrcSize = 1;
236  int domainIndex = 0, axisIndex = 0, scalarIndex = 0;
237  for (int idx = 0; idx < nbElement; ++idx)
238  {
239    nGlobSrc[idx] = globalSrcSize;
240    int elementDimension = axisDomainSrcOrder(idx);
241
242    // If this is a domain
243    if (2 == elementDimension)
244    {
245      globalSrcSize *= domainListSrcP[domainIndex]->nj_glo.getValue() * domainListSrcP[domainIndex]->ni_glo.getValue();
246      ++domainIndex;
247    }
248    else if (1 == elementDimension) // So it's an axis
249    {
250      globalSrcSize *= axisListSrcP[axisIndex]->n_glo.getValue();
251      ++axisIndex;
252    }
253    else
254    {
255      globalSrcSize *= 1;
256      ++scalarIndex;
257    }
258  }
259
260  std::vector<CDomain*> domainListDestP = gridDst->getDomains();
261  std::vector<CAxis*> axisListDestP = gridDst->getAxis();
262  std::vector<CScalar*> scalarListDestP = gridDst->getScalars();
263  CArray<int,1> axisDomainDstOrder = gridDst->axis_domain_order;
264
265  std::vector<size_t> nGlobDst(nbElement);
266  size_t globalDstSize = 1;
267  domainIndex = axisIndex = scalarIndex = 0;
268  for (int idx = 0; idx < nbElement; ++idx)
269  {
270    nGlobDst[idx] = globalDstSize;
271    int elementDimension = axisDomainSrcOrder(idx);
272
273    // If this is a domain
274    if (2 == elementDimension)
275    {
276      globalDstSize *= domainListDestP[domainIndex]->nj_glo.getValue() * domainListDestP[domainIndex]->ni_glo.getValue();
277      ++domainIndex;
278    }
279    else if (1 == elementDimension) // So it's an axis
280    {
281      globalDstSize *= axisListDestP[axisIndex]->n_glo.getValue();
282      ++axisIndex;
283    }
284    else
285    {
286      globalDstSize *= 1;
287      ++scalarIndex;
288    }
289  }
290
291  for (int i = 0; i < srcRank.size(); ++i)
292  {
293    size_t ssize = 1;
294    int rankSrc = srcRank[i];
295    for (int idx = 0; idx < nbElement; ++idx)
296    {
297      ssize *= (globalElementIndexOnProc[idx][rankSrc]).size();
298    }
299
300    std::vector<int> idxLoop(nbElement,0);
301    std::vector<int> currentIndexSrc(nbElement, 0);
302    std::vector<int> currentIndexDst(nbElement, 0);
303    int innnerLoopSize = (globalElementIndexOnProc[0])[rankSrc].size();
304    size_t idx = 0;
305    while (idx < ssize)
306    {
307      for (int ind = 0; ind < nbElement; ++ind)
308      {
309        if (idxLoop[ind] == (globalElementIndexOnProc[ind])[rankSrc].size())
310        {
311          idxLoop[ind] = 0;
312          ++idxLoop[ind+1];
313        }
314
315        currentIndexDst[ind] = currentIndexSrc[ind] = (globalElementIndexOnProc[ind])[rankSrc][idxLoop[ind]];
316      }
317
318      for (int ind = 0; ind < innnerLoopSize; ++ind)
319      {
320        currentIndexDst[0] = currentIndexSrc[0] = (globalElementIndexOnProc[0])[rankSrc][ind];
321        int globalElementDstIndexSize = 0;
322        if (1 == src2DstMap.count(currentIndexSrc[elementPositionInGrid]))
323        {
324          globalElementDstIndexSize = src2DstMap[currentIndexSrc[elementPositionInGrid]].size();
325        }
326
327        std::vector<size_t> globalDstVecIndex(globalElementDstIndexSize,0);
328        size_t globalSrcIndex = 0;
329        for (int idxElement = 0; idxElement < nbElement; ++idxElement)
330        {
331          if (idxElement == elementPositionInGrid)
332          {
333            for (int k = 0; k < globalElementDstIndexSize; ++k)
334            {
335              globalDstVecIndex[k] += src2DstMap[currentIndexSrc[elementPositionInGrid]][k].first * nGlobDst[idxElement];
336            }
337          }
338          else
339          {
340            for (int k = 0; k < globalElementDstIndexSize; ++k)
341            {
342              globalDstVecIndex[k] += currentIndexDst[idxElement] * nGlobDst[idxElement];
343            }
344          }
345          globalSrcIndex += currentIndexSrc[idxElement] * nGlobSrc[idxElement];
346        }
347
348        for (int k = 0; k < globalElementDstIndexSize; ++k)
349        {
350          globaIndexWeightFromSrcToDst[rankSrc][globalSrcIndex].push_back(make_pair(globalDstVecIndex[k],src2DstMap[currentIndexSrc[elementPositionInGrid]][k].second ));
351        }
352        ++idxLoop[0];
353      }
354      idx += innnerLoopSize;
355    }
356  }
357}
358
359/*!
360  Find out proc and global index of axis source which axis destination is on demande
361  \param[in] axisDst Axis destination
362  \param[in] axisSrc Axis source
363  \param[in] destGlobalIndexPositionInGrid Relative position of axis corresponds to other element of grid.
364  \param[out] globalAxisIndexOnProc Global index of axis source on different procs
365*/
366void CGenericAlgorithmTransformation::computeExchangeAxisIndex(CAxis* axisDst,
367                                                               CAxis* axisSrc,
368                                                               CArray<size_t,1>& destGlobalIndexPositionInGrid,
369                                                               boost::unordered_map<int,std::vector<size_t> >& globalAxisIndexOnProc)
370{
371  CContext* context = CContext::getCurrent();
372  CContextClient* client=context->client;
373  int clientRank = client->clientRank;
374  int clientSize = client->clientSize;
375
376  size_t globalIndex;
377  int nIndexSize = axisSrc->index.numElements();
378  CClientClientDHTInt::Index2VectorInfoTypeMap globalIndex2ProcRank;
379  globalIndex2ProcRank.rehash(std::ceil(nIndexSize/globalIndex2ProcRank.max_load_factor()));
380  for (int idx = 0; idx < nIndexSize; ++idx)
381  {
382    globalIndex = axisSrc->index(idx);
383    globalIndex2ProcRank[globalIndex].push_back(clientRank);
384  }
385
386  CClientClientDHTInt dhtIndexProcRank(globalIndex2ProcRank, client->intraComm);
387  CArray<size_t,1> globalAxisIndex(axisDst->index.numElements());
388  for (int idx = 0; idx < globalAxisIndex.numElements(); ++idx)
389  {
390    globalAxisIndex(idx) = axisDst->index(idx);
391  }
392  dhtIndexProcRank.computeIndexInfoMapping(globalAxisIndex);
393
394  std::vector<int> countIndex(clientSize,0);
395  const CClientClientDHTInt::Index2VectorInfoTypeMap& computedGlobalIndexOnProc = dhtIndexProcRank.getInfoIndexMap();
396  CClientClientDHTInt::Index2VectorInfoTypeMap::const_iterator itb = computedGlobalIndexOnProc.begin(), it,
397                                                               ite = computedGlobalIndexOnProc.end();
398  for (it = itb; it != ite; ++it)
399  {
400    const std::vector<int>& procList = it->second;
401    for (int idx = 0; idx < procList.size(); ++idx) ++countIndex[procList[idx]];
402  }
403
404  globalAxisIndexOnProc.rehash(std::ceil(clientSize/globalAxisIndexOnProc.max_load_factor()));
405  for (int idx = 0; idx < clientSize; ++idx)
406  {
407    if (0 != countIndex[idx])
408    {
409      globalAxisIndexOnProc[idx].resize(countIndex[idx]);
410      countIndex[idx] = 0;
411    }
412  }
413
414  for (it = itb; it != ite; ++it)
415  {
416    const std::vector<int>& procList = it->second;
417    for (int idx = 0; idx < procList.size(); ++idx)
418    {
419      globalAxisIndexOnProc[procList[idx]][countIndex[procList[idx]]] = it->first;
420      ++countIndex[procList[idx]];
421    }
422  }
423}
424
425/*!
426  Find out proc and global index of domain source which domain destination is on demande
427  \param[in] domainDst Domain destination
428  \param[in] domainSrc Domain source
429  \param[in] destGlobalIndexPositionInGrid Relative position of domain corresponds to other element of grid.
430  \param[out] globalDomainIndexOnProc Global index of domain source on different procs
431*/
432void CGenericAlgorithmTransformation::computeExchangeDomainIndex(CDomain* domainDst,
433                                                                 CDomain* domainSrc,
434                                                                 CArray<size_t,1>& destGlobalIndexPositionInGrid,
435                                                                 boost::unordered_map<int,std::vector<size_t> >& globalDomainIndexOnProc)
436{
437  CContext* context = CContext::getCurrent();
438  CContextClient* client=context->client;
439  int clientRank = client->clientRank;
440  int clientSize = client->clientSize;
441
442  int niGlobSrc = domainSrc->ni_glo.getValue();
443  size_t globalIndex;
444  int i_ind, j_ind;
445  int nIndexSize = (destGlobalIndexPositionInGrid.isEmpty()) ? domainSrc->i_index.numElements()
446                                                             : destGlobalIndexPositionInGrid.numElements();
447  CClientClientDHTInt::Index2VectorInfoTypeMap globalIndex2ProcRank;
448  globalIndex2ProcRank.rehash(std::ceil(nIndexSize/globalIndex2ProcRank.max_load_factor()));
449  if (destGlobalIndexPositionInGrid.isEmpty())
450  {
451    for (int idx = 0; idx < nIndexSize; ++idx)
452    {
453      i_ind=domainSrc->i_index(idx) ;
454      j_ind=domainSrc->j_index(idx) ;
455
456      globalIndex = i_ind + j_ind * niGlobSrc;
457      globalIndex2ProcRank[globalIndex].resize(1);
458      globalIndex2ProcRank[globalIndex][0] = clientRank;
459    }
460  }
461  else
462  {
463    for (int idx = 0; idx < nIndexSize; ++idx)
464    {
465      globalIndex2ProcRank[destGlobalIndexPositionInGrid(idx)].push_back(clientRank);
466    }
467  }
468
469  CArray<size_t,1> globalDomainIndex;
470  if (destGlobalIndexPositionInGrid.isEmpty())
471  {
472    int niGlobDst = domainDst->ni_glo.getValue();
473    globalDomainIndex.resize(domainDst->i_index.numElements());
474    nIndexSize = domainDst->i_index.numElements();
475
476    for (int idx = 0; idx < nIndexSize; ++idx)
477    {
478      i_ind=domainDst->i_index(idx) ;
479      j_ind=domainDst->j_index(idx) ;
480
481      globalDomainIndex(idx) = i_ind + j_ind * niGlobDst;
482    }
483  }
484  else
485  {
486    globalDomainIndex.reference(destGlobalIndexPositionInGrid);
487  }
488
489  CClientClientDHTInt dhtIndexProcRank(globalIndex2ProcRank, client->intraComm);
490  dhtIndexProcRank.computeIndexInfoMapping(globalDomainIndex);
491
492  std::vector<int> countIndex(clientSize,0);
493  const CClientClientDHTInt::Index2VectorInfoTypeMap& computedGlobalIndexOnProc = dhtIndexProcRank.getInfoIndexMap();
494  CClientClientDHTInt::Index2VectorInfoTypeMap::const_iterator itb = computedGlobalIndexOnProc.begin(), it,
495                                                               ite = computedGlobalIndexOnProc.end();
496  for (it = itb; it != ite; ++it)
497  {
498    const std::vector<int>& procList = it->second;
499    for (int idx = 0; idx < procList.size(); ++idx) ++countIndex[procList[idx]];
500  }
501
502  globalDomainIndexOnProc.rehash(std::ceil(clientSize/globalDomainIndexOnProc.max_load_factor()));
503  for (int idx = 0; idx < clientSize; ++idx)
504  {
505    if (0 != countIndex[idx])
506    {
507      globalDomainIndexOnProc[idx].resize(countIndex[idx]);
508      countIndex[idx] = 0;
509    }
510  }
511
512  for (it = itb; it != ite; ++it)
513  {
514    const std::vector<int>& procList = it->second;
515    for (int idx = 0; idx < procList.size(); ++idx)
516    {
517      globalDomainIndexOnProc[procList[idx]][countIndex[procList[idx]]] = it->first;
518      ++countIndex[procList[idx]];
519    }
520  }
521}
522
523/*!
524  Compute index mapping between element source and element destination with an auxiliary inputs which determine
525position of each mapped index in global index of grid destination.
526  \param [in] dataAuxInputs auxiliary inputs
527*/
528void CGenericAlgorithmTransformation::computeIndexSourceMapping(const std::vector<CArray<double,1>* >& dataAuxInputs)
529{
530  computeIndexSourceMapping_(dataAuxInputs);
531}
532
533std::vector<StdString> CGenericAlgorithmTransformation::getIdAuxInputs()
534{
535  return idAuxInputs_;
536}
537
538}
Note: See TracBrowser for help on using the repository browser.