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

Last change on this file since 846 was 846, checked in by ymipsl, 6 years ago

Management of masked cell for conservative interpolation.

YM

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