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

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

Fix problem using rectilinear grid from file, or self-generated

  • Defining boundaries at poles was incorrect

YM

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