source: XIOS/trunk/src/transformation/domain_algorithm_interpolate.cpp @ 978

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

Correcting various bugs relating to transformation

+) Fix the order of transformation selection
+) Correct domain transformation selection
+) Reorganize test_remap

Test
+) On Curie
+) All tests pass

File size: 24.9 KB
Line 
1/*!
2   \file domain_algorithm_interpolate_from_file.cpp
3   \author Ha NGUYEN
4   \since 09 Jul 2015
5   \date 15 Sep 2015
6
7   \brief Algorithm for interpolation on a domain.
8 */
9#include "domain_algorithm_interpolate.hpp"
10#include <boost/unordered_map.hpp>
11#include "context.hpp"
12#include "context_client.hpp"
13#include "distribution_client.hpp"
14#include "client_server_mapping_distributed.hpp"
15#include "netcdf.hpp"
16#include "mapper.hpp"
17#include "mpi_tag.hpp"
18#include "domain.hpp"
19#include "grid_transformation_factory_impl.hpp"
20#include "interpolate_domain.hpp"
21#include "grid.hpp"
22
23namespace xios {
24CGenericAlgorithmTransformation* CDomainAlgorithmInterpolate::create(CGrid* gridDst, CGrid* gridSrc,
25                                                                     CTransformation<CDomain>* transformation,
26                                                                     int elementPositionInGrid,
27                                                                     std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
28                                                                     std::map<int, int>& elementPositionInGridSrc2AxisPosition,
29                                                                     std::map<int, int>& elementPositionInGridSrc2DomainPosition,
30                                                                     std::map<int, int>& elementPositionInGridDst2ScalarPosition,
31                                                                     std::map<int, int>& elementPositionInGridDst2AxisPosition,
32                                                                     std::map<int, int>& elementPositionInGridDst2DomainPosition)
33{
34  std::vector<CDomain*> domainListDestP = gridDst->getDomains();
35  std::vector<CDomain*> domainListSrcP  = gridSrc->getDomains();
36
37  CInterpolateDomain* interpolateDomain = dynamic_cast<CInterpolateDomain*> (transformation);
38  int domainDstIndex = elementPositionInGridDst2DomainPosition[elementPositionInGrid];
39  int domainSrcIndex = elementPositionInGridSrc2DomainPosition[elementPositionInGrid];
40
41  return (new CDomainAlgorithmInterpolate(domainListDestP[domainDstIndex], domainListSrcP[domainSrcIndex], interpolateDomain));
42}
43
44bool CDomainAlgorithmInterpolate::registerTrans()
45{
46  CGridTransformationFactory<CDomain>::registerTransformation(TRANS_INTERPOLATE_DOMAIN, create);
47}
48
49CDomainAlgorithmInterpolate::CDomainAlgorithmInterpolate(CDomain* domainDestination, CDomain* domainSource, CInterpolateDomain* interpDomain)
50: CDomainAlgorithmTransformation(domainDestination, domainSource), interpDomain_(interpDomain)
51{
52  interpDomain_->checkValid(domainSource);
53}
54
55/*!
56  Compute remap with integrated remap calculation module
57*/
58void CDomainAlgorithmInterpolate::computeRemap()
59{
60  using namespace sphereRemap;
61
62  CContext* context = CContext::getCurrent();
63  CContextClient* client=context->client;
64  int clientRank = client->clientRank;
65  int i, j, k, idx;
66  std::vector<double> srcPole(3,0), dstPole(3,0);
67  int orderInterp = interpDomain_->order.getValue();
68  bool renormalize ;
69
70  if (interpDomain_->renormalize.isEmpty()) renormalize=true;
71  else renormalize = interpDomain_->renormalize;
72
73  const double poleValue = 90.0;
74  const int constNVertex = 4; // Value by default number of vertex for rectangular domain
75  int nVertexSrc, nVertexDest;
76  nVertexSrc = nVertexDest = constNVertex;
77
78  // First of all, try to retrieve the boundary values of domain source and domain destination
79  int localDomainSrcSize = domainSrc_->i_index.numElements();
80  int niSrc = domainSrc_->ni.getValue(), njSrc = domainSrc_->nj.getValue();
81  bool hasBoundSrc = domainSrc_->hasBounds;
82  if (hasBoundSrc) nVertexSrc = domainSrc_->nvertex.getValue();
83  CArray<double,2> boundsLonSrc(nVertexSrc,localDomainSrcSize);
84  CArray<double,2> boundsLatSrc(nVertexSrc,localDomainSrcSize);
85
86  if (domainSrc_->hasPole) srcPole[2] = 1;
87  if (hasBoundSrc)  // Suppose that domain source is curvilinear or unstructured
88  {
89    if (!domainSrc_->bounds_lon_2d.isEmpty())
90    {
91      for (j = 0; j < njSrc; ++j)
92        for (i = 0; i < niSrc; ++i)
93        {
94          k=j*niSrc+i;
95          for(int n=0;n<nVertexSrc;++n)
96          {
97            boundsLonSrc(n,k) = domainSrc_->bounds_lon_2d(n,i,j);
98            boundsLatSrc(n,k) = domainSrc_->bounds_lat_2d(n,i,j);
99          }
100        }
101    }
102    else
103    {
104      boundsLonSrc = domainSrc_->bounds_lon_1d;
105      boundsLatSrc = domainSrc_->bounds_lat_1d;
106    }
107  }
108  else // if domain source is rectilinear, not do anything now
109  {
110    CArray<double,1> lon_g ;
111    CArray<double,1> lat_g ;
112
113    if (!domainSrc_->lonvalue_1d.isEmpty() && !domainSrc_->latvalue_1d.isEmpty())
114    {
115      domainSrc_->AllgatherRectilinearLonLat(domainSrc_->lonvalue_1d,domainSrc_->latvalue_1d, lon_g,lat_g) ;
116    }
117    else if (! domainSrc_->latvalue_rectilinear_read_from_file.isEmpty() && ! domainSrc_->lonvalue_rectilinear_read_from_file.isEmpty() )
118    {
119      lat_g=domainSrc_->latvalue_rectilinear_read_from_file ;
120      lon_g=domainSrc_->lonvalue_rectilinear_read_from_file ;
121    }
122    else if (!domainSrc_->lon_start.isEmpty() && !domainSrc_->lon_end.isEmpty() &&
123             !domainSrc_->lat_start.isEmpty() && !domainSrc_->lat_end.isEmpty())
124    {
125      double step=(domainSrc_->lon_end-domainSrc_->lon_start)/domainSrc_->ni_glo ;
126      for (int i=0; i<domainSrc_->ni_glo; ++i) lon_g(i)=domainSrc_->lon_start+i*step ;
127      step=(domainSrc_->lat_end-domainSrc_->lat_start)/domainSrc_->nj_glo ;
128      for (int i=0; i<domainSrc_->ni_glo; ++i) lat_g(i)=domainSrc_->lat_start+i*step ;
129    }
130    else ERROR("void CDomainAlgorithmInterpolate::computeRemap()",<<"Cannot compute bounds for rectilinear domain") ;
131
132    nVertexSrc = constNVertex;
133    domainSrc_->fillInRectilinearBoundLonLat(lon_g,lat_g, boundsLonSrc, boundsLatSrc);
134  }
135
136  std::map<int,std::vector<std::pair<int,double> > > interpMapValueNorthPole;
137  std::map<int,std::vector<std::pair<int,double> > > interpMapValueSouthPole;
138
139  int localDomainDestSize = domainDest_->i_index.numElements();
140  int niDest = domainDest_->ni.getValue(), njDest = domainDest_->nj.getValue();
141  bool hasBoundDest = domainDest_->hasBounds;
142  if (hasBoundDest) nVertexDest = domainDest_->nvertex.getValue();
143  CArray<double,2> boundsLonDest(nVertexDest,localDomainDestSize);
144  CArray<double,2> boundsLatDest(nVertexDest,localDomainDestSize);
145
146  if (domainDest_->hasPole) dstPole[2] = 1;
147  if (hasBoundDest)
148  {
149    if (!domainDest_->bounds_lon_2d.isEmpty())
150    {
151      for (j = 0; j < njDest; ++j)
152        for (i = 0; i < niDest; ++i)
153        {
154          k=j*niDest+i;
155          for(int n=0;n<nVertexDest;++n)
156          {
157            boundsLonDest(n,k) = domainDest_->bounds_lon_2d(n,i,j);
158            boundsLatDest(n,k) = domainDest_->bounds_lat_2d(n,i,j);
159          }
160        }
161    }
162    else
163    {
164      boundsLonDest = domainDest_->bounds_lon_1d;
165      boundsLatDest = domainDest_->bounds_lat_1d;
166    }
167  }
168  else
169  {
170    bool isNorthPole = false;
171    bool isSouthPole = false;
172
173    CArray<double,1> lon_g ;
174    CArray<double,1> lat_g ;
175
176    if (!domainDest_->lonvalue_1d.isEmpty() && !domainDest_->latvalue_1d.isEmpty())
177    {
178      domainDest_->AllgatherRectilinearLonLat(domainDest_->lonvalue_1d,domainDest_->latvalue_1d, lon_g,lat_g) ;
179    }
180    else if (! domainDest_->latvalue_rectilinear_read_from_file.isEmpty() && ! domainDest_->lonvalue_rectilinear_read_from_file.isEmpty() )
181    {
182      lat_g=domainDest_->latvalue_rectilinear_read_from_file ;
183      lon_g=domainDest_->lonvalue_rectilinear_read_from_file ;
184    }
185    else if (!domainDest_->lon_start.isEmpty() && !domainDest_->lon_end.isEmpty() &&
186             !domainDest_->lat_start.isEmpty() && !domainDest_->lat_end.isEmpty())
187    {
188      double step=(domainDest_->lon_end-domainDest_->lon_start)/domainDest_->ni_glo ;
189      for(int i=0; i<domainDest_->ni_glo; ++i) lon_g(i)=domainDest_->lon_start+i*step ;
190      step=(domainDest_->lat_end-domainDest_->lat_start)/domainDest_->nj_glo ;
191      for(int i=0; i<domainDest_->ni_glo; ++i) lat_g(i)=domainDest_->lat_start+i*step ;
192    }
193    else ERROR("void CDomainAlgorithmInterpolate::computeRemap()",<<"Cannot compute bounds for rectilinear domain") ;
194   
195    if (std::abs(poleValue - std::abs(lat_g(0))) < NumTraits<double>::epsilon()) isNorthPole = true;
196    if (std::abs(poleValue - std::abs(lat_g(domainDest_->nj_glo-1))) < NumTraits<double>::epsilon()) isSouthPole = true;
197
198
199
200
201    if (isNorthPole && (0 == domainDest_->jbegin.getValue()))
202    {
203      int ibegin = domainDest_->ibegin.getValue();
204      for (i = 0; i < niDest; ++i)
205      {
206        interpMapValueNorthPole[i+ibegin];
207      }
208    }
209
210    if (isSouthPole && (domainDest_->nj_glo.getValue() == (domainDest_->jbegin.getValue() + njDest)))
211    {
212      int ibegin = domainDest_->ibegin.getValue();
213      int njGlo = domainDest_->nj_glo.getValue();
214      int niGlo = domainDest_->ni_glo.getValue();
215      for (i = 0; i < niDest; ++i)
216      {
217        k = (njGlo - 1)*niGlo + i + ibegin;
218        interpMapValueSouthPole[k];
219      }
220    }
221
222    // Ok, fill in boundary values for rectangular domain
223    nVertexDest = constNVertex;
224    domainDest_->fillInRectilinearBoundLonLat(lon_g,lat_g, boundsLonDest, boundsLatDest);
225  }
226
227
228
229  // Ok, now use mapper to calculate
230  int nSrcLocal = domainSrc_->i_index.numElements();
231  int nDstLocal = domainDest_->i_index.numElements();
232  long int * globalSrc = new long int [nSrcLocal];
233  long int * globalDst = new long int [nDstLocal];
234
235  long int globalIndex;
236  int i_ind, j_ind;
237  for (int idx = 0; idx < nSrcLocal; ++idx)
238  {
239    i_ind=domainSrc_->i_index(idx) ;
240    j_ind=domainSrc_->j_index(idx) ;
241
242    globalIndex = i_ind + j_ind * domainSrc_->ni_glo;
243    globalSrc[idx] = globalIndex;
244  }
245
246  for (int idx = 0; idx < nDstLocal; ++idx)
247  {
248    i_ind=domainDest_->i_index(idx) ;
249    j_ind=domainDest_->j_index(idx) ;
250
251    globalIndex = i_ind + j_ind * domainDest_->ni_glo;
252    globalDst[idx] = globalIndex;
253  }
254
255
256  // Calculate weight index
257  Mapper mapper(client->intraComm);
258  mapper.setVerbosity(PROGRESS) ;
259
260     
261  // supress masked data for the source
262  int nSrcLocalUnmasked = 0 ;
263  for (int idx=0 ; idx < nSrcLocal; idx++) if (domainSrc_->localMask(idx)) ++nSrcLocalUnmasked ;
264
265
266  CArray<double,2> boundsLonSrcUnmasked(nVertexSrc,nSrcLocalUnmasked);
267  CArray<double,2> boundsLatSrcUnmasked(nVertexSrc,nSrcLocalUnmasked);
268  long int * globalSrcUnmasked = new long int [nSrcLocalUnmasked];
269
270  nSrcLocalUnmasked=0 ;
271  for (int idx=0 ; idx < nSrcLocal; idx++)
272  {
273    if (domainSrc_->localMask(idx))
274    {
275      for(int n=0;n<nVertexSrc;++n)
276      {
277        boundsLonSrcUnmasked(n,nSrcLocalUnmasked) = boundsLonSrc(n,idx) ;
278        boundsLatSrcUnmasked(n,nSrcLocalUnmasked) = boundsLatSrc(n,idx) ;
279      }
280      globalSrcUnmasked[nSrcLocalUnmasked]=globalSrc[idx] ;
281      ++nSrcLocalUnmasked ;
282    }
283  }
284
285
286  int nDstLocalUnmasked = 0 ;
287  for (int idx=0 ; idx < nDstLocal; idx++) if (domainDest_->localMask(idx)) ++nDstLocalUnmasked ;
288
289  CArray<double,2> boundsLonDestUnmasked(nVertexDest,nDstLocalUnmasked);
290  CArray<double,2> boundsLatDestUnmasked(nVertexDest,nDstLocalUnmasked);
291  long int * globalDstUnmasked = new long int [nDstLocalUnmasked];
292
293  nDstLocalUnmasked=0 ;
294  for (int idx=0 ; idx < nDstLocal; idx++)
295  {
296    if (domainDest_->localMask(idx))
297    {
298      for(int n=0;n<nVertexDest;++n)
299      {
300        boundsLonDestUnmasked(n,nDstLocalUnmasked) = boundsLonDest(n,idx) ;
301        boundsLatDestUnmasked(n,nDstLocalUnmasked) = boundsLatDest(n,idx) ;
302      }
303      globalDstUnmasked[nDstLocalUnmasked]=globalDst[idx] ;
304      ++nDstLocalUnmasked ;
305    }
306  }
307
308  mapper.setSourceMesh(boundsLonSrcUnmasked.dataFirst(), boundsLatSrcUnmasked.dataFirst(), nVertexSrc, nSrcLocalUnmasked, &srcPole[0], globalSrcUnmasked);
309  mapper.setTargetMesh(boundsLonDestUnmasked.dataFirst(), boundsLatDestUnmasked.dataFirst(), nVertexDest, nDstLocalUnmasked, &dstPole[0], globalDstUnmasked);
310
311  std::vector<double> timings = mapper.computeWeights(orderInterp,renormalize);
312
313  std::map<int,std::vector<std::pair<int,double> > > interpMapValue;
314  std::map<int,std::vector<std::pair<int,double> > >::const_iterator iteNorthPole = interpMapValueNorthPole.end(),
315                                                                     iteSouthPole = interpMapValueSouthPole.end();
316  for (int idx = 0;  idx < mapper.nWeights; ++idx)
317  {
318    interpMapValue[mapper.targetWeightId[idx]].push_back(make_pair(mapper.sourceWeightId[idx],mapper.remapMatrix[idx]));
319    if (iteNorthPole != interpMapValueNorthPole.find(mapper.targetWeightId[idx]))
320    {
321      interpMapValueNorthPole[mapper.targetWeightId[idx]].push_back(make_pair(mapper.sourceWeightId[idx],mapper.remapMatrix[idx]));
322    }
323
324    if (iteSouthPole != interpMapValueSouthPole.find(mapper.targetWeightId[idx]))
325    {
326      interpMapValueSouthPole[mapper.targetWeightId[idx]].push_back(make_pair(mapper.sourceWeightId[idx],mapper.remapMatrix[idx]));
327    }
328  }
329  int niGloDst = domainDest_->ni_glo.getValue();
330  processPole(interpMapValueNorthPole, niGloDst);
331  processPole(interpMapValueSouthPole, niGloDst);
332
333  if (!interpMapValueNorthPole.empty())
334  {
335     std::map<int,std::vector<std::pair<int,double> > >::iterator itNorthPole = interpMapValueNorthPole.begin();
336     for (; itNorthPole != iteNorthPole; ++itNorthPole)
337     {
338       if (!(itNorthPole->second.empty()))
339        itNorthPole->second.swap(interpMapValue[itNorthPole->first]);
340     }
341  }
342
343  if (!interpMapValueSouthPole.empty())
344  {
345     std::map<int,std::vector<std::pair<int,double> > >::iterator itSouthPole = interpMapValueSouthPole.begin();
346     for (; itSouthPole != iteSouthPole; ++itSouthPole)
347     {
348       if (!(itSouthPole->second.empty()))
349        itSouthPole->second.swap(interpMapValue[itSouthPole->first]);
350     }
351  }
352
353  exchangeRemapInfo(interpMapValue);
354
355  delete [] globalSrc;
356  delete [] globalSrcUnmasked;
357  delete [] globalDst;
358  delete [] globalDstUnmasked;
359
360}
361
362void CDomainAlgorithmInterpolate::processPole(std::map<int,std::vector<std::pair<int,double> > >& interMapValuePole,
363                                              int nbGlobalPointOnPole)
364{
365  CContext* context = CContext::getCurrent();
366  CContextClient* client=context->client;
367
368  MPI_Comm poleComme(MPI_COMM_NULL);
369  MPI_Comm_split(client->intraComm, interMapValuePole.empty() ? MPI_UNDEFINED : 1, 0, &poleComme);
370  if (MPI_COMM_NULL != poleComme)
371  {
372    int nbClientPole;
373    MPI_Comm_size(poleComme, &nbClientPole);
374
375    std::map<int,std::vector<std::pair<int,double> > >::iterator itePole = interMapValuePole.end(), itPole,
376                                                                 itbPole = interMapValuePole.begin();
377
378    int nbWeight = 0;
379    for (itPole = itbPole; itPole != itePole; ++itPole)
380       nbWeight += itPole->second.size();
381
382    std::vector<int> recvCount(nbClientPole,0);
383    std::vector<int> displ(nbClientPole,0);
384    MPI_Allgather(&nbWeight,1,MPI_INT,&recvCount[0],1,MPI_INT,poleComme) ;
385
386    displ[0]=0;
387    for(int n=1;n<nbClientPole;++n) displ[n]=displ[n-1]+recvCount[n-1] ;
388    int recvSize=displ[nbClientPole-1]+recvCount[nbClientPole-1] ;
389
390    std::vector<int> sendSourceIndexBuff(nbWeight);
391    std::vector<double> sendSourceWeightBuff(nbWeight);
392    int k = 0;
393    for (itPole = itbPole; itPole != itePole; ++itPole)
394    {
395      for (int idx = 0; idx < itPole->second.size(); ++idx)
396      {
397        sendSourceIndexBuff[k] = (itPole->second)[idx].first;
398        sendSourceWeightBuff[k] = (itPole->second)[idx].second;
399        ++k;
400      }
401    }
402
403    std::vector<int> recvSourceIndexBuff(recvSize);
404    std::vector<double> recvSourceWeightBuff(recvSize);
405
406    // Gather all index and weight for pole
407    MPI_Allgatherv(&sendSourceIndexBuff[0],nbWeight,MPI_INT,&recvSourceIndexBuff[0],&recvCount[0],&displ[0],MPI_INT,poleComme);
408    MPI_Allgatherv(&sendSourceWeightBuff[0],nbWeight,MPI_DOUBLE,&recvSourceWeightBuff[0],&recvCount[0],&displ[0],MPI_DOUBLE,poleComme);
409
410    std::map<int,double> recvTemp;
411    for (int idx = 0; idx < recvSize; ++idx)
412    {
413      if (recvTemp.end() != recvTemp.find(recvSourceIndexBuff[idx]))
414        recvTemp[recvSourceIndexBuff[idx]] += recvSourceWeightBuff[idx]/nbGlobalPointOnPole;
415      else
416        recvTemp[recvSourceIndexBuff[idx]] = 0.0;
417    }
418
419    std::map<int,double>::const_iterator itRecvTemp, itbRecvTemp = recvTemp.begin(), iteRecvTemp = recvTemp.end();
420
421    for (itPole = itbPole; itPole != itePole; ++itPole)
422    {
423      itPole->second.clear();
424      for (itRecvTemp = itbRecvTemp; itRecvTemp != iteRecvTemp; ++itRecvTemp)
425          itPole->second.push_back(make_pair(itRecvTemp->first, itRecvTemp->second));
426    }
427  }
428
429}
430
431/*!
432  Compute the index mapping between domain on grid source and one on grid destination
433*/
434void CDomainAlgorithmInterpolate::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
435{
436  if (!interpDomain_->file.isEmpty())
437    readRemapInfo();
438  else
439    computeRemap();
440}
441
442void CDomainAlgorithmInterpolate::readRemapInfo()
443{
444  CContext* context = CContext::getCurrent();
445  CContextClient* client=context->client;
446  int clientRank = client->clientRank;
447
448  std::string filename = interpDomain_->file.getValue();
449  std::map<int,std::vector<std::pair<int,double> > > interpMapValue;
450  readInterpolationInfo(filename, interpMapValue);
451
452  exchangeRemapInfo(interpMapValue);
453}
454
455
456/*!
457  Read remap information from file then distribute it among clients
458*/
459void CDomainAlgorithmInterpolate::exchangeRemapInfo(std::map<int,std::vector<std::pair<int,double> > >& interpMapValue)
460{
461  CContext* context = CContext::getCurrent();
462  CContextClient* client=context->client;
463  int clientRank = client->clientRank;
464
465  this->transformationMapping_.resize(1);
466  this->transformationWeight_.resize(1);
467
468  TransformationIndexMap& transMap = this->transformationMapping_[0];
469  TransformationWeightMap& transWeight = this->transformationWeight_[0];
470
471  boost::unordered_map<size_t,int> globalIndexOfDomainDest;
472  int ni = domainDest_->ni.getValue();
473  int nj = domainDest_->nj.getValue();
474  int ni_glo = domainDest_->ni_glo.getValue();
475  size_t globalIndex;
476  int nIndexSize = domainDest_->i_index.numElements(), i_ind, j_ind;
477  for (int idx = 0; idx < nIndexSize; ++idx)
478  {
479    i_ind=domainDest_->i_index(idx) ;
480    j_ind=domainDest_->j_index(idx) ;
481
482    globalIndex = i_ind + j_ind * ni_glo;
483    globalIndexOfDomainDest[globalIndex] = clientRank;
484  }
485
486  CClientServerMappingDistributed domainIndexClientClientMapping(globalIndexOfDomainDest,
487                                                                 client->intraComm,
488                                                                 true);
489  CArray<size_t,1> globalIndexInterp(interpMapValue.size());
490  std::map<int,std::vector<std::pair<int,double> > >::const_iterator itb = interpMapValue.begin(), it,
491                                                                     ite = interpMapValue.end();
492  size_t globalIndexCount = 0;
493  for (it = itb; it != ite; ++it)
494  {
495    globalIndexInterp(globalIndexCount) = it->first;
496    ++globalIndexCount;
497  }
498
499  domainIndexClientClientMapping.computeServerIndexMapping(globalIndexInterp);
500  const CClientServerMapping::GlobalIndexMap& globalIndexInterpSendToClient = domainIndexClientClientMapping.getGlobalIndexOnServer();
501
502  //Inform each client number of index they will receive
503  int nbClient = client->clientSize;
504  int* sendBuff = new int[nbClient];
505  int* recvBuff = new int[nbClient];
506  for (int i = 0; i < nbClient; ++i)
507  {
508    sendBuff[i] = 0;
509    recvBuff[i] = 0;
510  }
511  int sendBuffSize = 0;
512  CClientServerMapping::GlobalIndexMap::const_iterator itbMap = globalIndexInterpSendToClient.begin(), itMap,
513                                                       iteMap = globalIndexInterpSendToClient.end();
514  for (itMap = itbMap; itMap != iteMap; ++itMap)
515  {
516    const std::vector<size_t>& tmp = itMap->second;
517    int sizeIndex = 0, mapSize = (itMap->second).size();
518    for (int idx = 0; idx < mapSize; ++idx)
519    {
520//      sizeIndex += interpMapValue.at((itMap->second)[idx]).size();
521      sizeIndex += (interpMapValue[(int)(itMap->second)[idx]]).size();
522    }
523    sendBuff[itMap->first] = sizeIndex;
524    sendBuffSize += sizeIndex;
525  }
526
527
528  MPI_Allreduce(sendBuff, recvBuff, nbClient, MPI_INT, MPI_SUM, client->intraComm);
529
530  int* sendIndexDestBuff = new int [sendBuffSize];
531  int* sendIndexSrcBuff  = new int [sendBuffSize];
532  double* sendWeightBuff = new double [sendBuffSize];
533
534  std::vector<MPI_Request> sendRequest;
535
536  int sendOffSet = 0, l = 0;
537  for (itMap = itbMap; itMap != iteMap; ++itMap)
538  {
539    const std::vector<size_t>& indexToSend = itMap->second;
540    int mapSize = indexToSend.size();
541    int k = 0;
542    for (int idx = 0; idx < mapSize; ++idx)
543    {
544      std::vector<std::pair<int,double> >& interpMap = interpMapValue[(int)indexToSend[idx]]; //interpMapValue.at(indexToSend[idx]);
545      for (int i = 0; i < interpMap.size(); ++i)
546      {
547        sendIndexDestBuff[l] = indexToSend[idx];
548        sendIndexSrcBuff[l]  = interpMap[i].first;
549        sendWeightBuff[l]    = interpMap[i].second;
550        ++k;
551        ++l;
552      }
553    }
554
555    sendRequest.push_back(MPI_Request());
556    MPI_Isend(sendIndexDestBuff + sendOffSet,
557             k,
558             MPI_INT,
559             itMap->first,
560             MPI_DOMAIN_INTERPOLATION_DEST_INDEX,
561             client->intraComm,
562             &sendRequest.back());
563    sendRequest.push_back(MPI_Request());
564    MPI_Isend(sendIndexSrcBuff + sendOffSet,
565             k,
566             MPI_INT,
567             itMap->first,
568             MPI_DOMAIN_INTERPOLATION_SRC_INDEX,
569             client->intraComm,
570             &sendRequest.back());
571    sendRequest.push_back(MPI_Request());
572    MPI_Isend(sendWeightBuff + sendOffSet,
573             k,
574             MPI_DOUBLE,
575             itMap->first,
576             MPI_DOMAIN_INTERPOLATION_WEIGHT,
577             client->intraComm,
578             &sendRequest.back());
579    sendOffSet += k;
580  }
581
582  int recvBuffSize = recvBuff[clientRank];
583  int* recvIndexDestBuff = new int [recvBuffSize];
584  int* recvIndexSrcBuff  = new int [recvBuffSize];
585  double* recvWeightBuff = new double [recvBuffSize];
586  int receivedSize = 0;
587  int clientSrcRank;
588  while (receivedSize < recvBuffSize)
589  {
590    MPI_Status recvStatus;
591    MPI_Recv((recvIndexDestBuff + receivedSize),
592             recvBuffSize,
593             MPI_INT,
594             MPI_ANY_SOURCE,
595             MPI_DOMAIN_INTERPOLATION_DEST_INDEX,
596             client->intraComm,
597             &recvStatus);
598
599    int countBuff = 0;
600    MPI_Get_count(&recvStatus, MPI_INT, &countBuff);
601    clientSrcRank = recvStatus.MPI_SOURCE;
602
603    MPI_Recv((recvIndexSrcBuff + receivedSize),
604             recvBuffSize,
605             MPI_INT,
606             clientSrcRank,
607             MPI_DOMAIN_INTERPOLATION_SRC_INDEX,
608             client->intraComm,
609             &recvStatus);
610
611    MPI_Recv((recvWeightBuff + receivedSize),
612             recvBuffSize,
613             MPI_DOUBLE,
614             clientSrcRank,
615             MPI_DOMAIN_INTERPOLATION_WEIGHT,
616             client->intraComm,
617             &recvStatus);
618
619    for (int idx = 0; idx < countBuff; ++idx)
620    {
621      transMap[*(recvIndexDestBuff + receivedSize + idx)].push_back(*(recvIndexSrcBuff + receivedSize + idx));
622      transWeight[*(recvIndexDestBuff + receivedSize + idx)].push_back(*(recvWeightBuff + receivedSize + idx));
623    }
624    receivedSize += countBuff;
625  }
626
627  std::vector<MPI_Status> requestStatus(sendRequest.size());
628  MPI_Waitall(sendRequest.size(), &sendRequest[0], MPI_STATUS_IGNORE);
629
630  delete [] sendIndexDestBuff;
631  delete [] sendIndexSrcBuff;
632  delete [] sendWeightBuff;
633  delete [] recvIndexDestBuff;
634  delete [] recvIndexSrcBuff;
635  delete [] recvWeightBuff;
636  delete [] sendBuff;
637  delete [] recvBuff;
638}
639
640/*!
641  Read interpolation information from a file
642  \param [in] filename interpolation file
643  \param [in/out] interpMapValue Mapping between (global) index of domain on grid destination and
644         corresponding global index of domain and associated weight value on grid source
645*/
646void CDomainAlgorithmInterpolate::readInterpolationInfo(std::string& filename,
647                                                        std::map<int,std::vector<std::pair<int,double> > >& interpMapValue)
648{
649  int ncid ;
650  int weightDimId ;
651  size_t nbWeightGlo ;
652
653  CContext* context = CContext::getCurrent();
654  CContextClient* client=context->client;
655  int clientRank = client->clientRank;
656  int clientSize = client->clientSize;
657
658  nc_open(filename.c_str(),NC_NOWRITE, &ncid) ;
659  nc_inq_dimid(ncid,"n_weight",&weightDimId) ;
660  nc_inq_dimlen(ncid,weightDimId,&nbWeightGlo) ;
661
662  size_t nbWeight ;
663  size_t start ;
664  size_t div = nbWeightGlo/clientSize ;
665  size_t mod = nbWeightGlo%clientSize ;
666  if (clientRank < mod)
667  {
668    nbWeight=div+1 ;
669    start=clientRank*(div+1) ;
670  }
671  else
672  {
673    nbWeight=div ;
674    start= mod * (div+1) + (clientRank-mod) * div ;
675  }
676
677  double* weight=new double[nbWeight] ;
678  int weightId ;
679  nc_inq_varid (ncid, "weight", &weightId) ;
680  nc_get_vara_double(ncid, weightId, &start, &nbWeight, weight) ;
681
682  long* srcIndex=new long[nbWeight] ;
683  int srcIndexId ;
684  nc_inq_varid (ncid, "src_idx", &srcIndexId) ;
685  nc_get_vara_long(ncid, srcIndexId, &start, &nbWeight, srcIndex) ;
686
687  long* dstIndex=new long[nbWeight] ;
688  int dstIndexId ;
689  nc_inq_varid (ncid, "dst_idx", &dstIndexId) ;
690  nc_get_vara_long(ncid, dstIndexId, &start, &nbWeight, dstIndex) ;
691
692  for(size_t ind=0; ind<nbWeight;++ind)
693    interpMapValue[dstIndex[ind]-1].push_back(make_pair(srcIndex[ind]-1,weight[ind]));
694}
695
696}
Note: See TracBrowser for help on using the repository browser.