source: XIOS/dev/XIOS_DEV_CMIP6/src/transformation/axis_algorithm_interpolate.cpp @ 1412

Last change on this file since 1412 was 1412, checked in by ymipsl, 4 years ago

More timer instrumentation...

YM

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