source: XIOS/trunk/src/transformation/domain_algorithm_interpolate_from_file.cpp @ 660

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

Adding interpolation test_remap

+) Add new test case for domain interpolation
+) Resolve circular dependence
+) Add function to read weigh file

Test
+) On Curie
+) test_remap can print out .nc

File size: 9.3 KB
Line 
1/*!
2   \file domain_algorithm_interpolate_from_file.cpp
3   \author Ha NGUYEN
4   \since 09 Jul 2015
5   \date 17 Jul 2015
6
7   \brief Algorithm for interpolation on a domain.
8 */
9#include "domain_algorithm_interpolate_from_file.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
17namespace xios {
18
19CDomainAlgorithmInterpolateFromFile::CDomainAlgorithmInterpolateFromFile(CDomain* domainDestination, CDomain* domainSource, CInterpolateFromFileDomain* interpDomain)
20: CDomainAlgorithmTransformation(domainDestination, domainSource), interpDomain_(interpDomain)
21{
22  interpDomain_->checkValid(domainSource);
23  computeIndexSourceMapping();
24}
25
26/*!
27  Compute the index mapping between domain on grid source and one on grid destination
28*/
29void CDomainAlgorithmInterpolateFromFile::computeIndexSourceMapping()
30{
31  CContext* context = CContext::getCurrent();
32  CContextClient* client=context->client;
33  int clientRank = client->clientRank;
34
35  std::string filename = interpDomain_->file.getValue();
36  std::map<int,std::vector<std::pair<int,double> > > interpMapValue;
37  readInterpolationInfo(filename, interpMapValue);
38
39  //randomizeInterpolationInfo(interpMapValue);
40  boost::unordered_map<size_t,int> globalIndexOfDomainDest;
41  int ni = domainDest_->ni.getValue();
42  int nj = domainDest_->nj.getValue();
43  int ni_glo = domainDest_->ni_glo.getValue();
44  size_t globalIndex;
45  int nIndexSize = domainDest_->i_index.numElements(), i_ind, j_ind;
46  for (int idx = 0; idx < nIndexSize; ++idx)
47  {
48    i_ind=domainDest_->i_index(idx) ;
49    j_ind=domainDest_->j_index(idx) ;
50
51    globalIndex = i_ind + j_ind * ni_glo;
52    globalIndexOfDomainDest[globalIndex] = clientRank;
53  }
54
55  CClientServerMappingDistributed domainIndexClientClientMapping(globalIndexOfDomainDest,
56                                                                 client->intraComm,
57                                                                 true);
58  CArray<size_t,1> globalIndexInterp(interpMapValue.size());
59  std::map<int,std::vector<std::pair<int,double> > >::const_iterator itb = interpMapValue.begin(), it,
60                                                                     ite = interpMapValue.end();
61  size_t globalIndexCount = 0;
62  for (it = itb; it != ite; ++it)
63  {
64    globalIndexInterp(globalIndexCount) = it->first;
65    ++globalIndexCount;
66  }
67
68  domainIndexClientClientMapping.computeServerIndexMapping(globalIndexInterp);
69  const std::map<int, std::vector<size_t> >& globalIndexInterpSendToClient = domainIndexClientClientMapping.getGlobalIndexOnServer();
70
71  //Inform each client number of index they will receive
72  int nbClient = client->clientSize;
73  int* sendBuff = new int[nbClient];
74  int* recvBuff = new int[nbClient];
75  for (int i = 0; i < nbClient; ++i) sendBuff[i] = 0;
76  int sendBuffSize = 0;
77  std::map<int, std::vector<size_t> >::const_iterator itbMap = globalIndexInterpSendToClient.begin(), itMap,
78                                                      iteMap = globalIndexInterpSendToClient.end();
79  for (itMap = itbMap; itMap != iteMap; ++itMap)
80  {
81    int sizeIndex = 0, mapSize = (itMap->second).size();
82    for (int idx = 0; idx < mapSize; ++idx)
83    {
84      sizeIndex += interpMapValue[(itMap->second)[idx]].size();
85    }
86    sendBuff[itMap->first] = sizeIndex;
87    sendBuffSize += sizeIndex;
88  }
89
90  MPI_Allreduce(sendBuff, recvBuff, nbClient, MPI_INT, MPI_SUM, client->intraComm);
91
92  int* sendIndexDestBuff = new int [sendBuffSize];
93  int* sendIndexSrcBuff  = new int [sendBuffSize];
94  double* sendWeightBuff = new double [sendBuffSize];
95
96  std::vector<MPI_Request> sendRequest;
97  // Now send index and weight
98  int sendOffSet = 0;
99  for (itMap = itbMap; itMap != iteMap; ++itMap)
100  {
101    int k = 0;
102    int mapSize = (itMap->second).size();
103
104    for (int idx = 0; idx < mapSize; ++idx)
105    {
106      std::vector<std::pair<int,double> >& interpMap = interpMapValue[(itMap->second)[idx]];
107      for (int i = 0; i < interpMap.size(); ++i)
108      {
109        sendIndexDestBuff[k] = (itMap->second)[idx];
110        sendIndexSrcBuff[k]  = interpMap[i].first;
111        sendWeightBuff[k]    = interpMap[i].second;
112        ++k;
113      }
114    }
115
116    sendRequest.push_back(MPI_Request());
117    MPI_Isend(sendIndexDestBuff + sendOffSet,
118             k,
119             MPI_INT,
120             itMap->first,
121             7,
122             client->intraComm,
123             &sendRequest.back());
124    sendRequest.push_back(MPI_Request());
125    MPI_Isend(sendIndexSrcBuff + sendOffSet,
126             k,
127             MPI_INT,
128             itMap->first,
129             8,
130             client->intraComm,
131             &sendRequest.back());
132    sendRequest.push_back(MPI_Request());
133    MPI_Isend(sendWeightBuff + sendOffSet,
134             k,
135             MPI_DOUBLE,
136             itMap->first,
137             9,
138             client->intraComm,
139             &sendRequest.back());
140    sendOffSet += k;
141  }
142
143  int recvBuffSize = recvBuff[clientRank];
144  int* recvIndexDestBuff = new int [recvBuffSize];
145  int* recvIndexSrcBuff  = new int [recvBuffSize];
146  double* recvWeightBuff = new double [recvBuffSize];
147  int receivedSize = 0;
148  int clientSrcRank;
149  while (receivedSize < recvBuffSize)
150  {
151    MPI_Status recvStatus;
152    MPI_Recv((recvIndexDestBuff + receivedSize),
153             recvBuffSize,
154             MPI_INT,
155             MPI_ANY_SOURCE,
156             7,
157             client->intraComm,
158             &recvStatus);
159
160    int countBuff = 0;
161    MPI_Get_count(&recvStatus, MPI_INT, &countBuff);
162    clientSrcRank = recvStatus.MPI_SOURCE;
163
164    MPI_Recv((recvIndexSrcBuff + receivedSize),
165             recvBuffSize,
166             MPI_INT,
167             clientSrcRank,
168             8,
169             client->intraComm,
170             &recvStatus);
171
172    MPI_Recv((recvWeightBuff + receivedSize),
173             recvBuffSize,
174             MPI_DOUBLE,
175             clientSrcRank,
176             9,
177             client->intraComm,
178             &recvStatus);
179
180    for (int idx = 0; idx < countBuff; ++idx)
181    {
182      transformationMapping_[*(recvIndexDestBuff + receivedSize + idx)].push_back(*(recvIndexSrcBuff + receivedSize + idx));
183      transformationWeight_[*(recvIndexDestBuff + receivedSize + idx)].push_back(*(recvWeightBuff + receivedSize + idx));
184    }
185    receivedSize += countBuff;
186  }
187
188  std::vector<MPI_Status> requestStatus(sendRequest.size());
189  MPI_Wait(&sendRequest[0], &requestStatus[0]);
190
191  delete [] sendIndexDestBuff;
192  delete [] sendIndexSrcBuff;
193  delete [] sendWeightBuff;
194  delete [] recvIndexDestBuff;
195  delete [] recvIndexSrcBuff;
196  delete [] recvWeightBuff;
197  delete [] sendBuff;
198  delete [] recvBuff;
199}
200
201void CDomainAlgorithmInterpolateFromFile::randomizeInterpolationInfo(std::map<int,std::vector<std::pair<int,double> > >& interpMapValue)
202{
203  int iDestBegin = 2;
204  int jDestBegin = 2;
205
206  int ni_glo_src = domainSrc_->ni_glo.getValue();
207  int ni_glo = domainDest_->ni_glo.getValue();
208  size_t globalIndex;
209  int nIndexSize = domainDest_->i_index.numElements(), i_ind, j_ind;
210  for (int idx = 0; idx < nIndexSize; ++idx)
211  {
212    i_ind=domainDest_->i_index(idx) ;
213    j_ind=domainDest_->j_index(idx) ;
214
215    globalIndex = i_ind + j_ind * ni_glo;
216
217    interpMapValue[globalIndex].push_back(make_pair((i_ind+iDestBegin) + (j_ind+jDestBegin-1)*ni_glo_src, 0.25));
218    interpMapValue[globalIndex].push_back(make_pair((i_ind+iDestBegin+1) + (j_ind+jDestBegin)*ni_glo_src, 0.25));
219    interpMapValue[globalIndex].push_back(make_pair((i_ind+iDestBegin) + (j_ind+jDestBegin+1)*ni_glo_src, 0.25));
220    interpMapValue[globalIndex].push_back(make_pair((i_ind+iDestBegin-1) + (j_ind+jDestBegin)*ni_glo_src, 0.25));
221
222
223  }
224}
225
226/*!
227  Read interpolation information from a file
228  \param [in] filename interpolation file
229  \param [in/out] interpMapValue Mapping between (global) index of domain on grid destination and
230         corresponding global index of domain and associated weight value on grid source
231*/
232void CDomainAlgorithmInterpolateFromFile::readInterpolationInfo(std::string& filename,
233                                                                std::map<int,std::vector<std::pair<int,double> > >& interpMapValue)
234{
235  int ncid ;
236  int weightDimId ;
237  size_t nbWeightGlo ;
238
239  CContext* context = CContext::getCurrent();
240  CContextClient* client=context->client;
241  int clientRank = client->clientRank;
242  int clientSize = client->clientSize;
243 
244  nc_open(filename.c_str(),NC_NOWRITE, &ncid) ;
245  nc_inq_dimid(ncid,"n_weight",&weightDimId) ;
246  nc_inq_dimlen(ncid,weightDimId,&nbWeightGlo) ;
247
248  size_t nbWeight ;
249  size_t start ;
250  size_t div = nbWeightGlo/clientSize ;
251  size_t mod = nbWeightGlo%clientSize ;
252  if (clientRank < mod)
253  {
254    nbWeight=div+1 ;
255    start=clientRank*(div+1) ;
256  }
257  else
258  {
259    nbWeight=div ;
260    start= mod * (div+1) + (clientRank-mod) * div ;
261  }
262 
263 
264
265  double* weight=new double[nbWeight] ;
266  int weightId ;
267  nc_inq_varid (ncid, "weight", &weightId) ;
268  nc_get_vara_double(ncid, weightId, &start, &nbWeight, weight) ;
269
270  long* srcIndex=new long[nbWeight] ; 
271  int srcIndexId ;
272  nc_inq_varid (ncid, "src_idx", &srcIndexId) ;
273  nc_get_vara_long(ncid, srcIndexId, &start, &nbWeight, srcIndex) ;
274
275  long* dstIndex=new long[nbWeight] ; 
276  int dstIndexId ;
277  nc_inq_varid (ncid, "dst_idx", &dstIndexId) ;
278  nc_get_vara_long(ncid, dstIndexId, &start, &nbWeight, dstIndex) ;
279       
280  for(size_t ind=0; ind<nbWeight;++ind)
281    interpMapValue[dstIndex[ind]].push_back(make_pair(srcIndex[ind],weight[ind]));
282}
283
284}
Note: See TracBrowser for help on using the repository browser.