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

Last change on this file since 743 was 743, checked in by mhnguyen, 8 years ago

Processing pole of rectangular grid with latitude = 90

+) Use average value of all points with latitude 90 (-90)

Test
+) On Curie
+) test_remap pass

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