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

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

Improvements for dht

+) Implement adaptive hierarchy for dht, level of hierarchy depends on number of processes
+) Remove some redundant codes

Test
+) On Curie
+) All tests are correct

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