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

Last change on this file since 1004 was 1004, checked in by mhnguyen, 7 years ago

Some modifications of reading/writing interpolation weight from/into a file
There are 3 attributes concerning read/write weights from/into a file

  • write_weight: true/false. Default value: false
  • weight_filename.
  • mode: Default value: compute + read: read from file whose name is weight_filename + compute: compute weights + read_or_compute: if weight_filename already exists, the weights will be read from it; otherwise the weights will be computed

In mode compute and read_or_compute, if there is no weight_filename, a filename whose format
xios_interpolation_weight_nameOfContext_nameOfDomainSource_nameOfDomainDestination.nc will be used for read/write.

Test
+) On Curie
+) Work

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