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

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

Update to new transformation.
Source filter working again + reading
YM

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