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

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

Implementing dynamic interpolation on axis

+) Change grid transformation to make it more flexible
+) Make some small improvements

Test
+) On Curie
+) All test pass

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