Changeset 2122


Ignore:
Timestamp:
04/01/21 11:13:01 (3 years ago)
Author:
ymipsl
Message:

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

  • level -> pressure
  • pressure -> pressure

YM

Location:
XIOS/dev/dev_ym/XIOS_COUPLING/src
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • XIOS/dev/dev_ym/XIOS_COUPLING/src/config/interpolate_axis_attribute.conf

    r827 r2122  
    33DECLARE_ATTRIBUTE(int, order) 
    44DECLARE_ATTRIBUTE(StdString, coordinate) 
     5DECLARE_ATTRIBUTE(StdString, coordinate_src) 
     6DECLARE_ATTRIBUTE(StdString, coordinate_dst) 

  • XIOS/dev/dev_ym/XIOS_COUPLING/src/node/interpolate_axis.cpp

    r2011 r2122  
    6060    } 
    6161 
     62    if (!this->coordinate.isEmpty() && (!this->coordinate_src.isEmpty() || !this->coordinate_dst.isEmpty())) 
     63    { 
     64       ERROR("CInterpolateAxis::checkValid(CAxis* axisSrc)", 
     65               << "<coordinate> is incompatible with <coordinate_src> or <coordinate_dst> attributes. <coordinate> attribute is present only for backward" 
     66               << "compatibility and is now declared obsolete") ; 
     67    } 
    6268 
    6369    if (!this->coordinate.isEmpty()) 
    6470    { 
    65       StdString coordinate = this->coordinate.getValue(); 
    6671      if (!CField::has(coordinate)) 
    6772        ERROR("CInterpolateAxis::checkValid(CAxis* axisSrc)", 
    6873               << "Coordinate field whose id " << coordinate << "does not exist " 
     74               << "Please define one"); 
     75    } 
     76 
     77    if (!this->coordinate_src.isEmpty()) 
     78    { 
     79      if (!CField::has(coordinate_src)) 
     80        ERROR("CInterpolateAxis::checkValid(CAxis* axisSrc)", 
     81               << "Coordinate field whose id " << coordinate_src << "does not exist " 
     82               << "Please define one"); 
     83    } 
     84 
     85    if (!this->coordinate_dst.isEmpty()) 
     86    { 
     87      if (!CField::has(coordinate_dst)) 
     88        ERROR("CInterpolateAxis::checkValid(CAxis* axisSrc)", 
     89               << "Coordinate field whose id " << coordinate_dst << "does not exist " 
    6990               << "Please define one"); 
    7091    } 
    … …  
    7495  { 
    7596    std::vector<StdString> auxInputs; 
     97     
     98    if (!this->coordinate.isEmpty() && (!this->coordinate_src.isEmpty() || !this->coordinate_dst.isEmpty())) 
     99    { 
     100       ERROR("std::vector<StdString> CInterpolateAxis::checkAuxInputs_()", 
     101               << "<coordinate> is incompatible with <coordinate_src> or <coordinate_dst> attributes. <coordinate> attribute is present only for backward" 
     102               << "compatibility and is now declared obsolete") ; 
     103    } 
     104 
    76105    if (!this->coordinate.isEmpty()) 
    77106    { 
    78107      StdString coordinate = this->coordinate.getValue(); 
    79108      if (!CField::has(coordinate)) 
    80         ERROR("CInterpolateAxis::checkValid(CAxis* axisSrc)", 
     109        ERROR("std::vector<StdString> CInterpolateAxis::checkAuxInputs_()", 
     110               << "Coordinate field whose id " << coordinate << "does not exist " 
     111               << "Please define one"); 
     112      auxInputs.push_back(coordinate); 
     113    } 
     114 
     115    if (!this->coordinate_src.isEmpty()) 
     116    { 
     117      StdString coordinate = this->coordinate_src.getValue(); 
     118      if (!CField::has(coordinate)) 
     119        ERROR("std::vector<StdString> CInterpolateAxis::checkAuxInputs_()", 
     120               << "Coordinate field whose id " << coordinate << "does not exist " 
     121               << "Please define one"); 
     122      auxInputs.push_back(coordinate); 
     123    } 
     124 
     125    if (!this->coordinate_dst.isEmpty()) 
     126    { 
     127      StdString coordinate = this->coordinate_dst.getValue(); 
     128      if (!CField::has(coordinate)) 
     129        ERROR("std::vector<StdString> CInterpolateAxis::checkAuxInputs_()", 
    81130               << "Coordinate field whose id " << coordinate << "does not exist " 
    82131               << "Please define one"); 
    … …  
    97146                                                        std::map<int, int>& elementPositionInGridDst2DomainPosition) 
    98147  { 
    99     if (!coordinate.isEmpty())  return CAxisAlgorithmInterpolateCoordinate::create(isSource, gridDst,  gridSrc, this, elementPositionInGrid,  
    100                                       elementPositionInGridSrc2ScalarPosition, elementPositionInGridSrc2AxisPosition, elementPositionInGridSrc2DomainPosition, 
    101                                       elementPositionInGridDst2ScalarPosition, elementPositionInGridDst2AxisPosition, elementPositionInGridDst2DomainPosition); 
     148    if (!coordinate.isEmpty() || !coordinate_src.isEmpty() || !coordinate_dst.isEmpty())  
     149       return CAxisAlgorithmInterpolateCoordinate::create(isSource, gridDst,  gridSrc, this, elementPositionInGrid,  
     150              elementPositionInGridSrc2ScalarPosition, elementPositionInGridSrc2AxisPosition, elementPositionInGridSrc2DomainPosition, 
     151              elementPositionInGridDst2ScalarPosition, elementPositionInGridDst2AxisPosition, elementPositionInGridDst2DomainPosition); 
     152     
    102153    else return CAxisAlgorithmInterpolate::create(isSource, gridDst,  gridSrc, this, elementPositionInGrid,  
    103154                elementPositionInGridSrc2ScalarPosition, elementPositionInGridSrc2AxisPosition, elementPositionInGridSrc2DomainPosition, 

  • XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_interpolate_coordinate.cpp

    r2011 r2122  
    5555  vector<string> CAxisAlgorithmInterpolateCoordinate::getAuxFieldId(void) 
    5656  { 
    57     if (hasCoordinate_) return vector<string>(1,coordinate_) ; 
     57    if (hasCoordinateSrc_ && hasCoordinateDest_) return {coordinateSrc_,coordinateDest_} ; 
     58    else if (hasCoordinateSrc_) return {coordinateSrc_} ; 
     59    else if (hasCoordinateDest_) return {coordinateDest_} ; 
    5860    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 ; 
    5971  } 
    6072 
    … …  
    6981    if (!interpAxis->coordinate.isEmpty()) 
    7082    { 
    71       coordinate_ = interpAxis->coordinate.getValue(); 
     83      coordinateSrc_ = interpAxis->coordinate.getValue(); 
    7284      hasCoordinate_=true ; 
    73     } 
    74    
    75     ngloSrc_=axisSource->n_glo ; 
     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 ;  
    76103    nDest_ =  axisDest_-> getLocalView(CElementView::WORKFLOW)->getSize() ; 
    77     CArray<double,1> coord ; 
    78     CLocalConnector destConnector(axisDest_->getLocalView(CElementView::FULL), axisDest_->getLocalView(CElementView::WORKFLOW)) ; 
    79     destConnector.computeConnector() ; 
    80     destConnector.transfer(axisDest_->value, coord) ; 
    81     destCoordinate_ = vector<double>(coord.dataFirst(), coord.dataFirst()+nDest_) ; 
    82104     
     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 
    83114    CArray<size_t,1> globalIndex(ngloSrc_) ; 
    84115    for(int i=0;i<ngloSrc_;i++)  
    … …  
    93124 
    94125    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    } 
    95137  } 
    96138  CATCH 
    … …  
    98140  CTransformFilter* CAxisAlgorithmInterpolateCoordinate::createTransformFilter(CGarbageCollector& gc, CGridAlgorithm* algo, bool detectMissingValues, double defaultValue) 
    99141  { 
    100     return new CTransformFilter(gc, 2, algo, detectMissingValues, defaultValue) ;   
     142    if (hasCoordinateSrc_ && hasCoordinateDest_) return new CTransformFilter(gc, 3, algo, detectMissingValues, defaultValue) ; 
     143    else return new CTransformFilter(gc, 2, algo, detectMissingValues, defaultValue) ;   
    101144  } 
    102145 
    … …  
    105148  { 
    106149    CArray<double,1> dataInTmp; 
    107     CArray<double,1> auxDataInTmp ; 
     150    CArray<double,1> auxDataInSrc ; 
    108151    transferTransformConnector_ -> transfer(dimBefore, dimAfter, dataIn, dataInTmp) ; 
    109     transferTransformConnector_ -> transfer(dimBefore, dimAfter, auxDataIn[0], auxDataInTmp) ; 
    110    
     152    if (hasCoordinateSrc_)  transferTransformConnector_ -> transfer(dimBefore, dimAfter, auxDataIn[0], auxDataInSrc) ; 
     153 
     154    
    111155    dataOut.resize(dimBefore*dimAfter*nDest_) ; 
    112     const double* pressure = auxDataInTmp.dataFirst() ; 
     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     
    113162    const double* in = dataInTmp.dataFirst() ; 
    114163    double* out = dataOut.dataFirst() ; 
    115164 
    116     size_t sliceSrc = dimAfter*ngloSrc_ ; 
     165    size_t sliceSrc  = dimAfter*ngloSrc_ ; 
    117166    size_t sliceDest = dimAfter*nDest_ ; 
    118167    vector<double> srcCoordinate(ngloSrc_) ; 
    … …  
    123172    std::vector<int> destIndex(nDest_); 
    124173 
    125     size_t nsrc ; 
     174    size_t nsrc, nDest ; 
     175 
    126176    for(size_t j=0, posJsrc=0, posJdest=0 ;  j<dimBefore ; j++, posJsrc+=sliceSrc, posJdest+=sliceDest ) 
    127177      for(size_t k=0, posKsrc=posJsrc, posKdest=posJdest ; k<dimAfter ; k++, posKsrc++,posKdest++) 
    128178      { 
    129         nsrc=0 ; 
    130         for(size_t i=0, posIsrc=posKsrc, posIdest=posKdest ; i<ngloSrc_ ; i++, posIsrc+=dimAfter, posIdest+=dimAfter) 
    131         { 
    132           if ( !( std::isnan(pressure[posIsrc]) || std::isnan(in[posIsrc]) ) ) 
    133           { 
    134             srcCoordinate[nsrc]=pressure[posIsrc] ; 
    135             srcValue[nsrc] = in[posIsrc] ; 
    136             nsrc++ ; 
    137           } 
    138         } 
    139         destCoordinate = destCoordinate_ ;     
    140         computeInterp(nsrc, srcCoordinate, srcValue, srcIndex, nDest_, destCoordinate, destValue, destIndex) ; 
    141  
    142         for(size_t i=0, posIdest=posKdest ; i<nDest_ ; i++, posIdest+=dimAfter)  out[posIdest] = destValue[i] ; 
     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        } 
    143259      } 
    144  
    145260  } 
    146261 

  • XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_interpolate_coordinate.hpp

    r2011 r2122  
    3232  virtual ~CAxisAlgorithmInterpolateCoordinate() {} 
    3333  virtual vector<string> getAuxFieldId(void)  ; 
     34  virtual bool transformAuxField(int pos) ; 
    3435  virtual void apply(int dimBefore, int dimAfter, const CArray<double,1>& dataIn,  
    3536                     const vector<CArray<double,1>>& auxDataIn, CArray<double,1>& dataOut) ; 
    … …  
    4344  // Interpolation order 
    4445  int order_; 
    45   StdString coordinate_; 
     46  StdString coordinateSrc_; // pressure src 
     47  StdString coordinateDest_; // pressure dst 
    4648  bool hasCoordinate_=false ; 
     49  bool hasCoordinateSrc_=false ; 
     50  bool hasCoordinateDest_=false ; 
    4751 
    4852  CAxis* axisSrc_=nullptr ; 
    … …  
    5054  size_t ngloSrc_ ; 
    5155  size_t nDest_ ; 
    52   vector<double> destCoordinate_ ; 
     56  vector<double> srcCoordinate_  ; // src axis value 
     57  vector<double> destCoordinate_ ; // dst axis value 
    5358 
    5459public: 

  • XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/generic_algorithm_transformation.hpp

    r2016 r2122  
    4141    virtual void apply(int dimBefore, int dimAfter, const CArray<double,1>& dataIn, const vector<CArray<double,1>>& auxData, CArray<double,1>& dataOut) { abort() ;} //=0 
    4242    virtual bool isGenerateTransformation(void) { return false ;} 
    43  
     43    virtual bool transformAuxField(int pos) { return true ;} 
    4444    virtual vector<string> getAuxFieldId(void) ; 
    4545  protected : 

  • XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/grid_algorithm_generic.cpp

    r2011 r2122  
    6464 
    6565    gridTransformConnector_->transfer(dataIn, dataOutTmp) ; 
    66     for (int i=0; i<auxData.size();i++)  gridTransformConnector_->transfer(auxData[i], auxDataOutTmp[i]) ; 
     66    for (int i=0; i<auxData.size();i++)   
     67    { 
     68      if (algorithm_->transformAuxField(i)) gridTransformConnector_->transfer(auxData[i], auxDataOutTmp[i]) ; 
     69      else auxDataOutTmp[i].reference(auxData[i]) ; 
     70    } 
    6771 
    6872    algorithm_->apply(dimBefore_, dimAfter_, dataOutTmp, auxDataOutTmp, dataOut) ; 
Note: See TracChangeset for help on using the changeset viewer.