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

Last change on this file since 2196 was 2196, checked in by jderouillat, 3 years ago

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