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

Last change on this file since 2004 was 2002, checked in by ymipsl, 3 years ago

Some cleaning of old transformation dead code

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