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

Last change on this file since 868 was 868, checked in by mhnguyen, 8 years ago

Fixing a stupid bug in transformation which sometimes makes wrong result

Test
+) On Curie
+) All results are correct
File size: 20.6 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<bool,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 (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 //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                                  gridDst,
204                                  gridSrc,
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<CAxis*> axisListDestP = gridDst->getAxis();
229  std::vector<CDomain*> domainListDestP = gridDst->getDomains();
230  CArray<bool,1> axisDomainDstOrder = gridDst->axis_domain_order;
231  std::vector<CAxis*> axisListSrcP = gridSrc->getAxis();
232  std::vector<CDomain*> domainListSrcP = gridSrc->getDomains();
233  CArray<bool,1> axisDomainSrcOrder = gridDst->axis_domain_order;
234  size_t nbElement = axisDomainSrcOrder.numElements();
235  std::vector<size_t> nGlobSrc(nbElement), nGlobDst(nbElement);
236  size_t globalSrcSize = 1, globalDstSize = 1;
237  int domainIndex = 0;
238  int axisIndex = 0;
239  for (int idx = 0; idx < nbElement; ++idx)
240  {
241    nGlobSrc[idx] = globalSrcSize;
242    nGlobDst[idx] = globalDstSize;
243    bool isDomain = axisDomainSrcOrder(idx);
244
245    // If this is a domain
246    if (isDomain)
247    {
248      globalSrcSize *= domainListSrcP[domainIndex]->nj_glo.getValue() * domainListSrcP[domainIndex]->ni_glo.getValue();
249      globalDstSize *= domainListDestP[domainIndex]->nj_glo.getValue() * domainListDestP[domainIndex]->ni_glo.getValue();
250      ++domainIndex;
251    }
252    else // So it's an axis
253    {
254      globalSrcSize *= axisListSrcP[axisIndex]->n_glo.getValue();
255      globalDstSize *= axisListDestP[axisIndex]->n_glo.getValue();
256      ++axisIndex;
257    }
258  }
259
260  for (int i = 0; i < srcRank.size(); ++i)
261  {
262    size_t ssize = 1;
263    int rankSrc = srcRank[i];
264    for (int idx = 0; idx < nbElement; ++idx)
265    {
266      ssize *= (globalElementIndexOnProc[idx][rankSrc]).size();
267    }
268
269    std::vector<int> idxLoop(nbElement,0);
270    std::vector<int> currentIndexSrc(nbElement, 0);
271    std::vector<int> currentIndexDst(nbElement, 0);
272    int innnerLoopSize = (globalElementIndexOnProc[0])[rankSrc].size();
273    size_t idx = 0;
274    while (idx < ssize)
275    {
276      for (int ind = 0; ind < nbElement; ++ind)
277      {
278        if (idxLoop[ind] == (globalElementIndexOnProc[ind])[rankSrc].size())
279        {
280          idxLoop[ind] = 0;
281          ++idxLoop[ind+1];
282        }
283
284        currentIndexDst[ind] = currentIndexSrc[ind] = (globalElementIndexOnProc[ind])[rankSrc][idxLoop[ind]];
285      }
286
287      for (int ind = 0; ind < innnerLoopSize; ++ind)
288      {
289        currentIndexDst[0] = currentIndexSrc[0] = (globalElementIndexOnProc[0])[rankSrc][ind];
290        int globalElementDstIndexSize = 0;
291        if (1 == src2DstMap.count(currentIndexSrc[elementPositionInGrid]))
292        {
293          globalElementDstIndexSize = src2DstMap[currentIndexSrc[elementPositionInGrid]].size();
294        }
295
296        std::vector<size_t> globalDstVecIndex(globalElementDstIndexSize,0);
297        size_t globalSrcIndex = 0;
298        for (int idxElement = 0; idxElement < nbElement; ++idxElement)
299        {
300          if (idxElement == elementPositionInGrid)
301          {
302            for (int k = 0; k < globalElementDstIndexSize; ++k)
303            {
304              globalDstVecIndex[k] += src2DstMap[currentIndexSrc[elementPositionInGrid]][k].first * nGlobDst[idxElement];
305            }
306          }
307          else
308          {
309            for (int k = 0; k < globalElementDstIndexSize; ++k)
310            {
311              globalDstVecIndex[k] += currentIndexDst[idxElement] * nGlobDst[idxElement];
312            }
313          }
314          globalSrcIndex += currentIndexSrc[idxElement] * nGlobSrc[idxElement];
315        }
316
317        for (int k = 0; k < globalElementDstIndexSize; ++k)
318        {
319          globaIndexWeightFromSrcToDst[rankSrc][globalSrcIndex].push_back(make_pair(globalDstVecIndex[k],src2DstMap[currentIndexSrc[elementPositionInGrid]][k].second ));
320        }
321        ++idxLoop[0];
322      }
323      idx += innnerLoopSize;
324    }
325  }
326}
327
328/*!
329  Find out proc and global index of axis source which axis destination is on demande
330  \param[in] axisDst Axis destination
331  \param[in] axisSrc Axis source
332  \param[in] destGlobalIndexPositionInGrid Relative position of axis corresponds to other element of grid.
333  \param[out] globalAxisIndexOnProc Global index of axis source on different procs
334*/
335void CGenericAlgorithmTransformation::computeExchangeAxisIndex(CAxis* axisDst,
336                                                               CAxis* axisSrc,
337                                                               CArray<size_t,1>& destGlobalIndexPositionInGrid,
338                                                               boost::unordered_map<int,std::vector<size_t> >& globalAxisIndexOnProc)
339{
340  CContext* context = CContext::getCurrent();
341  CContextClient* client=context->client;
342  int clientRank = client->clientRank;
343  int clientSize = client->clientSize;
344
345  size_t globalIndex;
346  int nIndexSize = axisSrc->index.numElements();
347  CClientClientDHTInt::Index2VectorInfoTypeMap globalIndex2ProcRank;
348  globalIndex2ProcRank.rehash(std::ceil(nIndexSize/globalIndex2ProcRank.max_load_factor()));
349  for (int idx = 0; idx < nIndexSize; ++idx)
350  {
351    globalIndex = axisSrc->index(idx);
352    globalIndex2ProcRank[globalIndex].push_back(clientRank);
353  }
354
355  CClientClientDHTInt dhtIndexProcRank(globalIndex2ProcRank, client->intraComm);
356  CArray<size_t,1> globalAxisIndex(axisDst->index.numElements());
357  for (int idx = 0; idx < globalAxisIndex.numElements(); ++idx)
358  {
359    globalAxisIndex(idx) = axisDst->index(idx);
360  }
361  dhtIndexProcRank.computeIndexInfoMapping(globalAxisIndex);
362
363  std::vector<int> countIndex(clientSize,0);
364  const CClientClientDHTInt::Index2VectorInfoTypeMap& computedGlobalIndexOnProc = dhtIndexProcRank.getInfoIndexMap();
365  CClientClientDHTInt::Index2VectorInfoTypeMap::const_iterator itb = computedGlobalIndexOnProc.begin(), it,
366                                                               ite = computedGlobalIndexOnProc.end();
367  for (it = itb; it != ite; ++it)
368  {
369    const std::vector<int>& procList = it->second;
370    for (int idx = 0; idx < procList.size(); ++idx) ++countIndex[procList[idx]];
371  }
372
373  globalAxisIndexOnProc.rehash(std::ceil(clientSize/globalAxisIndexOnProc.max_load_factor()));
374  for (int idx = 0; idx < clientSize; ++idx)
375  {
376    if (0 != countIndex[idx])
377    {
378      globalAxisIndexOnProc[idx].resize(countIndex[idx]);
379      countIndex[idx] = 0;
380    }
381  }
382
383  for (it = itb; it != ite; ++it)
384  {
385    const std::vector<int>& procList = it->second;
386    for (int idx = 0; idx < procList.size(); ++idx)
387    {
388      globalAxisIndexOnProc[procList[idx]][countIndex[procList[idx]]] = it->first;
389      ++countIndex[procList[idx]];
390    }
391  }
392}
393
394/*!
395  Find out proc and global index of domain source which domain destination is on demande
396  \param[in] domainDst Domain destination
397  \param[in] domainSrc Domain source
398  \param[in] destGlobalIndexPositionInGrid Relative position of domain corresponds to other element of grid.
399  \param[out] globalDomainIndexOnProc Global index of domain source on different procs
400*/
401void CGenericAlgorithmTransformation::computeExchangeDomainIndex(CDomain* domainDst,
402                                                                 CDomain* domainSrc,
403                                                                 CArray<size_t,1>& destGlobalIndexPositionInGrid,
404                                                                 boost::unordered_map<int,std::vector<size_t> >& globalDomainIndexOnProc)
405{
406  CContext* context = CContext::getCurrent();
407  CContextClient* client=context->client;
408  int clientRank = client->clientRank;
409  int clientSize = client->clientSize;
410
411  int niGlobSrc = domainSrc->ni_glo.getValue();
412  size_t globalIndex;
413  int i_ind, j_ind;
414  int nIndexSize = (destGlobalIndexPositionInGrid.isEmpty()) ? domainSrc->i_index.numElements()
415                                                             : destGlobalIndexPositionInGrid.numElements();
416  CClientClientDHTInt::Index2VectorInfoTypeMap globalIndex2ProcRank;
417  globalIndex2ProcRank.rehash(std::ceil(nIndexSize/globalIndex2ProcRank.max_load_factor()));
418  if (destGlobalIndexPositionInGrid.isEmpty())
419  {
420    for (int idx = 0; idx < nIndexSize; ++idx)
421    {
422      i_ind=domainSrc->i_index(idx) ;
423      j_ind=domainSrc->j_index(idx) ;
424
425      globalIndex = i_ind + j_ind * niGlobSrc;
426      globalIndex2ProcRank[globalIndex].resize(1);
427      globalIndex2ProcRank[globalIndex][0] = clientRank;
428    }
429  }
430  else
431  {
432    for (int idx = 0; idx < nIndexSize; ++idx)
433    {
434      globalIndex2ProcRank[destGlobalIndexPositionInGrid(idx)].push_back(clientRank);
435    }
436  }
437
438  CArray<size_t,1> globalDomainIndex;
439  if (destGlobalIndexPositionInGrid.isEmpty())
440  {
441    int niGlobDst = domainDst->ni_glo.getValue();
442    globalDomainIndex.resize(domainDst->i_index.numElements());
443    nIndexSize = domainDst->i_index.numElements();
444
445    for (int idx = 0; idx < nIndexSize; ++idx)
446    {
447      i_ind=domainDst->i_index(idx) ;
448      j_ind=domainDst->j_index(idx) ;
449
450      globalDomainIndex(idx) = i_ind + j_ind * niGlobDst;
451    }
452  }
453  else
454  {
455    globalDomainIndex.reference(destGlobalIndexPositionInGrid);
456  }
457
458  CClientClientDHTInt dhtIndexProcRank(globalIndex2ProcRank, client->intraComm);
459  dhtIndexProcRank.computeIndexInfoMapping(globalDomainIndex);
460
461  std::vector<int> countIndex(clientSize,0);
462  const CClientClientDHTInt::Index2VectorInfoTypeMap& computedGlobalIndexOnProc = dhtIndexProcRank.getInfoIndexMap();
463  CClientClientDHTInt::Index2VectorInfoTypeMap::const_iterator itb = computedGlobalIndexOnProc.begin(), it,
464                                                               ite = computedGlobalIndexOnProc.end();
465  for (it = itb; it != ite; ++it)
466  {
467    const std::vector<int>& procList = it->second;
468    for (int idx = 0; idx < procList.size(); ++idx) ++countIndex[procList[idx]];
469  }
470
471  globalDomainIndexOnProc.rehash(std::ceil(clientSize/globalDomainIndexOnProc.max_load_factor()));
472  for (int idx = 0; idx < clientSize; ++idx)
473  {
474    if (0 != countIndex[idx])
475    {
476      globalDomainIndexOnProc[idx].resize(countIndex[idx]);
477      countIndex[idx] = 0;
478    }
479  }
480
481  for (it = itb; it != ite; ++it)
482  {
483    const std::vector<int>& procList = it->second;
484    for (int idx = 0; idx < procList.size(); ++idx)
485    {
486      globalDomainIndexOnProc[procList[idx]][countIndex[procList[idx]]] = it->first;
487      ++countIndex[procList[idx]];
488    }
489  }
490}
491
492/*!
493  Compute index mapping between element source and element destination with an auxiliary inputs which determine
494position of each mapped index in global index of grid destination.
495  \param [in] dataAuxInputs auxiliary inputs
496*/
497void CGenericAlgorithmTransformation::computeIndexSourceMapping(const std::vector<CArray<double,1>* >& dataAuxInputs)
498{
499  computeIndexSourceMapping_(dataAuxInputs);
500}
501
502std::vector<StdString> CGenericAlgorithmTransformation::getIdAuxInputs()
503{
504  return idAuxInputs_;
505}
506
507}
Note: See TracBrowser for help on using the repository browser.