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

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

XIOS coupling branch
Some updates.

First coupling test is beginning to work...

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