source: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_interpolate_coordinate.cpp @ 2122

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

Improve axis interpolate feature :
introduce : coordinate_src and coordinate_dst feature :

  • level -> pressure
  • pressure -> pressure

YM
  • Property svn:eol-style set to native
  • Property svn:executable set to *
File size: 12.6 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_coordinate.hpp"
10#include "axis.hpp"
11#include "interpolate_axis.hpp"
12#include <algorithm>
13#include <numeric>
14#include "context.hpp"
15#include "context_client.hpp"
16#include "utils.hpp"
17#include "grid.hpp"
18#include "grid_transformation_factory_impl.hpp"
19#include "distribution_client.hpp"
20#include "transform_filter.hpp"
21#include "timer.hpp"
22
23namespace xios
24{
25  CGenericAlgorithmTransformation* CAxisAlgorithmInterpolateCoordinate::create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
26                                                                     CTransformation<CAxis>* transformation,
27                                                                     int elementPositionInGrid,
28                                                                     std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
29                                                                     std::map<int, int>& elementPositionInGridSrc2AxisPosition,
30                                                                     std::map<int, int>& elementPositionInGridSrc2DomainPosition,
31                                                                     std::map<int, int>& elementPositionInGridDst2ScalarPosition,
32                                                                     std::map<int, int>& elementPositionInGridDst2AxisPosition,
33                                                                     std::map<int, int>& elementPositionInGridDst2DomainPosition)
34  TRY
35  {
36    std::vector<CAxis*> axisListDestP = gridDst->getAxis();
37    std::vector<CAxis*> axisListSrcP  = gridSrc->getAxis();
38
39    CInterpolateAxis* interpolateAxis = dynamic_cast<CInterpolateAxis*> (transformation);
40    int axisDstIndex = elementPositionInGridDst2AxisPosition[elementPositionInGrid];
41    int axisSrcIndex = elementPositionInGridSrc2AxisPosition[elementPositionInGrid];
42
43    return (new CAxisAlgorithmInterpolateCoordinate(isSource, axisListDestP[axisDstIndex], axisListSrcP[axisSrcIndex], interpolateAxis));
44  }
45  CATCH
46
47  bool CAxisAlgorithmInterpolateCoordinate::dummyRegistered_ = CAxisAlgorithmInterpolateCoordinate::registerTrans();
48  bool CAxisAlgorithmInterpolateCoordinate::registerTrans()
49  TRY
50  {
51    return CGridTransformationFactory<CAxis>::registerTransformation(TRANS_INTERPOLATE_AXIS, create);
52  }
53  CATCH
54
55  vector<string> CAxisAlgorithmInterpolateCoordinate::getAuxFieldId(void)
56  {
57    if (hasCoordinateSrc_ && hasCoordinateDest_) return {coordinateSrc_,coordinateDest_} ;
58    else if (hasCoordinateSrc_) return {coordinateSrc_} ;
59    else if (hasCoordinateDest_) return {coordinateDest_} ;
60    else return vector<string>() ;
61  }
62 
63  bool CAxisAlgorithmInterpolateCoordinate::transformAuxField(int pos)
64  {
65    if (pos==0)
66    {
67      if (hasCoordinateSrc_) return true ;
68      else if(hasCoordinateDest_) return false ;
69    }
70    if (pos==1) return false ;
71  }
72
73  CAxisAlgorithmInterpolateCoordinate::CAxisAlgorithmInterpolateCoordinate(bool isSource, CAxis* axisDestination, CAxis* axisSource, CInterpolateAxis* interpAxis)
74  : CAlgorithmTransformationTransfer(isSource), axisSrc_(axisSource), axisDest_(axisDestination)
75  TRY
76  {
77    interpAxis->checkValid(axisSource);
78    axisDestination->checkAttributes() ;
79
80    order_ = interpAxis->order.getValue();
81    if (!interpAxis->coordinate.isEmpty())
82    {
83      coordinateSrc_ = interpAxis->coordinate.getValue();
84      hasCoordinate_=true ;
85      hasCoordinateSrc_=true ;
86    }
87
88    if (!interpAxis->coordinate_src.isEmpty())
89    {
90      coordinateSrc_ = interpAxis->coordinate_src.getValue();
91      hasCoordinate_=true ;
92      hasCoordinateSrc_=true ;
93    }
94
95    if (!interpAxis->coordinate_dst.isEmpty())
96    {
97      coordinateDest_ = interpAxis->coordinate_dst.getValue();
98      hasCoordinate_=true ;
99      hasCoordinateDest_=true ;
100    }
101
102    ngloSrc_=axisSource->n_glo ; 
103    nDest_ =  axisDest_-> getLocalView(CElementView::WORKFLOW)->getSize() ;
104   
105    if (!hasCoordinateDest_)
106    {
107      CArray<double,1> coord ;
108      CLocalConnector destConnector(axisDest_->getLocalView(CElementView::FULL), axisDest_->getLocalView(CElementView::WORKFLOW)) ;
109      destConnector.computeConnector() ;
110      destConnector.transfer(axisDest_->value, coord) ;
111      destCoordinate_ = vector<double>(coord.dataFirst(), coord.dataFirst()+nDest_) ;
112    }
113
114    CArray<size_t,1> globalIndex(ngloSrc_) ;
115    for(int i=0;i<ngloSrc_;i++) 
116    {
117      transformationMapping_[i] = i ;
118      globalIndex(i) = i ;
119    }
120   
121
122    CLocalElement axisSourceGlo(CContext::getCurrent()->getIntraCommRank(), ngloSrc_, globalIndex) ;
123    axisSourceGlo.addFullView() ; 
124
125    this->computeAlgorithm(axisSource->getLocalView(CElementView::WORKFLOW), axisSourceGlo.getView(CElementView::FULL)) ;
126
127    if (!hasCoordinateSrc_)
128    {
129      CArray<double,1> coord ;
130      CArray<double,1> coordGlo ;
131      CLocalConnector srcConnector(axisSrc_->getLocalView(CElementView::FULL), axisSrc_->getLocalView(CElementView::WORKFLOW)) ;
132      srcConnector.computeConnector() ;
133      srcConnector.transfer(axisSrc_->value, coord) ; // full view value -> workflow value
134      transferTransformConnector_ -> transfer(coord, coordGlo) ; // workflow view -> full global view
135      srcCoordinate_ = vector<double>(coordGlo.dataFirst(), coordGlo.dataFirst()+ngloSrc_) ;
136    }
137  }
138  CATCH
139
140  CTransformFilter* CAxisAlgorithmInterpolateCoordinate::createTransformFilter(CGarbageCollector& gc, CGridAlgorithm* algo, bool detectMissingValues, double defaultValue)
141  {
142    if (hasCoordinateSrc_ && hasCoordinateDest_) return new CTransformFilter(gc, 3, algo, detectMissingValues, defaultValue) ;
143    else return new CTransformFilter(gc, 2, algo, detectMissingValues, defaultValue) ; 
144  }
145
146  void CAxisAlgorithmInterpolateCoordinate::apply(int dimBefore, int dimAfter, const CArray<double,1>& dataIn, 
147                                                   const vector<CArray<double,1>>& auxDataIn, CArray<double,1>& dataOut)
148  {
149    CArray<double,1> dataInTmp;
150    CArray<double,1> auxDataInSrc ;
151    transferTransformConnector_ -> transfer(dimBefore, dimAfter, dataIn, dataInTmp) ;
152    if (hasCoordinateSrc_)  transferTransformConnector_ -> transfer(dimBefore, dimAfter, auxDataIn[0], auxDataInSrc) ;
153
154   
155    dataOut.resize(dimBefore*dimAfter*nDest_) ;
156    const double* pressureSrc ;
157    const double* pressureDest ;
158    if (hasCoordinateSrc_)  pressureSrc = auxDataInSrc.dataFirst() ;
159    if (hasCoordinateDest_ && hasCoordinateSrc_) pressureDest = auxDataIn[1].dataFirst() ;
160    else if (hasCoordinateDest_ && !hasCoordinateSrc_ ) pressureDest = auxDataIn[0].dataFirst() ;
161   
162    const double* in = dataInTmp.dataFirst() ;
163    double* out = dataOut.dataFirst() ;
164
165    size_t sliceSrc  = dimAfter*ngloSrc_ ;
166    size_t sliceDest = dimAfter*nDest_ ;
167    vector<double> srcCoordinate(ngloSrc_) ;
168    vector<double> destCoordinate(nDest_) ;
169    std::vector<int> srcIndex(ngloSrc_);
170    vector<double> srcValue(ngloSrc_) ;
171    vector<double> destValue(nDest_) ;
172    std::vector<int> destIndex(nDest_);
173
174    size_t nsrc, nDest ;
175
176    for(size_t j=0, posJsrc=0, posJdest=0 ;  j<dimBefore ; j++, posJsrc+=sliceSrc, posJdest+=sliceDest )
177      for(size_t k=0, posKsrc=posJsrc, posKdest=posJdest ; k<dimAfter ; k++, posKsrc++,posKdest++)
178      {
179        if (hasCoordinateSrc_)
180        {
181          nsrc=0 ;
182          for(size_t i=0, posIsrc=posKsrc, posIdest=posKdest ; i<ngloSrc_ ; i++, posIsrc+=dimAfter, posIdest+=dimAfter)
183          {
184            if ( !( std::isnan(pressureSrc[posIsrc]) || std::isnan(in[posIsrc]) ) )
185            {
186              srcCoordinate[nsrc]=pressureSrc[posIsrc] ;
187              srcValue[nsrc] = in[posIsrc] ;
188              nsrc++ ;
189            }
190          }
191        }
192        else
193        {
194          nsrc=0 ;
195          for(size_t i=0, posIsrc=posKsrc ; i<ngloSrc_ ; i++, posIsrc+=dimAfter)
196          {
197            if ( !( std::isnan(srcCoordinate_[i]) || std::isnan(in[posIsrc]) ) )
198            {
199              srcCoordinate[nsrc]=srcCoordinate_[i] ;
200              srcValue[nsrc] = in[posIsrc] ;
201              nsrc++ ;
202            }
203          }
204        } 
205
206        if (hasCoordinateDest_)
207        {
208          nDest=0 ;
209          for(size_t i=0, posIdest=posKdest ; i<nDest_ ; i++, posIdest+=dimAfter)
210          {
211            if ( !( std::isnan(pressureDest[posIdest]) ) )
212            {
213              destCoordinate[nDest]=pressureDest[posIdest] ;
214              nDest++ ;
215            }
216          }
217        }
218        else
219        {
220          nDest=0 ;
221          for(size_t i=0, posIdest=posKdest ; i<nDest_ ; i++, posIdest+=dimAfter)
222          {
223            if ( !( std::isnan(destCoordinate[i]) ) )
224            {
225              destCoordinate[nDest]=destCoordinate_[i] ;
226              nDest++ ;
227            }
228          }
229        }
230 
231        computeInterp(nsrc, srcCoordinate, srcValue, srcIndex, nDest, destCoordinate, destValue, destIndex) ;
232       
233        if (hasCoordinateDest_)
234        {
235          nDest=0 ;
236          for(size_t i=0, posIdest=posKdest ; i<nDest_ ; i++, posIdest+=dimAfter) 
237          {
238            if ( !( std::isnan(pressureDest[posIdest]) ) ) 
239            {
240              out[posIdest] = destValue[nDest] ;
241              nDest++ ;
242            }
243            else out[posIdest] = std::numeric_limits<double>::quiet_NaN() ;
244          }
245        }
246        else
247        {
248          nDest=0 ;
249          for(size_t i=0, posIdest=posKdest ; i<nDest_ ; i++, posIdest+=dimAfter) 
250          {
251            if ( !( std::isnan(destCoordinate[i]) ) ) 
252            {
253              out[posIdest] = destValue[nDest] ;
254              nDest++ ;
255            }
256            else out[posIdest] = std::numeric_limits<double>::quiet_NaN() ;
257          }
258        }
259      }
260  }
261
262  void CAxisAlgorithmInterpolateCoordinate::computeInterp(int nsrc, vector<double>& srcCoordinate, vector<double>& srcValue, vector<int>& srcIndex,
263                                                          int ndst, vector<double>& dstCoordinate, vector<double>& dstValue, vector<int>& dstIndex)
264  {
265    double x,y ;
266    double d ;
267 
268    iota(srcIndex.data(), srcIndex.data()+nsrc, 0); // sort array and retrive sorted index
269    stable_sort(srcIndex.data(), srcIndex.data()+nsrc, [&srcCoordinate](size_t i1, size_t i2) {return srcCoordinate[i1] < srcCoordinate[i2];});
270
271    iota(dstIndex.data(), dstIndex.data()+ndst, 0);
272    stable_sort(dstIndex.data(), dstIndex.data()+ndst, [&dstCoordinate](size_t i1, size_t i2) {return dstCoordinate[i1] < dstCoordinate[i2];});
273
274    if (order_==1 || nsrc<=2)
275    {
276      if (nsrc<=1) dstValue.assign(ndst,std::numeric_limits<double>::quiet_NaN());
277      else
278      {
279
280        double x0,x1 ;
281        double y0,y1 ;
282        int lastj=0;
283     
284        for(int i=0; i < ndst;i++)
285        {
286          x=dstCoordinate[dstIndex[i]] ;
287          if ( x<=srcCoordinate[srcIndex[0]]) lastj=0 ;
288          else if (x>=srcCoordinate[srcIndex[nsrc-1]]) lastj=nsrc-2 ;
289          else
290          {
291            for(int j=lastj; j<nsrc; j++)
292            { 
293              if (x >= srcCoordinate[srcIndex[j]] && x<srcCoordinate[srcIndex[j+1]]) ;
294              lastj=j ;
295              break ;
296            } 
297          }
298          x0=srcCoordinate[srcIndex[lastj]] ;
299          x1=srcCoordinate[srcIndex[lastj+1]] ;
300          y0=srcValue[srcIndex[lastj]] ;
301          y1=srcValue[srcIndex[lastj+1]] ;
302          y=((x-x1)*y0-(x-x0)*y1)/(x0-x1) ;
303          dstValue[dstIndex[i]]=y ;
304        }
305      }
306    }
307    else if (order_==2)
308    {
309      double x0,x1,x2 ;
310      double y0,y1,y2 ;
311      int lastj=0, cj ;
312     
313      for(int i=0; i < ndst;i++)
314      {
315        x=dstCoordinate[dstIndex[i]] ;
316        if ( x<=srcCoordinate[srcIndex[0]]) lastj=0 ;
317        else if (x>=srcCoordinate[srcIndex[nsrc-1]]) lastj=nsrc-2 ;
318        else
319        {
320          for(int j=lastj; j<nsrc; j++)
321          { 
322            if (x >= srcCoordinate[srcIndex[j]] && x<srcCoordinate[srcIndex[j+1]]) ;
323           
324            lastj=j ;
325            break ;
326          } 
327        }
328       
329        if (lastj==0) cj=1 ;
330        else if (lastj>=nsrc-2) cj=nsrc-2 ;
331        else
332        {
333          if ( (x-srcCoordinate[srcIndex[lastj-1]]) > (srcCoordinate[srcIndex[lastj+2]]-x) ) cj=lastj ;
334          else cj=lastj+1 ;
335        } 
336        x0=srcCoordinate[srcIndex[cj-1]] ;
337        x1=srcCoordinate[srcIndex[cj]] ;
338        x2=srcCoordinate[srcIndex[cj+1]] ;
339        y0=srcValue[srcIndex[cj-1]] ;
340        y1=srcValue[srcIndex[cj]] ;
341        y2=srcValue[srcIndex[cj+1]] ;
342           
343        y=y0*(x-x1)*(x-x2)/((x0-x1)*(x0-x2)) + y1*(x-x0)*(x-x2)/((x1-x0)*(x1-x2)) + y2*(x-x0)*(x-x1)/((x2-x0)*(x2-x1))  ;
344        dstValue[dstIndex[i]]=y ;
345      }
346    }
347  } 
348
349}
Note: See TracBrowser for help on using the repository browser.