source: XIOS2/trunk/src/transformation/domain_algorithm_interpolate.cpp @ 2442

Last change on this file since 2442 was 2442, checked in by jderouillat, 19 months ago

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