source: XIOS/dev/XIOS_DEV_CMIP6/src/transformation/generic_algorithm_transformation.cpp @ 1400

Last change on this file since 1400 was 1400, checked in by ymipsl, 4 years ago

2 bug fix for spatial transformation

  • Dynamic spatial transformation distributed :: missing information saved in the class for next time step
  • Non distributed transformation : solve local mask grid problem for intermediate grid
  • take into account some transformation that are always considered as distributed

YM

File size: 40.8 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#include "utils.hpp"
14#include "timer.hpp"
15#include "mpi.hpp"
16
17namespace xios {
18
19CGenericAlgorithmTransformation::CGenericAlgorithmTransformation()
20 : transformationMapping_(), transformationWeight_(), transformationPosition_(),
21   idAuxInputs_(), type_(ELEMENT_NO_MODIFICATION_WITH_DATA), indexElementSrc_(),
22   computedProcSrcNonTransformedElement_(false), eliminateRedondantSrc_(true), isDistributedComputed_(false)
23{
24}
25
26void CGenericAlgorithmTransformation::updateData(CArray<double,1>& dataOut)
27{
28
29}
30
31void CGenericAlgorithmTransformation::apply(const std::vector<std::pair<int,double> >& localIndex,
32                                            const double* dataInput,
33                                            CArray<double,1>& dataOut,
34                                            std::vector<bool>& flagInitial,
35                                            bool ignoreMissingValue, bool firstPass  )
36{
37  int nbLocalIndex = localIndex.size();   
38  double defaultValue = std::numeric_limits<double>::quiet_NaN();
39  if (ignoreMissingValue)
40  {
41    for (int idx = 0; idx < nbLocalIndex; ++idx)
42    {
43      if (NumTraits<double>::isnan(*(dataInput + idx)))
44      {
45        flagInitial[localIndex[idx].first] = false;
46      }
47      else
48      {
49        dataOut(localIndex[idx].first) += *(dataInput + idx) * localIndex[idx].second;
50        flagInitial[localIndex[idx].first] = true; // Reset flag to indicate not all data source are nan
51      }
52    }
53
54    // If all data source are nan then data destination must be nan
55    for (int idx = 0; idx < nbLocalIndex; ++idx)
56    {
57      if (!flagInitial[localIndex[idx].first])
58        dataOut(localIndex[idx].first) = defaultValue;
59    }
60  }
61  else
62  {
63    for (int idx = 0; idx < nbLocalIndex; ++idx)
64    {
65      dataOut(localIndex[idx].first) += *(dataInput + idx) * localIndex[idx].second;
66    }
67  }
68}
69
70void CGenericAlgorithmTransformation::computePositionElements(CGrid* dst, CGrid* src)
71{
72  int idxScalar = 0, idxAxis = 0, idxDomain = 0;
73  CArray<int,1> axisDomainOrderDst = dst->axis_domain_order;
74  for (int i = 0; i < axisDomainOrderDst.numElements(); ++i)
75  {
76    int dimElement = axisDomainOrderDst(i);
77    if (2 == dimElement)
78    {
79      elementPositionInGridDst2DomainPosition_[i] = idxDomain;
80      ++idxDomain;
81    }
82    else if (1 == dimElement)
83    {
84      elementPositionInGridDst2AxisPosition_[i] = idxAxis;
85      ++idxAxis;
86    }
87    else
88    {
89      elementPositionInGridDst2ScalarPosition_[i] = idxScalar;
90      ++idxScalar;
91    }
92  }
93
94  idxScalar = idxAxis = idxDomain = 0;
95  CArray<int,1> axisDomainOrderSrc = src->axis_domain_order;
96  for (int i = 0; i < axisDomainOrderSrc.numElements(); ++i)
97  {
98    int dimElement = axisDomainOrderSrc(i);
99    if (2 == dimElement)
100    {
101      elementPositionInGridSrc2DomainPosition_[i] = idxDomain;
102      ++idxDomain;
103    }
104    else if (1 == dimElement)
105    {
106      elementPositionInGridSrc2AxisPosition_[i] = idxAxis;
107      ++idxAxis;
108    }
109    else
110    {
111      elementPositionInGridSrc2ScalarPosition_[i] = idxScalar;
112      ++idxScalar;
113    }
114  }
115}
116
117bool CGenericAlgorithmTransformation::isDistributedTransformation(int elementPositionInGrid, CGrid* gridSrc, CGrid* gridDst)
118{
119
120  if (!isDistributedComputed_)
121  {
122    isDistributedComputed_=true ;
123    if (!eliminateRedondantSrc_) isDistributed_=true ;
124    else
125    {
126      CContext* context = CContext::getCurrent();
127      CContextClient* client = context->client;
128 
129      computePositionElements(gridSrc, gridDst);
130      std::vector<CScalar*> scalarListSrcP  = gridSrc->getScalars();
131      std::vector<CAxis*> axisListSrcP = gridSrc->getAxis();
132      std::vector<CDomain*> domainListSrcP = gridSrc->getDomains();
133      int distributed, distributed_glo ;
134 
135      CArray<int,1> axisDomainSrcOrder = gridSrc->axis_domain_order;
136      if (2 == axisDomainSrcOrder(elementPositionInGrid)) // It's domain
137      {
138        distributed=domainListSrcP[elementPositionInGridSrc2DomainPosition_[elementPositionInGrid]]->isDistributed() ;
139        MPI_Allreduce(&distributed,&distributed_glo, 1, MPI_INT, MPI_LOR, client->intraComm) ;
140   
141      }
142      else if (1 == axisDomainSrcOrder(elementPositionInGrid))//it's an axis
143      {
144        distributed=axisListSrcP[elementPositionInGridSrc2AxisPosition_[elementPositionInGrid]]->isDistributed() ;
145        MPI_Allreduce(&distributed,&distributed_glo, 1, MPI_INT, MPI_LOR, client->intraComm) ;
146      }
147      else //it's a scalar
148      {
149        distributed_glo=false ;
150      } 
151      isDistributed_=distributed_glo ;
152    }
153  }
154  return isDistributed_ ;
155} 
156/*!
157  This function computes the global indexes of grid source, which the grid destination is in demand.
158  \param[in] elementPositionInGrid position of an element in a grid .E.g: if grid is composed of domain and axis (in order),
159                then position of axis in grid is 1 and domain is positioned at 0.
160  \param[in] gridSrc Grid source
161  \param[in] gridDst Grid destination
162  \param[in\out] globaIndexWeightFromSrcToDst mapping of each global index source and weight to index destination
163*/
164void CGenericAlgorithmTransformation::computeGlobalSourceIndex(int elementPositionInGrid,
165                                                               CGrid* gridSrc,
166                                                               CGrid* gridDst,
167                                                               SourceDestinationIndexMap& globaIndexWeightFromSrcToDst)
168 {
169  CContext* context = CContext::getCurrent();
170  CContextClient* client = context->client;
171  int nbClient = client->clientSize;
172
173  typedef boost::unordered_map<int, std::vector<std::pair<int,double> > > SrcToDstMap;
174  int idx;
175
176  // compute position of elements on grids
177  computePositionElements(gridDst, gridSrc);
178  std::vector<CScalar*> scalarListDestP = gridDst->getScalars();
179  std::vector<CAxis*> axisListDestP = gridDst->getAxis();
180  std::vector<CDomain*> domainListDestP = gridDst->getDomains();
181  CArray<int,1> axisDomainDstOrder = gridDst->axis_domain_order;
182  std::vector<CScalar*> scalarListSrcP  = gridSrc->getScalars();
183  std::vector<CAxis*> axisListSrcP = gridSrc->getAxis();
184  std::vector<CDomain*> domainListSrcP = gridSrc->getDomains();
185  CArray<int,1> axisDomainSrcOrder = gridSrc->axis_domain_order; 
186 
187  bool isTransPosEmpty = transformationPosition_.empty();
188  CArray<size_t,1> transPos;
189  if (!isTransPosEmpty) transPos.resize(transformationMapping_.size());
190  std::set<size_t> allIndexSrc; // All index of source, which can be scattered among processes, need for doing transformation
191 
192  for (size_t idxTrans = 0; idxTrans < transformationMapping_.size(); ++idxTrans)
193  {
194    TransformationIndexMap::const_iterator itbTransMap = transformationMapping_[idxTrans].begin(), itTransMap,
195                                           iteTransMap = transformationMapping_[idxTrans].end();
196    TransformationWeightMap::const_iterator itbTransWeight = transformationWeight_[idxTrans].begin(), itTransWeight;
197
198    // Build mapping between global source element index and global destination element index.
199    itTransWeight = itbTransWeight;
200    for (itTransMap = itbTransMap; itTransMap != iteTransMap; ++itTransMap, ++itTransWeight)
201    {
202      const std::vector<int>& srcIndex = itTransMap->second;
203      for (idx = 0; idx < srcIndex.size(); ++idx)
204        allIndexSrc.insert(srcIndex[idx]);
205    }
206
207    if (!isTransPosEmpty)
208    {
209      TransformationPositionMap::const_iterator itPosMap = transformationPosition_[idxTrans].begin();
210      transPos(idxTrans) = itPosMap->second[0];
211    }
212  }
213
214  size_t indexSrcSize = 0;
215  CArray<size_t,1> indexSrc(allIndexSrc.size());
216  for (std::set<size_t>::iterator it = allIndexSrc.begin(); it != allIndexSrc.end(); ++it, ++indexSrcSize)
217    indexSrc(indexSrcSize) = *it;
218
219  // Flag to indicate whether we will recompute distribution of source global index  on processes
220  bool computeGlobalIndexOnProc = false; 
221  if (indexElementSrc_.size() != allIndexSrc.size())
222    computeGlobalIndexOnProc = true; 
223  else
224  {
225    for (std::set<size_t>::iterator it = allIndexSrc.begin(); it != allIndexSrc.end(); ++it)
226      if (0 == indexElementSrc_.count(*it))
227      {
228        computeGlobalIndexOnProc = true;
229        break;       
230      }
231  }
232
233  if (computeGlobalIndexOnProc)
234    indexElementSrc_.swap(allIndexSrc);
235     
236  int sendValue = (computeGlobalIndexOnProc) ? 1 : 0;
237  int recvValue = 0;
238  MPI_Allreduce(&sendValue, &recvValue, 1, MPI_INT, MPI_SUM, client->intraComm);
239  computeGlobalIndexOnProc = (0 < recvValue);
240
241//  CClientClientDHTInt::Index2VectorInfoTypeMap globalIndexOfTransformedElementOnProc;
242
243  if (computeGlobalIndexOnProc || !computedProcSrcNonTransformedElement_)
244  {   
245    {
246      CClientClientDHTInt::Index2VectorInfoTypeMap tmp ;
247      globalIndexOfTransformedElementOnProc_.swap(tmp) ;
248    }
249    // Find out global index source of transformed element on corresponding process.   
250    if (globalElementIndexOnProc_.empty())
251      globalElementIndexOnProc_.resize(axisDomainDstOrder.numElements());
252   
253    for (idx = 0; idx < axisDomainDstOrder.numElements(); ++idx)
254    {
255     
256      if (idx == elementPositionInGrid)
257        computeExchangeGlobalIndex(indexSrc, axisDomainSrcOrder(idx), globalIndexOfTransformedElementOnProc_); //globalElementIndexOnProc[idx]);
258      if (!computedProcSrcNonTransformedElement_)
259      {
260        if (2 == axisDomainDstOrder(idx)) // It's domain
261        {
262          if (idx != elementPositionInGrid)
263            computeExchangeDomainIndex(domainListDestP[elementPositionInGridDst2DomainPosition_[idx]],
264                                       domainListSrcP[elementPositionInGridSrc2DomainPosition_[idx]],
265                                       transPos,
266                                       globalElementIndexOnProc_[idx]);     
267
268        }
269        else if (1 == axisDomainDstOrder(idx))//it's an axis
270        {
271          if (idx != elementPositionInGrid)
272            computeExchangeAxisIndex(axisListDestP[elementPositionInGridDst2AxisPosition_[idx]],
273                                     axisListSrcP[elementPositionInGridSrc2AxisPosition_[idx]],
274                                     transPos,
275                                     globalElementIndexOnProc_[idx]);
276        }
277        else //it's a scalar
278        {
279          if (idx != elementPositionInGrid)
280            computeExchangeScalarIndex(scalarListDestP[elementPositionInGridDst2ScalarPosition_[idx]],
281                                       scalarListSrcP[elementPositionInGridSrc2ScalarPosition_[idx]],
282                                       transPos,
283                                       globalElementIndexOnProc_[idx]);
284
285        }
286      }
287    }
288
289    if (!isTransPosEmpty && !computedProcSrcNonTransformedElement_)
290    {
291      for (idx = 0; idx < globalElementIndexOnProc_.size(); ++idx)
292      {
293        if (idx != elementPositionInGrid)
294        {
295          boost::unordered_map<int,std::vector<size_t> >::iterator itb = globalElementIndexOnProc_[idx].begin(), it,
296                                                                   ite = globalElementIndexOnProc_[idx].end();
297          for (it = itb; it != ite; ++it) it->second.resize(1);
298        }
299      }
300    }
301
302/*     
303    if (!computedProcSrcNonTransformedElement_)
304    {
305      for (idx = 0; idx < globalElementIndexOnProc_.size(); ++idx)
306      {
307        if (idx != elementPositionInGrid)
308        {
309          boost::unordered_map<int,std::vector<size_t> >::iterator itb = globalElementIndexOnProc_[idx].begin(), it,
310                                                                   ite = globalElementIndexOnProc_[idx].end();
311          for (it = itb; it != ite; ++it) procOfNonTransformedElements_.insert(it->first);
312          if (procOfNonTransformedElements_.size() == nbClient)
313            break;
314        }
315      }
316    }
317   
318    // Processes contain the source index of transformed element
319    std::set<int> procOfTransformedElement;
320    CClientClientDHTInt::Index2VectorInfoTypeMap::iterator itIdxb = globalIndexOfTransformedElementOnProc.begin(),
321                                                           itIdxe = globalIndexOfTransformedElementOnProc.end(), itIdx;
322    for (itIdx = itIdxb; itIdx != itIdxe; ++itIdx)
323    {
324      std::vector<int>& tmp = itIdx->second;
325      for (int i = 0; i < tmp.size(); ++i)
326        procOfTransformedElement.insert(tmp[i]);
327      if (tmp.size() == nbClient)
328        break;
329    }                                                           
330   
331    std::set<int>& commonProc = (procOfTransformedElement.size() < procOfNonTransformedElements_.size()) ? procOfTransformedElement
332                              : (!procOfNonTransformedElements_.empty() ? procOfNonTransformedElements_ : procOfTransformedElement);
333   
334    std::vector<int> procContainSrcElementIdx(commonProc.size());
335    int count = 0;
336    for (std::set<int>::iterator it = commonProc.begin(); it != commonProc.end(); ++it)
337      procContainSrcElementIdx[count++] = *it;
338
339    procContainSrcElementIdx_.swap(procContainSrcElementIdx);   
340*/
341   
342    if (procElementList_.empty()) procElementList_.resize(axisDomainDstOrder.numElements()) ;
343    for (idx = 0; idx < axisDomainDstOrder.numElements(); ++idx)
344    {
345      std::set<int>& procList=procElementList_[idx] ;
346      std::set<int> commonTmp ;
347      if (idx == elementPositionInGrid)
348      {
349          set<int> tmpSet ; 
350          procList.swap(tmpSet) ;
351          CClientClientDHTInt::Index2VectorInfoTypeMap::iterator itIdxb = globalIndexOfTransformedElementOnProc_.begin(),
352                                                                 itIdxe = globalIndexOfTransformedElementOnProc_.end(), itIdx;
353          for (itIdx = itIdxb; itIdx != itIdxe; ++itIdx)
354          {
355             std::vector<int>& tmp = itIdx->second;
356             for (int i = 0; i < tmp.size(); ++i) procList.insert(tmp[i]);
357             if (tmp.size() == nbClient)
358             break;
359          }
360      }
361      else
362      {
363        if (!computedProcSrcNonTransformedElement_)
364        {
365          set<int> tmpSet ; 
366          procList.swap(tmpSet) ;
367          boost::unordered_map<int,std::vector<size_t> >::iterator itb = globalElementIndexOnProc_[idx].begin(), it,
368                                                                   ite = globalElementIndexOnProc_[idx].end();
369          for (it = itb; it != ite; ++it)
370          {
371            procList.insert(it->first);
372            if (procList.size() == nbClient)  break;
373          }
374        }
375      }
376
377      if (idx==0) commonProc_.swap(procList) ;
378      else
379      {
380        for (std::set<int>::iterator it = commonProc_.begin(); it != commonProc_.end(); ++it)
381          if (procList.count(*it)==1) commonTmp.insert(*it) ;
382        commonProc_.swap(commonTmp) ;
383      }
384    }
385    std::vector<int> procContainSrcElementIdx(commonProc_.size());
386    int count = 0;
387    for (std::set<int>::iterator it = commonProc_.begin(); it != commonProc_.end(); ++it) procContainSrcElementIdx[count++] = *it;
388    procContainSrcElementIdx_.swap(procContainSrcElementIdx);
389   
390        // For the first time, surely we calculate proc containing non transformed elements
391    if (!computedProcSrcNonTransformedElement_) computedProcSrcNonTransformedElement_ = true;
392  }
393 
394  for (size_t idxTrans = 0; idxTrans < transformationMapping_.size(); ++idxTrans)
395  {
396    TransformationIndexMap::const_iterator itbTransMap = transformationMapping_[idxTrans].begin(), itTransMap,
397                                           iteTransMap = transformationMapping_[idxTrans].end();
398    TransformationWeightMap::const_iterator itbTransWeight = transformationWeight_[idxTrans].begin(), itTransWeight;
399    SrcToDstMap src2DstMap;
400    src2DstMap.rehash(std::ceil(transformationMapping_[idxTrans].size()/src2DstMap.max_load_factor()));
401
402    // Build mapping between global source element index and global destination element index.
403    boost::unordered_map<int,std::vector<size_t> >().swap(globalElementIndexOnProc_[elementPositionInGrid]);
404    std::set<int> tmpCounter;
405    itTransWeight = itbTransWeight;
406    for (itTransMap = itbTransMap; itTransMap != iteTransMap; ++itTransMap, ++itTransWeight)
407    {
408      const std::vector<int>& srcIndex = itTransMap->second;
409      const std::vector<double>& weight = itTransWeight->second;
410      for (idx = 0; idx < srcIndex.size(); ++idx)
411      {
412        src2DstMap[srcIndex[idx]].push_back(make_pair(itTransMap->first, weight[idx]));
413        if (0 == tmpCounter.count(srcIndex[idx]))
414        {         
415          tmpCounter.insert(srcIndex[idx]);
416       
417          vector<int>& rankSrc = globalIndexOfTransformedElementOnProc_[srcIndex[idx]] ;
418          for (int n=0;n<rankSrc.size();++n)
419          {
420            if (commonProc_.count(rankSrc[n])==1) globalElementIndexOnProc_[elementPositionInGrid][rankSrc[n]].push_back(srcIndex[idx]);
421          }
422//          for (int j = 0; j < procContainSrcElementIdx_.size(); ++j)
423//            globalElementIndexOnProc_[elementPositionInGrid][procContainSrcElementIdx_[j]].push_back(srcIndex[idx]);
424        }
425      }
426    }
427 
428    if (!isTransPosEmpty)
429    {
430      for (idx = 0; idx < globalElementIndexOnProc_.size(); ++idx)
431      {
432        if (idx != elementPositionInGrid)
433        {
434          boost::unordered_map<int,std::vector<size_t> >::iterator itb = globalElementIndexOnProc_[idx].begin(), it,
435                                                                   ite = globalElementIndexOnProc_[idx].end();
436          for (it = itb; it != ite; ++it) it->second[0] = transPos(idxTrans);
437        }
438      }
439    }
440
441    // Ok, now compute global index of grid source and ones of grid destination
442    computeGlobalGridIndexMapping(elementPositionInGrid,
443                                  procContainSrcElementIdx_, //srcRank,
444                                  src2DstMap,
445                                  gridSrc,
446                                  gridDst,
447                                  globalElementIndexOnProc_,
448                                  globaIndexWeightFromSrcToDst);
449  } 
450 }
451
452/*!
453  Compute mapping of global index of grid source and grid destination
454  \param [in] elementPositionInGrid position of element in grid. E.x: grid composed of domain and axis, domain has position 0 and axis 1.
455  \param [in] srcRank rank of client from which we demand global index of element source
456  \param [in] src2DstMap mapping of global index of element source and global index of element destination
457  \param [in] gridSrc Grid source
458  \param [in] gridDst Grid destination
459  \param [in] globalElementIndexOnProc Global index of element source on different client rank
460  \param [out] globaIndexWeightFromSrcToDst Mapping of global index of grid source and grid destination
461*/
462void CGenericAlgorithmTransformation::computeGlobalGridIndexMapping(int elementPositionInGrid,
463                                                                   const std::vector<int>& srcRank,
464                                                                   boost::unordered_map<int, std::vector<std::pair<int,double> > >& src2DstMap,
465                                                                   CGrid* gridSrc,
466                                                                   CGrid* gridDst,
467                                                                   std::vector<boost::unordered_map<int,std::vector<size_t> > >& globalElementIndexOnProc,
468                                                                   SourceDestinationIndexMap& globaIndexWeightFromSrcToDst)
469{
470  SourceDestinationIndexMap globaIndexWeightFromSrcToDst_tmp ;
471 
472  CContext* context = CContext::getCurrent();
473  CContextClient* client=context->client;
474  int clientRank = client->clientRank;
475 
476  std::vector<CDomain*> domainListSrcP = gridSrc->getDomains();
477  std::vector<CAxis*> axisListSrcP = gridSrc->getAxis();
478  std::vector<CScalar*> scalarListSrcP = gridSrc->getScalars();
479  CArray<int,1> axisDomainSrcOrder = gridSrc->axis_domain_order;
480
481  size_t nbElement = axisDomainSrcOrder.numElements();
482  std::vector<size_t> nGlobSrc(nbElement);
483  size_t globalSrcSize = 1;
484  int domainIndex = 0, axisIndex = 0, scalarIndex = 0;
485  for (int idx = 0; idx < nbElement; ++idx)
486  {
487    nGlobSrc[idx] = globalSrcSize;
488    int elementDimension = axisDomainSrcOrder(idx);
489
490    // If this is a domain
491    if (2 == elementDimension)
492    {
493      globalSrcSize *= domainListSrcP[domainIndex]->nj_glo.getValue() * domainListSrcP[domainIndex]->ni_glo.getValue();
494      ++domainIndex;
495    }
496    else if (1 == elementDimension) // So it's an axis
497    {
498      globalSrcSize *= axisListSrcP[axisIndex]->n_glo.getValue();
499      ++axisIndex;
500    }
501    else
502    {
503      globalSrcSize *= 1;
504      ++scalarIndex;
505    }
506  }
507
508  std::vector<CDomain*> domainListDestP = gridDst->getDomains();
509  std::vector<CAxis*> axisListDestP = gridDst->getAxis();
510  std::vector<CScalar*> scalarListDestP = gridDst->getScalars();
511  CArray<int,1> axisDomainDstOrder = gridDst->axis_domain_order;
512
513  std::vector<size_t> nGlobDst(nbElement);
514  size_t globalDstSize = 1;
515  domainIndex = axisIndex = scalarIndex = 0;
516  set<size_t>  globalSrcStoreIndex ;
517 
518  for (int idx = 0; idx < nbElement; ++idx)
519  {
520    nGlobDst[idx] = globalDstSize;
521    int elementDimension = axisDomainDstOrder(idx);
522
523    // If this is a domain
524    if (2 == elementDimension)
525    {
526      globalDstSize *= domainListDestP[domainIndex]->nj_glo.getValue() * domainListDestP[domainIndex]->ni_glo.getValue();
527      ++domainIndex;
528    }
529    else if (1 == elementDimension) // So it's an axis
530    {
531      globalDstSize *= axisListDestP[axisIndex]->n_glo.getValue();
532      ++axisIndex;
533    }
534    else
535    {
536      globalDstSize *= 1;
537      ++scalarIndex;
538    }
539  }
540
541  std::map< std::pair<size_t,size_t>, int > rankMap ;
542  std::map< std::pair<size_t,size_t>, int >:: iterator rankMapIt ;
543 
544  for (int i = 0; i < srcRank.size(); ++i)
545  {
546    size_t ssize = 1;
547    int rankSrc = srcRank[i];
548    for (int idx = 0; idx < nbElement; ++idx)
549    {
550      ssize *= (globalElementIndexOnProc[idx][rankSrc]).size();
551    }
552
553    std::vector<int> idxLoop(nbElement,0);
554    std::vector<int> currentIndexSrc(nbElement, 0);
555    std::vector<int> currentIndexDst(nbElement, 0);
556    int innnerLoopSize = (globalElementIndexOnProc[0])[rankSrc].size();
557    size_t idx = 0;
558    while (idx < ssize)
559    {
560      for (int ind = 0; ind < nbElement; ++ind)
561      {
562        if (idxLoop[ind] == (globalElementIndexOnProc[ind])[rankSrc].size())
563        {
564          idxLoop[ind] = 0;
565          ++idxLoop[ind+1];
566        }
567
568        currentIndexDst[ind] = currentIndexSrc[ind] = (globalElementIndexOnProc[ind])[rankSrc][idxLoop[ind]];
569      }
570
571      for (int ind = 0; ind < innnerLoopSize; ++ind)
572      {
573        currentIndexDst[0] = currentIndexSrc[0] = (globalElementIndexOnProc[0])[rankSrc][ind];
574        int globalElementDstIndexSize = 0;
575        if (1 == src2DstMap.count(currentIndexSrc[elementPositionInGrid]))
576        {
577          globalElementDstIndexSize = src2DstMap[currentIndexSrc[elementPositionInGrid]].size();
578        }
579
580        std::vector<size_t> globalDstVecIndex(globalElementDstIndexSize,0);
581        size_t globalSrcIndex = 0;
582        for (int idxElement = 0; idxElement < nbElement; ++idxElement)
583        {
584          if (idxElement == elementPositionInGrid)
585          {
586            for (int k = 0; k < globalElementDstIndexSize; ++k)
587            {
588              globalDstVecIndex[k] += src2DstMap[currentIndexSrc[elementPositionInGrid]][k].first * nGlobDst[idxElement];
589            }
590          }
591          else
592          {
593            for (int k = 0; k < globalElementDstIndexSize; ++k)
594            {
595              globalDstVecIndex[k] += currentIndexDst[idxElement] * nGlobDst[idxElement];
596            }
597          }
598          globalSrcIndex += currentIndexSrc[idxElement] * nGlobSrc[idxElement];
599        }
600
601        for (int k = 0; k < globalElementDstIndexSize; ++k)
602        {
603         
604          globaIndexWeightFromSrcToDst_tmp[rankSrc][globalSrcIndex].push_back(make_pair(globalDstVecIndex[k],src2DstMap[currentIndexSrc[elementPositionInGrid]][k].second ));
605          rankMapIt=rankMap.find(make_pair(globalSrcIndex,globalDstVecIndex[k])) ;
606          if (rankMapIt==rankMap.end()) rankMap[make_pair(globalSrcIndex,globalDstVecIndex[k])] = rankSrc ;
607          else if (rankSrc==clientRank) rankMapIt->second = rankSrc ;
608        }
609        ++idxLoop[0];
610      }
611      idx += innnerLoopSize;
612    }
613  }
614
615  // eliminate redondant global src point owned by differrent processes.
616  // Avoid as possible to tranfer data from an other process if the src point is also owned by current process
617      int rankSrc ;
618      size_t globalSrcIndex ;
619      size_t globalDstIndex ;
620      double weight ;
621 
622      SourceDestinationIndexMap::iterator rankIt,rankIte ;
623      boost::unordered_map<size_t, std::vector<std::pair<size_t,double> > >::iterator globalSrcIndexIt, globalSrcIndexIte ;
624      std::vector<std::pair<size_t,double> >::iterator vectIt,vectIte ;
625   
626      rankIt=globaIndexWeightFromSrcToDst_tmp.begin() ; rankIte=globaIndexWeightFromSrcToDst_tmp.end() ;
627      for(;rankIt!=rankIte;++rankIt)
628      {
629        rankSrc = rankIt->first ;
630        globalSrcIndexIt = rankIt->second.begin() ; globalSrcIndexIte = rankIt->second.end() ;
631        for(;globalSrcIndexIt!=globalSrcIndexIte;++globalSrcIndexIt)
632        { 
633          globalSrcIndex = globalSrcIndexIt->first ;
634          vectIt = globalSrcIndexIt->second.begin() ; vectIte = globalSrcIndexIt->second.end() ;
635          for(vectIt; vectIt!=vectIte; vectIt++)
636          {
637            globalDstIndex = vectIt->first ;
638            weight = vectIt->second ;
639            if (eliminateRedondantSrc_)
640            {
641              if (rankMap[make_pair(globalSrcIndex,globalDstIndex)] == rankSrc) 
642                globaIndexWeightFromSrcToDst[rankSrc][globalSrcIndex].push_back(make_pair(globalDstIndex,weight)) ;
643            }
644            else globaIndexWeightFromSrcToDst[rankSrc][globalSrcIndex].push_back(make_pair(globalDstIndex,weight)) ;
645         }
646       }
647     }
648
649}
650
651/*!
652  Find out proc and global index of axis source which axis destination is on demande
653  \param[in] scalar Scalar destination
654  \param[in] scalar Scalar source
655  \param[in] destGlobalIndexPositionInGrid Relative position of axis corresponds to other element of grid.
656  \param[out] globalScalarIndexOnProc Global index of axis source on different procs
657*/
658void CGenericAlgorithmTransformation::computeExchangeScalarIndex(CScalar* scalarDst,
659                                                                 CScalar* scalarSrc,
660                                                                 CArray<size_t,1>& destGlobalIndexPositionInGrid,
661                                                                 boost::unordered_map<int,std::vector<size_t> >& globalScalarIndexOnProc)
662{
663  CContext* context = CContext::getCurrent();
664  CContextClient* client=context->client;
665  int clientRank = client->clientRank;
666  int clientSize = client->clientSize;
667
668  globalScalarIndexOnProc.rehash(std::ceil(clientSize/globalScalarIndexOnProc.max_load_factor()));
669  for (int idx = 0; idx < clientSize; ++idx)
670  {
671    globalScalarIndexOnProc[idx].push_back(0);
672  }
673}
674
675/*!
676  Find out proc and global index of axis source which axis destination is on demande
677  \param[in] axisDst Axis destination
678  \param[in] axisSrc Axis source
679  \param[in] destGlobalIndexPositionInGrid Relative position of axis corresponds to other element of grid.
680  \param[out] globalAxisIndexOnProc Global index of axis source on different procs
681*/
682void CGenericAlgorithmTransformation::computeExchangeAxisIndex(CAxis* axisDst,
683                                                               CAxis* axisSrc,
684                                                               CArray<size_t,1>& destGlobalIndexPositionInGrid,
685                                                               boost::unordered_map<int,std::vector<size_t> >& globalAxisIndexOnProc)
686{
687  CContext* context = CContext::getCurrent();
688  CContextClient* client=context->client;
689  int clientRank = client->clientRank;
690  int clientSize = client->clientSize;
691
692  size_t globalIndex;
693  int nIndexSize = axisSrc->index.numElements();
694  CClientClientDHTInt::Index2VectorInfoTypeMap globalIndex2ProcRank;
695  globalIndex2ProcRank.rehash(std::ceil(nIndexSize/globalIndex2ProcRank.max_load_factor()));
696  for (int idx = 0; idx < nIndexSize; ++idx)
697  {
698    globalIndex = axisSrc->index(idx);
699    globalIndex2ProcRank[globalIndex].push_back(clientRank);
700  }
701
702  CClientClientDHTInt dhtIndexProcRank(globalIndex2ProcRank, client->intraComm);
703  CArray<size_t,1> globalAxisIndex(axisDst->index.numElements());
704  for (int idx = 0; idx < globalAxisIndex.numElements(); ++idx)
705  {
706    globalAxisIndex(idx) = axisDst->index(idx);
707  }
708  dhtIndexProcRank.computeIndexInfoMapping(globalAxisIndex);
709
710  std::vector<int> countIndex(clientSize,0);
711  const CClientClientDHTInt::Index2VectorInfoTypeMap& computedGlobalIndexOnProc = dhtIndexProcRank.getInfoIndexMap();
712  CClientClientDHTInt::Index2VectorInfoTypeMap::const_iterator itb = computedGlobalIndexOnProc.begin(), it,
713                                                               ite = computedGlobalIndexOnProc.end();
714  for (it = itb; it != ite; ++it)
715  {
716    const std::vector<int>& procList = it->second;
717    for (int idx = 0; idx < procList.size(); ++idx) ++countIndex[procList[idx]];
718  }
719
720  globalAxisIndexOnProc.rehash(std::ceil(clientSize/globalAxisIndexOnProc.max_load_factor()));
721  for (int idx = 0; idx < clientSize; ++idx)
722  {
723    if (0 != countIndex[idx])
724    {
725      globalAxisIndexOnProc[idx].resize(countIndex[idx]);
726      countIndex[idx] = 0;
727    }
728  }
729
730  for (it = itb; it != ite; ++it)
731  {
732    const std::vector<int>& procList = it->second;
733    for (int idx = 0; idx < procList.size(); ++idx)
734    {
735      globalAxisIndexOnProc[procList[idx]][countIndex[procList[idx]]] = it->first;
736      ++countIndex[procList[idx]];
737    }
738  }
739}
740
741/*!
742  Find out proc and global index of domain source which domain destination is on demande
743  \param[in] domainDst Domain destination
744  \param[in] domainSrc Domain source
745  \param[in] destGlobalIndexPositionInGrid Relative position of domain corresponds to other element of grid.
746  \param[out] globalDomainIndexOnProc Global index of domain source on different procs
747*/
748void CGenericAlgorithmTransformation::computeExchangeDomainIndex(CDomain* domainDst,
749                                                                 CDomain* domainSrc,
750                                                                 CArray<size_t,1>& destGlobalIndexPositionInGrid,
751                                                                 boost::unordered_map<int,std::vector<size_t> >& globalDomainIndexOnProc)
752{
753  CContext* context = CContext::getCurrent();
754  CContextClient* client=context->client;
755  int clientRank = client->clientRank;
756  int clientSize = client->clientSize;
757
758  int niGlobSrc = domainSrc->ni_glo.getValue();
759  size_t globalIndex;
760  int i_ind, j_ind;
761  int nIndexSize = (destGlobalIndexPositionInGrid.isEmpty()) ? domainSrc->i_index.numElements()
762                                                             : destGlobalIndexPositionInGrid.numElements();
763  CClientClientDHTInt::Index2VectorInfoTypeMap globalIndex2ProcRank;
764  globalIndex2ProcRank.rehash(std::ceil(nIndexSize/globalIndex2ProcRank.max_load_factor()));
765  if (destGlobalIndexPositionInGrid.isEmpty())
766  {
767    for (int idx = 0; idx < nIndexSize; ++idx)
768    {
769      i_ind=domainSrc->i_index(idx) ;
770      j_ind=domainSrc->j_index(idx) ;
771
772      globalIndex = i_ind + j_ind * niGlobSrc;
773      globalIndex2ProcRank[globalIndex].resize(1);
774      globalIndex2ProcRank[globalIndex][0] = clientRank;
775    }
776  }
777  else
778  {
779    for (int idx = 0; idx < nIndexSize; ++idx)
780    {
781      globalIndex2ProcRank[destGlobalIndexPositionInGrid(idx)].push_back(clientRank);
782    }
783  }
784
785  CArray<size_t,1> globalDomainIndex;
786  if (destGlobalIndexPositionInGrid.isEmpty())
787  {
788    int niGlobDst = domainDst->ni_glo.getValue();
789    globalDomainIndex.resize(domainDst->i_index.numElements());
790    nIndexSize = domainDst->i_index.numElements();
791
792    for (int idx = 0; idx < nIndexSize; ++idx)
793    {
794      i_ind=domainDst->i_index(idx) ;
795      j_ind=domainDst->j_index(idx) ;
796
797      globalDomainIndex(idx) = i_ind + j_ind * niGlobDst;
798    }
799  }
800  else
801  {
802    globalDomainIndex.reference(destGlobalIndexPositionInGrid);
803  }
804
805  CClientClientDHTInt dhtIndexProcRank(globalIndex2ProcRank, client->intraComm);
806  dhtIndexProcRank.computeIndexInfoMapping(globalDomainIndex);
807
808  std::vector<int> countIndex(clientSize,0);
809  const CClientClientDHTInt::Index2VectorInfoTypeMap& computedGlobalIndexOnProc = dhtIndexProcRank.getInfoIndexMap();
810  CClientClientDHTInt::Index2VectorInfoTypeMap::const_iterator itb = computedGlobalIndexOnProc.begin(), it,
811                                                               ite = computedGlobalIndexOnProc.end();
812  for (it = itb; it != ite; ++it)
813  {
814    const std::vector<int>& procList = it->second;
815    for (int idx = 0; idx < procList.size(); ++idx) ++countIndex[procList[idx]];
816  }
817
818  globalDomainIndexOnProc.rehash(std::ceil(clientSize/globalDomainIndexOnProc.max_load_factor()));
819  for (int idx = 0; idx < clientSize; ++idx)
820  {
821    if (0 != countIndex[idx])
822    {
823      globalDomainIndexOnProc[idx].resize(countIndex[idx]);
824      countIndex[idx] = 0;
825    }
826  }
827
828  for (it = itb; it != ite; ++it)
829  {
830    const std::vector<int>& procList = it->second;
831    for (int idx = 0; idx < procList.size(); ++idx)
832    {
833      globalDomainIndexOnProc[procList[idx]][countIndex[procList[idx]]] = it->first;
834      ++countIndex[procList[idx]];
835    }
836  }
837}
838
839
840
841
842
843
844void CGenericAlgorithmTransformation::computeTransformationMappingNonDistributed(int elementPositionInGrid, CGrid* gridSrc, CGrid* gridDst,
845                                                                                 vector<int>& localSrc, vector<int>& localDst, vector<double>& weight, vector<bool>& localMaskOnGridDest)
846{
847
848  CContext* context = CContext::getCurrent();
849  CContextClient* client = context->client;
850  int nbClient = client->clientSize;
851
852  computePositionElements(gridDst, gridSrc);
853  std::vector<CScalar*> scalarListDstP = gridDst->getScalars();
854  std::vector<CAxis*> axisListDstP = gridDst->getAxis();
855  std::vector<CDomain*> domainListDstP = gridDst->getDomains();
856  CArray<int,1> axisDomainDstOrder = gridDst->axis_domain_order;
857  std::vector<CScalar*> scalarListSrcP  = gridSrc->getScalars();
858  std::vector<CAxis*> axisListSrcP = gridSrc->getAxis();
859  std::vector<CDomain*> domainListSrcP = gridSrc->getDomains();
860  CArray<int,1> axisDomainSrcOrder = gridSrc->axis_domain_order; 
861
862  int nElement=axisDomainSrcOrder.numElements() ;
863  int indSrc=1 ;
864  int indDst=1 ;
865  vector<int> nIndexSrc(nElement) ;
866  vector<int> nIndexDst(nElement) ;
867  int nlocalIndexSrc=1 ;
868  int nlocalIndexDest=1 ;
869 
870  for(int i=0 ; i<nElement; i++)
871  {
872    int dimElement = axisDomainSrcOrder(i);
873    if (2 == dimElement) //domain
874    {
875      CDomain* domain=domainListSrcP[elementPositionInGridSrc2DomainPosition_[i]] ;
876      nIndexSrc[i] = domain->i_index.numElements() ;
877    }
878    else if (1 == dimElement) //axis
879    {
880      CAxis* axis=axisListSrcP[elementPositionInGridSrc2DomainPosition_[i]] ;
881      nIndexSrc[i] = axis->index.numElements() ;
882    }
883    else  //scalar
884    {
885      nIndexSrc[i]=1 ;
886    }
887    nlocalIndexSrc=nlocalIndexSrc*nIndexSrc[i] ;
888  }
889
890  int offset=1 ;
891  for(int i=0 ; i<nElement; i++)
892  {
893    int dimElement = axisDomainDstOrder(i);
894    if (2 == dimElement) //domain
895    {
896      CDomain* domain=domainListDstP[elementPositionInGridDst2DomainPosition_[i]] ;
897      nIndexDst[i] = domain->i_index.numElements() ;
898    }
899    else if (1 == dimElement) //axis
900    {
901      CAxis* axis = axisListDstP[elementPositionInGridDst2DomainPosition_[i]] ;
902      nIndexDst[i] = axis->index.numElements() ;
903    }
904    else  //scalar
905    {
906      nIndexDst[i]=1 ;
907    }
908    if (i<elementPositionInGrid) offset=offset*nIndexDst[i] ;
909    nlocalIndexDest=nlocalIndexDest*nIndexDst[i] ;
910  }
911
912  vector<int> dstLocalInd ;
913  int dimElement = axisDomainDstOrder(elementPositionInGrid);
914  if (2 == dimElement) //domain
915  {
916    CDomain* domain = domainListSrcP[elementPositionInGridSrc2DomainPosition_[elementPositionInGrid]] ;
917    int ni_glo=domain->ni_glo ;
918    int nj_glo=domain->nj_glo ;
919    int nindex_glo=ni_glo*nj_glo ;
920    dstLocalInd.resize(nindex_glo,-1) ;
921    int nIndex=domain->i_index.numElements() ;
922    for(int i=0;i<nIndex;i++)
923    {
924      int globIndex=domain->j_index(i)*ni_glo+domain->i_index(i) ;
925      dstLocalInd[globIndex]=i ;
926    }
927  }
928  else if (1 == dimElement) //axis
929  {
930    CAxis* axis = axisListSrcP[elementPositionInGridSrc2AxisPosition_[elementPositionInGrid]] ;
931    int nindex_glo=axis->n_glo ;
932    dstLocalInd.resize(nindex_glo,-1) ;
933    int nIndex=axis->index.numElements() ;
934    for(int i=0;i<nIndex;i++)
935    {
936      dstLocalInd[axis->index(i)]=i ;
937    }
938  }
939  else  //scalar
940  {
941    dstLocalInd.resize(1) ; 
942    dstLocalInd[0]=0 ; 
943  }
944
945// just get the local src mask
946  CArray<bool,1> localMaskOnSrcGrid;
947  gridSrc->getLocalMask(localMaskOnSrcGrid) ;
948// intermediate grid, mask is not initialized => set up mask to true
949  if (localMaskOnSrcGrid.isEmpty()) localMaskOnSrcGrid.resize(nlocalIndexSrc) ;
950  localMaskOnSrcGrid=true ;
951 
952
953  localMaskOnGridDest.resize(nlocalIndexDest,false) ;
954
955  vector<vector<vector<pair<int,double> > > > dstIndWeight(transformationMapping_.size()) ;
956   
957  for(int t=0;t<transformationMapping_.size();++t)
958  {
959    TransformationIndexMap::const_iterator   itTransMap = transformationMapping_[t].begin(),
960                                             iteTransMap = transformationMapping_[t].end();
961    TransformationWeightMap::const_iterator itTransWeight = transformationWeight_[t].begin();
962    dstIndWeight[t].resize(nIndexSrc[elementPositionInGrid]) ;
963   
964    for(;itTransMap!=iteTransMap;++itTransMap,++itTransWeight)
965    {
966      int dst=dstLocalInd[itTransMap->first] ;
967      if (dst!=-1)
968      {
969        const vector<int>& srcs=itTransMap->second;
970        const vector<double>& weights=itTransWeight->second;
971        for(int i=0;i<srcs.size() ;i++) dstIndWeight[t][srcs[i]].push_back(pair<int,double>(dst*offset+t,weights[i])) ;
972      }
973    }
974  }
975  int srcInd=0 ; 
976  int currentInd ;
977  int t=0 ; 
978
979 
980  nonDistributedrecursiveFunct(nElement-1,elementPositionInGrid,srcInd, nIndexSrc, t, dstIndWeight, 
981                               currentInd,localSrc,localDst,weight, localMaskOnSrcGrid, localMaskOnGridDest );
982               
983}
984
985
986void CGenericAlgorithmTransformation::nonDistributedrecursiveFunct(int currentPos, int elementPositionInGrid, int& srcInd, vector<int>& nIndexSrc, int& t, vector<vector<vector<pair<int,double> > > >& dstIndWeight, int currentInd,
987                    vector<int>& localSrc, vector<int>& localDst, vector<double>& weight,  CArray<bool,1>& localMaskOnGridSrc, vector<bool>& localMaskOnGridDest )
988{
989  if (currentPos!=elementPositionInGrid)
990  {
991    if (currentPos!=0)
992    {
993      for(int i=0;i<nIndexSrc[currentPos];i++)
994      {
995        nonDistributedrecursiveFunct(currentPos-1, elementPositionInGrid, srcInd, nIndexSrc, t, dstIndWeight, currentInd, localSrc, localDst, weight, localMaskOnGridSrc, localMaskOnGridDest) ;
996      }
997    }
998    else
999    {
1000      for(int i=0;i<nIndexSrc[currentPos];i++)
1001      {
1002        if (dstIndWeight[t][currentInd].size()>0)
1003        {
1004          for(vector<pair<int,double> >::iterator it = dstIndWeight[t][currentInd].begin(); it!=dstIndWeight[t][currentInd].end(); ++it)
1005          {
1006            if (localMaskOnGridSrc(srcInd))
1007            {
1008              localSrc.push_back(srcInd) ;
1009              localDst.push_back(it->first) ;
1010              weight.push_back(it->second) ;
1011              localMaskOnGridDest[it->first]=true ;
1012            }
1013            (it->first)++ ;
1014          }
1015        }
1016        srcInd++ ;
1017        if (t < dstIndWeight.size()-1) t++ ;
1018      }
1019    }
1020  }
1021  else
1022  {
1023 
1024    if (currentPos!=0)
1025    {
1026
1027
1028      for(int i=0;i<nIndexSrc[currentPos];i++)
1029      {
1030        t=0 ;
1031        nonDistributedrecursiveFunct(currentPos-1, elementPositionInGrid, srcInd, nIndexSrc, t, dstIndWeight , i,  localSrc, localDst, weight, localMaskOnGridSrc, localMaskOnGridDest) ;
1032      }
1033    }
1034    else
1035    {
1036      for(int i=0;i<nIndexSrc[currentPos];i++)
1037      {
1038        t=0 ;       
1039        if (dstIndWeight[t][i].size()>0)
1040        {
1041          for(vector<pair<int,double> >::iterator it = dstIndWeight[t][i].begin(); it!=dstIndWeight[t][i].end(); ++it)
1042          {
1043            if (localMaskOnGridSrc(srcInd))
1044            {
1045              localSrc.push_back(srcInd) ;
1046              localDst.push_back(it->first) ;
1047              weight.push_back(it->second) ;
1048              localMaskOnGridDest[it->first]=true ;
1049            }
1050            (it->first)++ ;
1051          }
1052        }
1053        srcInd++ ;
1054        if (t < dstIndWeight.size()-1) t++ ;
1055      }
1056    }
1057  }
1058
1059}
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071/*!
1072  Compute index mapping between element source and element destination with an auxiliary inputs which determine
1073position of each mapped index in global index of grid destination.
1074  \param [in] dataAuxInputs auxiliary inputs
1075*/
1076void CGenericAlgorithmTransformation::computeIndexSourceMapping(const std::vector<CArray<double,1>* >& dataAuxInputs)
1077{
1078  computeIndexSourceMapping_(dataAuxInputs);
1079}
1080
1081std::vector<StdString> CGenericAlgorithmTransformation::getIdAuxInputs()
1082{
1083  return idAuxInputs_;
1084}
1085
1086CGenericAlgorithmTransformation::AlgoTransType CGenericAlgorithmTransformation::type()
1087{
1088  return type_;
1089}
1090
1091}
Note: See TracBrowser for help on using the repository browser.