Changeset 1078 for XIOS/trunk/src/transformation/domain_algorithm_expand.cpp

Timestamp:

03/17/17 16:02:40 (7 years ago)

Author:

mhnguyen

Message:

Adding rectilinear and curvilinear domain for expand_domain transformation
-) Rectilinear/curvilinear is expanded not only locally but also globally, its global size ni_glo, nj_glo become ni_glo+2 and nj_glo+2
-) Two attributes i_periodic, j_periodic are only used for rectilinear/curvilinear to process priodic condition

+) Do some minor modification

Test
+) Add test_connectivity_expand
+) On Curie
+) Work (but need more real tests)

File:

• XIOS/trunk/src/transformation/domain_algorithm_expand.cpp (modified) (9 diffs)

XIOS/trunk/src/transformation/domain_algorithm_expand.cpp

@@ -45 +45 @@
                                                CDomain* domainSource,
                                                CExpandDomain* expandDomain)
-: CDomainAlgorithmTransformation(domainDestination, domainSource)
+: CDomainAlgorithmTransformation(domainDestination, domainSource),
+  isXPeriodic_(false), isYPeriodic_(false)
 {
   if (domainDestination == domainSource)
@@ -55 +56 @@
   }
 
-//  if (!domainDestination->hasRefTo(domainSource))
-//  {
-//    ERROR("CDomainAlgorithmExpand::CDomainAlgorithmExpand(CDomain* domainDestination,CDomain* domainSource, CExpandDomain* expandDomain)",
-//           << "Domain domain destination must refer to domain source (directly or indirectly) by domain_ref " << std::endl
-//           << "Domain source " <<domainSource->getId() << std::endl
-//           << "Domain destination " <<domainDestination->getId() << std::endl);
-//  }
-
   this->type_ = (ELEMENT_MODIFICATION_WITH_DATA);
   expandDomain->checkValid(domainDestination);
+  if (!expandDomain->i_periodic.isEmpty()) isXPeriodic_ = expandDomain->i_periodic;
+  if (!expandDomain->j_periodic.isEmpty()) isYPeriodic_ = expandDomain->j_periodic;
 
   switch (expandDomain->type)
@@ -92 +87 @@
   CContextClient* client=context->client;
 
+  int type = 1; // For edge
+  CMesh mesh;
   CArray<double,2>& bounds_lon_src = domainSource->bounds_lon_1d;
   CArray<double,2>& bounds_lat_src = domainSource->bounds_lat_1d;
   CArray<int,2> neighborsSrc;
-
-  int type = 1; // For edge
-  CMesh mesh;
-  mesh.getGlobalNghbFaces(type, client->intraComm, domainSource->i_index, bounds_lon_src, bounds_lat_src, neighborsSrc);
-  updateDomainAttributes(domainDestination, domainSource, neighborsSrc);
+  switch (domainSource->type) {
+   case CDomain::type_attr::unstructured:      
+      mesh.getGlobalNghbFaces(type, client->intraComm, domainSource->i_index, bounds_lon_src, bounds_lat_src, neighborsSrc);
+      updateUnstructuredDomainAttributes(domainDestination, domainSource, neighborsSrc);
+      break;
+   default:
+      updateRectilinearDomainAttributes(domainDestination, domainSource, neighborsSrc);
+      break;      
+  }  
 }
 
@@ -113 +114 @@
   CContextClient* client=context->client;
 
+  int type = 1; // For edge
+  CMesh mesh;
   CArray<double,2>& bounds_lon_src = domainSource->bounds_lon_1d;
   CArray<double,2>& bounds_lat_src = domainSource->bounds_lat_1d;
   CArray<int,2> neighborsSrc;
-
-  int type = 0; // For node
-  CMesh mesh;
-  mesh.getGlobalNghbFaces(type, client->intraComm, domainSource->i_index, bounds_lon_src, bounds_lat_src, neighborsSrc);
-  updateDomainAttributes(domainDestination, domainSource, neighborsSrc);
+  switch (domainSource->type) {
+   case CDomain::type_attr::unstructured:      
+      mesh.getGlobalNghbFaces(type, client->intraComm, domainSource->i_index, bounds_lon_src, bounds_lat_src, neighborsSrc);
+      updateUnstructuredDomainAttributes(domainDestination, domainSource, neighborsSrc);
+      break;
+   default:
+      updateRectilinearDomainAttributes(domainDestination, domainSource, neighborsSrc);
+      break;      
+  }  
+}
+
+/*!
+ *  Extend rectilinear or curvilinear domain destination and update its attributes
+ *  Suppose that domain destination and domain source have the same values for all attributes (by inheritance)
+ *  \param [in/out] domainDestination domain destination
+ *  \param [in] domainSource domain source
+ *  \param [in] neighborsDomainSrc neighbor of domain source. For now, we don't need it for rectilinear
+ */
+void CDomainAlgorithmExpand::updateRectilinearDomainAttributes(CDomain* domainDestination,
+                                                               CDomain* domainSource,
+                                                               CArray<int,2>& neighborsDomainSrc)
+{
+  int index, globalIndex, idx;
+  int iindexDst, jindexDst, globIndexDst;
+  int iindexSrc, jindexSrc, globIndexSrc;
+  CContext* context = CContext::getCurrent();
+  CContextClient* client=context->client;
+
+  // First of all, "copy" all attributes of domain source to domain destination
+  StdString domainDstRef = (!domainDestination->domain_ref.isEmpty()) ? domainDestination->domain_ref.getValue()
+                                                                      : "";
+  if (domainDstRef != domainSource->getId())
+  {
+    domainDestination->domain_ref.setValue(domainSource->getId());
+    domainDestination->solveRefInheritance(true);
+  }
+
+  if (domainDstRef.empty()) domainDestination->domain_ref.reset();
+  else domainDestination->domain_ref.setValue(domainDstRef);
+
+  // Here are attributes of source need tranfering
+  int niGloSrc = domainSource->ni_glo;
+  int njGloSrc = domainSource->nj_glo;
+  int niSrc = domainSource->ni, ibegin = domainSource->ibegin;
+  int njSrc = domainSource->nj, jbegin = domainSource->jbegin;
+  CArray<bool,1>& mask_1d_src = domainSource->mask_1d;
+  CArray<int,1>& i_index_src = domainSource->i_index;
+  CArray<int,1>& j_index_src = domainSource->j_index;
+  CArray<int,1>& data_i_index_src = domainSource->data_i_index;
+  CArray<int,1>& data_j_index_src = domainSource->data_j_index;
+  int data_i_begin_src = domainSource->data_ibegin;
+  int data_j_begin_src = domainSource->data_jbegin;
+  CArray<double,1>& lon_src = domainSource->lonvalue_client;
+  CArray<double,1>& lat_src = domainSource->latvalue_client;
+
+  // We need to generate boundary for longitude and latitude
+  if (domainSource->bounds_lon_1d.isEmpty() || domainSource->bounds_lat_1d.isEmpty())
+  {
+    CArray<double,1> lon = lon_src(Range(0,niSrc-1));
+    CArray<double,1> lat = lat_src(Range(0,lat_src.numElements()-1,niSrc));
+    CArray<double,2>& bounds_lon_src = domainSource->bounds_lon_1d;
+    CArray<double,2>& bounds_lat_src = domainSource->bounds_lat_1d;
+    domainSource->fillInRectilinearBoundLonLat(lon_src, lat_src, bounds_lon_src, bounds_lat_src);
+  }
+
+
+  CArray<double,2>& bounds_lon_src = domainSource->bounds_lon_1d;
+  CArray<double,2>& bounds_lat_src = domainSource->bounds_lat_1d;
+
+  int nVertex       = bounds_lon_src.shape()[0];
+  int oldNbLocal = i_index_src.numElements();
+  int dataIindexBoundSrc = max(i_index_src) - min(i_index_src);
+  int dataJindexBoundSrc = max(j_index_src) - min(j_index_src);
+
+  // Uncompress data_i_index, data_j_index
+  CArray<int,1> data_i_index_src_full(oldNbLocal);
+  CArray<int,1> data_j_index_src_full(oldNbLocal);
+  int nbUnMaskedPointOnLocalDomain = 0;
+  data_i_index_src_full = -1; // Suppose all values are masked
+  data_j_index_src_full = -1; // Suppose all values are masked
+  for (idx = 0; idx < data_i_index_src.numElements(); ++idx)
+  {
+    int dataIidx = data_i_index_src(idx) + data_i_begin_src;
+    int dataJidx = data_j_index_src(idx) + data_j_begin_src;
+    if ((0 <= dataIidx) && (dataIidx <= dataIindexBoundSrc) &&
+        (0 <= dataJidx) && (dataJidx <= dataJindexBoundSrc))
+    {
+      data_i_index_src_full(nbUnMaskedPointOnLocalDomain) = dataIidx;
+      data_j_index_src_full(nbUnMaskedPointOnLocalDomain) = dataJidx;
+      ++nbUnMaskedPointOnLocalDomain;
+    }
+  }
+
+  // Expand domain destination, not only local but also global
+  int niGloDst = niGloSrc + 2;
+  int njGloDst = njGloSrc + 2;
+  int niDst = niSrc + 2;
+  int njDst = njSrc + 2;
+  domainDestination->ni_glo.setValue(niGloDst);
+  domainDestination->nj_glo.setValue(njGloDst);
+  domainDestination->ni.setValue(niDst);
+  domainDestination->nj.setValue(njDst);
+  domainDestination->global_zoom_ni.setValue(domainSource->global_zoom_ni+2);
+  domainDestination->global_zoom_nj.setValue(domainSource->global_zoom_nj+2);
+
+  CArray<bool,1>& mask_1d_dst = domainDestination->mask_1d;
+  CArray<int,1>& i_index_dst  = domainDestination->i_index;
+  CArray<int,1>& j_index_dst  = domainDestination->j_index;  
+  CArray<int,1>& data_i_index_dst  = domainDestination->data_i_index;
+  CArray<int,1>& data_j_index_dst  = domainDestination->data_j_index;
+  
+  // Make sure that we use only lonvalue_client, latvalue_client
+  if (!domainDestination->lonvalue_1d.isEmpty()) domainDestination->lonvalue_1d.reset();
+  if (!domainDestination->latvalue_1d.isEmpty()) domainDestination->latvalue_1d.reset();
+  if (!domainDestination->lonvalue_2d.isEmpty()) domainDestination->lonvalue_2d.reset();
+  if (!domainDestination->latvalue_2d.isEmpty()) domainDestination->latvalue_2d.reset();
+
+  // Recalculate i_index, j_index of extended domain
+  // Should be enough for common case, but if we have arbitrary distribution?
+  int newNbLocalDst = niDst * njDst;     
+
+  mask_1d_dst.resize(newNbLocalDst);
+  i_index_dst.resize(newNbLocalDst);
+  j_index_dst.resize(newNbLocalDst);
+  CArray<int,1> data_i_index_dst_full(newNbLocalDst, -1);
+  CArray<int,1> data_j_index_dst_full(newNbLocalDst, -1);
+
+  domainDestination->lonvalue_client.resizeAndPreserve(newNbLocalDst);  
+  domainDestination->latvalue_client.resizeAndPreserve(newNbLocalDst);
+  domainDestination->bounds_lon_1d.resizeAndPreserve(nVertex, newNbLocalDst);
+  domainDestination->bounds_lat_1d.resizeAndPreserve(nVertex, newNbLocalDst);
+  CArray<double,1>& lon_dst   = domainDestination->lonvalue_client;
+  CArray<double,1>& lat_dst   = domainDestination->latvalue_client;
+  CArray<double,2>& bounds_lon_dst = domainDestination->bounds_lon_1d;
+  CArray<double,2>& bounds_lat_dst = domainDestination->bounds_lat_1d;
+
+  // Update i_index, j_index
+  for (int j = 0; j < njDst; ++j)
+    for (int i = 0; i < niDst; ++i)
+    {
+      idx = j * niDst + i; 
+      i_index_dst(idx) = i + ibegin;
+      j_index_dst(idx) = j + jbegin;
+    }
+
+ 
+
+  // 1. Fill in array relating to global index (i_index, j_index, transmap, etc, ...)
+  // Global index mapping between destination and source
+  this->transformationMapping_.resize(1);
+  this->transformationWeight_.resize(1);
+  TransformationIndexMap& transMap = this->transformationMapping_[0];
+  TransformationWeightMap& transWeight = this->transformationWeight_[0];
+
+  transMap.rehash(std::ceil(newNbLocalDst/transMap.max_load_factor()));
+  transWeight.rehash(std::ceil(newNbLocalDst/transWeight.max_load_factor()));
+  
+  // Index mapping for local domain
+  // Mapping global index of expanded domain into original one 
+  // (Representing global index of expanded domain in form of global index of original one)
+  CArray<size_t,1> globalIndexSrcOnDstDomain(newNbLocalDst); 
+  for (idx = 0; idx < newNbLocalDst; ++idx)
+  {
+    iindexDst = i_index_dst(idx);
+    jindexDst = j_index_dst(idx);
+    globIndexDst = jindexDst * niGloDst + iindexDst;
+    globIndexSrc = (((jindexDst-1)+njGloSrc) % njGloSrc) * niGloSrc + (((iindexDst-1)+niGloSrc) % niGloSrc) ;
+    globalIndexSrcOnDstDomain(idx) = globIndexSrc;
+
+    transMap[globIndexDst].push_back(globIndexSrc);
+    transWeight[globIndexDst].push_back(1.0); 
+  }
+
+  // 2. Exchange local info among domains (lon,lat,bounds,mask,etc,...)
+  CClientClientDHTDouble::Index2VectorInfoTypeMap localData;
+  localData.rehash(std::ceil(oldNbLocal/localData.max_load_factor()));
+
+  // Information exchanged among domains (attention to their order), number in parentheses presents size of data
+  // lon(1) + lat(1) + bounds_lon(nVertex) + bounds_lat(nVertex) + mask(1) + data_i_index(1)
+  int dataPackageSize =  1 + 1 + // lon + lat
+                         nVertex + nVertex + //bounds_lon + bounds_lat
+                         1 + // mask_1d_dst;
+                         1 + 1; // data_i_index + data_j_index
+  // Initialize database
+  for (int idx = 0; idx < oldNbLocal; ++idx)
+  {
+    index = i_index_src(idx) + j_index_src(idx) * niGloSrc;
+    localData[index].resize(dataPackageSize);
+    std::vector<double>& data = localData[index];
+
+    //Pack data
+    int dataIdx = 0;
+    data[dataIdx] = lon_src(idx);++dataIdx;
+    data[dataIdx] = lat_src(idx);++dataIdx;
+    for (int i = 0; i < nVertex; ++i)
+    {
+      data[dataIdx] = bounds_lon_src(i,idx); ++dataIdx;
+    }
+    for (int i = 0; i < nVertex; ++i)
+    {
+      data[dataIdx] = bounds_lat_src(i,idx); ++dataIdx;
+    }
+    data[dataIdx] = mask_1d_src(idx) ? 1.0 : -1.0; ++dataIdx;
+    data[dataIdx] = data_i_index_src_full(idx);++dataIdx;
+    data[dataIdx] = data_j_index_src_full(idx);
+  }
+
+  CClientClientDHTDouble dhtData(localData,client->intraComm);
+  dhtData.computeIndexInfoMapping(globalIndexSrcOnDstDomain);
+  CClientClientDHTDouble::Index2VectorInfoTypeMap& neighborData = dhtData.getInfoIndexMap();
+  CClientClientDHTDouble::Index2VectorInfoTypeMap::iterator ite = neighborData.end(), it;
+
+  // Ok get all data for destination
+  // If domain is not periodic, then we mask all extended part.
+  int nbUnMaskedPointOnExtendedPart = 0, remainder = 0, dataIIndex, dataJIndex;
+  size_t nIdx;
+  double maskValue = 1.0;
+  for (index = 0; index < newNbLocalDst; ++index)
+  {
+     nIdx = globalIndexSrcOnDstDomain(index);
+     it = neighborData.find(nIdx);
+     if (ite != it)
+     {
+        std::vector<double>& data = it->second;
+        // Unpack data
+        int dataIdx = 0;
+        lon_dst(index) = data[dataIdx]; ++dataIdx;
+        lat_dst(index) = data[dataIdx]; ++dataIdx;
+        for (int i = 0; i < nVertex; ++i)
+        {
+          bounds_lon_dst(i,index) = data[dataIdx]; ++dataIdx;
+        }
+        for (int i = 0; i < nVertex; ++i)
+        {
+          bounds_lat_dst(i,index) = data[dataIdx]; ++dataIdx;
+        }
+        
+        // Check whether we have x periodic. If we don't, we should mask all point at 0 and niGloDst-1
+        maskValue = data[dataIdx];
+        if (!isXPeriodic_) 
+        {
+          remainder = i_index_dst(index) % (niGloDst-1);
+          if (0 == remainder) 
+          {
+            maskValue = -1.0;
+          }
+        }
+
+        if (!isYPeriodic_) 
+        {
+          remainder = j_index_dst(index) % (njGloDst-1);
+          if (0 == remainder) 
+          {
+            maskValue = -1.0;
+          }
+        }
+
+        mask_1d_dst(index) = (1.0 == maskValue) ? true : false; ++dataIdx;
+
+        dataIIndex = (int) data[dataIdx];
+        if (!isXPeriodic_) 
+        {
+          remainder = i_index_dst(index) % (niGloDst-1);
+          if (0 == remainder) 
+          {
+            dataIIndex = -1;
+          }
+        }
+        data_i_index_dst_full(index) = dataIIndex; ++dataIdx;
+        
+        dataJIndex = (int) data[dataIdx];
+        if (!isYPeriodic_) 
+        {
+          remainder = j_index_dst(index) % (njGloDst-1);
+          if (0 == remainder) 
+          {
+            dataJIndex = -1;
+          }
+        }        
+        data_j_index_dst_full(index) = dataJIndex;
+
+        if ((0 <= data_i_index_dst_full(index)) && (0 <= data_j_index_dst_full(index)))
+        {
+          ++nbUnMaskedPointOnExtendedPart;
+        }
+     }
+  }
+
+  
+  // Finally, update data_i_index, data_j_index
+  data_i_index_dst.resize(nbUnMaskedPointOnExtendedPart);
+  data_j_index_dst.resize(nbUnMaskedPointOnExtendedPart); 
+  int count = 0; 
+  for (idx = 0; idx < newNbLocalDst; ++idx)
+  {
+    dataIIndex = data_i_index_dst_full(idx);
+    dataJIndex = data_j_index_dst_full(idx);
+    if ((0 <= dataIIndex) && (0 <= dataJIndex))
+    {
+      data_i_index_dst(count) = i_index_dst(idx) - i_index_dst(0);
+      data_j_index_dst(count) = j_index_dst(idx) - j_index_dst(0);
+      ++count;
+    }
+  }
+
+  // Update data_ni, data_nj
+  domainDestination->data_ni.setValue(niDst);
+  domainDestination->data_nj.setValue(njDst);
+  domainDestination->data_ibegin.setValue(0);
+  domainDestination->data_jbegin.setValue(0);
+
+
+  // Update longitude and latitude 
+  if (niSrc == domainSource->lonvalue_1d.numElements() && njSrc == domainSource->latvalue_1d.numElements()) // Ok, we have rectilinear here
+  {
+     domainDestination->lonvalue_1d.resize(niDst);
+     domainDestination->lonvalue_1d = lon_dst(Range(0,niDst-1));
+     domainDestination->latvalue_1d.resize(njDst);
+     domainDestination->latvalue_1d = lat_dst(Range(0,lat_dst.numElements()-1,niDst));
+  }
+  else // It should be curvilinear
+  {
+     domainDestination->lonvalue_1d.resize(lon_dst.numElements());
+     domainDestination->lonvalue_1d = lon_dst;
+     domainDestination->latvalue_1d.resize(lat_dst.numElements());
+     domainDestination->latvalue_1d = (lat_dst);
+  }
+
 }
 
@@ -130 +456 @@
  *  \param [in] neighborsDomainSrc domain extended part
  */
-void CDomainAlgorithmExpand::updateDomainAttributes(CDomain* domainDestination,
-                                                    CDomain* domainSource,
-                                                    CArray<int,2>& neighborsDomainSrc)
-{
+void CDomainAlgorithmExpand::updateUnstructuredDomainAttributes(CDomain* domainDestination,
+                                                                CDomain* domainSource,
+                                                                CArray<int,2>& neighborsDomainSrc)
+{
+
   CContext* context = CContext::getCurrent();
   CContextClient* client=context->client;
@@ -255 +582 @@
       data[dataIdx] = bounds_lat_src(i,idx); ++dataIdx;
     }
-    data[dataIdx] = mask_1d_src(idx) ? 1.0 : -1; ++dataIdx;
+    data[dataIdx] = mask_1d_src(idx) ? 1.0 : -1.0; ++dataIdx;
     data[dataIdx] = data_i_index_src_full(idx);
   }
@@ -308 +635 @@
   {
     dataIdx = data_i_index_dst_full(idx);
-    if ((0 <= dataIdx) && (dataIdx < newNbLocalDst))
+    if ((0 <= dataIdx))
     {
       ++count;
@@ -322 +649 @@
   {
     dataIdx = data_i_index_dst_full(idx);
-    if ((0 <= dataIdx) && (dataIdx < newNbLocalDst))
+    if ((0 <= dataIdx))
     {
       data_i_index_dst(count) = dataIdx;
@@ -331 +658 @@
   // Update ni
   domainDestination->ni.setValue(newNbLocalDst);
-
 }