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

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

Smalll improvements for interpolation axis

+) Strictly check order of interpolation
+) Change default order to 1
+) Add precision check

Test
+) On Curie
+) Correct

File size: 14.2 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 "axis.hpp"
11#include "interpolate_axis.hpp"
12#include <algorithm>
13#include "context.hpp"
14#include "context_client.hpp"
15#include "utils.hpp"
16#include "grid.hpp"
17#include "grid_transformation_factory_impl.hpp"
18#include "distribution_client.hpp"
19
20namespace xios {
21CGenericAlgorithmTransformation* CAxisAlgorithmInterpolate::create(CGrid* gridDst, CGrid* gridSrc,
22                                                                   CTransformation<CAxis>* transformation,
23                                                                   int elementPositionInGrid,
24                                                                   std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
25                                                                   std::map<int, int>& elementPositionInGridSrc2AxisPosition,
26                                                                   std::map<int, int>& elementPositionInGridSrc2DomainPosition,
27                                                                   std::map<int, int>& elementPositionInGridDst2ScalarPosition,
28                                                                   std::map<int, int>& elementPositionInGridDst2AxisPosition,
29                                                                   std::map<int, int>& elementPositionInGridDst2DomainPosition)
30{
31  std::vector<CAxis*> axisListDestP = gridDst->getAxis();
32  std::vector<CAxis*> axisListSrcP  = gridSrc->getAxis();
33
34  CInterpolateAxis* interpolateAxis = dynamic_cast<CInterpolateAxis*> (transformation);
35  int axisDstIndex = elementPositionInGridDst2AxisPosition[elementPositionInGrid];
36  int axisSrcIndex = elementPositionInGridSrc2AxisPosition[elementPositionInGrid];
37
38  return (new CAxisAlgorithmInterpolate(axisListDestP[axisDstIndex], axisListSrcP[axisSrcIndex], interpolateAxis));
39}
40
41bool CAxisAlgorithmInterpolate::registerTrans()
42{
43  CGridTransformationFactory<CAxis>::registerTransformation(TRANS_INTERPOLATE_AXIS, create);
44}
45
46
47CAxisAlgorithmInterpolate::CAxisAlgorithmInterpolate(CAxis* axisDestination, CAxis* axisSource, CInterpolateAxis* interpAxis)
48: CAxisAlgorithmTransformation(axisDestination, axisSource), coordinate_(), transPosition_()
49{
50  interpAxis->checkValid(axisSource);
51  order_ = interpAxis->order.getValue();
52  if (!interpAxis->coordinate.isEmpty())
53  {
54    coordinate_ = interpAxis->coordinate.getValue();
55    this->idAuxInputs_.resize(1);
56    this->idAuxInputs_[0] = coordinate_;
57  }
58}
59
60/*!
61  Compute the index mapping between axis on grid source and one on grid destination
62*/
63void CAxisAlgorithmInterpolate::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
64{
65  CContext* context = CContext::getCurrent();
66  CContextClient* client=context->client;
67  int nbClient = client->clientSize;
68  CArray<bool,1>& axisMask = axisSrc_->mask;
69  int srcSize  = axisSrc_->n_glo.getValue();
70  std::vector<CArray<double,1> > vecAxisValue;
71
72  // Fill in axis value from coordinate
73  fillInAxisValue(vecAxisValue, dataAuxInputs);
74  std::vector<double> valueSrc(srcSize);
75  std::vector<double> recvBuff(srcSize);
76  std::vector<int> indexVec(srcSize);
77
78  for (int idx = 0; idx < vecAxisValue.size(); ++idx)
79  {
80    CArray<double,1>& axisValue = vecAxisValue[idx];
81    retrieveAllAxisValue(axisValue, axisMask, recvBuff, indexVec);
82    XIOSAlgorithms::sortWithIndex<double, CVectorStorage>(recvBuff, indexVec);
83    for (int i = 0; i < srcSize; ++i) valueSrc[i] = recvBuff[indexVec[i]];
84    computeInterpolantPoint(valueSrc, indexVec, idx);
85  }
86}
87
88/*!
89  Compute the interpolant points
90  Assume that we have all value of axis source, with these values, need to calculate weight (coeff) of Lagrange polynomial
91  \param [in] axisValue all value of axis source
92  \param [in] tranPos position of axis on a domain
93*/
94void CAxisAlgorithmInterpolate::computeInterpolantPoint(const std::vector<double>& axisValue,
95                                                        const std::vector<int>& indexVec,
96                                                        int transPos)
97{
98  std::vector<double>::const_iterator itb = axisValue.begin(), ite = axisValue.end();
99  std::vector<double>::const_iterator itLowerBound, itUpperBound, it, iteRange, itfirst, itsecond;
100  const double sfmax = NumTraits<double>::sfmax();
101  const double precision = NumTraits<double>::dummy_precision();
102
103  int ibegin = axisDest_->begin.getValue();
104  CArray<double,1>& axisDestValue = axisDest_->value;
105  int numValue = axisDestValue.numElements();
106  std::map<int, std::vector<std::pair<int,double> > > interpolatingIndexValues;
107
108  for (int idx = 0; idx < numValue; ++idx)
109  {
110    bool outOfRange = false;
111    double destValue = axisDestValue(idx);
112    if (destValue < *itb) outOfRange = true;
113
114    itLowerBound = std::lower_bound(itb, ite, destValue);
115    itUpperBound = std::upper_bound(itb, ite, destValue);
116    if ((ite != itUpperBound) && (sfmax == *itUpperBound)) itUpperBound = ite;
117
118    if ((ite == itLowerBound) || (ite == itUpperBound)) outOfRange = true;
119
120    // We don't do extrapolation FOR NOW, maybe in the future
121    if (!outOfRange)
122    {
123      if ((itLowerBound == itUpperBound) && (itb != itLowerBound)) --itLowerBound;
124      double distanceToLower = destValue - *itLowerBound;
125      double distanceToUpper = *itUpperBound - destValue;
126      int order = (order_ + 1) - 2;
127      bool down = (distanceToLower < distanceToUpper) ? true : false;
128      for (int k = 0; k < order; ++k)
129      {
130        if ((itb != itLowerBound) && down)
131        {
132          --itLowerBound;
133          distanceToLower = destValue - *itLowerBound;
134          down = (distanceToLower < distanceToUpper) ? true : false;
135          continue;
136        }
137        if ((ite != itUpperBound) && (sfmax != *itUpperBound))
138        {
139          ++itUpperBound;
140          distanceToUpper = *itUpperBound - destValue;
141          down = (distanceToLower < distanceToUpper) ? true : false;
142
143        }
144      }
145
146      iteRange = (ite == itUpperBound) ? itUpperBound : itUpperBound + 1;
147      itsecond = it = itLowerBound; ++itsecond;
148      while (it < iteRange)
149      {
150        while (((*itsecond -*it) < precision) && itsecond < ite)
151        { ++itsecond; ++it; }
152        int index = std::distance(itb, it);
153        interpolatingIndexValues[idx+ibegin].push_back(make_pair(indexVec[index],*it));
154        ++it; ++itsecond;
155      }
156
157    }
158  }
159  computeWeightedValueAndMapping(interpolatingIndexValues, transPos);
160}
161
162/*!
163  Compute weight (coeff) of Lagrange's polynomial
164  \param [in] interpolatingIndexValues the necessary axis value to calculate the coeffs
165*/
166void CAxisAlgorithmInterpolate::computeWeightedValueAndMapping(const std::map<int, std::vector<std::pair<int,double> > >& interpolatingIndexValues, int transPos)
167{
168  TransformationIndexMap& transMap = this->transformationMapping_[transPos];
169  TransformationWeightMap& transWeight = this->transformationWeight_[transPos];
170  std::map<int, std::vector<std::pair<int,double> > >::const_iterator itb = interpolatingIndexValues.begin(), it,
171                                                                      ite = interpolatingIndexValues.end();
172  int ibegin = axisDest_->begin.getValue();
173  for (it = itb; it != ite; ++it)
174  {
175    int globalIndexDest = it->first;
176    double localValue = axisDest_->value(globalIndexDest - ibegin);
177    const std::vector<std::pair<int,double> >& interpVal = it->second;
178    int interpSize = interpVal.size();
179    transMap[globalIndexDest].resize(interpSize);
180    transWeight[globalIndexDest].resize(interpSize);
181    for (int idx = 0; idx < interpSize; ++idx)
182    {
183      int index = interpVal[idx].first;
184      double weight = 1.0;
185
186      for (int k = 0; k < interpSize; ++k)
187      {
188        if (k == idx) continue;
189        weight *= (localValue - interpVal[k].second);
190        weight /= (interpVal[idx].second - interpVal[k].second);
191      }
192      transMap[globalIndexDest][idx] = index;
193      transWeight[globalIndexDest][idx] = weight;
194      if (!transPosition_.empty())
195      {
196        (this->transformationPosition_[transPos])[globalIndexDest] = transPosition_[transPos];
197      }
198    }
199  }
200  if (!transPosition_.empty() && this->transformationPosition_[transPos].empty())
201    (this->transformationPosition_[transPos])[0] = transPosition_[transPos];
202
203}
204
205/*!
206  Each client retrieves all values of an axis
207  \param [in/out] recvBuff buffer for receiving values (already allocated)
208  \param [in/out] indexVec mapping between values and global index of axis
209*/
210void CAxisAlgorithmInterpolate::retrieveAllAxisValue(const CArray<double,1>& axisValue, const CArray<bool,1>& axisMask,
211                                                     std::vector<double>& recvBuff, std::vector<int>& indexVec)
212{
213  CContext* context = CContext::getCurrent();
214  CContextClient* client=context->client;
215  int nbClient = client->clientSize;
216
217  int srcSize  = axisSrc_->n_glo.getValue();
218  int numValue = axisValue.numElements();
219
220  if (srcSize == numValue)  // Only one client or axis not distributed
221  {
222    for (int idx = 0; idx < srcSize; ++idx)
223    {
224      if (axisMask(idx))
225      {
226        recvBuff[idx] = axisValue(idx);
227        indexVec[idx] = idx;
228      }
229      else
230      {
231        recvBuff[idx] = NumTraits<double>::sfmax();
232        indexVec[idx] = -1;
233      }
234    }
235
236  }
237  else // Axis distributed
238  {
239    double* sendValueBuff = new double [numValue];
240    int* sendIndexBuff = new int [numValue];
241    int* recvIndexBuff = new int [srcSize];
242
243    int ibegin = axisSrc_->begin.getValue();
244    for (int idx = 0; idx < numValue; ++idx)
245    {
246      if (axisMask(idx))
247      {
248        sendValueBuff[idx] = axisValue(idx);
249        sendIndexBuff[idx] = idx + ibegin;
250      }
251      else
252      {
253        sendValueBuff[idx] = NumTraits<double>::sfmax();
254        sendIndexBuff[idx] = -1;
255      }
256    }
257
258    int* recvCount=new int[nbClient];
259    MPI_Allgather(&numValue,1,MPI_INT,recvCount,1,MPI_INT,client->intraComm);
260
261    int* displ=new int[nbClient];
262    displ[0]=0 ;
263    for(int n=1;n<nbClient;n++) displ[n]=displ[n-1]+recvCount[n-1];
264
265    // Each client have enough global info of axis
266    MPI_Allgatherv(sendIndexBuff,numValue,MPI_INT,recvIndexBuff,recvCount,displ,MPI_INT,client->intraComm);
267    MPI_Allgatherv(sendValueBuff,numValue,MPI_DOUBLE,&(recvBuff[0]),recvCount,displ,MPI_DOUBLE,client->intraComm);
268
269    for (int idx = 0; idx < srcSize; ++idx)
270    {
271      indexVec[idx] = recvIndexBuff[idx];
272    }
273
274    delete [] displ;
275    delete [] recvCount;
276    delete [] recvIndexBuff;
277    delete [] sendIndexBuff;
278    delete [] sendValueBuff;
279  }
280}
281
282/*!
283  Fill in axis value dynamically from a field whose grid is composed of a domain and an axis
284  \param [in/out] vecAxisValue vector axis value filled in from input field
285*/
286void CAxisAlgorithmInterpolate::fillInAxisValue(std::vector<CArray<double,1> >& vecAxisValue,
287                                                const std::vector<CArray<double,1>* >& dataAuxInputs)
288{
289  if (coordinate_.empty())
290  {
291    vecAxisValue.resize(1);
292    vecAxisValue[0].resize(axisSrc_->value.numElements());
293    vecAxisValue[0] = axisSrc_->value;
294    this->transformationMapping_.resize(1);
295    this->transformationWeight_.resize(1);
296  }
297  else
298  {
299    CField* field = CField::get(coordinate_);
300    CGrid* grid = field->grid;
301
302    std::vector<CDomain*> domListP = grid->getDomains();
303    std::vector<CAxis*> axisListP = grid->getAxis();
304    if (domListP.empty() || axisListP.empty() || (1 < domListP.size()) || (1 < axisListP.size()))
305      ERROR("CAxisAlgorithmInterpolate::fillInAxisValue(std::vector<CArray<double,1> >& vecAxisValue)",
306             << "XIOS only supports dynamic interpolation with coordinate (field) associated with grid composed of a domain and an axis"
307             << "Coordinate (field) id = " <<field->getId() << std::endl
308             << "Associated grid id = " << grid->getId());
309
310    CDomain* dom = domListP[0];
311    size_t vecAxisValueSize = dom->i_index.numElements();
312    size_t vecAxisValueSizeWithMask = 0;
313    for (size_t idx = 0; idx < vecAxisValueSize; ++idx)
314    {
315      if (dom->mask_1d(idx)) ++vecAxisValueSizeWithMask;
316    }
317
318    int niGlobDom = dom->ni_glo.getValue();
319    vecAxisValue.resize(vecAxisValueSizeWithMask);
320    if (transPosition_.empty())
321    {
322      size_t indexMask = 0;
323      transPosition_.resize(vecAxisValueSizeWithMask);
324      for (size_t idx = 0; idx < vecAxisValueSize; ++idx)
325      {
326        if (dom->mask_1d(idx))
327        {
328          transPosition_[indexMask].resize(1);
329          transPosition_[indexMask][0] = (dom->i_index)(idx) + niGlobDom * (dom->j_index)(idx);
330          ++indexMask;
331        }
332
333      }
334    }
335    this->transformationMapping_.resize(vecAxisValueSizeWithMask);
336    this->transformationWeight_.resize(vecAxisValueSizeWithMask);
337    this->transformationPosition_.resize(vecAxisValueSizeWithMask);
338
339    const CDistributionClient::GlobalLocalDataMap& globalLocalIndexSendToServer = grid->getDistributionClient()->getGlobalLocalDataSendToServer();
340    CDistributionClient::GlobalLocalDataMap::const_iterator itIndex, iteIndex = globalLocalIndexSendToServer.end();
341    size_t axisSrcSize = axisSrc_->index.numElements();
342    std::vector<int> globalDimension = grid->getGlobalDimension();
343
344    size_t indexMask = 0;
345    for (size_t idx = 0; idx < vecAxisValueSize; ++idx)
346    {
347      if (dom->mask_1d(idx))
348      {
349        size_t axisValueSize = 0;
350        for (size_t jdx = 0; jdx < axisSrcSize; ++jdx)
351        {
352          size_t globalIndex = ((dom->i_index)(idx) + (dom->j_index)(idx)*globalDimension[0]) + (axisSrc_->index)(jdx)*globalDimension[0]*globalDimension[1];
353          if (iteIndex != globalLocalIndexSendToServer.find(globalIndex))
354          {
355            ++axisValueSize;
356          }
357        }
358
359        vecAxisValue[indexMask].resize(axisValueSize);
360        axisValueSize = 0;
361        for (size_t jdx = 0; jdx < axisSrcSize; ++jdx)
362        {
363          size_t globalIndex = ((dom->i_index)(idx) + (dom->j_index)(idx)*globalDimension[0]) + (axisSrc_->index)(jdx)*globalDimension[0]*globalDimension[1];
364          itIndex = globalLocalIndexSendToServer.find(globalIndex);
365          if (iteIndex != itIndex)
366          {
367            vecAxisValue[indexMask](axisValueSize) = (*dataAuxInputs[0])(itIndex->second);
368            ++axisValueSize;
369          }
370        }
371        ++indexMask;
372      }
373    }
374  }
375}
376
377}
Note: See TracBrowser for help on using the repository browser.