Changeset 1980


Ignore:
Timestamp:
11/24/20 11:33:59 (4 years ago)
Author:
yushan
Message:

trunk : axis interpolate can have coordinate source (coordinate_src) and coordinate destination (coordinate_dst), while previous attribute coordinate compatible to source

Location:
XIOS/trunk/src
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • XIOS/trunk/src/config/interpolate_axis_attribute.conf

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

  • XIOS/trunk/src/node/interpolate_axis.cpp

    r937 r1980  
    4141  { 
    4242    if (this->order.isEmpty()) this->order.setValue(1); 
     43    if (this->coordinate.isEmpty() && !this->coordinate_src.isEmpty()) this->coordinate.setValue(this->coordinate_src.getValue()); 
     44    if (this->coordinate_src.isEmpty() && !this->coordinate.isEmpty()) this->coordinate_src.setValue(this->coordinate.getValue()); 
    4345    int order = this->order.getValue(); 
    4446    if (order >= axisSrc->n_glo.getValue()) 
    … …  
    6769               << "Please define one"); 
    6870    } 
     71     
     72    if (!this->coordinate_dst.isEmpty()) 
     73    { 
     74      StdString coordinate = this->coordinate_dst.getValue(); 
     75      if (!CField::has(coordinate)) 
     76        ERROR("CInterpolateAxis::checkValid(CAxis* axisSrc)", 
     77               << "Coordinate field whose id " << coordinate << "does not exist " 
     78               << "Please define one"); 
     79    } 
     80    
    6981  } 
    7082 
    … …  
    7284  { 
    7385    std::vector<StdString> auxInputs; 
    74     if (!this->coordinate.isEmpty()) 
     86    if (!this->coordinate.isEmpty() && this->coordinate_src.isEmpty()) 
    7587    { 
    7688      StdString coordinate = this->coordinate.getValue(); 
     89      this->coordinate_src.setValue(coordinate); 
     90      if (!CField::has(coordinate)) 
     91        ERROR("CInterpolateAxis::checkValid(CAxis* axisSrc)", 
     92               << "Coordinate field whose id " << coordinate << "does not exist " 
     93               << "Please define one"); 
     94      auxInputs.push_back(coordinate); 
     95    } 
     96    if (!this->coordinate_src.isEmpty() || !this->coordinate.isEmpty()) 
     97    { 
     98      StdString coordinate = !this->coordinate.isEmpty()? this->coordinate.getValue():this->coordinate_src.getValue(); 
     99      this->coordinate.setValue(coordinate); 
     100      if (!CField::has(coordinate)) 
     101        ERROR("CInterpolateAxis::checkValid(CAxis* axisSrc)", 
     102               << "Coordinate field whose id " << coordinate << "does not exist " 
     103               << "Please define one"); 
     104      auxInputs.push_back(coordinate); 
     105    } 
     106     
     107    if (!this->coordinate_dst.isEmpty()) 
     108    { 
     109      StdString coordinate = this->coordinate_dst.getValue(); 
    77110      if (!CField::has(coordinate)) 
    78111        ERROR("CInterpolateAxis::checkValid(CAxis* axisSrc)", 

  • XIOS/trunk/src/test/generic_testcase.f90

    r1976 r1980  
    1 PROGRAM generic_testcase 
     1PROGRAM generic_testcase  
    22  USE xios 
    33  USE mod_wait 
    … …  
    161161    DOUBLE PRECISION, POINTER :: field3D_sub(:,:), other_field3D_sub(:,:) 
    162162    DOUBLE PRECISION, POINTER :: field3D_recv(:,:), other_field3D_recv(:,:) 
    163     DOUBLE PRECISION, POINTER :: pressure(:,:), other_pressure(:,:) 
     163    DOUBLE PRECISION, POINTER :: pressure(:,:), other_pressure(:,:), pressure_shifted(:,:) 
    164164 
    165165    DOUBLE PRECISION, POINTER :: field2D_W(:,:), other_field2D_W(:,:) 
    … …  
    188188    INTEGER :: ierr       
    189189 
    190     LOGICAL :: ok_field2D, ok_field3D, ok_pressure, ok_field2D_sub, ok_field3D_sub,ok_field3D_recv, ok_field3D_send 
     190    LOGICAL :: ok_field2D, ok_field3D, ok_pressure_shifted, ok_pressure, ok_field2D_sub, ok_field3D_sub,ok_field3D_recv, ok_field3D_send 
    191191    LOGICAL :: ok_field_X, ok_field_Y, ok_field_XY, ok_field_Z, ok_field_XYZ, ok_field_XZ, ok_field_YZ 
    192192    LOGICAL :: ok_field2D_W, ok_field3D_W, ok_pressure_W, ok_field2D_sub_W, ok_field3D_sub_W,ok_field3D_recv_W, ok_field3D_send_W 
    … …  
    290290    w=size(ensemble_value) 
    291291 
     292 
    292293    ALLOCATE(field2D(0:xy-1)) 
    293294    ALLOCATE(field3D(0:xy-1,0:z-1)) 
     295    ALLOCATE(pressure_shifted(0:xy-1,0:z-1)) 
    294296    ALLOCATE(pressure(0:xy-1,0:z-1)) 
    295297    ALLOCATE(field3D_recv(0:xy-1,0:z-1)) 
    … …  
    364366            dist=sqrt((lat(k)/90.)**2+(lon(k)/180.)**2) ; 
    365367            pressure(i,j)=pressure(i,j)*(1+params%pressure_factor*exp(-(4*dist)**2)) 
     368            pressure_shifted(i,j)=pressure(i,j)+5000 
    366369          ENDIF 
    367370        ENDIF 
    … …  
    382385    field_YZW(:,:,0) = field_YZ(:,:)*(1+0.1*ensemble_value(0)) 
    383386     
    384      
    385387    ok_field2D=xios_is_valid_field("field2D") ; 
    386388    ok_field3D=xios_is_valid_field("field3D") ; 
     389    ok_pressure_shifted=xios_is_valid_field("pressure_shifted") ;     
    387390    ok_pressure=xios_is_valid_field("pressure") ; 
    388391    ok_field2D_sub=xios_is_valid_field("field2D_sub") ; 
    … …  
    628631      CALL xios_update_calendar(ts) 
    629632 
     633 
    630634      IF (ok_field2D) CALL xios_send_field("field2D",field2D) 
    631635      IF (ok_field3D) CALL xios_send_field("field3D",field3D) 
     636      IF (ok_pressure_shifted) CALL xios_send_field("pressure_shifted",pressure_shifted) 
    632637      IF (ok_pressure) CALL xios_send_field("pressure",pressure) 
    633638      IF (ok_field_X) CALL xios_send_field("field_X",field_X) 
    … …  
    23842389  END SUBROUTINE get_decomposition 
    23852390 
    2386 END PROGRAM generic_testcase 
    2387  
    2388  
     2391END PROGRAM generic_testcase  
     2392 
     2393 

  • XIOS/trunk/src/test/test_interpolate.f90

    r1979 r1980  
    162162    DOUBLE PRECISION, POINTER :: field3D_sub(:,:), other_field3D_sub(:,:) 
    163163    DOUBLE PRECISION, POINTER :: field3D_recv(:,:), other_field3D_recv(:,:) 
    164     DOUBLE PRECISION, POINTER :: pressure(:,:), other_pressure(:,:) 
     164    DOUBLE PRECISION, POINTER :: pressure(:,:), other_pressure(:,:), pressure_shifted(:,:) 
    165165 
    166166    DOUBLE PRECISION, POINTER :: field2D_W(:,:), other_field2D_W(:,:) 
    … …  
    190190 
    191191    LOGICAL :: ok_field_in 
    192     LOGICAL :: ok_field2D, ok_field3D, ok_pressure, ok_field2D_sub, ok_field3D_sub,ok_field3D_recv, ok_field3D_send 
     192    LOGICAL :: ok_field2D, ok_field3D, ok_pressure_shifted, ok_pressure, ok_field2D_sub, ok_field3D_sub,ok_field3D_recv, ok_field3D_send 
    193193    LOGICAL :: ok_field_X, ok_field_Y, ok_field_XY, ok_field_Z, ok_field_XYZ, ok_field_XZ, ok_field_YZ 
    194194    LOGICAL :: ok_field2D_W, ok_field3D_W, ok_pressure_W, ok_field2D_sub_W, ok_field3D_sub_W,ok_field3D_recv_W, ok_field3D_send_W 
    … …  
    296296    ALLOCATE(field2D(0:xy-1)) 
    297297    ALLOCATE(field3D(0:xy-1,0:z-1)) 
     298    ALLOCATE(pressure_shifted(0:xy-1,0:z-1)) 
    298299    ALLOCATE(pressure(0:xy-1,0:z-1)) 
    299300    ALLOCATE(field3D_recv(0:xy-1,0:z-1)) 
    … …  
    368369            dist=sqrt((lat(k)/90.)**2+(lon(k)/180.)**2) ; 
    369370            pressure(i,j)=pressure(i,j)*(1+params%pressure_factor*exp(-(4*dist)**2)) 
     371            pressure_shifted(i,j)=pressure(i,j)+5000 
    370372          ENDIF 
    371373        ENDIF 
    … …  
    386388    field_YZW(:,:,0) = field_YZ(:,:)*(1+0.1*ensemble_value(0)) 
    387389     
    388     ok_field_in=xios_is_valid_field("field_in") ; 
     390    ok_field_in=xios_is_valid_field("field_input") ; 
     391    ok_field_in=.FALSE. 
    389392    ok_field2D=xios_is_valid_field("field2D") ; 
    390393    ok_field3D=xios_is_valid_field("field3D") ; 
     394    ok_pressure_shifted=xios_is_valid_field("pressure_shifted") ;     
    391395    ok_pressure=xios_is_valid_field("pressure") ; 
    392396    ok_field2D_sub=xios_is_valid_field("field2D_sub") ; 
    … …  
    632636      CALL xios_update_calendar(ts) 
    633637 
    634       if (ok_field_in) CALL xios_recv_field("field_in", field_in) 
     638      if (ok_field_in) CALL xios_recv_field("field_input", field_in) 
    635639      if (ok_field_in) CALL xios_send_field("field_out", field_in) 
    636640 
    637641      IF (ok_field2D) CALL xios_send_field("field2D",field2D) 
    638642      IF (ok_field3D) CALL xios_send_field("field3D",field3D) 
     643      IF (ok_pressure_shifted) CALL xios_send_field("pressure_shifted",pressure_shifted) 
    639644      IF (ok_pressure) CALL xios_send_field("pressure",pressure) 
    640645      IF (ok_field_X) CALL xios_send_field("field_X",field_X) 

  • XIOS/trunk/src/transformation/axis_algorithm_interpolate.cpp

    r1852 r1980  
    5050 
    5151CAxisAlgorithmInterpolate::CAxisAlgorithmInterpolate(CAxis* axisDestination, CAxis* axisSource, CInterpolateAxis* interpAxis) 
    52 : CAxisAlgorithmTransformation(axisDestination, axisSource), coordinate_(), transPosition_() 
     52: CAxisAlgorithmTransformation(axisDestination, axisSource), coordinate_(), coordinateDST_(),transPosition_() 
    5353TRY 
    5454{ 
    5555  interpAxis->checkValid(axisSource); 
    5656  order_ = interpAxis->order.getValue(); 
     57  this->idAuxInputs_.clear(); 
    5758  if (!interpAxis->coordinate.isEmpty()) 
    5859  { 
    … …  
    6162    this->idAuxInputs_[0] = coordinate_; 
    6263  } 
     64  else if (!interpAxis->coordinate_src.isEmpty()) 
     65  { 
     66    coordinate_ = interpAxis->coordinate_src.getValue(); 
     67    this->idAuxInputs_.resize(1); 
     68    this->idAuxInputs_[0] = coordinate_; 
     69  } 
     70  if (!interpAxis->coordinate_dst.isEmpty()) 
     71  { 
     72    coordinateDST_ = interpAxis->coordinate_dst.getValue(); 
     73    this->idAuxInputs_.resize(this->idAuxInputs_.size()+1); 
     74    this->idAuxInputs_[this->idAuxInputs_.size()-1] = coordinateDST_; 
     75  } 
     76   
     77 
    6378} 
    6479CATCH 
    … …  
    7186{ 
    7287  CTimer::get("CAxisAlgorithmInterpolate::computeIndexSourceMapping_").resume() ; 
    73   CContext* context = CContext::getCurrent(); 
    74   CContextClient* client=context->client; 
    75   int nbClient = client->clientSize; 
     88   
    7689  CArray<bool,1>& axisMask = axisSrc_->mask; 
    7790  int srcSize  = axisSrc_->n_glo.getValue(); 
    … …  
    90103    XIOSAlgorithms::sortWithIndex<double, CVectorStorage>(recvBuff, indexVec); 
    91104    for (int i = 0; i < srcSize; ++i) valueSrc[i] = recvBuff[indexVec[i]]; 
    92     computeInterpolantPoint(valueSrc, indexVec, idx); 
     105    computeInterpolantPoint(valueSrc, indexVec, dataAuxInputs, idx); 
    93106  } 
    94107  CTimer::get("CAxisAlgorithmInterpolate::computeIndexSourceMapping_").suspend() ; 
    … …  
    97110 
    98111/*! 
    99   Compute the interpolant points 
     112  Compute the interpolant points  
    100113  Assume that we have all value of axis source, with these values, need to calculate weight (coeff) of Lagrange polynomial 
    101114  \param [in] axisValue all value of axis source 
     115  \param [in] dataAuxInputs data for setting values of axis destination 
    102116  \param [in] tranPos position of axis on a domain 
    103117*/ 
    104118void CAxisAlgorithmInterpolate::computeInterpolantPoint(const std::vector<double>& axisValue, 
    105119                                                        const std::vector<int>& indexVec, 
     120                                                        const std::vector<CArray<double,1>* >& dataAuxInputs, 
    106121                                                        int transPos) 
    107122TRY 
    … …  
    115130  CArray<double,1>& axisDestValue = axisDest_->value; 
    116131  int numValue = axisDestValue.numElements(); 
     132  if(!coordinateDST_.empty()) 
     133  { 
     134    int nDomPoint = (*dataAuxInputs[0]).numElements()/numValue ; 
     135    int dst_position_in_data = dataAuxInputs.size()-1; 
     136    for(int ii=0; ii<numValue; ii++) 
     137    { 
     138      axisDestValue(ii) = (*dataAuxInputs[dst_position_in_data])(ii*nDomPoint+transPos); 
     139    } 
     140  } 
     141   
    117142  std::map<int, std::vector<std::pair<int,double> > > interpolatingIndexValues; 
    118143 
    … …  
    172197CATCH 
    173198 
     199 
     200 
    174201/*! 
    175202  Compute weight (coeff) of Lagrange's polynomial 
    … …  
    233260  int numValue = axisValue.numElements(); 
    234261 
     262 
    235263  if (srcSize == numValue)  // Only one client or axis not distributed 
    236264  { 
    … …  
    248276      } 
    249277    } 
    250  
    251278  } 
    252279  else // Axis distributed 
    … …  
    304331TRY 
    305332{ 
    306   if (coordinate_.empty()) 
     333  bool has_src = !coordinate_.empty(); 
     334  bool has_dst = !coordinateDST_.empty(); 
     335   
     336  int nb_inputs=dataAuxInputs.size(); 
     337   
     338   
     339  if (!has_src) 
    307340  { 
    308341    vecAxisValue.resize(1); 
    … …  
    312345    this->transformationWeight_.resize(1); 
    313346  } 
    314   else 
     347  else  
    315348  { 
    316349    CField* field = CField::get(coordinate_); 

  • XIOS/trunk/src/transformation/axis_algorithm_interpolate.hpp

    r1704 r1980  
    3434 
    3535  virtual StdString getName() {return "Axis Trans. Filter \\n Interpolation";} 
     36  bool isInversed() {return isInversed_;} 
    3637 
    3738protected: 
    … …  
    4243                            std::vector<double>& recvBuff, std::vector<int>& indexVec); 
    4344  void computeInterpolantPoint(const std::vector<double>& recvBuff, const std::vector<int>&, int transPos = 0); 
     45  void computeInterpolantPoint(const std::vector<double>& recvBuff, const std::vector<int>&,  
     46                               const std::vector<CArray<double,1>* >& dataAuxInputs, int transPos = 0); 
    4447  void computeWeightedValueAndMapping(const std::map<int, std::vector<std::pair<int,double> > >& interpolatingIndexValues, int transPos = 0); 
    4548  void fillInAxisValue(std::vector<CArray<double,1> >& vecAxisValue, 
    … …  
    5053  int order_; 
    5154  StdString coordinate_; 
     55  StdString coordinateDST_; 
     56  bool isInversed_; 
    5257  std::vector<std::vector<int> > transPosition_; 
    5358 
Note: See TracChangeset for help on using the changeset viewer.