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

Last change on this file since 2313 was 2313, checked in by ymipsl, 2 years ago

Take into account detect_missing_value and renormalize attribute for interpolate class transformation.

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