source: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_interpolate.cpp @ 2009

Last change on this file since 2009 was 2009, checked in by ymipsl, 3 years ago

Implement axis interpolation, without pressure coordinate for now.

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