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

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

Improvements for dht

+) Implement adaptive hierarchy for dht, level of hierarchy depends on number of processes
+) Remove some redundant codes

Test
+) On Curie
+) All tests are correct

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