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

Last change on this file since 1982 was 1937, checked in by ymipsl, 4 years ago

Fix problem when using area for domain interpolation.

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