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

Last change on this file since 1163 was 1163, checked in by ymipsl, 4 years ago

Bug fix when using domain interpolation to a destination domain with overlapping (data indexed domain).

MH. N & YM

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