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

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

Improving transformation selection. Instead of modifying directly grid_transformation
we only need to register a new transformation with the framework

+) Update all transformations with this new method

Test
+) On Curie
+) Basic tests pass

File size: 25.2 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 = elementPositionInGridDst2AxisPosition[elementPositionInGrid];
39  int domainSrcIndex = elementPositionInGridSrc2AxisPosition[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 (CDomain::type_attr::rectilinear == domainSrc_->type) 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 (CDomain::type_attr::rectilinear == domainDest_->type) 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    if (std::abs(poleValue - std::abs(domainDest_->lat_start)) < NumTraits<double>::epsilon()) isNorthPole = true;
173    if (std::abs(poleValue - std::abs(domainDest_->lat_end)) < NumTraits<double>::epsilon()) isSouthPole = true;
174
175    CArray<double,1> lon_g ;
176    CArray<double,1> lat_g ;
177
178    if (!domainDest_->lonvalue_1d.isEmpty() && !domainDest_->latvalue_1d.isEmpty())
179    {
180                domainDest_->AllgatherRectilinearLonLat(domainDest_->lonvalue_1d,domainDest_->latvalue_1d, lon_g,lat_g) ;
181        }
182        else if (! domainDest_->latvalue_rectilinear_read_from_file.isEmpty() && ! domainDest_->lonvalue_rectilinear_read_from_file.isEmpty() )
183    {
184                lat_g=domainDest_->latvalue_rectilinear_read_from_file ;
185                lon_g=domainDest_->lonvalue_rectilinear_read_from_file ;
186        }
187        else if (!domainDest_->lon_start.isEmpty() && !domainDest_->lon_end.isEmpty() &&
188                 !domainDest_->lat_start.isEmpty() && !domainDest_->lat_end.isEmpty())
189        {
190          double step=(domainDest_->lon_end-domainDest_->lon_start)/domainDest_->ni_glo ;
191          for(int i=0; i<domainDest_->ni_glo; ++i) lon_g(i)=domainDest_->lon_start+i*step ;
192          step=(domainDest_->lat_end-domainDest_->lat_start)/domainDest_->nj_glo ;
193          for(int i=0; i<domainDest_->ni_glo; ++i) lat_g(i)=domainDest_->lat_start+i*step ;
194        }
195        else ERROR("void CDomainAlgorithmInterpolate::computeRemap()",<<"Cannot compute bounds for rectilinear domain") ;
196    if (std::abs(poleValue - std::abs(lat_g(0))) < NumTraits<double>::epsilon()) isNorthPole = true;
197    if (std::abs(poleValue - std::abs(lat_g(domainDest_->nj_glo-1))) < NumTraits<double>::epsilon()) isSouthPole = true;
198
199
200
201
202    if (isNorthPole && (0 == domainDest_->jbegin.getValue()))
203    {
204      int ibegin = domainDest_->ibegin.getValue();
205      for (i = 0; i < niDest; ++i)
206      {
207        interpMapValueNorthPole[i+ibegin];
208      }
209    }
210
211    if (isSouthPole && (domainDest_->nj_glo.getValue() == (domainDest_->jbegin.getValue() + njDest)))
212    {
213      int ibegin = domainDest_->ibegin.getValue();
214      int njGlo = domainDest_->nj_glo.getValue();
215      int niGlo = domainDest_->ni_glo.getValue();
216      for (i = 0; i < niDest; ++i)
217      {
218        k = (njGlo - 1)*niGlo + i + ibegin;
219        interpMapValueSouthPole[k];
220      }
221    }
222
223    // Ok, fill in boundary values for rectangular domain
224    nVertexDest = constNVertex;
225    domainDest_->fillInRectilinearBoundLonLat(lon_g,lat_g, boundsLonDest, boundsLatDest);
226  }
227
228
229
230  // Ok, now use mapper to calculate
231  int nSrcLocal = domainSrc_->i_index.numElements();
232  int nDstLocal = domainDest_->i_index.numElements();
233  long int * globalSrc = new long int [nSrcLocal];
234  long int * globalDst = new long int [nDstLocal];
235
236  long int globalIndex;
237  int i_ind, j_ind;
238  for (int idx = 0; idx < nSrcLocal; ++idx)
239  {
240    i_ind=domainSrc_->i_index(idx) ;
241    j_ind=domainSrc_->j_index(idx) ;
242
243    globalIndex = i_ind + j_ind * domainSrc_->ni_glo;
244    globalSrc[idx] = globalIndex;
245  }
246
247  for (int idx = 0; idx < nDstLocal; ++idx)
248  {
249    i_ind=domainDest_->i_index(idx) ;
250    j_ind=domainDest_->j_index(idx) ;
251
252    globalIndex = i_ind + j_ind * domainDest_->ni_glo;
253    globalDst[idx] = globalIndex;
254  }
255
256
257  // Calculate weight index
258  Mapper mapper(client->intraComm);
259  mapper.setVerbosity(PROGRESS) ;
260
261     
262  // supress masked data for the source
263  int nSrcLocalUnmasked = 0 ;
264  for (int idx=0 ; idx < nSrcLocal; idx++) if (domainSrc_->localMask(idx)) ++nSrcLocalUnmasked ;
265
266
267  CArray<double,2> boundsLonSrcUnmasked(nVertexSrc,nSrcLocalUnmasked);
268  CArray<double,2> boundsLatSrcUnmasked(nVertexSrc,nSrcLocalUnmasked);
269  long int * globalSrcUnmasked = new long int [nSrcLocalUnmasked];
270
271  nSrcLocalUnmasked=0 ;
272  for (int idx=0 ; idx < nSrcLocal; idx++)
273  {
274    if (domainSrc_->localMask(idx))
275    {
276      for(int n=0;n<nVertexSrc;++n)
277      {
278        boundsLonSrcUnmasked(n,nSrcLocalUnmasked) = boundsLonSrc(n,idx) ;
279        boundsLatSrcUnmasked(n,nSrcLocalUnmasked) = boundsLatSrc(n,idx) ;
280      }
281      globalSrcUnmasked[nSrcLocalUnmasked]=globalSrc[idx] ;
282      ++nSrcLocalUnmasked ;
283    }
284  }
285
286
287  int nDstLocalUnmasked = 0 ;
288  for (int idx=0 ; idx < nDstLocal; idx++) if (domainDest_->localMask(idx)) ++nDstLocalUnmasked ;
289
290  CArray<double,2> boundsLonDestUnmasked(nVertexDest,nDstLocalUnmasked);
291  CArray<double,2> boundsLatDestUnmasked(nVertexDest,nDstLocalUnmasked);
292  long int * globalDstUnmasked = new long int [nDstLocalUnmasked];
293
294  nDstLocalUnmasked=0 ;
295  for (int idx=0 ; idx < nDstLocal; idx++)
296  {
297    if (domainDest_->localMask(idx))
298    {
299      for(int n=0;n<nVertexDest;++n)
300      {
301        boundsLonDestUnmasked(n,nDstLocalUnmasked) = boundsLonDest(n,idx) ;
302        boundsLatDestUnmasked(n,nDstLocalUnmasked) = boundsLatDest(n,idx) ;
303      }
304      globalDstUnmasked[nDstLocalUnmasked]=globalDst[idx] ;
305      ++nDstLocalUnmasked ;
306    }
307  }
308
309  mapper.setSourceMesh(boundsLonSrcUnmasked.dataFirst(), boundsLatSrcUnmasked.dataFirst(), nVertexSrc, nSrcLocalUnmasked, &srcPole[0], globalSrcUnmasked);
310  mapper.setTargetMesh(boundsLonDestUnmasked.dataFirst(), boundsLatDestUnmasked.dataFirst(), nVertexDest, nDstLocalUnmasked, &dstPole[0], globalDstUnmasked);
311
312  std::vector<double> timings = mapper.computeWeights(orderInterp,renormalize);
313
314  std::map<int,std::vector<std::pair<int,double> > > interpMapValue;
315  std::map<int,std::vector<std::pair<int,double> > >::const_iterator iteNorthPole = interpMapValueNorthPole.end(),
316                                                                     iteSouthPole = interpMapValueSouthPole.end();
317  for (int idx = 0;  idx < mapper.nWeights; ++idx)
318  {
319    interpMapValue[mapper.targetWeightId[idx]].push_back(make_pair(mapper.sourceWeightId[idx],mapper.remapMatrix[idx]));
320    if (iteNorthPole != interpMapValueNorthPole.find(mapper.targetWeightId[idx]))
321    {
322      interpMapValueNorthPole[mapper.targetWeightId[idx]].push_back(make_pair(mapper.sourceWeightId[idx],mapper.remapMatrix[idx]));
323    }
324
325    if (iteSouthPole != interpMapValueSouthPole.find(mapper.targetWeightId[idx]))
326    {
327      interpMapValueSouthPole[mapper.targetWeightId[idx]].push_back(make_pair(mapper.sourceWeightId[idx],mapper.remapMatrix[idx]));
328    }
329  }
330  int niGloDst = domainDest_->ni_glo.getValue();
331  processPole(interpMapValueNorthPole, niGloDst);
332  processPole(interpMapValueSouthPole, niGloDst);
333
334  if (!interpMapValueNorthPole.empty())
335  {
336     std::map<int,std::vector<std::pair<int,double> > >::iterator itNorthPole = interpMapValueNorthPole.begin();
337     for (; itNorthPole != iteNorthPole; ++itNorthPole)
338     {
339       if (!(itNorthPole->second.empty()))
340        itNorthPole->second.swap(interpMapValue[itNorthPole->first]);
341     }
342  }
343
344  if (!interpMapValueSouthPole.empty())
345  {
346     std::map<int,std::vector<std::pair<int,double> > >::iterator itSouthPole = interpMapValueSouthPole.begin();
347     for (; itSouthPole != iteSouthPole; ++itSouthPole)
348     {
349       if (!(itSouthPole->second.empty()))
350        itSouthPole->second.swap(interpMapValue[itSouthPole->first]);
351     }
352  }
353
354  exchangeRemapInfo(interpMapValue);
355
356  delete [] globalSrc;
357  delete [] globalSrcUnmasked;
358  delete [] globalDst;
359  delete [] globalDstUnmasked;
360
361}
362
363void CDomainAlgorithmInterpolate::processPole(std::map<int,std::vector<std::pair<int,double> > >& interMapValuePole,
364                                              int nbGlobalPointOnPole)
365{
366  CContext* context = CContext::getCurrent();
367  CContextClient* client=context->client;
368
369  MPI_Comm poleComme(MPI_COMM_NULL);
370  MPI_Comm_split(client->intraComm, interMapValuePole.empty() ? MPI_UNDEFINED : 1, 0, &poleComme);
371  if (MPI_COMM_NULL != poleComme)
372  {
373    int nbClientPole;
374    MPI_Comm_size(poleComme, &nbClientPole);
375
376    std::map<int,std::vector<std::pair<int,double> > >::iterator itePole = interMapValuePole.end(), itPole,
377                                                                 itbPole = interMapValuePole.begin();
378
379    int nbWeight = 0;
380    for (itPole = itbPole; itPole != itePole; ++itPole)
381       nbWeight += itPole->second.size();
382
383    std::vector<int> recvCount(nbClientPole,0);
384    std::vector<int> displ(nbClientPole,0);
385    MPI_Allgather(&nbWeight,1,MPI_INT,&recvCount[0],1,MPI_INT,poleComme) ;
386
387    displ[0]=0;
388    for(int n=1;n<nbClientPole;++n) displ[n]=displ[n-1]+recvCount[n-1] ;
389    int recvSize=displ[nbClientPole-1]+recvCount[nbClientPole-1] ;
390
391    std::vector<int> sendSourceIndexBuff(nbWeight);
392    std::vector<double> sendSourceWeightBuff(nbWeight);
393    int k = 0;
394    for (itPole = itbPole; itPole != itePole; ++itPole)
395    {
396      for (int idx = 0; idx < itPole->second.size(); ++idx)
397      {
398        sendSourceIndexBuff[k] = (itPole->second)[idx].first;
399        sendSourceWeightBuff[k] = (itPole->second)[idx].second;
400        ++k;
401      }
402    }
403
404    std::vector<int> recvSourceIndexBuff(recvSize);
405    std::vector<double> recvSourceWeightBuff(recvSize);
406
407    // Gather all index and weight for pole
408    MPI_Allgatherv(&sendSourceIndexBuff[0],nbWeight,MPI_INT,&recvSourceIndexBuff[0],&recvCount[0],&displ[0],MPI_INT,poleComme);
409    MPI_Allgatherv(&sendSourceWeightBuff[0],nbWeight,MPI_DOUBLE,&recvSourceWeightBuff[0],&recvCount[0],&displ[0],MPI_DOUBLE,poleComme);
410
411    std::map<int,double> recvTemp;
412    for (int idx = 0; idx < recvSize; ++idx)
413    {
414      if (recvTemp.end() != recvTemp.find(recvSourceIndexBuff[idx]))
415        recvTemp[recvSourceIndexBuff[idx]] += recvSourceWeightBuff[idx]/nbGlobalPointOnPole;
416      else
417        recvTemp[recvSourceIndexBuff[idx]] = 0.0;
418    }
419
420    std::map<int,double>::const_iterator itRecvTemp, itbRecvTemp = recvTemp.begin(), iteRecvTemp = recvTemp.end();
421
422    for (itPole = itbPole; itPole != itePole; ++itPole)
423    {
424      itPole->second.clear();
425      for (itRecvTemp = itbRecvTemp; itRecvTemp != iteRecvTemp; ++itRecvTemp)
426          itPole->second.push_back(make_pair(itRecvTemp->first, itRecvTemp->second));
427    }
428  }
429
430}
431
432/*!
433  Compute the index mapping between domain on grid source and one on grid destination
434*/
435void CDomainAlgorithmInterpolate::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
436{
437  if (!interpDomain_->file.isEmpty())
438    readRemapInfo();
439  else
440    computeRemap();
441}
442
443void CDomainAlgorithmInterpolate::readRemapInfo()
444{
445  CContext* context = CContext::getCurrent();
446  CContextClient* client=context->client;
447  int clientRank = client->clientRank;
448
449  std::string filename = interpDomain_->file.getValue();
450  std::map<int,std::vector<std::pair<int,double> > > interpMapValue;
451  readInterpolationInfo(filename, interpMapValue);
452
453  exchangeRemapInfo(interpMapValue);
454}
455
456
457/*!
458  Read remap information from file then distribute it among clients
459*/
460void CDomainAlgorithmInterpolate::exchangeRemapInfo(std::map<int,std::vector<std::pair<int,double> > >& interpMapValue)
461{
462  CContext* context = CContext::getCurrent();
463  CContextClient* client=context->client;
464  int clientRank = client->clientRank;
465
466  this->transformationMapping_.resize(1);
467  this->transformationWeight_.resize(1);
468
469  TransformationIndexMap& transMap = this->transformationMapping_[0];
470  TransformationWeightMap& transWeight = this->transformationWeight_[0];
471
472  boost::unordered_map<size_t,int> globalIndexOfDomainDest;
473  int ni = domainDest_->ni.getValue();
474  int nj = domainDest_->nj.getValue();
475  int ni_glo = domainDest_->ni_glo.getValue();
476  size_t globalIndex;
477  int nIndexSize = domainDest_->i_index.numElements(), i_ind, j_ind;
478  for (int idx = 0; idx < nIndexSize; ++idx)
479  {
480    i_ind=domainDest_->i_index(idx) ;
481    j_ind=domainDest_->j_index(idx) ;
482
483    globalIndex = i_ind + j_ind * ni_glo;
484    globalIndexOfDomainDest[globalIndex] = clientRank;
485  }
486
487  CClientServerMappingDistributed domainIndexClientClientMapping(globalIndexOfDomainDest,
488                                                                 client->intraComm,
489                                                                 true);
490  CArray<size_t,1> globalIndexInterp(interpMapValue.size());
491  std::map<int,std::vector<std::pair<int,double> > >::const_iterator itb = interpMapValue.begin(), it,
492                                                                     ite = interpMapValue.end();
493  size_t globalIndexCount = 0;
494  for (it = itb; it != ite; ++it)
495  {
496    globalIndexInterp(globalIndexCount) = it->first;
497    ++globalIndexCount;
498  }
499
500  domainIndexClientClientMapping.computeServerIndexMapping(globalIndexInterp);
501  const CClientServerMapping::GlobalIndexMap& globalIndexInterpSendToClient = domainIndexClientClientMapping.getGlobalIndexOnServer();
502
503  //Inform each client number of index they will receive
504  int nbClient = client->clientSize;
505  int* sendBuff = new int[nbClient];
506  int* recvBuff = new int[nbClient];
507  for (int i = 0; i < nbClient; ++i)
508  {
509    sendBuff[i] = 0;
510    recvBuff[i] = 0;
511  }
512  int sendBuffSize = 0;
513  CClientServerMapping::GlobalIndexMap::const_iterator itbMap = globalIndexInterpSendToClient.begin(), itMap,
514                                                       iteMap = globalIndexInterpSendToClient.end();
515  for (itMap = itbMap; itMap != iteMap; ++itMap)
516  {
517    const std::vector<size_t>& tmp = itMap->second;
518    int sizeIndex = 0, mapSize = (itMap->second).size();
519    for (int idx = 0; idx < mapSize; ++idx)
520    {
521//      sizeIndex += interpMapValue.at((itMap->second)[idx]).size();
522      sizeIndex += (interpMapValue[(int)(itMap->second)[idx]]).size();
523    }
524    sendBuff[itMap->first] = sizeIndex;
525    sendBuffSize += sizeIndex;
526  }
527
528
529  MPI_Allreduce(sendBuff, recvBuff, nbClient, MPI_INT, MPI_SUM, client->intraComm);
530
531  int* sendIndexDestBuff = new int [sendBuffSize];
532  int* sendIndexSrcBuff  = new int [sendBuffSize];
533  double* sendWeightBuff = new double [sendBuffSize];
534
535  std::vector<MPI_Request> sendRequest;
536
537  int sendOffSet = 0, l = 0;
538  for (itMap = itbMap; itMap != iteMap; ++itMap)
539  {
540    const std::vector<size_t>& indexToSend = itMap->second;
541    int mapSize = indexToSend.size();
542    int k = 0;
543    for (int idx = 0; idx < mapSize; ++idx)
544    {
545      std::vector<std::pair<int,double> >& interpMap = interpMapValue[(int)indexToSend[idx]]; //interpMapValue.at(indexToSend[idx]);
546      for (int i = 0; i < interpMap.size(); ++i)
547      {
548        sendIndexDestBuff[l] = indexToSend[idx];
549        sendIndexSrcBuff[l]  = interpMap[i].first;
550        sendWeightBuff[l]    = interpMap[i].second;
551        ++k;
552        ++l;
553      }
554    }
555
556    sendRequest.push_back(MPI_Request());
557    MPI_Isend(sendIndexDestBuff + sendOffSet,
558             k,
559             MPI_INT,
560             itMap->first,
561             MPI_DOMAIN_INTERPOLATION_DEST_INDEX,
562             client->intraComm,
563             &sendRequest.back());
564    sendRequest.push_back(MPI_Request());
565    MPI_Isend(sendIndexSrcBuff + sendOffSet,
566             k,
567             MPI_INT,
568             itMap->first,
569             MPI_DOMAIN_INTERPOLATION_SRC_INDEX,
570             client->intraComm,
571             &sendRequest.back());
572    sendRequest.push_back(MPI_Request());
573    MPI_Isend(sendWeightBuff + sendOffSet,
574             k,
575             MPI_DOUBLE,
576             itMap->first,
577             MPI_DOMAIN_INTERPOLATION_WEIGHT,
578             client->intraComm,
579             &sendRequest.back());
580    sendOffSet += k;
581  }
582
583  int recvBuffSize = recvBuff[clientRank];
584  int* recvIndexDestBuff = new int [recvBuffSize];
585  int* recvIndexSrcBuff  = new int [recvBuffSize];
586  double* recvWeightBuff = new double [recvBuffSize];
587  int receivedSize = 0;
588  int clientSrcRank;
589  while (receivedSize < recvBuffSize)
590  {
591    MPI_Status recvStatus;
592    MPI_Recv((recvIndexDestBuff + receivedSize),
593             recvBuffSize,
594             MPI_INT,
595             MPI_ANY_SOURCE,
596             MPI_DOMAIN_INTERPOLATION_DEST_INDEX,
597             client->intraComm,
598             &recvStatus);
599
600    int countBuff = 0;
601    MPI_Get_count(&recvStatus, MPI_INT, &countBuff);
602    clientSrcRank = recvStatus.MPI_SOURCE;
603
604    MPI_Recv((recvIndexSrcBuff + receivedSize),
605             recvBuffSize,
606             MPI_INT,
607             clientSrcRank,
608             MPI_DOMAIN_INTERPOLATION_SRC_INDEX,
609             client->intraComm,
610             &recvStatus);
611
612    MPI_Recv((recvWeightBuff + receivedSize),
613             recvBuffSize,
614             MPI_DOUBLE,
615             clientSrcRank,
616             MPI_DOMAIN_INTERPOLATION_WEIGHT,
617             client->intraComm,
618             &recvStatus);
619
620    for (int idx = 0; idx < countBuff; ++idx)
621    {
622      transMap[*(recvIndexDestBuff + receivedSize + idx)].push_back(*(recvIndexSrcBuff + receivedSize + idx));
623      transWeight[*(recvIndexDestBuff + receivedSize + idx)].push_back(*(recvWeightBuff + receivedSize + idx));
624    }
625    receivedSize += countBuff;
626  }
627
628  std::vector<MPI_Status> requestStatus(sendRequest.size());
629  MPI_Waitall(sendRequest.size(), &sendRequest[0], MPI_STATUS_IGNORE);
630
631  delete [] sendIndexDestBuff;
632  delete [] sendIndexSrcBuff;
633  delete [] sendWeightBuff;
634  delete [] recvIndexDestBuff;
635  delete [] recvIndexSrcBuff;
636  delete [] recvWeightBuff;
637  delete [] sendBuff;
638  delete [] recvBuff;
639}
640
641/*!
642  Read interpolation information from a file
643  \param [in] filename interpolation file
644  \param [in/out] interpMapValue Mapping between (global) index of domain on grid destination and
645         corresponding global index of domain and associated weight value on grid source
646*/
647void CDomainAlgorithmInterpolate::readInterpolationInfo(std::string& filename,
648                                                        std::map<int,std::vector<std::pair<int,double> > >& interpMapValue)
649{
650  int ncid ;
651  int weightDimId ;
652  size_t nbWeightGlo ;
653
654  CContext* context = CContext::getCurrent();
655  CContextClient* client=context->client;
656  int clientRank = client->clientRank;
657  int clientSize = client->clientSize;
658
659  nc_open(filename.c_str(),NC_NOWRITE, &ncid) ;
660  nc_inq_dimid(ncid,"n_weight",&weightDimId) ;
661  nc_inq_dimlen(ncid,weightDimId,&nbWeightGlo) ;
662
663  size_t nbWeight ;
664  size_t start ;
665  size_t div = nbWeightGlo/clientSize ;
666  size_t mod = nbWeightGlo%clientSize ;
667  if (clientRank < mod)
668  {
669    nbWeight=div+1 ;
670    start=clientRank*(div+1) ;
671  }
672  else
673  {
674    nbWeight=div ;
675    start= mod * (div+1) + (clientRank-mod) * div ;
676  }
677
678  double* weight=new double[nbWeight] ;
679  int weightId ;
680  nc_inq_varid (ncid, "weight", &weightId) ;
681  nc_get_vara_double(ncid, weightId, &start, &nbWeight, weight) ;
682
683  long* srcIndex=new long[nbWeight] ;
684  int srcIndexId ;
685  nc_inq_varid (ncid, "src_idx", &srcIndexId) ;
686  nc_get_vara_long(ncid, srcIndexId, &start, &nbWeight, srcIndex) ;
687
688  long* dstIndex=new long[nbWeight] ;
689  int dstIndexId ;
690  nc_inq_varid (ncid, "dst_idx", &dstIndexId) ;
691  nc_get_vara_long(ncid, dstIndexId, &start, &nbWeight, dstIndex) ;
692
693  for(size_t ind=0; ind<nbWeight;++ind)
694    interpMapValue[dstIndex[ind]-1].push_back(make_pair(srcIndex[ind]-1,weight[ind]));
695}
696
697}
Note: See TracBrowser for help on using the repository browser.