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

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

Making changes in domain to make sure unstructed grid work with new method of index distribution

+) Change the way define i_index and j_index of a domain
+) Remove some redundant attributes of domain
+) Adjust the way to calculate index distribution on server side

Test
+) Make some minor change to test_unstruct_complete to work with new XIOS
+) On Curie
+) All test pass and correct

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