source: XIOS/dev/XIOS_DEV_CMIP6/src/transformation/domain_algorithm_interpolate.cpp @ 1264

Last change on this file since 1264 was 1264, checked in by ymipsl, 7 years ago
  • Add new attribute : detect_missing_value on "interpolate_domain" element. It will replace the standard detect_missing_value set on the field. It will be added progressively to all spatial transformation.
  • Now, when detecting missing value, horizontal interpolation do a correct renormalization.

YM

File size: 33.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#include "domain.hpp"
19#include "grid_transformation_factory_impl.hpp"
20#include "interpolate_domain.hpp"
21#include "grid.hpp"
22
23namespace xios {
24CGenericAlgorithmTransformation* CDomainAlgorithmInterpolate::create(CGrid* gridDst, CGrid* gridSrc,
25                                                                     CTransformation<CDomain>* transformation,
26                                                                     int elementPositionInGrid,
27                                                                     std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
28                                                                     std::map<int, int>& elementPositionInGridSrc2AxisPosition,
29                                                                     std::map<int, int>& elementPositionInGridSrc2DomainPosition,
30                                                                     std::map<int, int>& elementPositionInGridDst2ScalarPosition,
31                                                                     std::map<int, int>& elementPositionInGridDst2AxisPosition,
32                                                                     std::map<int, int>& elementPositionInGridDst2DomainPosition)
33{
34  std::vector<CDomain*> domainListDestP = gridDst->getDomains();
35  std::vector<CDomain*> domainListSrcP  = gridSrc->getDomains();
36
37  CInterpolateDomain* interpolateDomain = dynamic_cast<CInterpolateDomain*> (transformation);
38  int domainDstIndex = elementPositionInGridDst2DomainPosition[elementPositionInGrid];
39  int domainSrcIndex = elementPositionInGridSrc2DomainPosition[elementPositionInGrid];
40
41  return (new CDomainAlgorithmInterpolate(domainListDestP[domainDstIndex], domainListSrcP[domainSrcIndex], interpolateDomain));
42}
43
44bool CDomainAlgorithmInterpolate::registerTrans()
45{
46  CGridTransformationFactory<CDomain>::registerTransformation(TRANS_INTERPOLATE_DOMAIN, create);
47}
48
49CDomainAlgorithmInterpolate::CDomainAlgorithmInterpolate(CDomain* domainDestination, CDomain* domainSource, CInterpolateDomain* interpDomain)
50: CDomainAlgorithmTransformation(domainDestination, domainSource), interpDomain_(interpDomain), writeToFile_(false), readFromFile_(false)
51{
52  CContext* context = CContext::getCurrent();
53  interpDomain_->checkValid(domainSource);
54
55  detectMissingValue = interpDomain_->detect_missing_value ;
56  renormalize = interpDomain_->renormalize ;
57  quantity = interpDomain_->quantity ;
58 
59  fileToReadWrite_ = "xios_interpolation_weights_";
60
61  if (interpDomain_->weight_filename.isEmpty())
62  {
63    fileToReadWrite_ += context->getId() + "_" + 
64                    domainSource->getDomainOutputName() + "_" + 
65                    domainDestination->getDomainOutputName() + ".nc";   
66  }
67  else 
68    fileToReadWrite_ = interpDomain_->weight_filename;
69
70  ifstream f(fileToReadWrite_.c_str()); 
71  switch (interpDomain_->mode)
72  {
73    case CInterpolateDomain::mode_attr::read:
74      readFromFile_ = true;     
75      break;
76    case CInterpolateDomain::mode_attr::compute:
77      readFromFile_ = false;
78      break;
79    case CInterpolateDomain::mode_attr::read_or_compute:     
80      if (!f.good())
81        readFromFile_ = false;
82      else
83        readFromFile_ = true;
84      break;
85    default:
86      break;
87  } 
88
89  writeToFile_ = interpDomain_->write_weight; 
90   
91}
92
93/*!
94  Compute remap with integrated remap calculation module
95*/
96void CDomainAlgorithmInterpolate::computeRemap()
97{
98  using namespace sphereRemap;
99
100  CContext* context = CContext::getCurrent();
101  CContextClient* client=context->client;
102  int clientRank = client->clientRank;
103  int i, j, k, idx;
104  std::vector<double> srcPole(3,0), dstPole(3,0);
105  int orderInterp = interpDomain_->order.getValue();
106
107
108  const double poleValue = 90.0;
109  const int constNVertex = 4; // Value by default number of vertex for rectangular domain
110  int nVertexSrc, nVertexDest;
111  nVertexSrc = nVertexDest = constNVertex;
112
113  // First of all, try to retrieve the boundary values of domain source and domain destination
114  int localDomainSrcSize = domainSrc_->i_index.numElements();
115  int niSrc = domainSrc_->ni.getValue(), njSrc = domainSrc_->nj.getValue();
116  bool hasBoundSrc = domainSrc_->hasBounds;
117  if (hasBoundSrc) nVertexSrc = domainSrc_->nvertex.getValue();
118  CArray<double,2> boundsLonSrc(nVertexSrc,localDomainSrcSize);
119  CArray<double,2> boundsLatSrc(nVertexSrc,localDomainSrcSize);
120
121  if (domainSrc_->hasPole) srcPole[2] = 1;
122  if (hasBoundSrc)  // Suppose that domain source is curvilinear or unstructured
123  {
124    if (!domainSrc_->bounds_lon_2d.isEmpty())
125    {
126      for (j = 0; j < njSrc; ++j)
127        for (i = 0; i < niSrc; ++i)
128        {
129          k=j*niSrc+i;
130          for(int n=0;n<nVertexSrc;++n)
131          {
132            boundsLonSrc(n,k) = domainSrc_->bounds_lon_2d(n,i,j);
133            boundsLatSrc(n,k) = domainSrc_->bounds_lat_2d(n,i,j);
134          }
135        }
136    }
137    else
138    {
139      boundsLonSrc = domainSrc_->bounds_lon_1d;
140      boundsLatSrc = domainSrc_->bounds_lat_1d;
141    }
142  }
143  else // if domain source is rectilinear, not do anything now
144  {
145    CArray<double,1> lon_g ;
146    CArray<double,1> lat_g ;
147
148    if (!domainSrc_->lonvalue_1d.isEmpty() && !domainSrc_->latvalue_1d.isEmpty())
149    {
150      domainSrc_->AllgatherRectilinearLonLat(domainSrc_->lonvalue_1d,domainSrc_->latvalue_1d, lon_g,lat_g) ;
151    }
152    else if (! domainSrc_->latvalue_rectilinear_read_from_file.isEmpty() && ! domainSrc_->lonvalue_rectilinear_read_from_file.isEmpty() )
153    {
154      lat_g=domainSrc_->latvalue_rectilinear_read_from_file ;
155      lon_g=domainSrc_->lonvalue_rectilinear_read_from_file ;
156    }
157    else if (!domainSrc_->lon_start.isEmpty() && !domainSrc_->lon_end.isEmpty() &&
158             !domainSrc_->lat_start.isEmpty() && !domainSrc_->lat_end.isEmpty())
159    {
160      double step=(domainSrc_->lon_end-domainSrc_->lon_start)/domainSrc_->ni_glo ;
161      for (int i=0; i<domainSrc_->ni_glo; ++i) lon_g(i)=domainSrc_->lon_start+i*step ;
162      step=(domainSrc_->lat_end-domainSrc_->lat_start)/domainSrc_->nj_glo ;
163      for (int i=0; i<domainSrc_->ni_glo; ++i) lat_g(i)=domainSrc_->lat_start+i*step ;
164    }
165    else ERROR("void CDomainAlgorithmInterpolate::computeRemap()",<<"Cannot compute bounds for rectilinear domain") ;
166
167    nVertexSrc = constNVertex;
168    domainSrc_->fillInRectilinearBoundLonLat(lon_g,lat_g, boundsLonSrc, boundsLatSrc);
169  }
170
171  std::map<int,std::vector<std::pair<int,double> > > interpMapValueNorthPole;
172  std::map<int,std::vector<std::pair<int,double> > > interpMapValueSouthPole;
173
174  int localDomainDestSize = domainDest_->i_index.numElements();
175  int niDest = domainDest_->ni.getValue(), njDest = domainDest_->nj.getValue();
176  bool hasBoundDest = domainDest_->hasBounds;
177  if (hasBoundDest) nVertexDest = domainDest_->nvertex.getValue();
178  CArray<double,2> boundsLonDest(nVertexDest,localDomainDestSize);
179  CArray<double,2> boundsLatDest(nVertexDest,localDomainDestSize);
180
181  if (domainDest_->hasPole) dstPole[2] = 1;
182  if (hasBoundDest)
183  {
184    if (!domainDest_->bounds_lon_2d.isEmpty())
185    {
186      for (j = 0; j < njDest; ++j)
187        for (i = 0; i < niDest; ++i)
188        {
189          k=j*niDest+i;
190          for(int n=0;n<nVertexDest;++n)
191          {
192            boundsLonDest(n,k) = domainDest_->bounds_lon_2d(n,i,j);
193            boundsLatDest(n,k) = domainDest_->bounds_lat_2d(n,i,j);
194          }
195        }
196    }
197    else
198    {
199      boundsLonDest = domainDest_->bounds_lon_1d;
200      boundsLatDest = domainDest_->bounds_lat_1d;
201    }
202  }
203  else
204  {
205    bool isNorthPole = false;
206    bool isSouthPole = false;
207
208    CArray<double,1> lon_g ;
209    CArray<double,1> lat_g ;
210
211    if (!domainDest_->lonvalue_1d.isEmpty() && !domainDest_->latvalue_1d.isEmpty())
212    {
213      domainDest_->AllgatherRectilinearLonLat(domainDest_->lonvalue_1d,domainDest_->latvalue_1d, lon_g,lat_g) ;
214    }
215    else if (! domainDest_->latvalue_rectilinear_read_from_file.isEmpty() && ! domainDest_->lonvalue_rectilinear_read_from_file.isEmpty() )
216    {
217      lat_g=domainDest_->latvalue_rectilinear_read_from_file ;
218      lon_g=domainDest_->lonvalue_rectilinear_read_from_file ;
219    }
220    else if (!domainDest_->lon_start.isEmpty() && !domainDest_->lon_end.isEmpty() &&
221             !domainDest_->lat_start.isEmpty() && !domainDest_->lat_end.isEmpty())
222    {
223      double step=(domainDest_->lon_end-domainDest_->lon_start)/domainDest_->ni_glo ;
224      for(int i=0; i<domainDest_->ni_glo; ++i) lon_g(i)=domainDest_->lon_start+i*step ;
225      step=(domainDest_->lat_end-domainDest_->lat_start)/domainDest_->nj_glo ;
226      for(int i=0; i<domainDest_->ni_glo; ++i) lat_g(i)=domainDest_->lat_start+i*step ;
227    }
228    else ERROR("void CDomainAlgorithmInterpolate::computeRemap()",<<"Cannot compute bounds for rectilinear domain") ;
229   
230    if (std::abs(poleValue - std::abs(lat_g(0))) < NumTraits<double>::epsilon()) isNorthPole = true;
231    if (std::abs(poleValue - std::abs(lat_g(domainDest_->nj_glo-1))) < NumTraits<double>::epsilon()) isSouthPole = true;
232
233
234
235
236    if (isNorthPole && (0 == domainDest_->jbegin.getValue()))
237    {
238      int ibegin = domainDest_->ibegin.getValue();
239      for (i = 0; i < niDest; ++i)
240      {
241        interpMapValueNorthPole[i+ibegin];
242      }
243    }
244
245    if (isSouthPole && (domainDest_->nj_glo.getValue() == (domainDest_->jbegin.getValue() + njDest)))
246    {
247      int ibegin = domainDest_->ibegin.getValue();
248      int njGlo = domainDest_->nj_glo.getValue();
249      int niGlo = domainDest_->ni_glo.getValue();
250      for (i = 0; i < niDest; ++i)
251      {
252        k = (njGlo - 1)*niGlo + i + ibegin;
253        interpMapValueSouthPole[k];
254      }
255    }
256
257    // Ok, fill in boundary values for rectangular domain
258    nVertexDest = constNVertex;
259    domainDest_->fillInRectilinearBoundLonLat(lon_g,lat_g, boundsLonDest, boundsLatDest);
260  }
261
262
263
264  // Ok, now use mapper to calculate
265  int nSrcLocal = domainSrc_->i_index.numElements();
266  int nDstLocal = domainDest_->i_index.numElements();
267  long int * globalSrc = new long int [nSrcLocal];
268  long int * globalDst = new long int [nDstLocal];
269
270  long int globalIndex;
271  int i_ind, j_ind;
272  for (int idx = 0; idx < nSrcLocal; ++idx)
273  {
274    i_ind=domainSrc_->i_index(idx) ;
275    j_ind=domainSrc_->j_index(idx) ;
276
277    globalIndex = i_ind + j_ind * domainSrc_->ni_glo;
278    globalSrc[idx] = globalIndex;
279  }
280
281  for (int idx = 0; idx < nDstLocal; ++idx)
282  {
283    i_ind=domainDest_->i_index(idx) ;
284    j_ind=domainDest_->j_index(idx) ;
285
286    globalIndex = i_ind + j_ind * domainDest_->ni_glo;
287    globalDst[idx] = globalIndex;
288  }
289
290
291  // Calculate weight index
292  Mapper mapper(client->intraComm);
293  mapper.setVerbosity(PROGRESS) ;
294
295     
296  // supress masked data for the source
297  int nSrcLocalUnmasked = 0 ;
298  for (int idx=0 ; idx < nSrcLocal; idx++) if (domainSrc_->localMask(idx)) ++nSrcLocalUnmasked ;
299
300
301  CArray<double,2> boundsLonSrcUnmasked(nVertexSrc,nSrcLocalUnmasked);
302  CArray<double,2> boundsLatSrcUnmasked(nVertexSrc,nSrcLocalUnmasked);
303  long int * globalSrcUnmasked = new long int [nSrcLocalUnmasked];
304
305  nSrcLocalUnmasked=0 ;
306  for (int idx=0 ; idx < nSrcLocal; idx++)
307  {
308    if (domainSrc_->localMask(idx))
309    {
310      for(int n=0;n<nVertexSrc;++n)
311      {
312        boundsLonSrcUnmasked(n,nSrcLocalUnmasked) = boundsLonSrc(n,idx) ;
313        boundsLatSrcUnmasked(n,nSrcLocalUnmasked) = boundsLatSrc(n,idx) ;
314      }
315      globalSrcUnmasked[nSrcLocalUnmasked]=globalSrc[idx] ;
316      ++nSrcLocalUnmasked ;
317    }
318  }
319
320
321  int nDstLocalUnmasked = 0 ;
322  for (int idx=0 ; idx < nDstLocal; idx++) if (domainDest_->localMask(idx)) ++nDstLocalUnmasked ;
323
324  CArray<double,2> boundsLonDestUnmasked(nVertexDest,nDstLocalUnmasked);
325  CArray<double,2> boundsLatDestUnmasked(nVertexDest,nDstLocalUnmasked);
326  long int * globalDstUnmasked = new long int [nDstLocalUnmasked];
327
328  nDstLocalUnmasked=0 ;
329  for (int idx=0 ; idx < nDstLocal; idx++)
330  {
331    if (domainDest_->localMask(idx))
332    {
333      for(int n=0;n<nVertexDest;++n)
334      {
335        boundsLonDestUnmasked(n,nDstLocalUnmasked) = boundsLonDest(n,idx) ;
336        boundsLatDestUnmasked(n,nDstLocalUnmasked) = boundsLatDest(n,idx) ;
337      }
338      globalDstUnmasked[nDstLocalUnmasked]=globalDst[idx] ;
339      ++nDstLocalUnmasked ;
340    }
341  }
342
343  mapper.setSourceMesh(boundsLonSrcUnmasked.dataFirst(), boundsLatSrcUnmasked.dataFirst(), nVertexSrc, nSrcLocalUnmasked, &srcPole[0], globalSrcUnmasked);
344  mapper.setTargetMesh(boundsLonDestUnmasked.dataFirst(), boundsLatDestUnmasked.dataFirst(), nVertexDest, nDstLocalUnmasked, &dstPole[0], globalDstUnmasked);
345
346  std::vector<double> timings = mapper.computeWeights(orderInterp,renormalize,quantity);
347
348  std::map<int,std::vector<std::pair<int,double> > > interpMapValue;
349  std::map<int,std::vector<std::pair<int,double> > >::const_iterator iteNorthPole = interpMapValueNorthPole.end(),
350                                                                     iteSouthPole = interpMapValueSouthPole.end();
351  for (int idx = 0;  idx < mapper.nWeights; ++idx)
352  {
353    interpMapValue[mapper.targetWeightId[idx]].push_back(make_pair(mapper.sourceWeightId[idx],mapper.remapMatrix[idx]));
354    if (iteNorthPole != interpMapValueNorthPole.find(mapper.targetWeightId[idx]))
355    {
356      interpMapValueNorthPole[mapper.targetWeightId[idx]].push_back(make_pair(mapper.sourceWeightId[idx],mapper.remapMatrix[idx]));
357    }
358
359    if (iteSouthPole != interpMapValueSouthPole.find(mapper.targetWeightId[idx]))
360    {
361      interpMapValueSouthPole[mapper.targetWeightId[idx]].push_back(make_pair(mapper.sourceWeightId[idx],mapper.remapMatrix[idx]));
362    }
363  }
364  int niGloDst = domainDest_->ni_glo.getValue();
365  processPole(interpMapValueNorthPole, niGloDst);
366  processPole(interpMapValueSouthPole, niGloDst);
367
368  if (!interpMapValueNorthPole.empty())
369  {
370     std::map<int,std::vector<std::pair<int,double> > >::iterator itNorthPole = interpMapValueNorthPole.begin();
371     for (; itNorthPole != iteNorthPole; ++itNorthPole)
372     {
373       if (!(itNorthPole->second.empty()))
374        itNorthPole->second.swap(interpMapValue[itNorthPole->first]);
375     }
376  }
377
378  if (!interpMapValueSouthPole.empty())
379  {
380     std::map<int,std::vector<std::pair<int,double> > >::iterator itSouthPole = interpMapValueSouthPole.begin();
381     for (; itSouthPole != iteSouthPole; ++itSouthPole)
382     {
383       if (!(itSouthPole->second.empty()))
384        itSouthPole->second.swap(interpMapValue[itSouthPole->first]);
385     }
386  }
387
388  if (writeToFile_ && !readFromFile_) writeRemapInfo(interpMapValue);
389//  exchangeRemapInfo(interpMapValue);
390  convertRemapInfo(interpMapValue) ;
391
392  delete [] globalSrc;
393  delete [] globalSrcUnmasked;
394  delete [] globalDst;
395  delete [] globalDstUnmasked;
396
397}
398
399void CDomainAlgorithmInterpolate::processPole(std::map<int,std::vector<std::pair<int,double> > >& interMapValuePole,
400                                              int nbGlobalPointOnPole)
401{
402  CContext* context = CContext::getCurrent();
403  CContextClient* client=context->client;
404
405  MPI_Comm poleComme(MPI_COMM_NULL);
406  MPI_Comm_split(client->intraComm, interMapValuePole.empty() ? MPI_UNDEFINED : 1, 0, &poleComme);
407  if (MPI_COMM_NULL != poleComme)
408  {
409    int nbClientPole;
410    MPI_Comm_size(poleComme, &nbClientPole);
411
412    std::map<int,std::vector<std::pair<int,double> > >::iterator itePole = interMapValuePole.end(), itPole,
413                                                                 itbPole = interMapValuePole.begin();
414
415    int nbWeight = 0;
416    for (itPole = itbPole; itPole != itePole; ++itPole)
417       nbWeight += itPole->second.size();
418
419    std::vector<int> recvCount(nbClientPole,0);
420    std::vector<int> displ(nbClientPole,0);
421    MPI_Allgather(&nbWeight,1,MPI_INT,&recvCount[0],1,MPI_INT,poleComme) ;
422
423    displ[0]=0;
424    for(int n=1;n<nbClientPole;++n) displ[n]=displ[n-1]+recvCount[n-1] ;
425    int recvSize=displ[nbClientPole-1]+recvCount[nbClientPole-1] ;
426
427    std::vector<int> sendSourceIndexBuff(nbWeight);
428    std::vector<double> sendSourceWeightBuff(nbWeight);
429    int k = 0;
430    for (itPole = itbPole; itPole != itePole; ++itPole)
431    {
432      for (int idx = 0; idx < itPole->second.size(); ++idx)
433      {
434        sendSourceIndexBuff[k] = (itPole->second)[idx].first;
435        sendSourceWeightBuff[k] = (itPole->second)[idx].second;
436        ++k;
437      }
438    }
439
440    std::vector<int> recvSourceIndexBuff(recvSize);
441    std::vector<double> recvSourceWeightBuff(recvSize);
442
443    // Gather all index and weight for pole
444    MPI_Allgatherv(&sendSourceIndexBuff[0],nbWeight,MPI_INT,&recvSourceIndexBuff[0],&recvCount[0],&displ[0],MPI_INT,poleComme);
445    MPI_Allgatherv(&sendSourceWeightBuff[0],nbWeight,MPI_DOUBLE,&recvSourceWeightBuff[0],&recvCount[0],&displ[0],MPI_DOUBLE,poleComme);
446
447    std::map<int,double> recvTemp;
448    for (int idx = 0; idx < recvSize; ++idx)
449    {
450      if (recvTemp.end() != recvTemp.find(recvSourceIndexBuff[idx]))
451        recvTemp[recvSourceIndexBuff[idx]] += recvSourceWeightBuff[idx]/nbGlobalPointOnPole;
452      else
453        recvTemp[recvSourceIndexBuff[idx]] = 0.0;
454    }
455
456    std::map<int,double>::const_iterator itRecvTemp, itbRecvTemp = recvTemp.begin(), iteRecvTemp = recvTemp.end();
457
458    for (itPole = itbPole; itPole != itePole; ++itPole)
459    {
460      itPole->second.clear();
461      for (itRecvTemp = itbRecvTemp; itRecvTemp != iteRecvTemp; ++itRecvTemp)
462          itPole->second.push_back(make_pair(itRecvTemp->first, itRecvTemp->second));
463    }
464  }
465
466}
467
468/*!
469  Compute the index mapping between domain on grid source and one on grid destination
470*/
471void CDomainAlgorithmInterpolate::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
472{
473  if (readFromFile_ && !writeToFile_) 
474    readRemapInfo();
475  else
476  {
477    computeRemap(); 
478  }
479}
480
481void CDomainAlgorithmInterpolate::writeRemapInfo(std::map<int,std::vector<std::pair<int,double> > >& interpMapValue)
482{ 
483  writeInterpolationInfo(fileToReadWrite_, interpMapValue);
484}
485
486void CDomainAlgorithmInterpolate::readRemapInfo()
487{ 
488  std::map<int,std::vector<std::pair<int,double> > > interpMapValue;
489  readInterpolationInfo(fileToReadWrite_, interpMapValue);
490
491  exchangeRemapInfo(interpMapValue);
492}
493
494void CDomainAlgorithmInterpolate::convertRemapInfo(std::map<int,std::vector<std::pair<int,double> > >& interpMapValue)
495{
496  CContext* context = CContext::getCurrent();
497  CContextClient* client=context->client;
498  int clientRank = client->clientRank;
499
500  this->transformationMapping_.resize(1);
501  this->transformationWeight_.resize(1);
502
503  TransformationIndexMap& transMap = this->transformationMapping_[0];
504  TransformationWeightMap& transWeight = this->transformationWeight_[0];
505
506  std::map<int,std::vector<std::pair<int,double> > >::const_iterator itb = interpMapValue.begin(), it,
507                                                                     ite = interpMapValue.end();
508 
509  for (it = itb; it != ite; ++it)
510  {   
511    const std::vector<std::pair<int,double> >& tmp = it->second;
512    for (int i = 0; i < tmp.size(); ++i)
513    {
514      transMap[it->first].push_back(tmp[i].first);
515      transWeight[it->first].push_back(tmp[i].second);
516    }     
517  }     
518}
519
520/*!
521  Read remap information from file then distribute it among clients
522*/
523void CDomainAlgorithmInterpolate::exchangeRemapInfo(std::map<int,std::vector<std::pair<int,double> > >& interpMapValue)
524{
525  CContext* context = CContext::getCurrent();
526  CContextClient* client=context->client;
527  int clientRank = client->clientRank;
528
529  this->transformationMapping_.resize(1);
530  this->transformationWeight_.resize(1);
531
532  TransformationIndexMap& transMap = this->transformationMapping_[0];
533  TransformationWeightMap& transWeight = this->transformationWeight_[0];
534
535  boost::unordered_map<size_t,int> globalIndexOfDomainDest;
536  int ni = domainDest_->ni.getValue();
537  int nj = domainDest_->nj.getValue();
538  int ni_glo = domainDest_->ni_glo.getValue();
539  size_t globalIndex;
540  int nIndexSize = domainDest_->i_index.numElements(), i_ind, j_ind;
541  for (int idx = 0; idx < nIndexSize; ++idx)
542  {
543    i_ind=domainDest_->i_index(idx) ;
544    j_ind=domainDest_->j_index(idx) ;
545
546    globalIndex = i_ind + j_ind * ni_glo;
547    globalIndexOfDomainDest[globalIndex] = clientRank;
548  }
549
550  CClientServerMappingDistributed domainIndexClientClientMapping(globalIndexOfDomainDest,
551                                                                 client->intraComm,
552                                                                 true);
553  CArray<size_t,1> globalIndexInterp(interpMapValue.size());
554  std::map<int,std::vector<std::pair<int,double> > >::const_iterator itb = interpMapValue.begin(), it,
555                                                                     ite = interpMapValue.end();
556  size_t globalIndexCount = 0;
557  for (it = itb; it != ite; ++it)
558  {
559    globalIndexInterp(globalIndexCount) = it->first;
560    ++globalIndexCount;
561  }
562
563  domainIndexClientClientMapping.computeServerIndexMapping(globalIndexInterp);
564  const CClientServerMapping::GlobalIndexMap& globalIndexInterpSendToClient = domainIndexClientClientMapping.getGlobalIndexOnServer();
565
566  //Inform each client number of index they will receive
567  int nbClient = client->clientSize;
568  int* sendBuff = new int[nbClient];
569  int* recvBuff = new int[nbClient];
570  for (int i = 0; i < nbClient; ++i)
571  {
572    sendBuff[i] = 0;
573    recvBuff[i] = 0;
574  }
575  int sendBuffSize = 0;
576  CClientServerMapping::GlobalIndexMap::const_iterator itbMap = globalIndexInterpSendToClient.begin(), itMap,
577                                                       iteMap = globalIndexInterpSendToClient.end();
578  for (itMap = itbMap; itMap != iteMap; ++itMap)
579  {
580    const std::vector<size_t>& tmp = itMap->second;
581    int sizeIndex = 0, mapSize = (itMap->second).size();
582    for (int idx = 0; idx < mapSize; ++idx)
583    {
584//      sizeIndex += interpMapValue.at((itMap->second)[idx]).size();
585      sizeIndex += (interpMapValue[(int)(itMap->second)[idx]]).size();
586    }
587    sendBuff[itMap->first] = sizeIndex;
588    sendBuffSize += sizeIndex;
589  }
590
591
592  MPI_Allreduce(sendBuff, recvBuff, nbClient, MPI_INT, MPI_SUM, client->intraComm);
593
594  int* sendIndexDestBuff = new int [sendBuffSize];
595  int* sendIndexSrcBuff  = new int [sendBuffSize];
596  double* sendWeightBuff = new double [sendBuffSize];
597
598  std::vector<MPI_Request> sendRequest;
599
600  int sendOffSet = 0, l = 0;
601  for (itMap = itbMap; itMap != iteMap; ++itMap)
602  {
603    const std::vector<size_t>& indexToSend = itMap->second;
604    int mapSize = indexToSend.size();
605    int k = 0;
606    for (int idx = 0; idx < mapSize; ++idx)
607    {
608      std::vector<std::pair<int,double> >& interpMap = interpMapValue[(int)indexToSend[idx]]; //interpMapValue.at(indexToSend[idx]);
609      for (int i = 0; i < interpMap.size(); ++i)
610      {
611        sendIndexDestBuff[l] = indexToSend[idx];
612        sendIndexSrcBuff[l]  = interpMap[i].first;
613        sendWeightBuff[l]    = interpMap[i].second;
614        ++k;
615        ++l;
616      }
617    }
618
619    sendRequest.push_back(MPI_Request());
620    MPI_Isend(sendIndexDestBuff + sendOffSet,
621             k,
622             MPI_INT,
623             itMap->first,
624             MPI_DOMAIN_INTERPOLATION_DEST_INDEX,
625             client->intraComm,
626             &sendRequest.back());
627    sendRequest.push_back(MPI_Request());
628    MPI_Isend(sendIndexSrcBuff + sendOffSet,
629             k,
630             MPI_INT,
631             itMap->first,
632             MPI_DOMAIN_INTERPOLATION_SRC_INDEX,
633             client->intraComm,
634             &sendRequest.back());
635    sendRequest.push_back(MPI_Request());
636    MPI_Isend(sendWeightBuff + sendOffSet,
637             k,
638             MPI_DOUBLE,
639             itMap->first,
640             MPI_DOMAIN_INTERPOLATION_WEIGHT,
641             client->intraComm,
642             &sendRequest.back());
643    sendOffSet += k;
644  }
645
646  int recvBuffSize = recvBuff[clientRank];
647  int* recvIndexDestBuff = new int [recvBuffSize];
648  int* recvIndexSrcBuff  = new int [recvBuffSize];
649  double* recvWeightBuff = new double [recvBuffSize];
650  int receivedSize = 0;
651  int clientSrcRank;
652  while (receivedSize < recvBuffSize)
653  {
654    MPI_Status recvStatus;
655    MPI_Recv((recvIndexDestBuff + receivedSize),
656             recvBuffSize,
657             MPI_INT,
658             MPI_ANY_SOURCE,
659             MPI_DOMAIN_INTERPOLATION_DEST_INDEX,
660             client->intraComm,
661             &recvStatus);
662
663    int countBuff = 0;
664    MPI_Get_count(&recvStatus, MPI_INT, &countBuff);
665    clientSrcRank = recvStatus.MPI_SOURCE;
666
667    MPI_Recv((recvIndexSrcBuff + receivedSize),
668             recvBuffSize,
669             MPI_INT,
670             clientSrcRank,
671             MPI_DOMAIN_INTERPOLATION_SRC_INDEX,
672             client->intraComm,
673             &recvStatus);
674
675    MPI_Recv((recvWeightBuff + receivedSize),
676             recvBuffSize,
677             MPI_DOUBLE,
678             clientSrcRank,
679             MPI_DOMAIN_INTERPOLATION_WEIGHT,
680             client->intraComm,
681             &recvStatus);
682
683    for (int idx = 0; idx < countBuff; ++idx)
684    {
685      transMap[*(recvIndexDestBuff + receivedSize + idx)].push_back(*(recvIndexSrcBuff + receivedSize + idx));
686      transWeight[*(recvIndexDestBuff + receivedSize + idx)].push_back(*(recvWeightBuff + receivedSize + idx));
687    }
688    receivedSize += countBuff;
689  }
690
691  std::vector<MPI_Status> requestStatus(sendRequest.size());
692  MPI_Waitall(sendRequest.size(), &sendRequest[0], MPI_STATUS_IGNORE);
693
694  delete [] sendIndexDestBuff;
695  delete [] sendIndexSrcBuff;
696  delete [] sendWeightBuff;
697  delete [] recvIndexDestBuff;
698  delete [] recvIndexSrcBuff;
699  delete [] recvWeightBuff;
700  delete [] sendBuff;
701  delete [] recvBuff;
702}
703 
704/*! Redefined some functions of CONetCDF4 to make use of them */
705CDomainAlgorithmInterpolate::WriteNetCdf::WriteNetCdf(const StdString& filename, const MPI_Comm comm)
706  : CNc4DataOutput(NULL, filename, false, false, true, comm, false, true) {}
707int CDomainAlgorithmInterpolate::WriteNetCdf::addDimensionWrite(const StdString& name, 
708                                                                const StdSize size)
709{
710  return CONetCDF4::addDimension(name, size); 
711}
712
713int CDomainAlgorithmInterpolate::WriteNetCdf::addVariableWrite(const StdString& name, nc_type type,
714                                                               const std::vector<StdString>& dim)
715{
716  return CONetCDF4::addVariable(name, type, dim);
717}
718
719void CDomainAlgorithmInterpolate::WriteNetCdf::endDefinition()
720{
721  CONetCDF4::definition_end();
722}
723
724void CDomainAlgorithmInterpolate::WriteNetCdf::writeDataIndex(const CArray<int,1>& data, const StdString& name,
725                                                              bool collective, StdSize record,
726                                                              const std::vector<StdSize>* start,
727                                                              const std::vector<StdSize>* count)
728{
729  CONetCDF4::writeData<int,1>(data, name, collective, record, start, count);
730}
731
732void CDomainAlgorithmInterpolate::WriteNetCdf::writeDataIndex(const CArray<double,1>& data, const StdString& name,
733                                                              bool collective, StdSize record,
734                                                              const std::vector<StdSize>* start,
735                                                              const std::vector<StdSize>* count)
736{
737  CONetCDF4::writeData<double,1>(data, name, collective, record, start, count);
738}
739
740/*
741   Write interpolation weights into a file
742   \param [in] filename name of output file
743   \param interpMapValue mapping of global index of domain destination and domain source as well as the corresponding weight
744*/
745void CDomainAlgorithmInterpolate::writeInterpolationInfo(std::string& filename,
746                                                         std::map<int,std::vector<std::pair<int,double> > >& interpMapValue)
747{
748  CContext* context = CContext::getCurrent();
749  CContextClient* client=context->client;
750
751  size_t n_src = domainSrc_->ni_glo * domainSrc_->nj_glo;
752  size_t n_dst = domainDest_->ni_glo * domainDest_->nj_glo;
753
754  long localNbWeight = 0;
755  long globalNbWeight;
756  long startIndex;
757  typedef std::map<int,std::vector<std::pair<int,double> > > IndexRemap;
758  IndexRemap::iterator itb = interpMapValue.begin(), it,
759                       ite = interpMapValue.end();
760  for (it = itb; it!=ite; ++it)
761  {
762    localNbWeight += (it->second).size();
763  }
764
765  CArray<int,1> src_idx(localNbWeight);
766  CArray<int,1> dst_idx(localNbWeight);
767  CArray<double,1> weights(localNbWeight);
768
769  int index = 0;
770  for (it = itb; it !=ite; ++it)
771  {
772    std::vector<std::pair<int,double> >& tmp = it->second;
773    for (int idx = 0; idx < tmp.size(); ++idx)
774    {
775      dst_idx(index) = it->first + 1;
776      src_idx(index) = tmp[idx].first + 1;
777      weights(index) = tmp[idx].second;
778      ++index;
779    }   
780  }
781
782  MPI_Allreduce(&localNbWeight, &globalNbWeight, 1, MPI_LONG, MPI_SUM, client->intraComm);
783  MPI_Scan(&localNbWeight, &startIndex, 1, MPI_LONG, MPI_SUM, client->intraComm);
784 
785  if (0 == globalNbWeight)
786  {
787    info << "There is no interpolation weights calculated between "
788         << "domain source: " << domainSrc_->getDomainOutputName()
789         << " and domain destination: " << domainDest_->getDomainOutputName()
790         << std::endl;
791    return;
792  }
793
794  std::vector<StdSize> start(1, startIndex - localNbWeight);
795  std::vector<StdSize> count(1, localNbWeight);
796 
797  WriteNetCdf netCdfWriter(filename, client->intraComm); 
798
799  // Define some dimensions
800  netCdfWriter.addDimensionWrite("n_src", n_src);
801  netCdfWriter.addDimensionWrite("n_dst", n_dst);
802  netCdfWriter.addDimensionWrite("n_weight", globalNbWeight);
803 
804  std::vector<StdString> dims(1,"n_weight");
805
806  // Add some variables
807  netCdfWriter.addVariableWrite("src_idx", NC_INT, dims);
808  netCdfWriter.addVariableWrite("dst_idx", NC_INT, dims);
809  netCdfWriter.addVariableWrite("weight", NC_DOUBLE, dims);
810
811  // End of definition
812  netCdfWriter.endDefinition();
813
814  // // Write variables
815  if (0 != localNbWeight)
816  {
817    netCdfWriter.writeDataIndex(src_idx, "src_idx", false, 0, &start, &count);
818    netCdfWriter.writeDataIndex(dst_idx, "dst_idx", false, 0, &start, &count);
819    netCdfWriter.writeDataIndex(weights, "weight", false, 0, &start, &count);
820  }
821
822  netCdfWriter.closeFile();
823}
824
825/*!
826  Read interpolation information from a file
827  \param [in] filename interpolation file
828  \param [in/out] interpMapValue Mapping between (global) index of domain on grid destination and
829         corresponding global index of domain and associated weight value on grid source
830*/
831void CDomainAlgorithmInterpolate::readInterpolationInfo(std::string& filename,
832                                                        std::map<int,std::vector<std::pair<int,double> > >& interpMapValue)
833{
834  int ncid ;
835  int weightDimId ;
836  size_t nbWeightGlo ;
837
838  CContext* context = CContext::getCurrent();
839  CContextClient* client=context->client;
840  int clientRank = client->clientRank;
841  int clientSize = client->clientSize;
842
843  nc_open(filename.c_str(),NC_NOWRITE, &ncid) ;
844  nc_inq_dimid(ncid,"n_weight",&weightDimId) ;
845  nc_inq_dimlen(ncid,weightDimId,&nbWeightGlo) ;
846
847  size_t nbWeight ;
848  size_t start ;
849  size_t div = nbWeightGlo/clientSize ;
850  size_t mod = nbWeightGlo%clientSize ;
851  if (clientRank < mod)
852  {
853    nbWeight=div+1 ;
854    start=clientRank*(div+1) ;
855  }
856  else
857  {
858    nbWeight=div ;
859    start= mod * (div+1) + (clientRank-mod) * div ;
860  }
861
862  double* weight=new double[nbWeight] ;
863  int weightId ;
864  nc_inq_varid (ncid, "weight", &weightId) ;
865  nc_get_vara_double(ncid, weightId, &start, &nbWeight, weight) ;
866
867  long* srcIndex=new long[nbWeight] ;
868  int srcIndexId ;
869  nc_inq_varid (ncid, "src_idx", &srcIndexId) ;
870  nc_get_vara_long(ncid, srcIndexId, &start, &nbWeight, srcIndex) ;
871
872  long* dstIndex=new long[nbWeight] ;
873  int dstIndexId ;
874  nc_inq_varid (ncid, "dst_idx", &dstIndexId) ;
875  nc_get_vara_long(ncid, dstIndexId, &start, &nbWeight, dstIndex) ;
876
877  for(size_t ind=0; ind<nbWeight;++ind)
878    interpMapValue[dstIndex[ind]-1].push_back(make_pair(srcIndex[ind]-1,weight[ind]));
879}
880
881void CDomainAlgorithmInterpolate::apply(const std::vector<std::pair<int,double> >& localIndex,
882                                            const double* dataInput,
883                                            CArray<double,1>& dataOut,
884                                            std::vector<bool>& flagInitial,
885                                            bool ignoreMissingValue, bool firstPass  )
886{
887  int nbLocalIndex = localIndex.size();   
888  double defaultValue = std::numeric_limits<double>::quiet_NaN();
889   
890  if (detectMissingValue)
891  {
892     if (firstPass && renormalize)
893     {
894       renormalizationFactor.resize(dataOut.numElements()) ;
895       renormalizationFactor=1 ;
896     }
897
898    for (int idx = 0; idx < nbLocalIndex; ++idx)
899    {
900      if (NumTraits<double>::isnan(*(dataInput + idx)))
901      {
902        flagInitial[localIndex[idx].first] = false;
903        if (renormalize) renormalizationFactor(localIndex[idx].first)-=localIndex[idx].second ;
904      }
905      else
906      {
907        dataOut(localIndex[idx].first) += *(dataInput + idx) * localIndex[idx].second;
908        flagInitial[localIndex[idx].first] = true; // Reset flag to indicate not all data source are nan
909      }
910    }
911
912    // If all data source are nan then data destination must be nan
913    for (int idx = 0; idx < nbLocalIndex; ++idx)
914    {
915      if (!flagInitial[localIndex[idx].first])
916        dataOut(localIndex[idx].first) = defaultValue;
917    }
918  }
919  else
920  {
921    for (int idx = 0; idx < nbLocalIndex; ++idx)
922    {
923      dataOut(localIndex[idx].first) += *(dataInput + idx) * localIndex[idx].second;
924    }
925  }
926}
927
928void CDomainAlgorithmInterpolate::updateData(CArray<double,1>& dataOut)
929{
930  if (detectMissingValue && renormalize) dataOut/=renormalizationFactor ;
931}
932
933}
Note: See TracBrowser for help on using the repository browser.