source: XIOS/trunk/src/transformation/axis_algorithm_interpolate.cpp @ 827

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

Implementing dynamic interpolation on axis

+) Change grid transformation to make it more flexible
+) Make some small improvements

Test
+) On Curie
+) All test pass

File size: 11.4 KB
Line 
1/*!
2   \file axis_algorithm_interpolate.cpp
3   \author Ha NGUYEN
4   \since 23 June 2015
5   \date 02 Jul 2015
6
7   \brief Algorithm for interpolation on an axis.
8 */
9#include "axis_algorithm_interpolate.hpp"
10#include <algorithm>
11#include "context.hpp"
12#include "context_client.hpp"
13#include "utils.hpp"
14#include "grid.hpp"
15#include "distribution_client.hpp"
16
17namespace xios {
18
19CAxisAlgorithmInterpolate::CAxisAlgorithmInterpolate(CAxis* axisDestination, CAxis* axisSource, CInterpolateAxis* interpAxis)
20: CAxisAlgorithmTransformation(axisDestination, axisSource), coordinate_(), transPosition_()
21{
22  interpAxis->checkValid(axisSource);
23  order_ = interpAxis->order.getValue();
24  if (!interpAxis->coordinate.isEmpty())
25  {
26    coordinate_ = interpAxis->coordinate.getValue();
27    this->idAuxInputs_.resize(1);
28    this->idAuxInputs_[0] = coordinate_;
29  }
30
31//  computeIndexSourceMapping();
32}
33
34/*!
35  Compute the index mapping between axis on grid source and one on grid destination
36*/
37void CAxisAlgorithmInterpolate::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
38{
39  CContext* context = CContext::getCurrent();
40  CContextClient* client=context->client;
41  int nbClient = client->clientSize;
42  CArray<bool,1>& axisMask = axisSrc_->mask;
43  int srcSize  = axisSrc_->n_glo.getValue();
44  std::vector<CArray<double,1> > vecAxisValue;
45
46  // Fill in axis value from coordinate
47  fillInAxisValue(vecAxisValue, dataAuxInputs);
48
49  for (int idx = 0; idx < vecAxisValue.size(); ++idx)
50  {
51    CArray<double,1>& axisValue = vecAxisValue[idx];
52    std::vector<double> recvBuff(srcSize);
53    std::vector<int> indexVec(srcSize);
54    retrieveAllAxisValue(axisValue, axisMask, recvBuff, indexVec);
55    XIOSAlgorithms::sortWithIndex<double, CVectorStorage>(recvBuff, indexVec);
56    computeInterpolantPoint(recvBuff, indexVec, idx);
57  }
58}
59
60/*!
61  Compute the interpolant points
62  Assume that we have all value of axis source, with these values, need to calculate weight (coeff) of Lagrange polynomial
63  \param [in] axisValue all value of axis source
64  \param [in] indexVec permutation index of axisValue
65*/
66void CAxisAlgorithmInterpolate::computeInterpolantPoint(const std::vector<double>& axisValue, const std::vector<int>& indexVec, int transPos)
67{
68  std::vector<double>::const_iterator itb = axisValue.begin(), ite = axisValue.end();
69  std::vector<double>::const_iterator itLowerBound, itUpperBound, it;
70  std::vector<int>::const_iterator itbVec = indexVec.begin(), itVec;
71  const double sfmax = NumTraits<double>::sfmax();
72
73  int ibegin = axisDest_->begin.getValue();
74  CArray<double,1>& axisDestValue = axisDest_->value;
75  int numValue = axisDestValue.numElements();
76  std::map<int, std::vector<std::pair<int,double> > > interpolatingIndexValues;
77
78  for (int idx = 0; idx < numValue; ++idx)
79  {
80    double destValue = axisDestValue(idx);
81    itLowerBound = std::lower_bound(itb, ite, destValue);
82    itUpperBound = std::upper_bound(itb, ite, destValue);
83    if ((ite != itUpperBound) && (sfmax == *itUpperBound)) itUpperBound = ite;
84
85    // If the value is not in the range, that means we'll do extra-polation
86    if (ite == itLowerBound) // extra-polation
87    {
88      itLowerBound = itb;
89      itUpperBound = itb + order_+1;
90    }
91    else if (ite == itUpperBound) // extra-polation
92    {
93      itLowerBound = itUpperBound - order_-1;
94    }
95    else
96    {
97      if (itb != itLowerBound) --itLowerBound;
98      if (ite != itUpperBound) ++itUpperBound;
99      int order = (order_ + 1) - 2;
100      bool down = true;
101      for (int k = 0; k < order; ++k)
102      {
103        if ((itb != itLowerBound) && down)
104        {
105          --itLowerBound;
106          down = false;
107          continue;
108        }
109        if ((ite != itUpperBound) && (sfmax != *itUpperBound))
110        {
111          ++itUpperBound;
112          down = true;
113        }
114      }
115    }
116
117    for (it = itLowerBound; it != itUpperBound; ++it)
118    {
119      int index = std::distance(itb, it);
120      interpolatingIndexValues[idx+ibegin].push_back(make_pair(indexVec[index],*it));
121    }
122  }
123  computeWeightedValueAndMapping(interpolatingIndexValues, transPos);
124}
125
126/*!
127  Compute weight (coeff) of Lagrange's polynomial
128  \param [in] interpolatingIndexValues the necessary axis value to calculate the coeffs
129*/
130void CAxisAlgorithmInterpolate::computeWeightedValueAndMapping(const std::map<int, std::vector<std::pair<int,double> > >& interpolatingIndexValues, int transPos)
131{
132  std::map<int, std::vector<int> >& transMap = this->transformationMapping_[transPos];
133  std::map<int, std::vector<double> >& transWeight = this->transformationWeight_[transPos];
134  std::map<int, std::vector<std::pair<int,double> > >::const_iterator itb = interpolatingIndexValues.begin(), it,
135                                                                      ite = interpolatingIndexValues.end();
136  int ibegin = axisDest_->begin.getValue();
137  for (it = itb; it != ite; ++it)
138  {
139    int globalIndexDest = it->first;
140    double localValue = axisDest_->value(globalIndexDest - ibegin);
141    const std::vector<std::pair<int,double> >& interpVal = it->second;
142    int interpSize = interpVal.size();
143    transMap[globalIndexDest].resize(interpSize);
144    transWeight[globalIndexDest].resize(interpSize);
145    for (int idx = 0; idx < interpSize; ++idx)
146    {
147      int index = interpVal[idx].first;
148      double weight = 1.0;
149
150      for (int k = 0; k < interpSize; ++k)
151      {
152        if (k == idx) continue;
153        weight *= (localValue - interpVal[k].second);
154        weight /= (interpVal[idx].second - interpVal[k].second);
155      }
156      transMap[globalIndexDest][idx] = index;
157      transWeight[globalIndexDest][idx] = weight;
158      if (!transPosition_.empty())
159      {
160        (this->transformationPosition_[transPos])[globalIndexDest] = transPosition_[transPos];
161      }
162    }
163  }
164}
165
166/*!
167  Each client retrieves all values of an axis
168  \param [in/out] recvBuff buffer for receiving values (already allocated)
169  \param [in/out] indexVec mapping between values and global index of axis
170*/
171void CAxisAlgorithmInterpolate::retrieveAllAxisValue(const CArray<double,1>& axisValue, const CArray<bool,1>& axisMask,
172                                                     std::vector<double>& recvBuff, std::vector<int>& indexVec)
173{
174  CContext* context = CContext::getCurrent();
175  CContextClient* client=context->client;
176  int nbClient = client->clientSize;
177
178  int srcSize  = axisSrc_->n_glo.getValue();
179  int numValue = axisValue.numElements();
180
181  if (srcSize == numValue)  // Only one client or axis not distributed
182  {
183    for (int idx = 0; idx < srcSize; ++idx)
184    {
185      if (axisMask(idx))
186      {
187        recvBuff[idx] = axisValue(idx);
188        indexVec[idx] = idx;
189      }
190      else recvBuff[idx] = NumTraits<double>::sfmax();
191    }
192
193  }
194  else // Axis distributed
195  {
196    double* sendValueBuff = new double [numValue];
197    int* sendIndexBuff = new int [numValue];
198    int* recvIndexBuff = new int [srcSize];
199
200    int ibegin = axisSrc_->begin.getValue();
201    for (int idx = 0; idx < numValue; ++idx)
202    {
203      if (axisMask(idx))
204      {
205        sendValueBuff[idx] = axisValue(idx);
206        sendIndexBuff[idx] = idx + ibegin;
207      }
208      else
209      {
210        sendValueBuff[idx] = NumTraits<double>::sfmax();
211        sendIndexBuff[idx] = -1;
212      }
213    }
214
215    int* recvCount=new int[nbClient];
216    MPI_Allgather(&numValue,1,MPI_INT,recvCount,1,MPI_INT,client->intraComm);
217
218    int* displ=new int[nbClient];
219    displ[0]=0 ;
220    for(int n=1;n<nbClient;n++) displ[n]=displ[n-1]+recvCount[n-1];
221
222    // Each client have enough global info of axis
223    MPI_Allgatherv(sendIndexBuff,numValue,MPI_INT,recvIndexBuff,recvCount,displ,MPI_INT,client->intraComm);
224    MPI_Allgatherv(sendValueBuff,numValue,MPI_DOUBLE,&(recvBuff[0]),recvCount,displ,MPI_DOUBLE,client->intraComm);
225
226    for (int idx = 0; idx < srcSize; ++idx)
227    {
228      indexVec[idx] = recvIndexBuff[idx];
229    }
230
231    delete [] displ;
232    delete [] recvCount;
233    delete [] recvIndexBuff;
234    delete [] sendIndexBuff;
235    delete [] sendValueBuff;
236  }
237}
238
239/*!
240  Fill in axis value dynamically from a field whose grid is composed of a domain and an axis
241  \param [in/out] vecAxisValue vector axis value filled in from input field
242*/
243void CAxisAlgorithmInterpolate::fillInAxisValue(std::vector<CArray<double,1> >& vecAxisValue,
244                                                const std::vector<CArray<double,1>* >& dataAuxInputs)
245{
246  if (coordinate_.empty())
247  {
248    vecAxisValue.resize(1);
249    vecAxisValue[0].resize(axisSrc_->value.numElements());
250    vecAxisValue[0] = axisSrc_->value;
251    this->transformationMapping_.resize(1);
252    this->transformationWeight_.resize(1);
253  }
254  else
255  {
256    CField* field = CField::get(coordinate_);
257    CGrid* grid = field->grid;
258
259    std::vector<CDomain*> domListP = grid->getDomains();
260    std::vector<CAxis*> axisListP = grid->getAxis();
261    if (domListP.empty() || axisListP.empty() || (1 < domListP.size()) || (1 < axisListP.size()))
262      ERROR("CAxisAlgorithmInterpolate::fillInAxisValue(std::vector<CArray<double,1> >& vecAxisValue)",
263             << "XIOS only supports dynamic interpolation with coordinate (field) associated with grid composed of a domain and an axis"
264             << "Coordinate (field) id = " <<field->getId() << std::endl
265             << "Associated grid id = " << grid->getId());
266
267    CDomain* dom = domListP[0];
268    size_t vecAxisValueSize = dom->i_index.numElements();
269    vecAxisValue.resize(vecAxisValueSize);
270    if (transPosition_.empty())
271    {
272      transPosition_.resize(vecAxisValueSize);
273      for (size_t idx = 0; idx < vecAxisValueSize; ++idx)
274      {
275        transPosition_[idx].resize(2);
276        transPosition_[idx][0] = (dom->i_index)(idx);
277        transPosition_[idx][1] = (dom->j_index)(idx);
278      }
279    }
280    this->transformationMapping_.resize(vecAxisValueSize);
281    this->transformationWeight_.resize(vecAxisValueSize);
282    this->transformationPosition_.resize(vecAxisValueSize);
283
284    const std::vector<size_t>& globalIndexSendToServer = grid->getDistributionClient()->getGlobalDataIndexSendToServer();
285    const std::vector<int>& localDataSendToServer = grid->getDistributionClient()->getLocalDataIndexSendToServer();
286    std::vector<int> permutIndex(globalIndexSendToServer.size());
287    std::vector<int>::iterator itVec;
288    XIOSAlgorithms::fillInIndex(globalIndexSendToServer.size(), permutIndex);
289    XIOSAlgorithms::sortWithIndex<size_t, CVectorStorage>(globalIndexSendToServer, permutIndex);
290    size_t axisSrcSize = axisSrc_->index.numElements();
291    std::vector<int> globalDimension = grid->getGlobalDimension();
292    XIOSBinarySearchWithIndex<size_t> binSearch(globalIndexSendToServer);
293    for (size_t idx = 0; idx < vecAxisValueSize; ++idx)
294    {
295      size_t axisValueSize = 0;
296      for (size_t jdx = 0; jdx < axisSrcSize; ++jdx)
297      {
298        size_t globalIndex = ((dom->i_index)(idx) + (dom->j_index)(idx)*globalDimension[0]) + (axisSrc_->index)(jdx)*globalDimension[0]*globalDimension[1];
299        if (binSearch.search(permutIndex.begin(), permutIndex.end(), globalIndex, itVec))
300        {
301          ++axisValueSize;
302        }
303      }
304
305      vecAxisValue[idx].resize(axisValueSize);
306      axisValueSize = 0;
307      for (size_t jdx = 0; jdx < axisSrcSize; ++jdx)
308      {
309        size_t globalIndex = ((dom->i_index)(idx) + (dom->j_index)(idx)*globalDimension[0]) + (axisSrc_->index)(jdx)*globalDimension[0]*globalDimension[1];
310        if (binSearch.search(permutIndex.begin(), permutIndex.end(), globalIndex, itVec))
311        {
312          vecAxisValue[idx](axisValueSize) = (*dataAuxInputs[0])(localDataSendToServer[*itVec]);
313          ++axisValueSize;
314        }
315      }
316    }
317  }
318}
319
320}
Note: See TracBrowser for help on using the repository browser.