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

Last change on this file since 809 was 809, checked in by ymipsl, 7 years ago

Fix problems for interpolation onto regular domain.

YM

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