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

Last change on this file since 1638 was 1638, checked in by yushan, 5 years ago

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