Changeset 821 for XIOS/trunk/src

Timestamp:

02/10/16 14:35:03 (8 years ago)

Author:

mhnguyen

Message:

Making some improvements of transformation algorithm

+) Correct the way to enlisting transformations in an element (domain, axis)
+) Optimize generic transformation to make sure temporary grid to be created on demand
+) Update some mpi tag to prevent conflict
+) Correct some minor stuffs
+) Update documents

Test
+) On Curie
+) all test pass

Location:

XIOS/trunk/src

Files:

• config/axis_attribute_private.conf (added)
• config/domain_attribute_private.conf (modified) (1 diff)
• distribution_client.cpp (modified) (1 diff)
• io/nc4_data_output.cpp (modified) (1 diff)
• mpi_tag.hpp (modified) (1 diff)
• node/axis.cpp (modified) (7 diffs)
• node/axis.hpp (modified) (2 diffs)
• node/domain.cpp (modified) (5 diffs)
• node/domain.hpp (modified) (3 diffs)
• node/grid.cpp (modified) (1 diff)
• node/grid.hpp (modified) (4 diffs)
• test/test_new_features.f90 (modified) (4 diffs)
• transformation/axis_algorithm_zoom.cpp (modified) (1 diff)
• transformation/domain_algorithm_interpolate.cpp (modified) (10 diffs)
• transformation/grid_transformation.cpp (modified) (16 diffs)
• transformation/grid_transformation.hpp (modified) (3 diffs)

XIOS/trunk/src/config/domain_attribute_private.conf

@@ -13 +13 @@
 DECLARE_ARRAY_PRIVATE(double, 1, lonvalue_rectilinear_read_from_file)
 DECLARE_ARRAY_PRIVATE(double, 1, latvalue_rectilinear_read_from_file)
+
+DECLARE_ATTRIBUTE_PRIVATE(int,  global_zoom_ibegin)
+DECLARE_ATTRIBUTE_PRIVATE(int,  global_zoom_ni)
+DECLARE_ATTRIBUTE_PRIVATE(int,  global_zoom_jbegin)
+DECLARE_ATTRIBUTE_PRIVATE(int,  global_zoom_nj)

XIOS/trunk/src/distribution_client.cpp

@@ -235 +235 @@
       nBeginGlobal_.at(indexMap_[idx]) = axisList[axisIndex]->begin.getValue();
       nZoomBegin_.at((indexMap_[idx])) = axisList[axisIndex]->global_zoom_begin;
-      nZoomEnd_.at((indexMap_[idx])) = axisList[axisIndex]->global_zoom_begin + axisList[axisIndex]->global_zoom_size-1;
+      nZoomEnd_.at((indexMap_[idx])) = axisList[axisIndex]->global_zoom_begin + axisList[axisIndex]->global_zoom_n-1;
 
       dataBegin_.at(indexMap_[idx]) = axisList[axisIndex]->data_begin.getValue();

XIOS/trunk/src/io/nc4_data_output.cpp

@@ -625 +625 @@
         int zoom_begin_srv = axis->zoom_begin_srv;
         int zoom_size  = (MULTI_FILE == SuperClass::type) ? zoom_size_srv
-                                                              : axis->global_zoom_size;
+                                                              : axis->global_zoom_n;
         int zoom_begin = (MULTI_FILE == SuperClass::type) ? zoom_begin_srv
                                                               : axis->global_zoom_begin;

XIOS/trunk/src/mpi_tag.hpp

@@ -14 +14 @@
 #define MPI_DHT_INFO_0 22
 
+/* Tag for mpi communication to send and receive info of current grid source in grid transformation*/
+#define MPI_GRID_TRANSFORMATION_CURRENT_GRID_INDEX 31
+
+/* Tag for mpi communication to send and receive info of current grid source in grid transformation*/
+#define MPI_GRID_TRANSFORMATION_ORIGINAL_GRID_INDEX 32
+
+/* Tag for mpi communication to send and receive info of current grid source in grid transformation*/
+#define MPI_GRID_TRANSFORMATION_ORIGINAL_GRID_WEIGHT 33
+
+/* Tag for mpi communication to send and receive info of DOMAIN in domain interpolation*/
+#define MPI_DOMAIN_INTERPOLATION_DEST_INDEX 7
+
+/* Tag for mpi communication to send and receive info of DOMAIN in domain interpolation*/
+#define MPI_DOMAIN_INTERPOLATION_SRC_INDEX 8
+
+/* Tag for mpi communication to send and receive info of DOMAIN in domain interpolation*/
+#define MPI_DOMAIN_INTERPOLATION_WEIGHT 9
+
 #endif

XIOS/trunk/src/node/axis.cpp

@@ -26 +26 @@
       , isDistributed_(false), hasBounds_(false), isCompressible_(false)
       , numberWrittenIndexes_(0), totalNumberWrittenIndexes_(0), offsetWrittenIndexes_(0)
-      , transformationMap_(), global_zoom_begin(0), global_zoom_size(0), hasValue(false)
+      , transformationMap_(), hasValue(false)
    {
    }
@@ -35 +35 @@
       , isDistributed_(false), hasBounds_(false), isCompressible_(false)
       , numberWrittenIndexes_(0), totalNumberWrittenIndexes_(0), offsetWrittenIndexes_(0)
-      , transformationMap_(), global_zoom_begin(0), global_zoom_size(0), hasValue(false)
+      , transformationMap_(), hasValue(false)
    {
    }
@@ -266 +266 @@
    void CAxis::checkZoom(void)
    {
-     if (0 == global_zoom_size) global_zoom_size = this->n_glo.getValue();
+     if (global_zoom_begin.isEmpty()) global_zoom_begin.setValue(0);
+     if (global_zoom_n.isEmpty()) global_zoom_n.setValue(n_glo.getValue());
    }
 
@@ -391 +392 @@
     size_t ni = this->n.getValue();
     size_t ibegin = this->begin.getValue();
-    size_t zoom_end = global_zoom_begin+global_zoom_size-1;
+    size_t zoom_end = global_zoom_begin+global_zoom_n-1;
     size_t nZoomCount = 0;
     size_t nbIndex = index.numElements();
@@ -517 +518 @@
     CEventClient event(getType(), EVENT_ID_NON_DISTRIBUTED_VALUE);
 
-    int zoom_end = global_zoom_begin + global_zoom_size - 1;
+    int zoom_end = global_zoom_begin + global_zoom_n - 1;
     int nb = 0;
     for (size_t idx = 0; idx < n; ++idx)
@@ -830 +831 @@
         msg << this->getId();
         msg << ni << begin << end;
-        msg << global_zoom_begin << global_zoom_size;
+        msg << global_zoom_begin.getValue() << global_zoom_n.getValue();
         msg << isCompressible_;
 
@@ -850 +851 @@
   void CAxis::recvServerAttribut(CBufferIn& buffer)
   {
-    int ni_srv, begin_srv, end_srv, global_zoom_begin_tmp, global_zoom_size_tmp;
+    int ni_srv, begin_srv, end_srv, global_zoom_begin_tmp, global_zoom_n_tmp;
 
     buffer >> ni_srv >> begin_srv >> end_srv;
-    buffer >> global_zoom_begin_tmp >> global_zoom_size_tmp;
+    buffer >> global_zoom_begin_tmp >> global_zoom_n_tmp;
     buffer >> isCompressible_;
     global_zoom_begin = global_zoom_begin_tmp;
-    global_zoom_size  = global_zoom_size_tmp;
-    int global_zoom_end = global_zoom_begin + global_zoom_size - 1;
+    global_zoom_n  = global_zoom_n_tmp;
+    int global_zoom_end = global_zoom_begin + global_zoom_n - 1;
 
     zoom_begin_srv = global_zoom_begin > begin_srv ? global_zoom_begin : begin_srv ;

XIOS/trunk/src/node/axis.hpp

@@ -29 +29 @@
    BEGIN_DECLARE_ATTRIBUTE_MAP(CAxis)
 #  include "axis_attribute.conf"
+#  include "axis_attribute_private.conf"
    END_DECLARE_ATTRIBUTE_MAP(CAxis)
 
@@ -116 +117 @@
         int zoom_begin_srv, zoom_end_srv, zoom_size_srv;
         int ni_srv, begin_srv, end_srv;
-        int global_zoom_begin, global_zoom_size;  // The real global zoom begin and zoom size after axis is transformed (zoomed)
         CArray<double,1> value_srv;
         CArray<double,2> bound_srv;

XIOS/trunk/src/node/domain.cpp

@@ -32 +32 @@
       , isChecked(false), relFiles(), isClientChecked(false), nbConnectedClients_(), indSrv_(), connectedServerRank_()
       , hasBounds(false), hasArea(false), isDistributed_(false), nGlobDomain_(), isCompressible_(false), isUnstructed_(false)
-      , global_zoom_ni(0), global_zoom_ibegin(0), global_zoom_nj(0), global_zoom_jbegin(0)
       , isClientAfterTransformationChecked(false), hasLonLat(false)
       , lonvalue_client(), latvalue_client(), bounds_lon_client(), bounds_lat_client()
@@ -42 +41 @@
       , isChecked(false), relFiles(), isClientChecked(false), nbConnectedClients_(), indSrv_(), connectedServerRank_()
       , hasBounds(false), hasArea(false), isDistributed_(false), nGlobDomain_(), isCompressible_(false), isUnstructed_(false)
-      , global_zoom_ni(0), global_zoom_ibegin(0), global_zoom_nj(0), global_zoom_jbegin(0)
       , isClientAfterTransformationChecked(false), hasLonLat(false)
       , lonvalue_client(), latvalue_client(), bounds_lon_client(), bounds_lat_client()
@@ -663 +661 @@
      }
      computeNGlobDomain();
-
-     if (0 == global_zoom_ni) global_zoom_ni = ni_glo;
-     if (0 == global_zoom_nj) global_zoom_nj = nj_glo;
+     checkZoom();
+   }
+
+   void CDomain::checkZoom(void)
+   {
+     if (global_zoom_ibegin.isEmpty())
+      global_zoom_ibegin.setValue(0);
+     if (global_zoom_ni.isEmpty())
+      global_zoom_ni.setValue(ni_glo);
+     if (global_zoom_jbegin.isEmpty())
+      global_zoom_jbegin.setValue(0);
+     if (global_zoom_nj.isEmpty())
+      global_zoom_nj.setValue(nj_glo);
    }
 
@@ -1253 +1261 @@
         msg << this->getId() ;
         msg << ni_srv << ibegin_srv << iend_srv << nj_srv << jbegin_srv << jend_srv;
-        msg << global_zoom_ni << global_zoom_ibegin << global_zoom_nj << global_zoom_jbegin;
+        msg << global_zoom_ni.getValue() << global_zoom_ibegin.getValue() << global_zoom_nj.getValue() << global_zoom_jbegin.getValue();
         msg << isCompressible_;
 
@@ -1680 +1688 @@
   void CDomain::recvServerAttribut(CBufferIn& buffer)
   {
+    int global_zoom_ni_tmp, global_zoom_ibegin_tmp, global_zoom_nj_tmp, global_zoom_jbegin_tmp;
     buffer >> ni_srv >> ibegin_srv >> iend_srv >> nj_srv >> jbegin_srv >> jend_srv
-           >> global_zoom_ni >> global_zoom_ibegin >> global_zoom_nj >> global_zoom_jbegin
+           >> global_zoom_ni_tmp >> global_zoom_ibegin_tmp >> global_zoom_nj_tmp >> global_zoom_jbegin_tmp
            >> isCompressible_;
+
+    global_zoom_ni.setValue(global_zoom_ni_tmp);
+    global_zoom_ibegin.setValue(global_zoom_ibegin_tmp);
+    global_zoom_nj.setValue(global_zoom_nj_tmp);
+    global_zoom_jbegin.setValue(global_zoom_jbegin_tmp);
 
     int zoom_iend = global_zoom_ibegin + global_zoom_ni - 1;

XIOS/trunk/src/node/domain.hpp

@@ -95 +95 @@
          bool isCompressible(void) const;
 
-         int global_zoom_ibegin, global_zoom_ni;
-         int global_zoom_jbegin, global_zoom_nj;
-
          int ni_srv,ibegin_srv,iend_srv ;
          int zoom_ni_srv,zoom_ibegin_srv,zoom_iend_srv ;
@@ -124 +121 @@
          void sendLonLatArea(void);
          void computeConnectedServer(void) ;
-         
-         void AllgatherRectilinearLonLat(CArray<double,1>& lon, CArray<double,1>& lat, 
+
+         void AllgatherRectilinearLonLat(CArray<double,1>& lon, CArray<double,1>& lat,
                                          CArray<double,1>& lon_g, CArray<double,1>& lat_g);
-         
-         void fillInRectilinearBoundLonLat(CArray<double,1>& lon, CArray<double,1>& lat, 
+
+         void fillInRectilinearBoundLonLat(CArray<double,1>& lon, CArray<double,1>& lat,
                                            CArray<double,2>& boundsLon, CArray<double,2>& boundsLat);
          void fillInRectilinearLonLat();
@@ -170 +167 @@
          void checkArea(void);
          void checkLonLat();
+         void checkZoom(void);
 
          void checkTransformations();

XIOS/trunk/src/node/grid.cpp

@@ -262 +262 @@
      }
      this->isChecked = true;
+   }
+
+   void CGrid::createMask(void)
+   {
+      using namespace std;
+      std::vector<CDomain*> domainP = this->getDomains();
+      std::vector<CAxis*> axisP = this->getAxis();
+      int dim = domainP.size() * 2 + axisP.size();
+
+      std::vector<CArray<bool,1>* > domainMasks(domainP.size());
+      for (int i = 0; i < domainMasks.size(); ++i) domainMasks[i] = &(domainP[i]->mask_1d);
+      std::vector<CArray<bool,1>* > axisMasks(axisP.size());
+      for (int i = 0; i < axisMasks.size(); ++i) axisMasks[i] = &(axisP[i]->mask);
+
+      switch (dim) {
+        case 1:
+          checkGridMask(mask_1d, domainMasks, axisMasks, axis_domain_order, true);
+          break;
+        case 2:
+          checkGridMask(mask_2d, domainMasks, axisMasks, axis_domain_order, true);
+          break;
+        case 3:
+          checkGridMask(mask_3d, domainMasks, axisMasks, axis_domain_order, true);
+          break;
+        default:
+          break;
+      }
    }
 

XIOS/trunk/src/node/grid.hpp

@@ -206 +206 @@
          map<int, CArray<size_t, 1> > outIndexFromClient, compressedOutIndexFromClient;
          void checkMask(void);
+         void createMask(void);
          void modifyMask(const CArray<int,1>& indexToModify);
 
@@ -217 +218 @@
                           const std::vector<CArray<bool,1>* >& domainMasks,
                           const std::vector<CArray<bool,1>* >& axisMasks,
-                          const CArray<bool,1>& axisDomainOrder);
+                          const CArray<bool,1>& axisDomainOrder,
+                          bool createMask = false);
         template<int N>
         void modifyGridMask(CArray<bool,N>& gridMask, const CArray<int,1>& indexToModify);
@@ -295 +297 @@
                              const std::vector<CArray<bool,1>* >& domainMasks,
                              const std::vector<CArray<bool,1>* >& axisMasks,
-                             const CArray<bool,1>& axisDomainOrder)
+                             const CArray<bool,1>& axisDomainOrder,
+                             bool createMask)
    {
      int idx = 0;
@@ -319 +322 @@
     }
 
-    if (!gridMask.isEmpty())
+    if (!gridMask.isEmpty() && !createMask)
     {
       for (int i = 0; i < dim; ++i)

XIOS/trunk/src/test/test_new_features.f90

@@ -17 +17 @@
   INTEGER,PARAMETER :: ni_glo=100
   INTEGER,PARAMETER :: nj_glo=100
-  INTEGER,PARAMETER :: llm=2
-  INTEGER,PARAMETER :: llmInterPolated=5
+  INTEGER,PARAMETER :: llm=5
+  INTEGER,PARAMETER :: llmInterPolated=4
   DOUBLE PRECISION  :: lval(llm)=1, tsTemp, lvalInterPolated(llmInterPolated)=1
   TYPE(xios_field) :: field_hdl
@@ -45 +45 @@
 
 !!! XIOS Initialization (get the local communicator)
-
   CALL xios_initialize(id,return_comm=comm)
 
@@ -147 +146 @@
   CALL xios_set_axis_attr("axis_C", n_glo=llm, value=lval)
   CALL xios_set_axis_attr("axis_D", n_glo=llm, begin=axisBegin, n=nAxis, value=lvaln)
-  CALL xios_set_axis_attr("axis_G", n_glo=llmInterPolated, value=lvalnInterp, begin=axisterpBegin, n=nAxisinterp)
+  CALL xios_set_axis_attr("axis_E", n_glo=llmInterPolated, value=lvalnInterp, begin=axisterpBegin, n=nAxisinterp)
   CALL xios_set_domain_attr("domain_A",ni_glo=ni_glo, nj_glo=nj_glo, ibegin=ibegin, ni=ni,jbegin=jbegin,nj=nj)
   CALL xios_set_domain_attr("domain_A",data_dim=2, data_ibegin=-1, data_ni=ni+2, data_jbegin=-2, data_nj=nj+4)
@@ -206 +205 @@
     CALL xios_send_field("field_A",field_A)
     CALL xios_send_field("field_Axis",field_Axis)
-
     CALL xios_send_field("field_Two_Axis",field_Two_Axis)
     CALL xios_send_field("field_All_Axis",field_All_Axis)

XIOS/trunk/src/transformation/axis_algorithm_zoom.cpp

@@ -62 +62 @@
 {
   axisDest_->global_zoom_begin = zoomBegin_;
-  axisDest_->global_zoom_size  = zoomSize_;
+  axisDest_->global_zoom_n  = zoomSize_;
 }
 

XIOS/trunk/src/transformation/domain_algorithm_interpolate.cpp

@@ -15 +15 @@
 #include "netcdf.hpp"
 #include "mapper.hpp"
+#include "mpi_tag.hpp"
 
 namespace xios {
@@ -95 +96 @@
           double step=(domainSrc_->lon_end-domainSrc_->lon_start)/domainSrc_->ni_glo ;
           for(int i=0; i<domainSrc_->ni_glo; ++i) lon_g(i)=domainSrc_->lon_start+i*step ;
-          step=(domainSrc_->lat_end-domainSrc_->lat_start)/domainSrc_->nj_glo ; 
+          step=(domainSrc_->lat_end-domainSrc_->lat_start)/domainSrc_->nj_glo ;
           for(int i=0; i<domainSrc_->ni_glo; ++i) lat_g(i)=domainSrc_->lat_start+i*step ;
         }
@@ -160 +161 @@
           double step=(domainDest_->lon_end-domainDest_->lon_start)/domainDest_->ni_glo ;
           for(int i=0; i<domainDest_->ni_glo; ++i) lon_g(i)=domainDest_->lon_start+i*step ;
-          step=(domainDest_->lat_end-domainDest_->lat_start)/domainDest_->nj_glo ; 
+          step=(domainDest_->lat_end-domainDest_->lat_start)/domainDest_->nj_glo ;
           for(int i=0; i<domainDest_->ni_glo; ++i) lat_g(i)=domainDest_->lat_start+i*step ;
         }
@@ -167 +168 @@
     if (std::abs(poleValue - std::abs(lat_g(domainDest_->nj_glo-1))) < NumTraits<double>::epsilon()) isSouthPole = true;
 
-    
-    
-    
+
+
+
     if (isNorthPole && (0 == domainDest_->jbegin.getValue()))
     {
@@ -469 +470 @@
              MPI_INT,
              itMap->first,
-             7,
+             MPI_DOMAIN_INTERPOLATION_DEST_INDEX,
              client->intraComm,
              &sendRequest.back());
@@ -477 +478 @@
              MPI_INT,
              itMap->first,
-             8,
+             MPI_DOMAIN_INTERPOLATION_SRC_INDEX,
              client->intraComm,
              &sendRequest.back());
@@ -485 +486 @@
              MPI_DOUBLE,
              itMap->first,
-             9,
+             MPI_DOMAIN_INTERPOLATION_WEIGHT,
              client->intraComm,
              &sendRequest.back());
@@ -504 +505 @@
              MPI_INT,
              MPI_ANY_SOURCE,
-             7,
+             MPI_DOMAIN_INTERPOLATION_DEST_INDEX,
              client->intraComm,
              &recvStatus);
@@ -516 +517 @@
              MPI_INT,
              clientSrcRank,
-             8,
+             MPI_DOMAIN_INTERPOLATION_SRC_INDEX,
              client->intraComm,
              &recvStatus);
@@ -524 +525 @@
              MPI_DOUBLE,
              clientSrcRank,
-             9,
+             MPI_DOMAIN_INTERPOLATION_WEIGHT,
              client->intraComm,
              &recvStatus);

XIOS/trunk/src/transformation/grid_transformation.cpp

@@ -18 +18 @@
 #include "axis_algorithm_transformation.hpp"
 #include "distribution_client.hpp"
+#include "mpi_tag.hpp"
 
 namespace xios {
 CGridTransformation::CGridTransformation(CGrid* destination, CGrid* source)
 : gridSource_(source), gridDestination_(destination), originalGridSource_(source),
-  globalIndexOfCurrentGridSource_(), globalIndexOfOriginalGridSource_(), weightOfGlobalIndexOfOriginalGridSource_(0), algoTypes_()
+  algoTypes_(), nbAlgos_(0), currentGridIndexToOriginalGridIndex_()
+
 {
   //Verify the compatibity between two grids
@@ -41 +43 @@
   }
 
-  std::vector<CAxis*> axisSrcTmp = gridSource_->getAxis(), axisSrc;
-  std::vector<CDomain*> domainSrcTmp = gridSource_->getDomains(), domainSrc;
-  for (int idx = 0; idx < axisSrcTmp.size(); ++idx)
-  {
-    CAxis* axis = CAxis::createAxis();
-    axis->axis_ref.setValue(axisSrcTmp[idx]->getId());
-    axis->solveRefInheritance(true);
-    axis->solveInheritanceTransformation();
-    axis->checkAttributesOnClient();
-    axisSrc.push_back(axis);
-  }
-
-  for (int idx = 0; idx < domainSrcTmp.size(); ++idx)
-  {
-    CDomain* domain = CDomain::createDomain();
-    domain->domain_ref.setValue(domainSrcTmp[idx]->getId());
-    domain->solveRefInheritance(true);
-    domain->solveInheritanceTransformation();
-    domain->checkAttributesOnClient();
-    domainSrc.push_back(domain);
-  }
-
-  gridSource_ = CGrid::createGrid(domainSrc, axisSrc, gridDestination_->axis_domain_order);
-  gridSource_->computeGridGlobalDimension(domainSrc, axisSrc, gridDestination_->axis_domain_order);
+//  std::vector<CAxis*> axisSrcTmp = gridSource_->getAxis(), axisSrc;
+//  std::vector<CDomain*> domainSrcTmp = gridSource_->getDomains(), domainSrc;
+//  for (int idx = 0; idx < axisSrcTmp.size(); ++idx)
+//  {
+//    CAxis* axis = CAxis::createAxis();
+//    axis->axis_ref.setValue(axisSrcTmp[idx]->getId());
+//    axis->solveRefInheritance(true);
+//    axis->solveInheritanceTransformation();
+//    axis->checkAttributesOnClient();
+//    axisSrc.push_back(axis);
+//  }
+//
+//  for (int idx = 0; idx < domainSrcTmp.size(); ++idx)
+//  {
+//    CDomain* domain = CDomain::createDomain();
+//    domain->domain_ref.setValue(domainSrcTmp[idx]->getId());
+//    domain->solveRefInheritance(true);
+//    domain->solveInheritanceTransformation();
+//    domain->checkAttributesOnClient();
+//    domainSrc.push_back(domain);
+//  }
+//
+//  gridSource_ = CGrid::createGrid(domainSrc, axisSrc, gridDestination_->axis_domain_order);
+//  gridSource_->computeGridGlobalDimension(domainSrc, axisSrc, gridDestination_->axis_domain_order);
 
   initializeMappingOfOriginalGridSource();
-  initializeAlgorithms();
+//  initializeAlgorithms();
 }
 
@@ -81 +83 @@
   CContextClient* client = context->client;
 
-  CDistributionClient distribution(client->clientRank, originalGridSource_);
-  const std::vector<size_t>& globalIndexGridSrcSendToServer = distribution.getGlobalDataIndexSendToServer();
-
-  weightOfGlobalIndexOfOriginalGridSource_.resize(globalIndexGridSrcSendToServer.size());
-  globalIndexOfCurrentGridSource_  = globalIndexGridSrcSendToServer;
-  globalIndexOfOriginalGridSource_ = globalIndexGridSrcSendToServer;
-  weightOfGlobalIndexOfOriginalGridSource_ = 1.0;
+  // Initialize algorithms
+  initializeAlgorithms();
+
+  ListAlgoType::const_iterator itb = listAlgos_.begin(),
+                               ite = listAlgos_.end(), it;
+
+  for (it = itb; it != ite; ++it)
+  {
+    ETranformationType transType = (it->second).first;
+    if (!isSpecialTransformation(transType)) ++nbAlgos_;
+  }
+
+  if (1<nbAlgos_)  // Only when there are more than 1 algorithm, will we create temporary grid source
+  {
+    std::vector<CAxis*> axisSrcTmp = gridSource_->getAxis(), axisSrc;
+    std::vector<CDomain*> domainSrcTmp = gridSource_->getDomains(), domainSrc;
+    for (int idx = 0; idx < axisSrcTmp.size(); ++idx)
+    {
+      CAxis* axis = CAxis::createAxis();
+      axis->axis_ref.setValue(axisSrcTmp[idx]->getId());
+      axis->solveRefInheritance(true);
+      axis->solveInheritanceTransformation();
+      axis->checkAttributesOnClient();
+      axisSrc.push_back(axis);
+    }
+
+    for (int idx = 0; idx < domainSrcTmp.size(); ++idx)
+    {
+      CDomain* domain = CDomain::createDomain();
+      domain->domain_ref.setValue(domainSrcTmp[idx]->getId());
+      domain->solveRefInheritance(true);
+      domain->solveInheritanceTransformation();
+      domain->checkAttributesOnClient();
+      domainSrc.push_back(domain);
+    }
+
+    gridSource_ = CGrid::createGrid(domainSrc, axisSrc, gridDestination_->axis_domain_order);
+    gridSource_->computeGridGlobalDimension(domainSrc, axisSrc, gridDestination_->axis_domain_order);
+  }
 }
 
@@ -149 +183 @@
   }
 }
-
-
 
 /*!
@@ -335 +367 @@
       break;
   }
+  gridSource_->createMask();
+  gridSource_->computeGridGlobalDimension(domListSrcP, axisListSrcP, gridSource_->axis_domain_order);
 }
 
@@ -352 +386 @@
   ListAlgoType::const_iterator itb = listAlgos_.begin(),
                                ite = listAlgos_.end(), it;
+
   CGenericAlgorithmTransformation* algo = 0;
   int nbAgloTransformation = 0; // Only count for executed transformation. Generate domain is a special one, not executed in the list
@@ -381 +416 @@
                                      globaIndexWeightFromDestToSource);
 
-      // Compute transformation of global indexes among grids
-      computeTransformationFromOriginalGridSource(globaIndexWeightFromDestToSource);
-
-      // Now grid destination becomes grid source in a new transformation
-      setUpGrid(elementPositionInGrid, transType);
+      if (1 < nbAlgos_)
+      {
+        // Compute transformation of global indexes among grids
+        computeTransformationFromOriginalGridSource(globaIndexWeightFromDestToSource);
+
+        // Now grid destination becomes grid source in a new transformation
+        if (nbAgloTransformation != (nbAlgos_-1)) setUpGrid(elementPositionInGrid, transType);
+      }
+      else
+      {
+        currentGridIndexToOriginalGridIndex_.swap(globaIndexWeightFromDestToSource);
+      }
+
       ++nbAgloTransformation;
     }
@@ -416 +459 @@
   iteArr = globalIndexOnClientDest.end();
 
-  // Then find out which index became invalid (become masked after being applied the algorithms, or demande some masked points from grid source)
-  int num = globalIndexOfOriginalGridSource_.size();
+  GlobalIndexMap::const_iterator iteGlobalMap = currentGridIndexToOriginalGridIndex_.end();
   const size_t sfmax = NumTraits<unsigned long>::sfmax();
   int maskIndexNum = 0;
-  for (int idx = 0; idx < num; ++idx)
-  {
-    if (sfmax == globalIndexOfOriginalGridSource_[idx])
-    {
-      size_t maskedGlobalIndex = globalIndexOfCurrentGridSource_[idx];
-      itArr = std::find(itbArr, iteArr, maskedGlobalIndex);
-      if (iteArr != itArr) ++maskIndexNum;
+  for (itArr = itbArr; itArr != iteArr; ++itArr)
+  {
+    if (iteGlobalMap != currentGridIndexToOriginalGridIndex_.find(*itArr))
+    {
+      const std::vector<std::pair<size_t,double> >& vecIndex = currentGridIndexToOriginalGridIndex_[*itArr];
+      for (int idx = 0; idx < vecIndex.size(); ++idx)
+      {
+        if (sfmax == vecIndex[idx].first)
+        {
+          ++maskIndexNum;
+          break;
+        }
+      }
     }
   }
@@ -432 +480 @@
   CArray<int,1> maskIndexToModify(maskIndexNum);
   maskIndexNum = 0;
-  for (int idx = 0; idx < num; ++idx)
-  {
-    if (sfmax == globalIndexOfOriginalGridSource_[idx])
-    {
-      size_t maskedGlobalIndex = globalIndexOfCurrentGridSource_[idx];
-      itArr = std::find(itbArr, iteArr, maskedGlobalIndex);
-      if (iteArr != itArr)
+  for (itArr = itbArr; itArr != iteArr; ++itArr)
+  {
+    if (iteGlobalMap != currentGridIndexToOriginalGridIndex_.find(*itArr))
+    {
+      const std::vector<std::pair<size_t,double> >& vecIndex = currentGridIndexToOriginalGridIndex_[*itArr];
+      for (int idx = 0; idx < vecIndex.size(); ++idx)
       {
-        int localIdx = std::distance(itbArr, itArr);
-        maskIndexToModify(maskIndexNum) = localMaskIndexOnClientDest[localIdx];
-        ++maskIndexNum;
+        if (sfmax == vecIndex[idx].first)
+        {
+          int localIdx = std::distance(itbArr, itArr);
+          maskIndexToModify(maskIndexNum) = localMaskIndexOnClientDest[localIdx];
+          ++maskIndexNum;
+          break;
+        }
       }
     }
@@ -460 +511 @@
   CContextClient* client = context->client;
 
+  if (currentGridIndexToOriginalGridIndex_.empty())
+  {
+    currentGridIndexToOriginalGridIndex_ = globaIndexMapFromDestToSource;
+    return;
+  }
+
   CTransformationMapping transformationMap(gridDestination_, gridSource_);
 
@@ -471 +528 @@
   std::map<int,std::vector<size_t> >::const_iterator itbSend = globalIndexToSend.begin(), itSend,
                                                      iteSend = globalIndexToSend.end();
-  std::vector<size_t>::const_iterator itbArr = globalIndexOfCurrentGridSource_.begin(), itArr,
-                                      iteArr = globalIndexOfCurrentGridSource_.end();
  int sendBuffSize = 0;
  for (itSend = itbSend; itSend != iteSend; ++itSend) sendBuffSize += (itSend->second).size();
+ // We use the first element of each block to send number of element in this block
+ sendBuffSize += globalIndexToSend.size();
+
 
  typedef unsigned long Scalar;
- unsigned long* sendBuff, *currentSendBuff;
- if (0 != sendBuffSize) sendBuff = new unsigned long [sendBuffSize];
+ Scalar* sendBuff, *currentSendBuff;
+ if (0 != sendBuffSize) sendBuff = new Scalar [sendBuffSize];
  for (StdSize idx = 0; idx < sendBuffSize; ++idx) sendBuff[idx] = NumTraits<Scalar>::sfmax();
 
- std::map<int, MPI_Request> requests;
-
- std::vector<int> permutIndex(globalIndexOfCurrentGridSource_.size());
- typedef XIOSBinarySearchWithIndex<size_t> BinarySearch;
- XIOSAlgorithms::fillInIndex(globalIndexOfCurrentGridSource_.size(), permutIndex);
- XIOSAlgorithms::sortWithIndex<size_t, CVectorStorage>(globalIndexOfCurrentGridSource_, permutIndex);
- BinarySearch searchCurrentSrc(globalIndexOfCurrentGridSource_);
- std::vector<int>::iterator itbIndex = permutIndex.begin(), itIndex,
-                            iteIndex = permutIndex.end();
-
-  // Find out local index on grid destination (received)
+ std::map<int, MPI_Request> requestsCurrentGrid, requestsOriginalGrid, requestsWeightGrid;
+ GlobalIndexMap::const_iterator iteGlobalIndex = currentGridIndexToOriginalGridIndex_.end();
+
+  // Only send global index of original source corresponding to non-masked index
+  // Use first position of each block to specify the number of elemnt in this block
+ int globalIndexOriginalSrcSendBuffSize = 0;
  int currentBuffPosition = 0;
  for (itSend = itbSend; itSend != iteSend; ++itSend)
@@ -497 +550 @@
    int destRank = itSend->first;
    const std::vector<size_t>& globalIndexOfCurrentGridSourceToSend = itSend->second;
-   int countSize = globalIndexOfCurrentGridSourceToSend.size();
-   for (int idx = 0; idx < (countSize); ++idx)
+   int countSize  = globalIndexOfCurrentGridSourceToSend.size();
+   size_t countBlock = 0;
+   for (int idx = 0; idx < countSize; ++idx)
    {
-     if (searchCurrentSrc.search(itbIndex, iteIndex, globalIndexOfCurrentGridSourceToSend[idx], itIndex))
+     size_t index = globalIndexOfCurrentGridSourceToSend[idx];
+     if (iteGlobalIndex != currentGridIndexToOriginalGridIndex_.find(index))// searchCurrentSrc.search(itbIndex, iteIndex, globalIndexOfCurrentGridSourceToSend[idx], itIndex))
      {
-       sendBuff[idx+currentBuffPosition] = globalIndexOfOriginalGridSource_[*itIndex];
+       globalIndexOriginalSrcSendBuffSize += currentGridIndexToOriginalGridIndex_[index].size() + 1; // 1 for number of elements in this block
+       sendBuff[idx+currentBuffPosition+1] = index;
+       countBlock += currentGridIndexToOriginalGridIndex_[index].size() + 1;
      }
    }
+   sendBuff[currentBuffPosition] = countBlock;
    currentSendBuff = sendBuff + currentBuffPosition;
-   MPI_Isend(currentSendBuff, countSize, MPI_UNSIGNED_LONG, destRank, 14, client->intraComm, &requests[destRank]);
-   currentBuffPosition += countSize;
+   MPI_Isend(currentSendBuff, countSize +1, MPI_UNSIGNED_LONG, destRank, MPI_GRID_TRANSFORMATION_CURRENT_GRID_INDEX, client->intraComm, &requestsCurrentGrid[destRank]);
+   currentBuffPosition += countSize + 1;
  }
+
+ Scalar* sendOriginalIndexBuff, *currentOriginalIndexSendBuff;
+ if (0 != globalIndexOriginalSrcSendBuffSize) sendOriginalIndexBuff = new Scalar [globalIndexOriginalSrcSendBuffSize];
+ double* sendOriginalWeightBuff, *currentOriginalWeightSendBuff;
+ if (0 != globalIndexOriginalSrcSendBuffSize) sendOriginalWeightBuff = new double [globalIndexOriginalSrcSendBuffSize];
+
+ currentBuffPosition = 0;
+ for (itSend = itbSend; itSend != iteSend; ++itSend)
+ {
+   int destRank = itSend->first;
+   const std::vector<size_t>& globalIndexOfCurrentGridSourceToSend = itSend->second;
+   int countSize = globalIndexOfCurrentGridSourceToSend.size();
+   int increaseStep = 0;
+   for (int idx = 0; idx < countSize; ++idx)
+   {
+     size_t index = globalIndexOfCurrentGridSourceToSend[idx];
+     if (iteGlobalIndex != currentGridIndexToOriginalGridIndex_.find(index))
+     {
+       size_t vectorSize = currentGridIndexToOriginalGridIndex_[index].size();
+       sendOriginalIndexBuff[currentBuffPosition+increaseStep]  = vectorSize;
+       sendOriginalWeightBuff[currentBuffPosition+increaseStep] = (double)vectorSize;
+       const std::vector<std::pair<size_t,double> >& indexWeightPair = currentGridIndexToOriginalGridIndex_[index];
+       for (size_t i = 0; i < vectorSize; ++i)
+       {
+         ++increaseStep;
+         sendOriginalIndexBuff[currentBuffPosition+increaseStep]  = indexWeightPair[i].first;
+         sendOriginalWeightBuff[currentBuffPosition+increaseStep] = indexWeightPair[i].second;
+       }
+       ++increaseStep;
+     }
+   }
+
+   currentOriginalIndexSendBuff = sendOriginalIndexBuff + currentBuffPosition;
+   currentOriginalWeightSendBuff = sendOriginalWeightBuff + currentBuffPosition;
+   if (0 != increaseStep)
+   {
+     MPI_Isend(currentOriginalIndexSendBuff, increaseStep, MPI_UNSIGNED_LONG, destRank, MPI_GRID_TRANSFORMATION_ORIGINAL_GRID_INDEX, client->intraComm, &requestsOriginalGrid[destRank]);
+     MPI_Isend(currentOriginalWeightSendBuff, increaseStep, MPI_DOUBLE, destRank, MPI_GRID_TRANSFORMATION_ORIGINAL_GRID_WEIGHT, client->intraComm, &requestsWeightGrid[destRank]);
+   }
+   currentBuffPosition += increaseStep;
+ }
+
 
  // Receiving global index of grid source sending from current grid source
@@ -515 +615 @@
  int recvBuffSize = 0;
  for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv) recvBuffSize += (itRecv->second).size();
-
- unsigned long* recvBuff, *currentRecvBuff;
- if (0 != recvBuffSize) recvBuff = new unsigned long [recvBuffSize];
+ recvBuffSize += globalIndexToReceive.size();
+
+ Scalar* recvBuff, *currentRecvBuff;
+ if (0 != recvBuffSize) recvBuff = new Scalar [recvBuffSize];
  for (StdSize idx = 0; idx < recvBuffSize; ++idx) recvBuff[idx] = NumTraits<Scalar>::sfmax();
 
+ std::map<int,int> countBlockMap;
+ int globalIndexOriginalSrcRecvBuffSize = 0;
  int currentRecvBuffPosition = 0;
  for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv)
@@ -527 +630 @@
    int countSize = (itRecv->second).size();
    currentRecvBuff = recvBuff + currentRecvBuffPosition;
-   MPI_Recv(currentRecvBuff, countSize, MPI_UNSIGNED_LONG, srcRank, 14, client->intraComm, &status);
-   currentRecvBuffPosition += countSize;
+   MPI_Recv(currentRecvBuff, countSize +1, MPI_UNSIGNED_LONG, srcRank, MPI_GRID_TRANSFORMATION_CURRENT_GRID_INDEX, client->intraComm, &status);
+   globalIndexOriginalSrcRecvBuffSize += *currentRecvBuff;
+   countBlockMap[srcRank] = *currentRecvBuff;
+   currentRecvBuffPosition += countSize +1;
  }
 
- int nbCurrentGridSource = 0;
+ Scalar* recvOriginalIndexBuff, *currentOriginalIndexRecvBuff;
+ if (0 != globalIndexOriginalSrcRecvBuffSize) recvOriginalIndexBuff = new Scalar [globalIndexOriginalSrcRecvBuffSize];
+ double* recvOriginalWeightBuff, *currentOriginalWeightRecvBuff;
+ if (0 != globalIndexOriginalSrcRecvBuffSize) recvOriginalWeightBuff = new double [globalIndexOriginalSrcRecvBuffSize];
+
+ int countBlock = 0;
+ currentRecvBuffPosition = 0;
+ currentBuffPosition = 0;
  for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv)
  {
-   int ssize = (itRecv->second).size();
-   for (int idx = 0; idx < ssize; ++idx)
+   MPI_Status statusIndex, statusWeight;
+   int srcRank = itRecv->first;
+   countBlock = countBlockMap[srcRank];
+   currentOriginalIndexRecvBuff = recvOriginalIndexBuff + currentBuffPosition;
+   currentOriginalWeightRecvBuff = recvOriginalWeightBuff + currentBuffPosition;
+   if (0 != countBlock)
    {
-     nbCurrentGridSource += (itRecv->second)[idx].size();
+     MPI_Recv(currentOriginalIndexRecvBuff, countBlock, MPI_UNSIGNED_LONG, srcRank, MPI_GRID_TRANSFORMATION_ORIGINAL_GRID_INDEX, client->intraComm, &statusIndex);
+     MPI_Recv(currentOriginalWeightRecvBuff, countBlock, MPI_DOUBLE, srcRank, MPI_GRID_TRANSFORMATION_ORIGINAL_GRID_WEIGHT, client->intraComm, &statusWeight);
+   }
+   currentBuffPosition += countBlock;
+ }
+
+ // We process everything in here, even case of masked index
+ // The way to process masked index needs discussing
+ const size_t sfmax = NumTraits<unsigned long>::sfmax();
+ GlobalIndexMap currentToOriginalTmp;
+
+ currentRecvBuffPosition = 0;
+ currentRecvBuff = recvBuff;
+ currentOriginalIndexRecvBuff  = recvOriginalIndexBuff;
+ currentOriginalWeightRecvBuff = recvOriginalWeightBuff;
+ for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv)
+ {
+   int countBlockRank = countBlockMap[itRecv->first];
+
+   ++currentRecvBuff;  // it's very subtle here, pay attention
+   int countSize = (itRecv->second).size();
+   for (int idx = 0; idx < countSize; ++idx)
+   {
+      ++currentRecvBuff;
+     int ssize = (itRecv->second)[idx].size();
+     if (sfmax != *currentRecvBuff)
+     {
+       if (0 != countBlockRank)
+       {
+         countBlock = *(currentOriginalIndexRecvBuff+currentRecvBuffPosition);
+         for (int i = 0; i < ssize; ++i)
+         {
+           for (int j = 0; j < countBlock; ++j)
+           {
+             size_t globalOriginalIndex = *(currentOriginalIndexRecvBuff+currentRecvBuffPosition+j+1);
+             double weightGlobal = *(currentOriginalWeightRecvBuff+currentRecvBuffPosition+j+1) * (itRecv->second)[idx][i].second;
+             currentToOriginalTmp[(itRecv->second)[idx][i].first].push_back(make_pair(globalOriginalIndex,weightGlobal));
+           }
+         }
+         currentRecvBuffPosition += countBlock+1;
+       }
+     }
+//     else
+//     {
+//       for (int i = 0; i < ssize; ++i)
+//       {
+//         currentToOriginalTmp[(itRecv->second)[idx][i].first].push_back(make_pair(sfmax,1.0));
+//       }
+//     }
    }
  }
 
- if (globalIndexOfCurrentGridSource_.size() != nbCurrentGridSource)
- {
-   globalIndexOfCurrentGridSource_.resize(nbCurrentGridSource);
-   globalIndexOfOriginalGridSource_.resize(nbCurrentGridSource);
-   weightOfGlobalIndexOfOriginalGridSource_.resize(nbCurrentGridSource);
- }
-
- int k = 0;
- currentRecvBuff = recvBuff;
- for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv)
- {
-   int countSize = (itRecv->second).size();
-   for (int idx = 0; idx < countSize; ++idx, ++currentRecvBuff)
-   {
-     int ssize = (itRecv->second)[idx].size();
-     for (int i = 0; i < ssize; ++i)
-     {
-       globalIndexOfCurrentGridSource_[k] = ((itRecv->second)[idx][i]).first;
-       weightOfGlobalIndexOfOriginalGridSource_(k) = ((itRecv->second)[idx][i]).second;
-       globalIndexOfOriginalGridSource_[k] = *currentRecvBuff;
-       ++k;
-     }
-   }
- }
+ currentGridIndexToOriginalGridIndex_.swap(currentToOriginalTmp);
 
  std::map<int, MPI_Request>::iterator itRequest;
- for (itRequest = requests.begin(); itRequest != requests.end(); ++itRequest)
+ for (itRequest = requestsCurrentGrid.begin(); itRequest != requestsCurrentGrid.end(); ++itRequest)
+   MPI_Wait(&itRequest->second, MPI_STATUS_IGNORE);
+ for (itRequest = requestsOriginalGrid.begin(); itRequest != requestsOriginalGrid.end(); ++itRequest)
+   MPI_Wait(&itRequest->second, MPI_STATUS_IGNORE);
+ for (itRequest = requestsWeightGrid.begin(); itRequest != requestsWeightGrid.end(); ++itRequest)
    MPI_Wait(&itRequest->second, MPI_STATUS_IGNORE);
 
  if (0 != sendBuffSize) delete [] sendBuff;
  if (0 != recvBuffSize) delete [] recvBuff;
+ if (0 != globalIndexOriginalSrcSendBuffSize) delete [] sendOriginalIndexBuff;
+ if (0 != globalIndexOriginalSrcSendBuffSize) delete [] sendOriginalWeightBuff;
+ if (0 != globalIndexOriginalSrcRecvBuffSize) delete [] recvOriginalIndexBuff;
+ if (0 != globalIndexOriginalSrcRecvBuffSize) delete [] recvOriginalWeightBuff;
 }
 
@@ -586 +735 @@
   CTransformationMapping transformationMap(gridDestination_, originalGridSource_);
 
-  std::map<size_t, std::vector<std::pair<size_t,double> > > globaIndexWeightFromDestToSource;
-  int nb = globalIndexOfCurrentGridSource_.size();
-  const size_t sfmax = NumTraits<unsigned long>::sfmax();
-  for (int idx = 0; idx < nb; ++idx)
-  {
-    if (sfmax != globalIndexOfOriginalGridSource_[idx])
-      globaIndexWeightFromDestToSource[globalIndexOfCurrentGridSource_[idx]].push_back(make_pair(globalIndexOfOriginalGridSource_[idx], weightOfGlobalIndexOfOriginalGridSource_(idx))) ;
-  }
-
-  // Then compute transformation mapping among clients
-  transformationMap.computeTransformationMapping(globaIndexWeightFromDestToSource);
+  transformationMap.computeTransformationMapping(currentGridIndexToOriginalGridIndex_);
 
   const std::map<int,std::vector<std::vector<std::pair<size_t,double> > > >& globalIndexToReceive = transformationMap.getGlobalIndexReceivedOnGridDestMapping();
@@ -667 +806 @@
 }
 
+bool CGridTransformation::isSpecialTransformation(ETranformationType transType)
+{
+  bool res;
+  switch (transType)
+  {
+    case TRANS_GENERATE_RECTILINEAR_DOMAIN:
+     res = true;
+     break;
+    default:
+     res = false;
+     break;
+  }
+
+  return res;
+}
+
 /*!
   Local index of data which need sending from the grid source

XIOS/trunk/src/transformation/grid_transformation.hpp

@@ -63 +63 @@
   void computeTransformationFromOriginalGridSource(const std::map<size_t, std::vector<std::pair<size_t,double> > >& globaIndexMapFromDestToSource);
   void updateFinalGridDestination();
+  bool isSpecialTransformation(ETranformationType transType);
 
 protected:
@@ -77 +78 @@
   //! List of algorithm types and their order
   ListAlgoType listAlgos_;
+
+  //! Number of algorithm
+  int nbAlgos_;
+
+  typedef std::map<size_t, std::vector<std::pair<size_t,double> > > GlobalIndexMap;
+
   // true if domain algorithm and false if axis algorithm (can be replaced by tuple with listAlgos_
   std::vector<bool> algoTypes_;
@@ -96 +103 @@
 
   //! (Grid) Global index of grid source
-  std::vector<size_t> globalIndexOfCurrentGridSource_;
-  std::vector<size_t> globalIndexOfOriginalGridSource_;
-  CArray<double,1> weightOfGlobalIndexOfOriginalGridSource_;
+  GlobalIndexMap currentGridIndexToOriginalGridIndex_;
 };