source: XIOS/dev/dev_trunk_omp/src/transformation/domain_algorithm_interpolate.cpp @ 1770

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

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