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

Last change on this file since 982 was 982, checked in by mhnguyen, 5 years ago

Implementing domain interpolation weights exportation
The option 'file' locates file to read if it exists, otherwise a new file with name = 'file' will be created

+) Write interpolation weights into a file

Test
+) On Curie
+) Work

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