Changeset 630


Ignore:
Timestamp:
07/07/15 10:46:25 (9 years ago)
Author:
mhnguyen
Message:

Implementing interpolation (polynomial) and correct some bugs

+) Implement interpolation (polynomial)
+) Correct some minor bugs relating to memory allocation
+) Clear some redundant codes

Test
+) On Curie
+) test_client and test_complete pass

Location:
XIOS/trunk
Files:
5 added
30 edited

Legend:

Unmodified
Added
Removed

  • XIOS/trunk/inputs/Version2/iodef.xml

    r624 r630  
    55   <calendar type="Gregorian" start_date="2012-03-01 15:00:00" time_origin="2012-02-29 15:00:00" /> 
    66 
    7  
    87   <field_definition level="1" enabled=".FALSE." default_value="9.96921e+36"> 
    9      <field id="field_AA"  operation="average" freq_op="3600s" domain_ref="domain_A"  axis_ref="axis_A" /> 
     8     <field id="field_Domain_Axis"  operation="average" freq_op="3600s" domain_ref="domain_A"  axis_ref="axis_A" /> 
    109     <field id="field_A"  operation="average" freq_op="3600s" grid_ref="grid_A" /> 
    1110     <field id="field_Axis"  operation="average" freq_op="3600s" grid_ref="grid_Axis" /> 
    1211     <field id="field_Two_Axis"  operation="average" freq_op="3600s" grid_ref="grid_Two_Axis" /> 
     12     <field id="field_All_Axis" operation="average" freq_op="3600s" grid_ref="grid_All_Axis" /> 
    1313     <field id="field_Axis_transformed"  operation="average" freq_op="3600s" field_ref="field_A" grid_ref="grid_Axis_tranformed" /> 
    14      <field id="field_All_Axis" operation="average" freq_op="3600s" grid_ref="grid_All_Axis" /> 
     14     <field id="field_Axis_transformed_Interpolated"  operation="average" freq_op="3600s" field_ref="field_Axis" grid_ref="grid_Axis_tranformed_Interpolated" /> 
    1515     <field id="field_Scalar" operation="average" freq_op="3600s" grid_ref="ScalarGrid" /> 
    1616   </field_definition> 
    1717 
    18  
    1918   <file_definition type="multiple_file" par_access="collective" output_freq="6h" output_level="10" enabled=".TRUE."> 
    2019     <file id="output" name="output"> 
    21 <!--        <field field_ref="field_A" />--> 
     20        <field field_ref="field_A" /> 
    2221     </file> 
    2322     <file id="output_Axis" name="output_Axis" type="one_file"> 
    24 <!--        <field field_ref="field_Axis" />--> 
     23        <field field_ref="field_Axis" /> 
    2524     </file> 
    26      <file id="output_Axis_transformed" name="output_Axis_transformed" type="one_file"> 
     25     <file id="output_All_Axis" name="output_All_Axis" type="one_file"> 
     26        <field field_ref="field_All_Axis" /> 
     27     </file> 
     28     <file id="output_Axis_transformed" name="output_Axis_transformed"> 
    2729        <field field_ref="field_Axis_transformed" /> 
    2830     </file> 
    29      <file id="output_All_Axis" name="output_All_Axis" type="one_file"> 
    30 <!--        <field field_ref="field_All_Axis" />--> 
     31     <file id="output_Axis_transformed_interpolated" name="output_Axis_transformed_interpolated"> 
     32        <field field_ref="field_Axis_transformed_Interpolated" /> 
    3133     </file> 
    3234     <file id="output_Scalar" name="output_Scalar" type="one_file"> 
    … …  
    4042     <axis id="axis_C" /> 
    4143     <axis id="axis_D" /> 
    42      <axis id="axis_E" axis_ref="axis_C"> 
    43 <!--        <inverse_axis />--> 
    44  
    45         <zoom_axis zoom_begin="1" zoom_size="2" /> 
    46 <!--        <inverse_axis />--> 
     44     <axis id="axis_E" axis_ref="axis_A"> 
     45       <inverse_axis /> 
    4746     </axis> 
    48      <axis id="axis_F" axis_ref="axis_A"> 
    49        <inverse_axis /> 
     47     <axis id="axis_F" axis_ref="axis_B"> 
     48        <inverse_axis /> 
     49        <zoom_axis zoom_begin="0" zoom_size="4" /> 
     50        <inverse_axis /> 
     51     </axis> 
     52     <axis id="axis_G" axis_ref="axis_D"> 
     53       <interpolate_axis type="polynomial" order="3"/> 
     54     </axis> 
     55     <axis id="axis_H" axis_ref="axis_C"> 
     56<!--       <inverse_axis />--> 
     57       <zoom_axis zoom_begin="0" zoom_size="3" /> 
    5058     </axis> 
    5159   </axis_definition> 
    … …  
    6371       <grid id="grid_Axis"> 
    6472         <axis axis_ref="axis_D" /> 
    65  
    6673       </grid> 
    6774       <grid id="grid_Two_Axis"> 
    … …  
    6976         <axis axis_ref="axis_B" /> 
    7077       </grid> 
    71        <grid id="grid_Axis_tranformed"> 
    72          <domain domain_ref="domain_A" /> 
    73             <axis axis_ref="axis_E" /> 
    74 <!--         <axis axis_ref="axis_F" />--> 
    75 <!--         <axis axis_ref="axis_E" />--> 
    76        </grid> 
    7778       <grid id="grid_All_Axis"> 
    7879         <axis axis_ref="axis_A" /> 
    7980         <axis axis_ref="axis_B" /> 
    8081         <axis axis_ref="axis_C" /> 
     82       </grid> 
     83       <grid id="grid_Axis_tranformed"> 
     84         <domain domain_ref="domain_A" /> 
     85            <axis axis_ref="axis_H" /> 
     86<!--         <axis axis_ref="axis_E" />--> 
     87<!--         <axis axis_ref="axis_F" />--> 
     88       </grid> 
     89       <grid id="grid_Axis_tranformed_Interpolated"> 
     90         <axis axis_ref="axis_G" /> 
    8191       </grid> 
    8292       <grid id="ScalarGrid"> 

  • XIOS/trunk/src/client_server_mapping_distributed.cpp

    r624 r630  
    5959    { 
    6060      int indexClient = std::distance(itbClientHash, itClientHash)-1; 
    61  
    6261      { 
    6362        client2ClientIndexGlobal[indexClient].push_back(globalIndexClient); 
    … …  
    6968  int* sendBuff = new int[nbClient_]; 
    7069  for (int i = 0; i < nbClient_; ++i) sendBuff[i] = 0; 
    71   std::map<int, std::vector<size_t> >::iterator it  = client2ClientIndexGlobal.begin(), 
    72                                                 ite = client2ClientIndexGlobal.end(); 
    73   for (; it != ite; ++it) sendBuff[it->first] = 1; 
     70  std::map<int, std::vector<size_t> >::iterator itb  = client2ClientIndexGlobal.begin(), it, 
     71                                                ite  = client2ClientIndexGlobal.end(); 
     72  for (it = itb; it != ite; ++it) sendBuff[it->first] = 1; 
    7473  int* recvBuff = new int[nbClient_]; 
    7574  MPI_Allreduce(sendBuff, recvBuff, nbClient_, MPI_INT, MPI_SUM, clientIntraComm_); 
    … …  
    7776  std::list<MPI_Request> sendRequest; 
    7877  if (0 != nbIndexSendToOthers) 
    79       for (it = client2ClientIndexGlobal.begin(); it != ite; ++it) 
     78      for (it = itb; it != ite; ++it) 
    8079         sendIndexGlobalToClients(it->first, it->second, clientIntraComm_, sendRequest); 
    8180 
    … …  
    8584  // The demand message contains global index; meanwhile the responds have server index information 
    8685  // Buffer to receive demand from other clients, it can be allocated or not depending whether it has demand(s) 
     86    // There are some cases we demand duplicate index so need to determine maxium size of demanding buffer 
     87  for (it = itb; it != ite; ++it) sendBuff[it->first] = (it->second).size(); 
     88  MPI_Allreduce(sendBuff, recvBuff, nbClient_, MPI_INT, MPI_SUM, clientIntraComm_); 
     89 
    8790  unsigned long* recvBuffIndexGlobal = 0; 
    88 //  int maxNbIndexDemandedFromOthers = (nbIndexAlreadyOnClient >= globalIndexToServerMapping_.size()) 
    89 //                                   ? 0 : (globalIndexToServerMapping_.size() - nbIndexAlreadyOnClient); 
    90   int maxNbIndexDemandedFromOthers = (globalIndexToServerMapping_.size() > nbIndexSendToOthers) 
    91                                       ? globalIndexToServerMapping_.size() : nbIndexSendToOthers; 
    92  
     91  int maxNbIndexDemandedFromOthers = recvBuff[clientRank_]; 
    9392  if (!isDataDistributed_) maxNbIndexDemandedFromOthers = nbDemandingClient * nbIndexSendToOthers; //globalIndexToServerMapping_.size(); // Not very optimal but it's general 
    9493 

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

    r621 r630  
    3939#endif //__XIOS_CZoomAxis__ 
    4040 
     41#ifdef __XIOS_CInterpolateAxis__ 
     42   DECLARE_NODE(InterpolateAxis, interpolate_axis) 
     43#endif //__XIOS_CInterpolateAxis__ 
     44 
    4145#ifdef __XIOS_CContext__ 
    4246   DECLARE_NODE_PAR(Context, context) 

  • XIOS/trunk/src/group_factory_decl.cpp

    r621 r630  
    2727  macro(CInverseAxisGroup) 
    2828  macro(CZoomAxisGroup) 
     29  macro(CInterpolateAxisGroup) 
    2930} 

  • XIOS/trunk/src/group_template_decl.cpp

    r621 r630  
    1616  macro(InverseAxis) 
    1717  macro(ZoomAxis) 
    18  
     18  macro(InterpolateAxis) 
    1919 
    2020} 

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

    r624 r630  
    99#include "context_client.hpp" 
    1010#include "xios_spl.hpp" 
     11#include "inverse_axis.hpp" 
     12#include "zoom_axis.hpp" 
     13#include "interpolate_axis.hpp" 
    1114 
    1215namespace xios { 
    … …  
    379382      StdString zoomAxisDefRoot("zoom_axis_definition"); 
    380383      StdString zoom("zoom_axis"); 
     384      StdString interpAxisDefRoot("interpolate_axis_definition"); 
     385      StdString interp("interpolate_axis"); 
    381386      do 
    382387      { 
    … …  
    390395          transformationMap_.push_back(std::make_pair(TRANS_ZOOM_AXIS,tmp)); 
    391396        } 
     397        else if (node.getElementName() == interp) { 
     398          CInterpolateAxis* tmp = (CInterpolateAxisGroup::get(interpAxisDefRoot))->createChild(); 
     399          tmp->parse(node); 
     400          transformationMap_.push_back(std::make_pair(TRANS_INTERPOLATE_AXIS,tmp)); 
     401        } 
    392402      } while (node.goToNextElement()) ; 
    393403      node.goToParentElement(); 

  • XIOS/trunk/src/node/axis.hpp

    r624 r630  
    1616#include "transformation.hpp" 
    1717#include "transformation_enum.hpp" 
    18 #include "inverse_axis.hpp" 
    19 #include "zoom_axis.hpp" 
    2018 
    2119namespace xios { 

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

    r624 r630  
    814814     { 
    815815        itFilterSrc = filterSources_.begin(); iteFilterSrc = filterSources_.end(); 
     816        dataToReceive = 0.0; // Reset all data destination 
    816817        for (; itFilterSrc != iteFilterSrc; ++itFilterSrc) 
    817818        { 
    … …  
    819820          { 
    820821             const std::map<int, CArray<int,1>* >& localIndexToSend = (*itFilterSrc)->grid->getTransformations()->getLocalIndexToSendFromGridSource(); 
    821              const std::map<int, std::vector<CArray<int,1>* > >& localIndexToReceive = (*itFilterSrc)->grid->getTransformations()->getLocalIndexToReceiveOnGridDest(); 
     822             const std::map<int,std::vector<std::vector<std::pair<int,double> > > >& localIndexToReceive = (*itFilterSrc)->grid->getTransformations()->getLocalIndexToReceiveOnGridDest(); 
    822823 
    823824             sendAndReceiveTransformedData(localIndexToSend, dataToSend, 
    … …  
    831832   void CField::sendAndReceiveTransformedData(const std::map<int, CArray<int,1>* >& localIndexToSend, 
    832833                                              const CArray<double, 1>& dataSrc, 
    833                                               const std::map<int, std::vector<CArray<int,1>* > >& localIndexToReceive, 
     834                                              const std::map<int,std::vector<std::vector<std::pair<int,double> > > >& localIndexToReceive, 
    834835                                              CArray<double,1>& dataDest) 
    835836   { 
    … …  
    858859 
    859860     // Receiving data on destination fields 
    860      std::map<int, std::vector<CArray<int,1>* > >::const_iterator itbRecv = localIndexToReceive.begin(), itRecv, 
     861     std::map<int,std::vector<std::vector<std::pair<int,double> > > >::const_iterator itbRecv = localIndexToReceive.begin(), itRecv, 
    861862                                                                  iteRecv = localIndexToReceive.end(); 
    862863     int recvBuffSize = 0; 
    … …  
    873874       for (int idx = 0; idx < countSize; ++idx) 
    874875       { 
    875          CArray<int,1>* localIndex_p = (itRecv->second)[idx]; 
    876          int numIndex = localIndex_p->numElements(); 
     876         const std::vector<std::pair<int,double> >& localIndex_p = (itRecv->second)[idx]; 
     877         int numIndex = localIndex_p.size(); 
    877878         for (int i = 0; i < numIndex; ++i) 
    878879         { 
    879            dataDest((*localIndex_p)(i)) = recvBuff[idx]; 
     880//           if ((localIndex_p)[i].first >= dataDest.numElements() ) 
     881           dataDest((localIndex_p)[i].first) += recvBuff[idx] * ((localIndex_p)[i].second); 
    880882         } 
    881883       } 

  • XIOS/trunk/src/node/field.hpp

    r624 r630  
    183183        void sendAndReceiveTransformedData(const std::map<int, CArray<int,1>* >& localIndexToSend, 
    184184                                           const CArray<double, 1>& dataSrc, 
    185                                            const std::map<int, std::vector<CArray<int,1>* > >& localIndexToReceive, 
     185                                           const std::map<int,std::vector<std::vector<std::pair<int,double> > > >& localIndexToReceive, 
    186186                                           CArray<double,1>& dataDest); 
    187187      public : 

  • XIOS/trunk/src/node/node_enum.hpp

    r621 r630  
    2121         eTransformation, 
    2222         eInverseAxis, 
    23          eZoomAxis 
    24 //#include "node_type.conf" 
     23         eZoomAxis, 
     24         eInterpolateAxis 
    2525 
    2626      } ENodeType; 

  • XIOS/trunk/src/node/node_type.hpp

    r621 r630  
    1212#include "inverse_axis.hpp" 
    1313#include "zoom_axis.hpp" 
     14#include "interpolate_axis.hpp" 
    1415 
    1516#endif // __XIOS_NODE_TYPE__ 

  • XIOS/trunk/src/node/transformation_enum.hpp

    r621 r630  
    1 #ifndef __XMLIO_TRANSFORMATION_ENUM__ 
    2 #define __XMLIO_TRANSFORMATION_ENUM__ 
    3  
    4 //#define DECLARE_NODE(Name_, name_)     ,e##Name_, g##Name_ 
    5 //#define DECLARE_NODE_PAR(Name_, name_) ,e##Name_, g##Name_ 
     1#ifndef __XIOS_TRANSFORMATION_ENUM__ 
     2#define __XIOS_TRANSFORMATION_ENUM__ 
    63 
    74namespace xios 
    … …  
    118      { 
    129        TRANS_ZOOM_AXIS, 
    13         TRANS_INVERSE_AXIS 
     10        TRANS_INVERSE_AXIS, 
     11        TRANS_INTERPOLATE_AXIS 
    1412      } ETranformationType; 
    1513 
    1614} // namespace xios 
    1715 
    18 #endif // __XMLIO_TRANSFORMATION_ENUM__ 
     16#endif // __XIOS_TRANSFORMATION_ENUM__ 

  • XIOS/trunk/src/object_factory_decl.cpp

    r621 r630  
    2626  macro(CInverseAxis) 
    2727  macro(CZoomAxis) 
     28  macro(CInterpolateAxis) 
    2829 
    2930  macro(CFieldGroup) 
    … …  
    3637  macro(CInverseAxisGroup) 
    3738  macro(CZoomAxisGroup) 
     39  macro(CInterpolateAxisGroup) 
    3840} 

  • XIOS/trunk/src/object_template_decl.cpp

    r621 r630  
    1111#include "inverse_axis.hpp" 
    1212#include "zoom_axis.hpp" 
     13#include "interpolate_axis.hpp" 
    1314 
    1415namespace xios 
    … …  
    2425  template class CObjectTemplate<CInverseAxis>; 
    2526  template class CObjectTemplate<CZoomAxis>; 
     27  template class CObjectTemplate<CInterpolateAxis>; 
    2628 
    2729  template class CObjectTemplate<CContextGroup>; 
    … …  
    3436  template class CObjectTemplate<CInverseAxisGroup>; 
    3537  template class CObjectTemplate<CZoomAxisGroup>; 
     38  template class CObjectTemplate<CInterpolateAxisGroup>; 
    3639} 

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

    r627 r630  
    1515  CHARACTER(len=15) :: calendar_type 
    1616  TYPE(xios_context) :: ctx_hdl 
    17   INTEGER,PARAMETER :: ni_glo=100 
    18   INTEGER,PARAMETER :: nj_glo=100 
    19   INTEGER,PARAMETER :: llm=5 
    20   DOUBLE PRECISION  :: lval(llm)=1, tsTemp 
     17  INTEGER,PARAMETER :: ni_glo=5 
     18  INTEGER,PARAMETER :: nj_glo=5 
     19  INTEGER,PARAMETER :: llm=10 
     20  INTEGER,PARAMETER :: llmInterPolated=5 
     21  DOUBLE PRECISION  :: lval(llm)=1, tsTemp, lvalInterPolated(llmInterPolated)=1 
    2122  TYPE(xios_field) :: field_hdl 
    2223  TYPE(xios_fieldgroup) :: fieldgroup_hdl 
    … …  
    2728  DOUBLE PRECISION,DIMENSION(ni_glo,nj_glo) :: lon_glo,lat_glo 
    2829  DOUBLE PRECISION :: field_A_glo(ni_glo,nj_glo,llm), lval_ni_glo(ni_glo), lval_nj_glo(nj_glo) 
    29   DOUBLE PRECISION,ALLOCATABLE :: lon(:,:), lat(:,:), field_A(:,:,:), field_All_Axis(:,:,:), lonvalue(:), & 
    30                                   field_Axis(:), lvaln(:), lval_ni(:), lval_nj(:), field_Two_Axis(:,:) 
    31   INTEGER :: ni,ibegin,iend,nj,jbegin,jend, nAxis, axisBegin, axisEnd 
     30  DOUBLE PRECISION,ALLOCATABLE :: lon(:,:),lat(:,:),field_A(:,:,:), field_All_Axis(:,:,:), lonvalue(:) , & 
     31                                  field_Axis(:), lvaln(:), lval_ni(:), lval_nj(:), field_Two_Axis(:,:), lvalnInterp(:) 
     32  INTEGER :: ni,ibegin,iend,nj,jbegin,jend, nAxis, axisBegin, axisEnd, axisterpBegin, nAxisinterp, axisinterpEnd 
    3233  INTEGER :: i,j,l,ts,n 
    3334 
    … …  
    5960 
    6061  jbegin=0 
    61   DO n=0,size-1 
    62     nj=nj_glo/size 
    63     IF (n<MOD(nj_glo,size)) nj=nj+1 
    64     IF (n==rank) exit 
    65     jbegin=jbegin+nj 
     62  CALL Distribute_index(jbegin, jend, nj, nj_glo, rank, size) 
     63 
     64  DO j=1,llm 
     65    lval(j) = j *10 
    6666  ENDDO 
    67   iend=ibegin+ni-1 ; jend=jbegin+nj-1 
     67  axisBegin = 0 
     68  CALL Distribute_index(axisBegin, axisEnd, nAxis, llm, rank, size) 
    6869 
    69   axisBegin = 0 
    70   nAxis = llm 
    71 !  DO n=0, size -1 
    72 !    nAxis = llm/size 
    73 !    IF (n<MOD(llm,size)) nAxis=nAxis+1 
    74 !    IF (n==rank) exit 
    75 !    axisBegin=axisBegin+nAxis 
    76 !  ENDDO 
    77   axisEnd=axisBegin+nAxis-1 
     70  DO j=1,llmInterPolated 
     71    lvalInterPolated(j) = j * 12 
     72  ENDDO 
     73  axisterpBegin = 0 
     74  CALL Distribute_index(axisterpBegin, axisinterpEnd, nAxisinterp, llmInterPolated, rank, size) 
    7875 
    79   DO i=1,llm 
    80     lval(i) = i 
    81   ENDDO 
     76  ALLOCATE(field_A(0:ni+1,-1:nj+2,llm), field_Two_Axis(ni_glo,1:nj), field_Axis(nAxis), field_All_Axis(1:ni,1:nj,llm), & 
     77          lon(ni,nj),lat(ni,nj), lonvalue(ni*nj), & 
     78          lvaln(nAxis), lval_ni(ni), lval_nj(nj), lvalnInterp(nAxisinterp)) 
    8279 
    83  
    84   ALLOCATE(lon(ni,nj), lat(ni,nj), field_A(0:ni+1,-1:nj+2,llm), lonvalue(ni*nj), field_Axis(nAxis), & 
    85            field_All_Axis(1:ni,1:nj,llm), lvaln(nAxis), lval_ni(ni), lval_nj(nj), field_Two_Axis(ni_glo,1:nj)) 
    8680  ALLOCATE(mask(nj)) 
    8781  DO i = 1, nj 
    … …  
    9690  field_A(1:ni,1:nj,:) = field_A_glo(ibegin+1:iend+1,jbegin+1:jend+1,:) 
    9791  field_Axis(1:nAxis)  = field_A_glo(1,1,axisBegin+1:axisEnd+1) 
     92  field_Two_Axis(:,1:nj)  = field_A_glo(:,jbegin+1:jend+1,1) 
    9893  field_All_Axis(1:ni,1:nj,:) = field_A_glo(ibegin+1:iend+1,jbegin+1:jend+1,:) 
    99   field_Two_Axis(:,1:nj)  = field_A_glo(:,jbegin+1:jend+1,1) 
     94 
    10095  lvaln(1:nAxis) = lval(axisBegin+1:axisEnd+1) 
    10196  lval_nj(1:nj) = lval_nj_glo(jbegin+1:jend+1); 
    10297  lval_ni(1:ni) = lval_ni_glo(ibegin+1:iend+1); 
     98  lvalnInterp(1:nAxisinterp) = lvalInterPolated(axisterpBegin+1:axisinterpEnd+1) 
    10399 
    104100  CALL xios_context_initialize("test",comm) 
    … …  
    113109  CALL xios_set_axis_attr("axis_C", size=llm, value=lval) 
    114110  CALL xios_set_axis_attr("axis_D", size=llm, ibegin=axisBegin, ni=nAxis, value=lvaln) 
     111!  CALL xios_set_axis_attr("axis_D", size=llm, value=lval) 
     112  CALL xios_set_axis_attr("axis_G", size=llmInterPolated, value=lvalnInterp, ibegin=axisterpBegin, ni=nAxisinterp) 
    115113  CALL xios_set_domain_attr("domain_A",ni_glo=ni_glo, nj_glo=nj_glo, ibegin=ibegin, ni=ni,jbegin=jbegin,nj=nj) 
    116114  CALL xios_set_domain_attr("domain_A",data_dim=2, data_ibegin=-1, data_ni=ni+2, data_jbegin=-2, data_nj=nj+4) 
    … …  
    119117  CALL xios_set_fieldgroup_attr("field_definition",enabled=.TRUE.) 
    120118 
    121   CALL xios_get_handle("field_definition",fieldgroup_hdl) 
    122   CALL xios_add_child(fieldgroup_hdl,field_hdl,"field_B") 
    123   CALL xios_set_attr(field_hdl,field_ref="field_A",name="field_B") 
    124  
    125   CALL xios_get_handle("output",file_hdl) 
    126   CALL xios_add_child(file_hdl,field_hdl) 
    127   CALL xios_set_attr(field_hdl,field_ref="field_A",name="field_C") 
     119!  CALL xios_get_handle("field_definition",fieldgroup_hdl) 
     120!  CALL xios_add_child(fieldgroup_hdl,field_hdl,"field_B") 
     121!  CALL xios_set_attr(field_hdl,field_ref="field_A",name="field_B") 
     122! 
     123!  CALL xios_get_handle("output",file_hdl) 
     124!  CALL xios_add_child(file_hdl,field_hdl) 
     125!  CALL xios_set_attr(field_hdl,field_ref="field_A",name="field_C") 
    128126! 
    129127!  CALL xios_get_handle("output_All_Axis",file_hdl) 
    … …  
    171169    CALL xios_send_field("field_Two_Axis",field_Two_Axis) 
    172170    CALL xios_send_field("field_All_Axis",field_All_Axis) 
    173 !    tsTemp = ts 
    174 !    CALL xios_send_scalar("field_Scalar", tsTemp) 
     171    tsTemp = ts 
     172    CALL xios_send_scalar("field_Scalar", tsTemp) 
    175173    CALL wait_us(5000) ; 
    176174  ENDDO 
    … …  
    181179  CALL MPI_FINALIZE(ierr) 
    182180 
     181CONTAINS 
     182  SUBROUTINE  Test_Interpolate 
     183 
     184 
     185  END SUBROUTINE Test_Interpolate 
     186 
     187  SUBROUTINE Distribute_index(ibegin, iend, ni, nglob, rank, size) 
     188    INTEGER, INTENT(INOUT) :: ibegin, iend, ni 
     189    INTEGER, INTENT(IN)    :: nglob, rank, size 
     190    DO n=0,size-1 
     191      ni=nglob/size 
     192      IF (n<MOD(nglob,size)) ni=ni+1 
     193      IF (n==rank) exit 
     194      ibegin=ibegin+ni 
     195    ENDDO 
     196    iend=ibegin+ni-1 
     197  END SUBROUTINE Distribute_index 
    183198END PROGRAM test_new_features 

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

    r624 r630  
    33   \author Ha NGUYEN 
    44   \since 14 May 2015 
    5    \date 09 June 2015 
     5   \date 29 June 2015 
    66 
    77   \brief Algorithm for inversing an axis.. 
    88 */ 
    99#include "axis_algorithm_inverse.hpp" 
     10#include "transformation_mapping.hpp" 
     11#include "context.hpp" 
     12#include "context_client.hpp" 
    1013 
    1114namespace xios { 
    … …  
    2225  } 
    2326 
    24  
    25   axisDestGlobalSize_ = axisDestination->size.getValue(); 
    26   int niDest = axisDestination->ni.getValue(); 
    27   int ibeginDest = axisDestination->ibegin.getValue(); 
    28  
    29   for (int idx = 0; idx < niDest; ++idx) 
    30     if ((axisDestination->mask)(idx)) axisDestGlobalIndex_.push_back(ibeginDest+idx); 
    3127  this->computeIndexSourceMapping(); 
     28  int niSrc   = axisSrc_->ni.getValue(); 
     29  int sizeSrc = axisSrc_->size.getValue(); 
     30  if (niSrc != sizeSrc) updateAxisValue(); 
     31  else 
     32  { 
     33    for (int idx = 0; idx < sizeSrc; ++idx) 
     34    { 
     35      axisDest_->value(idx) = axisSrc_->value(sizeSrc-idx-1); 
     36    } 
     37  } 
    3238} 
    3339 
    … …  
    3541{ 
    3642  std::map<int, std::vector<int> >& transMap = this->transformationMapping_; 
     43  std::map<int, std::vector<double> >& transWeight = this->transformationWeight_; 
    3744 
    3845  int globalIndexSize = axisDestGlobalIndex_.size(); 
    3946  for (int idx = 0; idx < globalIndexSize; ++idx) 
     47  { 
    4048    transMap[axisDestGlobalIndex_[idx]].push_back(axisDestGlobalSize_-axisDestGlobalIndex_[idx]-1); 
     49    transWeight[axisDestGlobalIndex_[idx]].push_back(1.0); 
     50  } 
     51} 
     52 
     53/*! 
     54  Update value on axis after inversing 
     55  After an axis is inversed, not only the data on it must be inversed but also the value 
     56*/ 
     57void CAxisAlgorithmInverse::updateAxisValue() 
     58{ 
     59  CContext* context = CContext::getCurrent(); 
     60  CContextClient* client=context->client; 
     61 
     62  int niSrc     = axisSrc_->ni.getValue(); 
     63  int ibeginSrc = axisSrc_->ibegin.getValue(); 
     64 
     65  CTransformationMapping transformationMap(axisDest_, axisSrc_); 
     66 
     67  std::map<int, std::vector<int> >& transMap = this->transformationMapping_; 
     68  std::map<int, std::vector<double> >& transWeight = this->transformationWeight_; 
     69 
     70  std::map<size_t, std::vector<std::pair<size_t,double> > > globaIndexMapFromDestToSource; 
     71  std::map<int, std::vector<int> >::const_iterator it = transMap.begin(), ite = transMap.end(); 
     72  for (; it != ite; ++it) 
     73  { 
     74    globaIndexMapFromDestToSource[it->first].push_back(make_pair((it->second)[0], (transWeight[it->first])[0])); 
     75  } 
     76 
     77  transformationMap.computeTransformationMapping(globaIndexMapFromDestToSource); 
     78 
     79  const std::map<int,std::vector<std::vector<std::pair<size_t,double> > > >& globalIndexToReceive = transformationMap.getGlobalIndexReceivedOnGridDestMapping(); 
     80  const std::map<int,std::vector<size_t> >& globalIndexToSend = transformationMap.getGlobalIndexSendToGridDestMapping(); 
     81 
     82 // Sending global index of original grid source 
     83  std::map<int,std::vector<size_t> >::const_iterator itbSend = globalIndexToSend.begin(), itSend, 
     84                                                     iteSend = globalIndexToSend.end(); 
     85 int sendBuffSize = 0; 
     86 for (itSend = itbSend; itSend != iteSend; ++itSend) sendBuffSize += (itSend->second).size(); 
     87 
     88 typedef double Scalar; 
     89 Scalar* sendBuff, *currentSendBuff; 
     90 if (0 != sendBuffSize) sendBuff = new Scalar[sendBuffSize]; 
     91 for (StdSize idx = 0; idx < sendBuffSize; ++idx) sendBuff[idx] = NumTraits<Scalar>::sfmax(); 
     92 
     93 int currentBuffPosition = 0; 
     94 for (itSend = itbSend; itSend != iteSend; ++itSend) 
     95 { 
     96   int destRank = itSend->first; 
     97   const std::vector<size_t>& globalIndexOfCurrentGridSourceToSend = itSend->second; 
     98   int countSize = globalIndexOfCurrentGridSourceToSend.size(); 
     99   for (int idx = 0; idx < (countSize); ++idx) 
     100   { 
     101     int index = globalIndexOfCurrentGridSourceToSend[idx] - ibeginSrc; 
     102     sendBuff[idx+currentBuffPosition] = (axisSrc_->value)(index); 
     103   } 
     104   currentSendBuff = sendBuff + currentBuffPosition; 
     105   MPI_Send(currentSendBuff, countSize, MPI_DOUBLE, destRank, 14, client->intraComm); 
     106   currentBuffPosition += countSize; 
     107 } 
     108 
     109 // Receiving global index of grid source sending from current grid source 
     110 std::map<int,std::vector<std::vector<std::pair<size_t,double> > > >::const_iterator itbRecv = globalIndexToReceive.begin(), itRecv, 
     111                                                                                     iteRecv = globalIndexToReceive.end(); 
     112 int recvBuffSize = 0; 
     113 for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv) recvBuffSize += (itRecv->second).size(); 
     114 
     115 Scalar* recvBuff, *currentRecvBuff; 
     116 if (0 != recvBuffSize) recvBuff = new Scalar [recvBuffSize]; 
     117 for (StdSize idx = 0; idx < recvBuffSize; ++idx) recvBuff[idx] = NumTraits<Scalar>::sfmax(); 
     118 
     119 int currentRecvBuffPosition = 0; 
     120 for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv) 
     121 { 
     122   MPI_Status status; 
     123   int srcRank = itRecv->first; 
     124   int countSize = (itRecv->second).size(); 
     125   currentRecvBuff = recvBuff + currentRecvBuffPosition; 
     126   MPI_Recv(currentRecvBuff, countSize, MPI_DOUBLE, srcRank, 14, client->intraComm, &status); 
     127   currentRecvBuffPosition += countSize; 
     128 } 
     129 
     130 int ibeginDest = axisDest_->ibegin.getValue(); 
     131 currentRecvBuff = recvBuff; 
     132 for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv) 
     133 { 
     134   int countSize = (itRecv->second).size(); 
     135   for (int idx = 0; idx < countSize; ++idx, ++currentRecvBuff) 
     136   { 
     137     int ssize = (itRecv->second)[idx].size(); 
     138     for (int i = 0; i < ssize; ++i) 
     139     { 
     140       int index = ((itRecv->second)[idx][i]).first - ibeginDest; 
     141       (axisDest_->value)(index) = *currentRecvBuff; 
     142     } 
     143   } 
     144 } 
     145 
     146 if (0 != sendBuffSize) delete [] sendBuff; 
     147 if (0 != recvBuffSize) delete [] recvBuff; 
    41148} 
    42149 

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

    r624 r630  
    2626 
    2727  virtual void computeIndexSourceMapping(); 
     28 
     29private: 
     30  void updateAxisValue(); 
    2831}; 
    2932 

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

    r624 r630  
    33   \author Ha NGUYEN 
    44   \since 14 May 2015 
    5    \date 09 June 2015 
     5   \date 29 June 2015 
    66 
    77   \brief Interface for all axis transformation algorithms. 
    … …  
    1515 
    1616CAxisAlgorithmTransformation::CAxisAlgorithmTransformation(CAxis* axisDestination, CAxis* axisSource) 
    17  : CGenericAlgorithmTransformation() 
     17 : CGenericAlgorithmTransformation(), axisDest_(axisDestination), axisSrc_(axisSource) 
    1818{ 
     19  axisDestGlobalSize_ = axisDestination->size.getValue(); 
     20  int niDest = axisDestination->ni.getValue(); 
     21  int ibeginDest = axisDestination->ibegin.getValue(); 
     22 
     23  for (int idx = 0; idx < niDest; ++idx) 
     24    if ((axisDestination->mask)(idx)) axisDestGlobalIndex_.push_back(ibeginDest+idx); 
    1925} 
    2026 
    … …  
    3743  \param[in/out] globalIndexSrcGrid array of global index of source grid (for 2d grid, this array is a line, for 3d, this array represents a plan). It should be preallocated 
    3844*/ 
    39 void CAxisAlgorithmTransformation::computeGlobalIndexFromGlobalIndexElement(int axisDestGlobalIndex, 
     45void CAxisAlgorithmTransformation::computeGlobalGridIndexFromGlobalIndexElement(int axisDestGlobalIndex, 
    4046                                                                          const std::vector<int>& axisSrcGlobalIndex, 
    4147                                                                          int axisPositionInGrid, 
    … …  
    4349                                                                          const CArray<size_t,1>& globalIndexGridDestSendToServer, 
    4450                                                                          CArray<size_t,1>& globalIndexDestGrid, 
    45                                                                           std::vector<CArray<size_t,1> >& globalIndexSrcGrid) 
     51                                                                          std::vector<std::vector<size_t> >& globalIndexSrcGrid) 
    4652{ 
    4753  int globalDim = gridDestGlobalDim.size(); 
    … …  
    9399    globalIndexDestGrid.resize(realGlobalIndexSize); 
    94100 
    95   if (axisSrcGlobalIndex.size() != globalIndexSrcGrid.size()) globalIndexSrcGrid.resize(axisSrcGlobalIndex.size()); 
     101  if (realGlobalIndexSize != globalIndexSrcGrid.size()) globalIndexSrcGrid.resize(realGlobalIndexSize); 
    96102  for (int i = 0; i < globalIndexSrcGrid.size(); ++i) 
    97     if (globalIndexSrcGrid[i].numElements() != realGlobalIndexSize) 
    98       globalIndexSrcGrid[i].resize(realGlobalIndexSize); 
    99  
     103    if (globalIndexSrcGrid[i].size() != axisSrcGlobalIndex.size()) 
     104      globalIndexSrcGrid[i].resize(axisSrcGlobalIndex.size()); 
    100105 
    101106  size_t realGlobalIndex = 0; 
    … …  
    128133    { 
    129134      globalIndexDestGrid(realGlobalIndex) = globIndex; 
    130       for (int i = 0; i < globalIndexSrcGrid.size(); ++i) 
     135      for (int i = 0; i < globalIndexSrcGrid[realGlobalIndex].size(); ++i) 
    131136      { 
    132137        currentIndex[axisPositionInGrid] = axisSrcGlobalIndex[i]; 
    … …  
    138143          globIndex += (currentIndex[k])*mulDim; 
    139144        } 
    140         (globalIndexSrcGrid[i])(realGlobalIndex) = globIndex; 
     145        (globalIndexSrcGrid[realGlobalIndex])[i] = globIndex; 
    141146      } 
    142147      ++realGlobalIndex; 

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

    r624 r630  
    33   \author Ha NGUYEN 
    44   \since 14 May 2015 
    5    \date 09 June 2015 
     5   \date 29 June 2015 
    66 
    77   \brief Interface for all axis transformation algorithms. 
    … …  
    2727 
    2828protected: 
    29   virtual void computeGlobalIndexFromGlobalIndexElement(int axisDestGlobalIndex, 
     29  virtual void computeGlobalGridIndexFromGlobalIndexElement(int axisDestGlobalIndex, 
    3030                                                        const std::vector<int>& axisSrcGlobalIndex, 
    3131                                                        int axisPositionInGrid, 
    … …  
    3333                                                        const CArray<size_t,1>& globalIndexGridDestSendToServer, 
    3434                                                        CArray<size_t,1>& globalIndexDestGrid, 
    35                                                         std::vector<CArray<size_t,1> >& globalIndexSrcGrid); 
     35                                                        std::vector<std::vector<size_t> >& globalIndexSrcGrid); 
    3636  void computeIndexSourceMapping(); 
    3737 
    … …  
    4343  int axisDestGlobalSize_; 
    4444 
     45    //! Axis on grid destination 
     46  CAxis* axisDest_; 
     47 
     48  //! Axis on grid source 
     49  CAxis* axisSrc_; 
    4550}; 
    4651 

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

    r624 r630  
    1212 
    1313CAxisAlgorithmZoom::CAxisAlgorithmZoom(CAxis* axisDestination, CAxis* axisSource, CZoomAxis* zoomAxis) 
    14 : CAxisAlgorithmTransformation(axisDestination, axisSource), axisDest_(axisDestination), axisSrc_(axisSource) 
     14: CAxisAlgorithmTransformation(axisDestination, axisSource) 
    1515{ 
    1616  zoomAxis->checkValid(axisSource); 
    … …  
    4545 
    4646  std::map<int, std::vector<int> >& transMap = this->transformationMapping_; 
     47  std::map<int, std::vector<double> >& transWeight = this->transformationWeight_; 
    4748  for (StdSize idx = 0; idx < ni; ++idx) 
    4849  { 
    4950    transMap[ibegin+idx].push_back(ibegin+idx); 
     51    transWeight[ibegin+idx].push_back(1.0); 
    5052  } 
    5153 

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

    r624 r630  
    4343  //! Global zoom size on axis 
    4444  StdSize zoomSize_; 
    45  
    46 private: 
    47   //! Axis on grid destination 
    48   CAxis* axisDest_; 
    49  
    50   //! Axis on grid source 
    51   CAxis* axisSrc_; 
    5245}; 
    5346 

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

    r624 r630  
    33   \author Ha NGUYEN 
    44   \since 14 May 2015 
    5    \date 09 June 2015 
     5   \date 29 June 2015 
    66 
    77   \brief Interface for all transformation algorithms. 
    … …  
    1212 
    1313CGenericAlgorithmTransformation::CGenericAlgorithmTransformation() 
    14  : transformationMapping_() 
     14 : transformationMapping_(), transformationWeight_() 
    1515{ 
    1616} 
    … …  
    2222  \param[in] gridDestGlobalDim global size of each dimension of grid source (all dimension must have the same size except of the one on which transformation is performed) 
    2323  \param[in] globalIndexGridDestSendToServer global index of grid destination on the current client to send to server 
    24   \param[in/out] globaIndexMapFromDestToSource mapping between transformed global index of grid destination 
    25                     and the demanded global index of grid source 
     24  \param[in/out] globaIndexWeightFromDestToSource mapping between transformed global index of grid destination 
     25             and the weighted value as long as global index from grid index source 
    2626*/ 
    2727void CGenericAlgorithmTransformation::computeGlobalSourceIndex(int elementPositionInGrid, 
    2828                                                             const std::vector<int>& gridDestGlobalDim, 
    2929                                                             const CArray<size_t,1>& globalIndexGridDestSendToServer, 
    30                                                              std::map<size_t, std::set<size_t> >& globaIndexMapFromDestToSource) 
     30                                                             std::map<size_t, std::vector<std::pair<size_t,double> > >& globaIndexWeightFromDestToSource) 
    3131{ 
    3232  std::map<int, std::vector<int> >::const_iterator itbTransMap = transformationMapping_.begin(), 
    3333                                                   itTransMap = itbTransMap, 
    3434                                                   iteTransMap = transformationMapping_.end(); 
     35  std::map<int, std::vector<double> >::const_iterator itTransWeight = transformationWeight_.begin(); 
     36  std::map<size_t, std::vector<std::pair<size_t,double> > >::iterator iteWeight, itWeight; 
    3537  std::vector<int>::const_iterator itbVec, itVec, iteVec; 
    36   std::vector<CArray<size_t,1> > globalIndexSrcGrid; //((itTransMap->second).size()); 
     38  std::vector<std::vector<size_t> > globalIndexSrcGrid; 
    3739  CArray<size_t,1> globalIndexDestGrid; 
    38   for (itTransMap = itbTransMap; itTransMap != iteTransMap; ++itTransMap) 
     40 
     41  for (itTransMap = itbTransMap; itTransMap != iteTransMap; ++itTransMap, ++itTransWeight) 
    3942  { 
    40  
    41     this->computeGlobalIndexFromGlobalIndexElement(itTransMap->first, 
     43    this->computeGlobalGridIndexFromGlobalIndexElement(itTransMap->first, 
    4244                                                   itTransMap->second, 
    4345                                                   elementPositionInGrid, 
    … …  
    4749                                                   globalIndexSrcGrid); 
    4850    size_t globalIndexSize = globalIndexDestGrid.numElements(); 
     51    std::vector<double> currentVecWeight = itTransWeight->second; 
    4952    for (size_t idx = 0; idx < globalIndexSize; ++idx) 
    5053    { 
    51       for (int i = 0; i < globalIndexSrcGrid.size(); ++i) 
     54      size_t globalIndex = globalIndexDestGrid(idx); 
     55      for (int i = 0; i < globalIndexSrcGrid[idx].size(); ++i) 
    5256      { 
    53         globaIndexMapFromDestToSource[globalIndexDestGrid(idx)].insert(globalIndexSrcGrid[i](idx)); 
     57        globaIndexWeightFromDestToSource[globalIndex].push_back(make_pair(globalIndexSrcGrid[idx][i], currentVecWeight[i])); 
    5458      } 
    5559    } 

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

    r624 r630  
    33   \author Ha NGUYEN 
    44   \since 14 May 2015 
    5    \date 09 June 2015 
     5   \date 29 June 2015 
    66 
    77   \brief Interface for all transformation algorithms. 
    … …  
    1111 
    1212#include <map> 
    13 #include <vector> 
     13#include <set> 
    1414#include "array_new.hpp" 
    1515 
    … …  
    2929                                const std::vector<int>& gridDestGlobalDim, 
    3030                                const CArray<size_t,1>& globalIndexGridDestSendToServer, 
    31                                 std::map<size_t, std::set<size_t> >& globaIndexMapFromDestToSource); 
     31                                std::map<size_t, std::vector<std::pair<size_t,double> > >& globaIndexWeightFromDestToSource); 
    3232 
    3333  /*! 
    … …  
    4747    \param[in/out] globalIndexSrcGrid array of global index of source grid (for 2d grid, this array is a line, for 3d, this array represents a plan). It should be preallocated 
    4848  */ 
    49   virtual void computeGlobalIndexFromGlobalIndexElement(int destGlobalIndex, 
     49  virtual void computeGlobalGridIndexFromGlobalIndexElement(int destGlobalIndex, 
    5050                                                        const std::vector<int>& srcGlobalIndex, 
    5151                                                        int elementPositionInGrid, 
    … …  
    5353                                                        const CArray<size_t,1>& globalIndexGridDestSendToServer, 
    5454                                                        CArray<size_t,1>& globalIndexDestGrid, 
    55                                                         std::vector<CArray<size_t,1> >& globalIndexSrcGrid) = 0; 
     55                                                        std::vector<std::vector<size_t> >& globalIndexSrcGrid) = 0; 
    5656 
    5757 
    … …  
    5959protected: 
    6060  std::map<int, std::vector<int> > transformationMapping_; 
     61  std::map<int, std::vector<double> > transformationWeight_; 
    6162}; 
    6263 

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

    r624 r630  
    33   \author Ha NGUYEN 
    44   \since 14 May 2015 
    5    \date 18 June 2015 
     5   \date 02 Jul 2015 
    66 
    77   \brief Interface for all transformations. 
    … …  
    1010#include "axis_algorithm_inverse.hpp" 
    1111#include "axis_algorithm_zoom.hpp" 
     12#include "axis_algorithm_interpolate.hpp" 
    1213#include "context.hpp" 
    1314#include "context_client.hpp" 
    1415#include "transformation_mapping.hpp" 
    15  
    1616#include "axis_algorithm_transformation.hpp" 
    1717 
    … …  
    1919CGridTransformation::CGridTransformation(CGrid* destination, CGrid* source) 
    2020: gridSource_(source), gridDestination_(destination), originalGridSource_(source), 
    21   globalIndexOfCurrentGridSource_(0), globalIndexOfOriginalGridSource_(0) 
     21  globalIndexOfCurrentGridSource_(0), globalIndexOfOriginalGridSource_(0), weightOfGlobalIndexOfOriginalGridSource_(0) 
    2222{ 
    2323  //Verify the compatibity between two grids 
    … …  
    7878  globalIndexOfCurrentGridSource_   = new CArray<size_t,1>(globalIndexGridDestSendToServer.numElements()); 
    7979  globalIndexOfOriginalGridSource_  = new CArray<size_t,1>(globalIndexGridDestSendToServer.numElements()); 
     80  weightOfGlobalIndexOfOriginalGridSource_  = new CArray<double,1>(globalIndexGridDestSendToServer.numElements()); 
    8081  *globalIndexOfCurrentGridSource_  = globalIndexGridDestSendToServer; 
    8182  *globalIndexOfOriginalGridSource_ = globalIndexGridDestSendToServer; 
     83  *weightOfGlobalIndexOfOriginalGridSource_ = 1.0; 
    8284} 
    8385 
    … …  
    8789                                                              ite = algoTransformation_.end(); 
    8890  for (it = itb; it != ite; ++it) delete (*it); 
    89  
    90   std::map<int, std::vector<CArray<int,1>* > >::const_iterator itMapRecv, iteMapRecv; 
    91   itMapRecv = localIndexToReceiveOnGridDest_.begin(); 
    92   iteMapRecv = localIndexToReceiveOnGridDest_.end(); 
    93   for (; itMapRecv != iteMapRecv; ++itMapRecv) 
    94   { 
    95     int numVec = (itMapRecv->second).size(); 
    96     for (int idx = 0; idx < numVec; ++idx) delete (itMapRecv->second)[idx]; 
    97   } 
    9891 
    9992  std::map<int, CArray<int,1>* >::const_iterator itMap, iteMap; 
    … …  
    10497  if (0 != globalIndexOfCurrentGridSource_) delete globalIndexOfCurrentGridSource_; 
    10598  if (0 != globalIndexOfOriginalGridSource_) delete globalIndexOfOriginalGridSource_; 
     99  if (0 != weightOfGlobalIndexOfOriginalGridSource_) delete weightOfGlobalIndexOfOriginalGridSource_; 
    106100} 
    107101 
    … …  
    194188 
    195189  CZoomAxis* zoomAxis = 0; 
     190  CInterpolateAxis* interpAxis = 0; 
    196191  CGenericAlgorithmTransformation* algo = 0; 
    197192  switch (transType) 
    198193  { 
     194    case TRANS_INTERPOLATE_AXIS: 
     195      interpAxis = dynamic_cast<CInterpolateAxis*> (it->second); 
     196      algo = new CAxisAlgorithmInterpolate(axisListDestP[axisIndex], axisListSrcP[axisIndex], interpAxis); 
     197      break; 
    199198    case TRANS_ZOOM_AXIS: 
    200199      zoomAxis = dynamic_cast<CZoomAxis*> (it->second); 
    … …  
    237236  switch (transType) 
    238237  { 
     238    case TRANS_INTERPOLATE_AXIS: 
    239239    case TRANS_ZOOM_AXIS: 
    240240    case TRANS_INVERSE_AXIS: 
    … …  
    263263                               ite = listAlgos_.end(), it; 
    264264  CGenericAlgorithmTransformation* algo = 0; 
     265 
    265266  for (it = itb; it != ite; ++it) 
    266267  { 
    267     std::map<size_t, std::set<size_t> > globaIndexMapFromDestToSource; 
    268268    int elementPositionInGrid = it->first; 
    269269    ETranformationType transType = (it->second).first; 
    270270    int transformationOrder = (it->second).second; 
     271    std::map<size_t, std::vector<std::pair<size_t,double> > > globaIndexWeightFromDestToSource; 
    271272 
    272273    // First of all, select an algorithm 
    … …  
    284285                                   gridDestinationDimensionSize, 
    285286                                   globalIndexGridDestSendToServer, 
    286                                    globaIndexMapFromDestToSource); 
     287                                   globaIndexWeightFromDestToSource); 
    287288 
    288289    // Compute transformation of global indexes among grids 
    289     computeTransformationFromOriginalGridSource(globaIndexMapFromDestToSource); 
     290    computeTransformationFromOriginalGridSource(globaIndexWeightFromDestToSource); 
    290291 
    291292    // Now grid destination becomes grid source in a new transformation 
    … …  
    357358the final grid destination 
    358359*/ 
    359 void CGridTransformation::computeTransformationFromOriginalGridSource(const std::map<size_t, std::set<size_t> >& globaIndexMapFromDestToSource) 
     360void CGridTransformation::computeTransformationFromOriginalGridSource(const std::map<size_t, std::vector<std::pair<size_t,double> > >& globaIndexMapFromDestToSource) 
    360361{ 
    361362  CContext* context = CContext::getCurrent(); 
    … …  
    367368  transformationMap.computeTransformationMapping(globaIndexMapFromDestToSource); 
    368369 
    369   const std::map<int,std::vector<std::vector<size_t> > >& globalIndexToReceive = transformationMap.getGlobalIndexReceivedOnGridDestMapping(); 
     370  const std::map<int,std::vector<std::vector<std::pair<size_t,double> > > >& globalIndexToReceive = transformationMap.getGlobalIndexReceivedOnGridDestMapping(); 
    370371  const std::map<int,std::vector<size_t> >& globalIndexToSend = transformationMap.getGlobalIndexSendToGridDestMapping(); 
    371372 
    … …  
    404405 
    405406 // Receiving global index of grid source sending from current grid source 
    406  std::map<int,std::vector<std::vector<size_t> > >::const_iterator itbRecv = globalIndexToReceive.begin(), itRecv, 
    407                                                                   iteRecv = globalIndexToReceive.end(); 
     407 std::map<int,std::vector<std::vector<std::pair<size_t,double> > > >::const_iterator itbRecv = globalIndexToReceive.begin(), itRecv, 
     408                                                                                     iteRecv = globalIndexToReceive.end(); 
    408409 int recvBuffSize = 0; 
    409410 for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv) recvBuffSize += (itRecv->second).size(); 
    … …  
    438439   delete globalIndexOfCurrentGridSource_; 
    439440   globalIndexOfCurrentGridSource_ = new CArray<size_t,1>(nbCurrentGridSource); 
    440  } 
    441  
    442  if (globalIndexOfOriginalGridSource_->numElements() != nbCurrentGridSource) 
    443  { 
    444441   delete globalIndexOfOriginalGridSource_; 
    445442   globalIndexOfOriginalGridSource_ = new CArray<size_t,1>(nbCurrentGridSource); 
     443   delete weightOfGlobalIndexOfOriginalGridSource_; 
     444   weightOfGlobalIndexOfOriginalGridSource_ = new CArray<double,1>(nbCurrentGridSource); 
    446445 } 
    447446 
    … …  
    456455     for (int i = 0; i < ssize; ++i) 
    457456     { 
    458        (*globalIndexOfCurrentGridSource_)(k) = (itRecv->second)[idx][i]; 
     457       (*globalIndexOfCurrentGridSource_)(k) = ((itRecv->second)[idx][i]).first; 
     458       (*weightOfGlobalIndexOfOriginalGridSource_)(k) = ((itRecv->second)[idx][i]).second; 
    459459       (*globalIndexOfOriginalGridSource_)(k) = *currentRecvBuff; 
    460460       ++k; 
    … …  
    479479  CTransformationMapping transformationMap(gridDestination_, originalGridSource_); 
    480480 
    481   std::map<size_t, std::set<size_t> > globaIndexMapFromDestToSource; 
    482  
     481  std::map<size_t, std::vector<std::pair<size_t,double> > > globaIndexWeightFromDestToSource; 
    483482  int nb = globalIndexOfCurrentGridSource_->numElements(); 
    484483  const size_t sfmax = NumTraits<unsigned long>::sfmax(); 
    … …  
    486485  { 
    487486    if (sfmax != (*globalIndexOfOriginalGridSource_)(idx)) 
    488       globaIndexMapFromDestToSource[(*globalIndexOfCurrentGridSource_)(idx)].insert((*globalIndexOfOriginalGridSource_)(idx)); 
     487      globaIndexWeightFromDestToSource[(*globalIndexOfCurrentGridSource_)(idx)].push_back(make_pair((*globalIndexOfOriginalGridSource_)(idx),(*weightOfGlobalIndexOfOriginalGridSource_)(idx))) ; 
    489488  } 
    490489 
    491490  // Then compute transformation mapping among clients 
    492   transformationMap.computeTransformationMapping(globaIndexMapFromDestToSource); 
    493  
    494   const std::map<int,std::vector<std::vector<size_t> > >& globalIndexToReceive = transformationMap.getGlobalIndexReceivedOnGridDestMapping(); 
     491  transformationMap.computeTransformationMapping(globaIndexWeightFromDestToSource); 
     492 
     493  const std::map<int,std::vector<std::vector<std::pair<size_t,double> > > >& globalIndexToReceive = transformationMap.getGlobalIndexReceivedOnGridDestMapping(); 
    495494  const std::map<int,std::vector<size_t> >& globalIndexToSend = transformationMap.getGlobalIndexSendToGridDestMapping(); 
    496495 
    … …  
    498497  CDistributionClient distributionClientSrc(client->clientRank, originalGridSource_); 
    499498 
    500 //  const CArray<int, 1>& localIndexOnClientDest = distributionClientDest.getLocalDataIndexOnClient(); //gridDestination_->getDistributionClient()->getLocalDataIndexOnClient(); 
    501499  const CArray<int, 1>& localIndexOnClientDest = distributionClientDest.getLocalDataIndexSendToServer(); 
    502   const CArray<size_t,1>& globalIndexOnClientDest = distributionClientDest.getGlobalDataIndexSendToServer(); //gridDestination_->getDistributionClient()->getGlobalDataIndexSendToServer(); 
    503  
    504   const CArray<int, 1>& localIndexOnClientSrc = distributionClientSrc.getLocalDataIndexOnClient(); //gridSource_->getDistributionClient()->getLocalDataIndexOnClient(); 
    505   const CArray<size_t,1>& globalIndexOnClientSrc = distributionClientSrc.getGlobalDataIndexSendToServer(); //gridSource_->getDistributionClient()->getGlobalDataIndexSendToServer(); 
     500  const CArray<size_t,1>& globalIndexOnClientDest = distributionClientDest.getGlobalDataIndexSendToServer(); 
     501 
     502  const CArray<int, 1>& localIndexOnClientSrc = distributionClientSrc.getLocalDataIndexOnClient(); 
     503  const CArray<size_t,1>& globalIndexOnClientSrc = distributionClientSrc.getGlobalDataIndexSendToServer(); 
    506504 
    507505  std::vector<size_t>::const_iterator itbVec, itVec, iteVec; 
    508506  CArray<size_t, 1>::const_iterator itbArr, itArr, iteArr; 
    509507 
    510   std::map<int,std::vector<std::vector<size_t> > >::const_iterator itbMapRecv, itMapRecv, iteMapRecv; 
     508  std::map<int,std::vector<std::vector<std::pair<size_t,double> > > >::const_iterator itbMapRecv, itMapRecv, iteMapRecv; 
    511509 
    512510  // Find out local index on grid destination (received) 
    … …  
    522520    { 
    523521      int vecSize = ((itMapRecv->second)[i]).size(); 
    524       CArray<int,1>* ptr = new CArray<int,1>(vecSize); 
    525       localIndexToReceiveOnGridDest_[sourceRank].push_back(ptr); 
     522      std::vector<std::pair<int,double> > tmpVec; 
    526523      for (int idx = 0; idx < vecSize; ++idx) 
    527524      { 
    528         itArr = std::find(itbArr, iteArr, (itMapRecv->second)[i][idx]); 
     525        size_t globalIndex = (itMapRecv->second)[i][idx].first; 
     526        double weight = (itMapRecv->second)[i][idx].second; 
     527        itArr = std::find(itbArr, iteArr, globalIndex); 
    529528        if (iteArr != itArr) 
    530529        { 
    531530          int localIdx = std::distance(itbArr, itArr); 
    532 //          (*localIndexToReceiveOnGridDest_[sourceRank][i])(idx) = localIndexOnClientDest(localIdx); // Local index of un-extracted data (only domain) 
    533           (*localIndexToReceiveOnGridDest_[sourceRank][i])(idx) = (localIdx); // Local index of extracted data 
     531          tmpVec.push_back(make_pair(localIdx, weight)); 
    534532        } 
    535533      } 
     534      localIndexToReceiveOnGridDest_[sourceRank].push_back(tmpVec); 
    536535    } 
    537536  } 
    538537 
     538  // Find out local index on grid source (to send) 
    539539  std::map<int,std::vector<size_t> >::const_iterator itbMap, itMap, iteMap; 
    540   // Find out local index on grid source (to send) 
    541540  itbMap = globalIndexToSend.begin(); 
    542541  iteMap = globalIndexToSend.end(); 
    … …  
    555554      { 
    556555        int localIdx = std::distance(itbArr, itArr); 
    557 //        (*localIndexToSendFromGridSource_[destRank])(idx) = localIndexOnClientSrc(localIdx); 
    558556        (*localIndexToSendFromGridSource_[destRank])(idx) = (localIdx); 
    559557      } 
    … …  
    575573  \return local index of data 
    576574*/ 
    577 const std::map<int, std::vector<CArray<int,1>* > >& CGridTransformation::getLocalIndexToReceiveOnGridDest() const 
     575const std::map<int,std::vector<std::vector<std::pair<int,double> > > >& CGridTransformation::getLocalIndexToReceiveOnGridDest() const 
    578576{ 
    579577  return localIndexToReceiveOnGridDest_; 

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

    r624 r630  
    2424  This class is an interface for all transformations to interact with the rest of XIOS. 
    2525The class, firstly, tries to get all information relating to requested transformations by retrieving directly from grid. 
    26 Then with all these information, all necessary transformations will be be created by generic class \class CGenericAlgorithmTransformation. 
     26Then with all these information, all necessary transformations will be created by generic class \class CGenericAlgorithmTransformation. 
    2727Because there are information exchange among clients to accomplish the transformations (e.g: some index need retrieving from other clients), 
    2828this class uses class \class CTransformationMapping to fulfill this demand. 
    … …  
    4141 
    4242  const std::map<int, CArray<int,1>* >& getLocalIndexToSendFromGridSource() const; 
    43   const std::map<int, std::vector<CArray<int,1>* > >& getLocalIndexToReceiveOnGridDest() const; 
     43  const std::map<int, std::vector<std::vector<std::pair<int,double> > > >& getLocalIndexToReceiveOnGridDest() const; 
    4444 
    4545private: 
    … …  
    5555  void setUpGrid(int elementPositionInGrid, ETranformationType transType); 
    5656  void computeFinalTransformationMapping(); 
    57   void computeTransformationFromOriginalGridSource(const std::map<size_t, std::set<size_t> >& globaIndexMapFromDestToSource); 
     57  void computeTransformationFromOriginalGridSource(const std::map<size_t, std::vector<std::pair<size_t,double> > >& globaIndexMapFromDestToSource); 
    5858  void updateFinalGridDestination(); 
    5959 
    … …  
    8989 
    9090  //! Local index of data to receive on grid destination 
    91   std::map<int, std::vector<CArray<int,1>* > > localIndexToReceiveOnGridDest_; 
     91  std::map<int,std::vector<std::vector<std::pair<int,double> > > > localIndexToReceiveOnGridDest_; 
    9292 
    9393  //! Position of axis and domain in grid 
    … …  
    9797  CArray<size_t,1>* globalIndexOfCurrentGridSource_; 
    9898  CArray<size_t,1>* globalIndexOfOriginalGridSource_; 
     99  CArray<double,1>* weightOfGlobalIndexOfOriginalGridSource_; 
    99100}; 
    100101 

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

    r624 r630  
    3838} 
    3939 
     40CTransformationMapping::CTransformationMapping(CAxis* destination, CAxis* source) 
     41  : gridSource_(0), gridDestination_(0) 
     42{ 
     43  CContext* context = CContext::getCurrent(); 
     44  CContextClient* client=context->client; 
     45  int clientRank = client->clientRank; 
     46 
     47  int niSrc     = source->ni.getValue(); 
     48  int ibeginSrc = source->ibegin.getValue(); 
     49 
     50  boost::unordered_map<size_t,int> globalIndexOfAxisSource; 
     51  for (int idx = 0; idx < niSrc; ++idx) 
     52  { 
     53    globalIndexOfAxisSource[idx+ibeginSrc] = clientRank; 
     54  } 
     55 
     56  gridIndexClientClientMapping_ = new CClientServerMappingDistributed(globalIndexOfAxisSource, 
     57                                                                      client->intraComm, 
     58                                                                      true); 
     59} 
     60 
    4061CTransformationMapping::~CTransformationMapping() 
    4162{ 
    … …  
    4667  Suppose that we have transformations between two grids, which are represented in form of mapping between global indexes of these two grids, 
    4768this function tries to find out which clients a client needs to send and receive these global indexes to accomplish the transformations. 
    48   The grid destination is the grid whose global indexes demande global indexes from the grid source 
     69  The grid destination is the grid whose global indexes demande global indexes from the other grid 
    4970  Grid destination and grid source are also distributed among clients but in different manners. 
    50   \param [in] globaIndexMapFromDestToSource mapping representing the transformations 
     71  \param [in] globaIndexWeightFromDestToSource mapping representing the transformations 
    5172*/ 
    52 void CTransformationMapping::computeTransformationMapping(const std::map<size_t, std::set<size_t> >& globaIndexMapFromDestToSource) 
     73void CTransformationMapping::computeTransformationMapping(const std::map<size_t, std::vector<std::pair<size_t,double> > >& globaIndexWeightFromDestToSource) 
    5374{ 
    5475  CContext* context = CContext::getCurrent(); 
    5576  CContextClient* client=context->client; 
    5677 
    57   int numMappingPoints = 0; 
    58   std::map<size_t, std::set<size_t> >::const_iterator itbMap = globaIndexMapFromDestToSource.begin(), itMap, 
    59                                                       iteMap = globaIndexMapFromDestToSource.end(); 
     78  std::map<size_t, std::vector<std::pair<size_t,double> > >::const_iterator itbMap = globaIndexWeightFromDestToSource.begin(), itMap, 
     79                                                                            iteMap = globaIndexWeightFromDestToSource.end(); 
     80 
     81  // Not only one index on grid destination can demande two indexes from grid source 
     82  // but an index on grid source has to be sent to two indexes of grid destination 
     83  std::map<size_t, std::vector<std::pair<size_t,double> > > globalIndexMapFromSrcToDest; 
     84  std::vector<std::pair<size_t,double> >::const_iterator itbVecPair, itVecPair, iteVecPair; 
    6085  for (itMap = itbMap; itMap != iteMap; ++itMap) 
    6186  { 
    62     numMappingPoints += (itMap->second).size(); 
     87    itbVecPair = (itMap->second).begin(); 
     88    iteVecPair = (itMap->second).end(); 
     89    for (itVecPair = itbVecPair; itVecPair != iteVecPair; ++itVecPair) 
     90    { 
     91      globalIndexMapFromSrcToDest[itVecPair->first].push_back(std::make_pair(itMap->first, itVecPair->second)); 
     92    } 
    6393  } 
    6494 
    6595  // All global indexes of a client on grid destination 
    66   CArray<size_t,1> globalIndexMap(numMappingPoints); 
    67   // Not only one index on grid destination can demande two indexes from grid source 
    68   // but an index on grid destination have to be sent to two indexes of grid destination 
    69   std::map<size_t, std::vector<size_t> > globalIndexMapFromSrcToDest; 
    70   std::set<size_t>::const_iterator itbSet, itSet, iteSet; 
     96  CArray<size_t,1> globalIndexMap(globalIndexMapFromSrcToDest.size()); 
     97  itbMap = globalIndexMapFromSrcToDest.begin(); 
     98  iteMap = globalIndexMapFromSrcToDest.end(); 
    7199  int idx = 0; 
    72100  for (itMap = itbMap; itMap != iteMap; ++itMap) 
    73101  { 
    74     itbSet = (itMap->second).begin(); 
    75     iteSet = (itMap->second).end(); 
    76     for (itSet = itbSet; itSet != iteSet; ++itSet) 
    77     { 
    78       globalIndexMap(idx) = *itSet; 
    79       globalIndexMapFromSrcToDest[*itSet].push_back(itMap->first); 
    80       ++idx; 
    81     } 
     102    globalIndexMap(idx) = itMap->first; 
     103    ++idx; 
    82104  } 
    83105 
    … …  
    165187  \return global index mapping to receive on grid destination 
    166188*/ 
    167 const std::map<int,std::vector<std::vector<size_t> > >& CTransformationMapping::getGlobalIndexReceivedOnGridDestMapping() const 
     189const std::map<int,std::vector<std::vector<std::pair<size_t,double> > > >& CTransformationMapping::getGlobalIndexReceivedOnGridDestMapping() const 
    168190{ 
    169191  return globalIndexReceivedOnGridDestMapping_; 

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

    r624 r630  
    1313#include <set> 
    1414#include "grid.hpp" 
     15#include "axis.hpp" 
    1516#include "array_new.hpp" 
    1617#include "client_server_mapping_distributed.hpp" 
    … …  
    2930  /** Default constructor */ 
    3031  CTransformationMapping(CGrid* destination, CGrid* source); 
     32  CTransformationMapping(CAxis* destination, CAxis* source); 
     33 
    3134  ~CTransformationMapping(); 
    3235 
    33   void computeTransformationMapping(const std::map<size_t, std::set<size_t> >& globaIndexMapFromDestToSource); 
    34   const std::map<int,std::vector<std::vector<size_t> > >& getGlobalIndexReceivedOnGridDestMapping() const; 
     36  void computeTransformationMapping(const std::map<size_t, std::vector<std::pair<size_t,double> > >& globaIndexMapFromDestToSource); 
     37  const std::map<int,std::vector<std::vector<std::pair<size_t,double> > > >& getGlobalIndexReceivedOnGridDestMapping() const; 
    3538  const std::map<int,std::vector<size_t> >& getGlobalIndexSendToGridDestMapping() const; 
    3639 
    … …  
    4346 
    4447  //! Mapping of client rank of grid source and global index received in grid destination 
    45   std::map<int,std::vector<std::vector<size_t> > > globalIndexReceivedOnGridDestMapping_; 
     48  std::map<int,std::vector<std::vector<std::pair<size_t,double> > > > globalIndexReceivedOnGridDestMapping_; 
    4649 
    4750  //! Mapping of client rank of grid destination and global index to send from grid source 

  • XIOS/trunk/src/type/type_util.hpp

    r621 r630  
    2424    class CZoomAxis; 
    2525    class CZoomAxisGroup; 
     26    class CInterpolateAxis; 
     27    class CInterpolateAxisGroup; 
    2628 
    2729  template <typename T> inline string getStrType(void); 
    … …  
    6466  macro(CZoomAxis) 
    6567  macro(CZoomAxisGroup) 
     68  macro(CInterpolateAxis) 
     69  macro(CInterpolateAxisGroup) 
    6670#undef macro 
    6771} 

  • XIOS/trunk/src/utils.hpp

    r623 r630  
    207207  } 
    208208}; 
     209 
     210template<> 
     211struct NumTraits<double> 
     212{ 
     213  typedef double Scalar; 
     214  typedef Scalar magnitudeType; 
     215 
     216  static inline Scalar sfmin() { 
     217    return std::numeric_limits<Scalar>::min(); 
     218  } 
     219  static inline Scalar sfmax() { 
     220    return std::numeric_limits<Scalar>::max(); 
     221  } 
     222}; 
     223 
     224template<class T> 
     225class sorter 
     226{ 
     227  const std::vector<T>& values; 
     228public: 
     229  sorter(const std::vector<T> &v) : values(v) {} 
     230  bool operator()(int a, int b) { return values[a] < values[b]; } 
     231}; 
     232 
     233template<class T> 
     234void order(const std::vector<T>& values, std::vector<int>& rv) 
     235{ 
     236  std::sort(rv.begin(), rv.end(), sorter<T>(values)); 
    209237} 
    210238 
     239} 
     240 
    211241#endif // __UTILS_HPP__ 

  • XIOS/trunk/src/xml_parser_decl.cpp

    r621 r630  
    77#include "file.hpp" 
    88#include "variable.hpp" 
    9 #include "transformation.hpp" 
     9//#include "transformation.hpp" 
     10#include "inverse_axis.hpp" 
     11#include "zoom_axis.hpp" 
     12#include "interpolate_axis.hpp" 
    1013 
    1114namespace xios 
    … …  
    2730    macro( InverseAxis ) 
    2831    macro( ZoomAxis ) 
     32    macro( InterpolateAxis ) 
    2933  } 
    3034} 
Note: See TracChangeset for help on using the changeset viewer.