source: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm_interpolate.cpp @ 1787

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

More cleaning : replace contextClient rank and size by context rank and size.
Because for couling it will have more than one contextClient.

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