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

Last change on this file since 896 was 896, checked in by mhnguyen, 6 years ago

Correcting a bug in axis interpolation

+) Ascending or Descending values of source axis wont affect on the results

Test
+) On Curie
+) Test is correct

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