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

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

Modifying vertical interpolation

+) Only process interpolation, the extrapolation will be ignored (value will be unidentified (masked))
+) Make sure even unidenfitified from one transformation can be known in the next transformation

Test
+) On Curie
+) Vertical interpolation: Correct
+) Vertical interpolation + horizontal interpolation: Correct

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