source: XIOS/dev/dev_trunk_omp/src/transformation/domain_algorithm_interpolate.cpp @ 1619

Last change on this file since 1619 was 1619, checked in by yushan, 4 years ago

branch_openmp merged and NOT tested with trunk r1618

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