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

Last change on this file since 862 was 862, checked in by mhnguyen, 8 years ago

Chaning the way to process transformation to improve the performance.
Instead of exchanging global index and weights on full GRID, each process only
sends and receives the global index and weights on each ELEMENT, which can reduce
the message size of DHT.

+) Domain and axis now have their own exchange function to transfer global index and weight
+) Generic transformation now plays the role of "synthesizer" for all elements
+) Grid transformation now plays the role of transformation mapping, e.x: exchange final global index and weight
among processes.

Test
+) On Curie
+) Pass on all basic tests
+) Dynamic interpolation on axis hasn't been tested (and it seems to need more change to make it rework)

File size: 11.1 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    int niGlobDom = dom->ni_glo.getValue();
268    vecAxisValue.resize(vecAxisValueSize);
269    if (transPosition_.empty())
270    {
271      transPosition_.resize(vecAxisValueSize);
272      for (size_t idx = 0; idx < vecAxisValueSize; ++idx)
273      {
274        transPosition_[idx].resize(1);
275        transPosition_[idx][0] = (dom->i_index)(idx) + niGlobDom * (dom->j_index)(idx);
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 CDistributionClient::GlobalLocalDataMap& globalLocalIndexSendToServer = grid->getDistributionClient()->getGlobalLocalDataSendToServer();
285    CDistributionClient::GlobalLocalDataMap::const_iterator itIndex, iteIndex = globalLocalIndexSendToServer.end();
286    size_t axisSrcSize = axisSrc_->index.numElements();
287    std::vector<int> globalDimension = grid->getGlobalDimension();
288
289    for (size_t idx = 0; idx < vecAxisValueSize; ++idx)
290    {
291      size_t axisValueSize = 0;
292      for (size_t jdx = 0; jdx < axisSrcSize; ++jdx)
293      {
294        size_t globalIndex = ((dom->i_index)(idx) + (dom->j_index)(idx)*globalDimension[0]) + (axisSrc_->index)(jdx)*globalDimension[0]*globalDimension[1];
295        if (iteIndex != globalLocalIndexSendToServer.find(globalIndex))
296        {
297          ++axisValueSize;
298        }
299      }
300
301      vecAxisValue[idx].resize(axisValueSize);
302      axisValueSize = 0;
303      for (size_t jdx = 0; jdx < axisSrcSize; ++jdx)
304      {
305        size_t globalIndex = ((dom->i_index)(idx) + (dom->j_index)(idx)*globalDimension[0]) + (axisSrc_->index)(jdx)*globalDimension[0]*globalDimension[1];
306        itIndex = globalLocalIndexSendToServer.find(globalIndex);
307        if (iteIndex != itIndex)
308        {
309          vecAxisValue[idx](axisValueSize) = (*dataAuxInputs[0])(itIndex->second);
310          ++axisValueSize;
311        }
312      }
313    }
314  }
315}
316
317}
Note: See TracBrowser for help on using the repository browser.