source: XIOS/dev/XIOS_DEV_CMIP6/src/transformation/domain_algorithm_interpolate.cpp @ 1480

Last change on this file since 1480 was 1480, checked in by ymipsl, 6 years ago

Fix for transformation & interpolation : better detection of missing values.

YM

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