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

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

Making some improvements of transformation algorithm

+) Correct the way to enlisting transformations in an element (domain, axis)
+) Optimize generic transformation to make sure temporary grid to be created on demand
+) Update some mpi tag to prevent conflict
+) Correct some minor stuffs
+) Update documents

Test
+) On Curie
+) all test pass

File size: 21.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
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()
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  boost::unordered_map<size_t,int> globalIndexOfDomainDest;
385  int ni = domainDest_->ni.getValue();
386  int nj = domainDest_->nj.getValue();
387  int ni_glo = domainDest_->ni_glo.getValue();
388  size_t globalIndex;
389  int nIndexSize = domainDest_->i_index.numElements(), i_ind, j_ind;
390  for (int idx = 0; idx < nIndexSize; ++idx)
391  {
392    i_ind=domainDest_->i_index(idx) ;
393    j_ind=domainDest_->j_index(idx) ;
394
395    globalIndex = i_ind + j_ind * ni_glo;
396    globalIndexOfDomainDest[globalIndex] = clientRank;
397  }
398
399  CClientServerMappingDistributed domainIndexClientClientMapping(globalIndexOfDomainDest,
400                                                                 client->intraComm,
401                                                                 true);
402  CArray<size_t,1> globalIndexInterp(interpMapValue.size());
403  std::map<int,std::vector<std::pair<int,double> > >::const_iterator itb = interpMapValue.begin(), it,
404                                                                     ite = interpMapValue.end();
405  size_t globalIndexCount = 0;
406  for (it = itb; it != ite; ++it)
407  {
408    globalIndexInterp(globalIndexCount) = it->first;
409    ++globalIndexCount;
410  }
411
412  domainIndexClientClientMapping.computeServerIndexMapping(globalIndexInterp);
413  const std::map<int, std::vector<size_t> >& globalIndexInterpSendToClient = domainIndexClientClientMapping.getGlobalIndexOnServer();
414
415  //Inform each client number of index they will receive
416  int nbClient = client->clientSize;
417  int* sendBuff = new int[nbClient];
418  int* recvBuff = new int[nbClient];
419  for (int i = 0; i < nbClient; ++i)
420  {
421    sendBuff[i] = 0;
422    recvBuff[i] = 0;
423  }
424  int sendBuffSize = 0;
425  std::map<int, std::vector<size_t> >::const_iterator itbMap = globalIndexInterpSendToClient.begin(), itMap,
426                                                      iteMap = globalIndexInterpSendToClient.end();
427  for (itMap = itbMap; itMap != iteMap; ++itMap)
428  {
429    const std::vector<size_t>& tmp = itMap->second;
430    int sizeIndex = 0, mapSize = (itMap->second).size();
431    for (int idx = 0; idx < mapSize; ++idx)
432    {
433      sizeIndex += interpMapValue.at((itMap->second)[idx]).size();
434    }
435    sendBuff[itMap->first] = sizeIndex;
436    sendBuffSize += sizeIndex;
437  }
438
439
440  MPI_Allreduce(sendBuff, recvBuff, nbClient, MPI_INT, MPI_SUM, client->intraComm);
441
442  int* sendIndexDestBuff = new int [sendBuffSize];
443  int* sendIndexSrcBuff  = new int [sendBuffSize];
444  double* sendWeightBuff = new double [sendBuffSize];
445
446  std::vector<MPI_Request> sendRequest;
447
448  int sendOffSet = 0, l = 0;
449  for (itMap = itbMap; itMap != iteMap; ++itMap)
450  {
451    const std::vector<size_t>& indexToSend = itMap->second;
452    int mapSize = indexToSend.size();
453    int k = 0;
454    for (int idx = 0; idx < mapSize; ++idx)
455    {
456      const std::vector<std::pair<int,double> >& interpMap = interpMapValue.at(indexToSend[idx]);
457      for (int i = 0; i < interpMap.size(); ++i)
458      {
459        sendIndexDestBuff[l] = indexToSend[idx];
460        sendIndexSrcBuff[l]  = interpMap[i].first;
461        sendWeightBuff[l]    = interpMap[i].second;
462        ++k;
463        ++l;
464      }
465    }
466
467    sendRequest.push_back(MPI_Request());
468    MPI_Isend(sendIndexDestBuff + sendOffSet,
469             k,
470             MPI_INT,
471             itMap->first,
472             MPI_DOMAIN_INTERPOLATION_DEST_INDEX,
473             client->intraComm,
474             &sendRequest.back());
475    sendRequest.push_back(MPI_Request());
476    MPI_Isend(sendIndexSrcBuff + sendOffSet,
477             k,
478             MPI_INT,
479             itMap->first,
480             MPI_DOMAIN_INTERPOLATION_SRC_INDEX,
481             client->intraComm,
482             &sendRequest.back());
483    sendRequest.push_back(MPI_Request());
484    MPI_Isend(sendWeightBuff + sendOffSet,
485             k,
486             MPI_DOUBLE,
487             itMap->first,
488             MPI_DOMAIN_INTERPOLATION_WEIGHT,
489             client->intraComm,
490             &sendRequest.back());
491    sendOffSet += k;
492  }
493
494  int recvBuffSize = recvBuff[clientRank];
495  int* recvIndexDestBuff = new int [recvBuffSize];
496  int* recvIndexSrcBuff  = new int [recvBuffSize];
497  double* recvWeightBuff = new double [recvBuffSize];
498  int receivedSize = 0;
499  int clientSrcRank;
500  while (receivedSize < recvBuffSize)
501  {
502    MPI_Status recvStatus;
503    MPI_Recv((recvIndexDestBuff + receivedSize),
504             recvBuffSize,
505             MPI_INT,
506             MPI_ANY_SOURCE,
507             MPI_DOMAIN_INTERPOLATION_DEST_INDEX,
508             client->intraComm,
509             &recvStatus);
510
511    int countBuff = 0;
512    MPI_Get_count(&recvStatus, MPI_INT, &countBuff);
513    clientSrcRank = recvStatus.MPI_SOURCE;
514
515    MPI_Recv((recvIndexSrcBuff + receivedSize),
516             recvBuffSize,
517             MPI_INT,
518             clientSrcRank,
519             MPI_DOMAIN_INTERPOLATION_SRC_INDEX,
520             client->intraComm,
521             &recvStatus);
522
523    MPI_Recv((recvWeightBuff + receivedSize),
524             recvBuffSize,
525             MPI_DOUBLE,
526             clientSrcRank,
527             MPI_DOMAIN_INTERPOLATION_WEIGHT,
528             client->intraComm,
529             &recvStatus);
530
531    for (int idx = 0; idx < countBuff; ++idx)
532    {
533      transformationMapping_[*(recvIndexDestBuff + receivedSize + idx)].push_back(*(recvIndexSrcBuff + receivedSize + idx));
534      transformationWeight_[*(recvIndexDestBuff + receivedSize + idx)].push_back(*(recvWeightBuff + receivedSize + idx));
535    }
536    receivedSize += countBuff;
537  }
538
539  std::vector<MPI_Status> requestStatus(sendRequest.size());
540  MPI_Waitall(sendRequest.size(), &sendRequest[0], MPI_STATUS_IGNORE);
541
542  delete [] sendIndexDestBuff;
543  delete [] sendIndexSrcBuff;
544  delete [] sendWeightBuff;
545  delete [] recvIndexDestBuff;
546  delete [] recvIndexSrcBuff;
547  delete [] recvWeightBuff;
548  delete [] sendBuff;
549  delete [] recvBuff;
550}
551
552/*!
553  Read interpolation information from a file
554  \param [in] filename interpolation file
555  \param [in/out] interpMapValue Mapping between (global) index of domain on grid destination and
556         corresponding global index of domain and associated weight value on grid source
557*/
558void CDomainAlgorithmInterpolate::readInterpolationInfo(std::string& filename,
559                                                        std::map<int,std::vector<std::pair<int,double> > >& interpMapValue)
560{
561  int ncid ;
562  int weightDimId ;
563  size_t nbWeightGlo ;
564
565  CContext* context = CContext::getCurrent();
566  CContextClient* client=context->client;
567  int clientRank = client->clientRank;
568  int clientSize = client->clientSize;
569
570  nc_open(filename.c_str(),NC_NOWRITE, &ncid) ;
571  nc_inq_dimid(ncid,"n_weight",&weightDimId) ;
572  nc_inq_dimlen(ncid,weightDimId,&nbWeightGlo) ;
573
574  size_t nbWeight ;
575  size_t start ;
576  size_t div = nbWeightGlo/clientSize ;
577  size_t mod = nbWeightGlo%clientSize ;
578  if (clientRank < mod)
579  {
580    nbWeight=div+1 ;
581    start=clientRank*(div+1) ;
582  }
583  else
584  {
585    nbWeight=div ;
586    start= mod * (div+1) + (clientRank-mod) * div ;
587  }
588
589  double* weight=new double[nbWeight] ;
590  int weightId ;
591  nc_inq_varid (ncid, "weight", &weightId) ;
592  nc_get_vara_double(ncid, weightId, &start, &nbWeight, weight) ;
593
594  long* srcIndex=new long[nbWeight] ;
595  int srcIndexId ;
596  nc_inq_varid (ncid, "src_idx", &srcIndexId) ;
597  nc_get_vara_long(ncid, srcIndexId, &start, &nbWeight, srcIndex) ;
598
599  long* dstIndex=new long[nbWeight] ;
600  int dstIndexId ;
601  nc_inq_varid (ncid, "dst_idx", &dstIndexId) ;
602  nc_get_vara_long(ncid, dstIndexId, &start, &nbWeight, dstIndex) ;
603
604  for(size_t ind=0; ind<nbWeight;++ind)
605    interpMapValue[dstIndex[ind]-1].push_back(make_pair(srcIndex[ind]-1,weight[ind]));
606}
607
608}
Note: See TracBrowser for help on using the repository browser.