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

Last change on this file since 1982 was 1980, checked in by yushan, 3 years ago

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