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

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

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