Changeset 1399 for XIOS/dev/XIOS_DEV_CMIP6

Timestamp:

01/10/18 16:17:57 (6 years ago)

Author:

ymipsl

Message:

2 update for spatial tranformations :

• One bug fix introduce in r1389 concerning domain interpolation or other distributed spatial distransfomation

• Specific optimisation for transformation on element wich is not distributed across processes

This optimisation is experimental and will be activated only when using a variable in xios context :

<variable id="activate_non_distributed_transformation" type="bool">true</variable>

Later when more succesfull test, the variable use will removed...

YM

Location:

XIOS/dev/XIOS_DEV_CMIP6/src/transformation

Files:

• generic_algorithm_transformation.cpp (modified) (8 diffs)
• generic_algorithm_transformation.hpp (modified) (4 diffs)
• grid_transformation.cpp (modified) (2 diffs)

XIOS/dev/XIOS_DEV_CMIP6/src/transformation/generic_algorithm_transformation.cpp

@@ -13 +13 @@
 #include "utils.hpp"
 #include "timer.hpp"
+#include "mpi.hpp"
 
 namespace xios {
@@ -19 +20 @@
  : transformationMapping_(), transformationWeight_(), transformationPosition_(),
    idAuxInputs_(), type_(ELEMENT_NO_MODIFICATION_WITH_DATA), indexElementSrc_(),
-   computedProcSrcNonTransformedElement_(false), eliminateRedondantSrc_(true)
+   computedProcSrcNonTransformedElement_(false), eliminateRedondantSrc_(true), isDistributedComputed_(false)
 {
 }
@@ -114 +115 @@
 }
 
+bool CGenericAlgorithmTransformation::isDistributedTransformation(int elementPositionInGrid, CGrid* gridSrc, CGrid* gridDst)
+{
+
+  if (!isDistributedComputed_)
+  {
+    CContext* context = CContext::getCurrent();
+    CContextClient* client = context->client;
+  
+    computePositionElements(gridSrc, gridDst);
+    std::vector<CScalar*> scalarListSrcP  = gridSrc->getScalars();
+    std::vector<CAxis*> axisListSrcP = gridSrc->getAxis();
+    std::vector<CDomain*> domainListSrcP = gridSrc->getDomains();
+    int distributed, distributed_glo ;
+  
+    CArray<int,1> axisDomainSrcOrder = gridSrc->axis_domain_order;
+    if (2 == axisDomainSrcOrder(elementPositionInGrid)) // It's domain
+    {
+      distributed=domainListSrcP[elementPositionInGridSrc2DomainPosition_[elementPositionInGrid]]->isDistributed() ;
+      MPI_Allreduce(&distributed,&distributed_glo, 1, MPI_INT, MPI_LOR, client->intraComm) ;
+    
+    }
+    else if (1 == axisDomainSrcOrder(elementPositionInGrid))//it's an axis
+    {
+      distributed=axisListSrcP[elementPositionInGridSrc2AxisPosition_[elementPositionInGrid]]->isDistributed() ;
+      MPI_Allreduce(&distributed,&distributed_glo, 1, MPI_INT, MPI_LOR, client->intraComm) ;
+    }
+    else //it's a scalar
+    {
+      distributed_glo=false ;
+    } 
+    isDistributed_=distributed_glo ;
+    isDistributedComputed_=true ;
+  }
+  return isDistributed_ ;
+}  
 /*!
   This function computes the global indexes of grid source, which the grid destination is in demand.
@@ -200 +236 @@
 
   CClientClientDHTInt::Index2VectorInfoTypeMap globalIndexOfTransformedElementOnProc;
-  
+
   if (computeGlobalIndexOnProc || !computedProcSrcNonTransformedElement_)
   {    
@@ -294 +330 @@
     procContainSrcElementIdx_.swap(procContainSrcElementIdx);    
 */
-    std::set<int> commonProc ;
+    
     if (procElementList_.empty()) procElementList_.resize(axisDomainDstOrder.numElements()) ;
     for (idx = 0; idx < axisDomainDstOrder.numElements(); ++idx)
@@ -330 +366 @@
       }
 
-      if (idx==0) commonProc.swap(procList) ;
+      if (idx==0) commonProc_.swap(procList) ;
       else
       {
-        for (std::set<int>::iterator it = commonProc.begin(); it != commonProc.end(); ++it)
+        for (std::set<int>::iterator it = commonProc_.begin(); it != commonProc_.end(); ++it)
           if (procList.count(*it)==1) commonTmp.insert(*it) ;
-        commonProc.swap(commonTmp) ;
-      }
-    }
-    std::vector<int> procContainSrcElementIdx(commonProc.size());
+        commonProc_.swap(commonTmp) ;
+      }
+    }
+    std::vector<int> procContainSrcElementIdx(commonProc_.size());
     int count = 0;
-    for (std::set<int>::iterator it = commonProc.begin(); it != commonProc.end(); ++it) procContainSrcElementIdx[count++] = *it;
+    for (std::set<int>::iterator it = commonProc_.begin(); it != commonProc_.end(); ++it) procContainSrcElementIdx[count++] = *it;
     procContainSrcElementIdx_.swap(procContainSrcElementIdx);
     
@@ -370 +406 @@
           tmpCounter.insert(srcIndex[idx]);
        
-          for (int j = 0; j < procContainSrcElementIdx_.size(); ++j)
-            globalElementIndexOnProc_[elementPositionInGrid][procContainSrcElementIdx_[j]].push_back(srcIndex[idx]);
+          vector<int>& rankSrc = globalIndexOfTransformedElementOnProc[srcIndex[idx]] ;
+          for (int n=0;n<rankSrc.size();++n)
+          {
+            if (commonProc_.count(rankSrc[n])==1) globalElementIndexOnProc_[elementPositionInGrid][rankSrc[n]].push_back(srcIndex[idx]);
+          }
+//          for (int j = 0; j < procContainSrcElementIdx_.size(); ++j)
+//            globalElementIndexOnProc_[elementPositionInGrid][procContainSrcElementIdx_[j]].push_back(srcIndex[idx]);
         }
       }
@@ -787 +828 @@
 }
 
+
+
+
+
+
+void CGenericAlgorithmTransformation::computeTransformationMappingNonDistributed(int elementPositionInGrid, CGrid* gridSrc, CGrid* gridDst,
+                                                                                 vector<int>& localSrc, vector<int>& localDst, vector<double>& weight, vector<bool>& localMaskOnGridDest)
+{
+
+  CContext* context = CContext::getCurrent();
+  CContextClient* client = context->client;
+  int nbClient = client->clientSize;
+
+  computePositionElements(gridDst, gridSrc);
+  std::vector<CScalar*> scalarListDstP = gridDst->getScalars();
+  std::vector<CAxis*> axisListDstP = gridDst->getAxis();
+  std::vector<CDomain*> domainListDstP = gridDst->getDomains();
+  CArray<int,1> axisDomainDstOrder = gridDst->axis_domain_order;
+  std::vector<CScalar*> scalarListSrcP  = gridSrc->getScalars();
+  std::vector<CAxis*> axisListSrcP = gridSrc->getAxis();
+  std::vector<CDomain*> domainListSrcP = gridSrc->getDomains();
+  CArray<int,1> axisDomainSrcOrder = gridSrc->axis_domain_order;  
+
+  int nElement=axisDomainSrcOrder.numElements() ;
+  int indSrc=1 ;
+  int indDst=1 ;
+  vector<int> nIndexSrc(nElement) ;
+  vector<int> nIndexDst(nElement) ;
+  int nlocalIndexDest=1 ;
+  
+  for(int i=0 ; i<nElement; i++)
+  {
+    int dimElement = axisDomainSrcOrder(i);
+    if (2 == dimElement) //domain
+    {
+      CDomain* domain=domainListSrcP[elementPositionInGridSrc2DomainPosition_[i]] ;
+      nIndexSrc[i] = domain->i_index.numElements() ;
+    }
+    else if (1 == dimElement) //axis
+    {
+      CAxis* axis=axisListSrcP[elementPositionInGridSrc2DomainPosition_[i]] ;
+      nIndexSrc[i] = axis->index.numElements() ;
+    }
+    else  //scalar
+    {
+      nIndexSrc[i]=1 ;
+    }
+  }
+
+  int offset=1 ;
+  for(int i=0 ; i<nElement; i++)
+  {
+    int dimElement = axisDomainDstOrder(i);
+    if (2 == dimElement) //domain
+    {
+      CDomain* domain=domainListDstP[elementPositionInGridDst2DomainPosition_[i]] ;
+      nIndexDst[i] = domain->i_index.numElements() ;
+    }
+    else if (1 == dimElement) //axis
+    {
+      CAxis* axis = axisListDstP[elementPositionInGridDst2DomainPosition_[i]] ;
+      nIndexDst[i] = axis->index.numElements() ;
+    }
+    else  //scalar
+    {
+      nIndexDst[i]=1 ;
+    }
+    if (i<elementPositionInGrid) offset=offset*nIndexDst[i] ;
+    nlocalIndexDest=nlocalIndexDest*nIndexDst[i] ;
+  }
+
+  vector<int> dstLocalInd ;
+  int dimElement = axisDomainDstOrder(elementPositionInGrid);
+  if (2 == dimElement) //domain
+  {
+    CDomain* domain = domainListSrcP[elementPositionInGridSrc2DomainPosition_[elementPositionInGrid]] ;
+    int ni_glo=domain->ni_glo ;
+    int nj_glo=domain->nj_glo ;
+    int nindex_glo=ni_glo*nj_glo ;
+    dstLocalInd.resize(nindex_glo,-1) ;
+    int nIndex=domain->i_index.numElements() ;
+    for(int i=0;i<nIndex;i++)
+    {
+      int globIndex=domain->j_index(i)*ni_glo+domain->i_index(i) ;
+      dstLocalInd[globIndex]=i ;
+    }
+  }
+  else if (1 == dimElement) //axis
+  {
+    CAxis* axis = axisListSrcP[elementPositionInGridSrc2AxisPosition_[elementPositionInGrid]] ;
+    int nindex_glo=axis->n_glo ;
+    dstLocalInd.resize(nindex_glo,-1) ;
+    int nIndex=axis->index.numElements() ;
+    for(int i=0;i<nIndex;i++)
+    {
+      dstLocalInd[axis->index(i)]=i ;
+    }
+  }
+  else  //scalar
+  {
+    dstLocalInd.resize(1) ; 
+    dstLocalInd[0]=0 ; 
+  }
+
+// just get the local src mask
+  CArray<bool,1> localMaskOnSrcGrid;
+  gridSrc->getLocalMask(localMaskOnSrcGrid) ; 
+  
+
+  localMaskOnGridDest.resize(nlocalIndexDest,false) ;
+
+  vector<vector<vector<pair<int,double> > > > dstIndWeight(transformationMapping_.size()) ;
+   
+  for(int t=0;t<transformationMapping_.size();++t)
+  {
+    TransformationIndexMap::const_iterator   itTransMap = transformationMapping_[t].begin(),
+                                             iteTransMap = transformationMapping_[t].end();
+    TransformationWeightMap::const_iterator itTransWeight = transformationWeight_[t].begin();
+    dstIndWeight[t].resize(nIndexSrc[elementPositionInGrid]) ;
+    
+    for(;itTransMap!=iteTransMap;++itTransMap,++itTransWeight)
+    {
+      int dst=dstLocalInd[itTransMap->first] ;
+      if (dst!=-1)
+      {
+        const vector<int>& srcs=itTransMap->second;
+        const vector<double>& weights=itTransWeight->second;
+        for(int i=0;i<srcs.size() ;i++) dstIndWeight[t][srcs[i]].push_back(pair<int,double>(dst*offset+t,weights[i])) ;
+      }
+    }
+  }
+  int srcInd=0 ;  
+  int currentInd ;
+  int t=0 ;  
+
+  
+  nonDistributedrecursiveFunct(nElement-1,elementPositionInGrid,srcInd, nIndexSrc, t, dstIndWeight,  
+                               currentInd,localSrc,localDst,weight, localMaskOnSrcGrid, localMaskOnGridDest );
+               
+}
+
+
+void CGenericAlgorithmTransformation::nonDistributedrecursiveFunct(int currentPos, int elementPositionInGrid, int& srcInd, vector<int>& nIndexSrc, int& t, vector<vector<vector<pair<int,double> > > >& dstIndWeight, int currentInd,
+                    vector<int>& localSrc, vector<int>& localDst, vector<double>& weight,  CArray<bool,1>& localMaskOnGridSrc, vector<bool>& localMaskOnGridDest )
+{
+  if (currentPos!=elementPositionInGrid)
+  {
+    if (currentPos!=0)
+    {
+      for(int i=0;i<nIndexSrc[currentPos];i++)
+      {
+        nonDistributedrecursiveFunct(currentPos-1, elementPositionInGrid, srcInd, nIndexSrc, t, dstIndWeight, currentInd, localSrc, localDst, weight, localMaskOnGridSrc, localMaskOnGridDest) ;
+      }
+    }
+    else
+    {
+      for(int i=0;i<nIndexSrc[currentPos];i++)
+      {
+        if (dstIndWeight[t][currentInd].size()>0)
+        {
+          for(vector<pair<int,double> >::iterator it = dstIndWeight[t][currentInd].begin(); it!=dstIndWeight[t][currentInd].end(); ++it)
+          {
+            if (localMaskOnGridSrc(srcInd))
+            {
+              localSrc.push_back(srcInd) ;
+              localDst.push_back(it->first) ;
+              weight.push_back(it->second) ;
+              localMaskOnGridDest[it->first]=true ;
+            }
+            (it->first)++ ;
+          }
+        }
+        srcInd++ ;
+        if (t < dstIndWeight.size()-1) t++ ;
+      }
+    }
+  }
+  else
+  {
+ 
+    if (currentPos!=0)
+    {
+
+
+      for(int i=0;i<nIndexSrc[currentPos];i++)
+      {
+        t=0 ;
+        nonDistributedrecursiveFunct(currentPos-1, elementPositionInGrid, srcInd, nIndexSrc, t, dstIndWeight , i,  localSrc, localDst, weight, localMaskOnGridSrc, localMaskOnGridDest) ;
+      }
+    }
+    else
+    {
+      for(int i=0;i<nIndexSrc[currentPos];i++)
+      {
+        t=0 ;        
+        if (dstIndWeight[t][i].size()>0)
+        {
+          for(vector<pair<int,double> >::iterator it = dstIndWeight[t][i].begin(); it!=dstIndWeight[t][i].end(); ++it)
+          {
+            if (localMaskOnGridSrc(srcInd))
+            {
+              localSrc.push_back(srcInd) ;
+              localDst.push_back(it->first) ;
+              weight.push_back(it->second) ;
+              localMaskOnGridDest[it->first]=true ;
+            }
+            (it->first)++ ;
+          }
+        }
+        srcInd++ ;
+        if (t < dstIndWeight.size()-1) t++ ;
+      }
+    }
+  }
+
+}
+
+
+
+
+
+
+
+
+
+
+
 /*!
   Compute index mapping between element source and element destination with an auxiliary inputs which determine

XIOS/dev/XIOS_DEV_CMIP6/src/transformation/generic_algorithm_transformation.hpp

@@ -52 +52 @@
   virtual ~CGenericAlgorithmTransformation() {}
 
+  bool isDistributedTransformation(int elementPositionInGrid, CGrid* gridSrc, CGrid* gridDst) ;
+
   void computeGlobalSourceIndex(int elementPositionInGrid,
                                CGrid* gridSrc,
@@ -85 +87 @@
   */
   void computeIndexSourceMapping(const std::vector<CArray<double,1>* >& dataAuxInputs = std::vector<CArray<double,1>* >());
+  void computeTransformationMappingNonDistributed(int elementPositionInGrid, CGrid* gridSrc, CGrid* gridDst,
+                                                  vector<int>& localSrc, vector<int>& localDst, vector<double>& weight, vector<bool>& localMaskOnGridDest);
+  void nonDistributedrecursiveFunct(int currentPos, int elementPositionInGrid, int& srcInd, vector<int>& nIndexSrc, int& t, vector<vector<vector<pair<int,double> > > >& dstIndWeight, int currentInd,
+                                     vector<int>& localSrc, vector<int>& localDst, vector<double>& weight, CArray<bool,1>& localMaskOnGridSrc, vector<bool>& localMaskOnGridDest) ;
 
 protected:
@@ -126 +132 @@
 
 protected:
+  //! indicate if the transformation is performed on a distributed element
+  bool isDistributed_ ;
+  //! indicate if the method  isDistributedTransformation has been called before
+  bool isDistributedComputed_ ;
+  
   //! Map between global index of destination element and source element
   std::vector<TransformationIndexMap> transformationMapping_;
@@ -144 +155 @@
   std::vector<int> procContainSrcElementIdx_;  // List of processes containing source index of transformed elements
 //  std::set<int> procOfNonTransformedElements_; // Processes contain the source index of non-transformed elements
+  std::set<int> commonProc_;
   std::vector< set<int> > procElementList_ ; // List of processes containing source index of elements
 

XIOS/dev/XIOS_DEV_CMIP6/src/transformation/grid_transformation.cpp

@@ -16 +16 @@
 #include "grid.hpp"
 #include <boost/unordered_map.hpp>
+#include "timer.hpp"
 
 namespace xios {
@@ -389 +390 @@
       // ComputeTransformation of global index of each element
       int elementPosition = it->first;
-      algo->computeGlobalSourceIndex(elementPosition,
-                                     gridSource_,
-                                     tmpGridDestination_,
-                                     globaIndexWeightFromSrcToDst);
-
-      // Compute transformation of global indexes among grids
-      computeTransformationMapping(globaIndexWeightFromSrcToDst);
-
+      bool nonDistributedActivated = CXios::getin<bool>("activate_non_distributed_transformation",false);
+      
+      if (nonDistributedActivated && !algo->isDistributedTransformation(elementPositionInGrid, gridSource_, tmpGridDestination_) )
+      {
+        vector<int> localSrc ;
+        vector<int> localDst ;
+        vector<double> weight ;
+        localMaskOnGridDest_.push_back(vector<bool>()) ;
+        CTimer::get("computeTransformationMappingNonDistributed").resume();  
+        algo->computeTransformationMappingNonDistributed(elementPosition, gridSource_, tmpGridDestination_, 
+                                                         localSrc, localDst, weight, localMaskOnGridDest_.back()) ;
+        CTimer::get("computeTransformationMappingNonDistributed").suspend();  
+
+        CTimer::get("computeTransformationMappingConvert").resume();  
+        nbLocalIndexOnGridDest_.push_back(localMaskOnGridDest_.back().size()) ;
+        int clientRank=client->clientRank ;
+        {
+          SendingIndexGridSourceMap tmp;
+          localIndexToSendFromGridSource_.push_back(tmp) ;
+          SendingIndexGridSourceMap& src=localIndexToSendFromGridSource_.back() ;
+          CArray<int,1> arrayTmp ;
+          src.insert( pair<int,CArray<int,1> >(clientRank,arrayTmp)) ;
+          CArray<int,1>& array=src[clientRank] ;
+          array.resize(localSrc.size()) ;
+          for(int i=0;i< localSrc.size();++i) array(i)=localSrc[i]  ;    
+        }
+        {
+          RecvIndexGridDestinationMap tmp;
+          localIndexToReceiveOnGridDest_.push_back(tmp) ;
+          RecvIndexGridDestinationMap& dst=localIndexToReceiveOnGridDest_.back() ;
+          vector<pair<int,double> > vectTmp ;
+          dst.insert( pair<int,vector<pair<int,double> >  >(clientRank,vectTmp)) ;
+          vector<pair<int,double> >& vect=dst[clientRank] ;
+          vect.resize(localDst.size()) ;
+          for(int i=0;i< localDst.size();++i) vect[i]=pair<int,double>(localDst[i],weight[i])  ;      
+        }
+        CTimer::get("computeTransformationMappingConvert").suspend();  
+      }
+      else
+      {
+        CTimer::get("computeGlobalSourceIndex").resume();           
+        algo->computeGlobalSourceIndex(elementPosition,
+                                       gridSource_,
+                                       tmpGridDestination_,
+                                       globaIndexWeightFromSrcToDst);
+                                     
+        CTimer::get("computeGlobalSourceIndex").suspend();           
+        CTimer::get("computeTransformationMapping").resume();     
+        // Compute transformation of global indexes among grids
+        computeTransformationMapping(globaIndexWeightFromSrcToDst);
+        CTimer::get("computeTransformationMapping").suspend(); 
+      } 
       if (1 < nbNormalAlgos_)
       {