Changeset 941

Timestamp:

09/22/16 10:59:17 (8 years ago)

Author:

mhnguyen

Message:

Finishing the implementation of expand domain transformation

+) Make use of updated new functions of class Mesh to compute neighboring cells
+) Make change to some minor stuffs

Test
+) On Curie
+) The transformation works correctly

Location:

XIOS/trunk/src

Files:

• client_client_dht_decl.cpp (modified) (1 diff)
• client_client_dht_template.hpp (modified) (1 diff)
• config/expand_domain_attribute.conf (modified) (1 diff)
• declare_ref_func.hpp (modified) (3 diffs)
• node/expand_domain.cpp (modified) (1 diff)
• node/mesh.cpp (modified) (1 diff, 1 prop)
• transformation/algo_types.hpp (modified) (1 diff)
• transformation/domain_algorithm_expand.cpp (modified) (3 diffs)
• transformation/domain_algorithm_expand.hpp (modified) (1 diff)
• transformation/grid_transformation.cpp (modified) (1 diff)
• transformation/grid_transformation_selector.cpp (modified) (1 diff)

XIOS/trunk/src/client_client_dht_decl.cpp

@@ -14 +14 @@
 template class CClientClientDHTTemplate<int>;
 template class CClientClientDHTTemplate<size_t>;
+template class CClientClientDHTTemplate<double>;
 template class CClientClientDHTTemplate<PairIntInt>;
 

XIOS/trunk/src/client_client_dht_template.hpp

@@ -126 +126 @@
 typedef CClientClientDHTTemplate<int> CClientClientDHTInt;
 typedef CClientClientDHTTemplate<size_t> CClientClientDHTSizet;
+typedef CClientClientDHTTemplate<double> CClientClientDHTDouble;
 typedef CClientClientDHTTemplate<PairIntInt> CClientClientDHTPairIntInt;
 

XIOS/trunk/src/config/expand_domain_attribute.conf

@@ -1 +1 @@
 DECLARE_ATTRIBUTE(int, order)
+
+DECLARE_ENUM2(type,node,edge)

XIOS/trunk/src/declare_ref_func.hpp

@@ -18 +18 @@
   bool hasDirect##type##Reference(void) const;              \
   C##type* getDirect##type##Reference(void) const;          \
-  const StdString& get##type##OutputName(void) const;        \
+  const StdString& get##type##OutputName(void) const;       \
   void setAttributesReference(bool apply = true);           \
+  bool hasRefTo(C##type* ref) const;                        \
                                                                        \
 private:                                                               \
-  std::vector<C##type*> refObjects;                                    \
-  StdString outputName_;
+  std::vector<C##type*> refObjects;
 
 // Definitions
@@ -32 +32 @@
   std::set<C##type*> tmpRefObjects;                                    \
   C##type* refer_ptr = this;                                           \
+  std::vector<C##type*>().swap(refObjects);                            \
   refObjects.push_back(this);                                          \
                                                                        \
@@ -107 +108 @@
     return getId();                                                    \
 }                                                                      \
+                                                                       \
+bool C##type::hasRefTo(C##type* ref) const                             \
+{                                                                      \
+  bool found = false;                                                  \
+  for (int idx = 0; idx < refObjects.size(); ++idx)                    \
+    if (ref == refObjects[idx]) { found = true; break; }                \
+  return found;                                                        \
+}
 
 #endif // __XIOS_DECLARE_REF_FUNC_HPP__

XIOS/trunk/src/node/expand_domain.cpp

@@ -38 +38 @@
   void CExpandDomain::checkValid(CDomain* domainDst)
   {
+    if (CDomain::type_attr::unstructured != domainDst->type)
+    {
+      ERROR("CExpandDomain::checkValid(CDomain* domainDst)",
+            << "Domain extension is only supported for unstructured" << std::endl
+            << "Check type of domain destination, id = " << domainDst->getId());
+    }
+
+    if (this->type.isEmpty()) this->type.setValue(CExpandDomain::type_attr::edge);
   }
 

XIOS/trunk/src/node/mesh.cpp

@@ -1923 +1923 @@
     nbNghbFaces.resize(nbFaces);
     nbNghbFaces = 0;
-
-    // nodeToFaces connectivity
-    CClientClientDHTSizet::Index2VectorInfoTypeMap nodeToFaces;
+    int maxNb = 10;  // assumption on max number of edges sharing the same node
+
+    // faceToNodes connectivity
+    CArray<int, 2> nodeToFaces(maxNb, nbFaces*nvertex);
+    nodeToFaces = -1;
+    boost::unordered_map <pair<double,double>, int> mapNodes;
+    vector<int> countNodes(nbFaces*nvertex);
+    countNodes.assign(nbFaces*nvertex, 0);
+
     for (int nf = 0; nf < nbFaces; ++nf)
       for (int nv = 0; nv < nvertex; ++nv)
       {
-        size_t nodeHash = (CMesh::createHashes(bounds_lon(nv, nf), bounds_lat(nv ,nf)))[0];
-        nodeToFaces[nodeHash].push_back(face_idx(nf));
+        if (mapNodes.find(make_pair (bounds_lon(nv, nf), bounds_lat(nv ,nf))) == mapNodes.end())
+        {
+          mapNodes[make_pair (bounds_lon(nv, nf), bounds_lat(nv, nf))] = nbNodes;
+          nodeToFaces(countNodes[nbNodes], nbNodes) = face_idx(nf);
+          ++(countNodes[nbNodes]);
+          ++nbNodes;
+        }
+        else
+        {
+          int node = mapNodes[make_pair (bounds_lon(nv, nf), bounds_lat(nv ,nf))];
+          nodeToFaces(countNodes[node], node) = face_idx(nf);
+          ++(countNodes[node]);
+        }
       }
+    nodeToFaces.resizeAndPreserve(maxNb, nbNodes);
+
+    vector<int> tmpVec (nbNodes);
+    for (int i = 0; i < nbNodes; ++i)
+    {
+      for (int j = 0; j < maxNb; ++j)
+        tmpVec[j] = nodeToFaces(j,i);
+    }
+
 
     // faceToFaces connectivity
     boost::unordered_map <int, int> mapFaces;  // mapFaces = < hash(face1, face2), hash> (the mapped value is irrelevant)
-    int maxNb = 20;                            // some assumption on the max possible number of neighboring cells
-    faceToFaces.resize(maxNb, nbFaces);
-    CClientClientDHTSizet::Index2VectorInfoTypeMap::iterator it;
-    for (it = nodeToFaces.begin(); it != nodeToFaces.end(); ++it)
-    {
-      int size = it->second.size();
-      for (int i = 0; i < (size-1); ++i)
-      {
-        int face1 = it->second[i];
-        for (int j = i+1; j < size; ++j)
-        {
-          int face2 = it->second[j];
-          if (face2 != face1)
-          {
-            int hashFace = hashPairOrdered(face1, face2);
-            if (mapFaces.count(hashFace) == 0)
-            {
-              faceToFaces(nbNghbFaces(face1), face1) = face2;
-              faceToFaces(nbNghbFaces(face2), face2) = face1;
-              ++nbNghbFaces(face1);
-              ++nbNghbFaces(face2);
-              mapFaces[hashFace] = hashFace;
-            }
-          }
-        }
+    faceToFaces.resize(nvertex, nbFaces);
+    for (int nn = 0; nn < nbNodes; ++nn)
+    {
+      int nf1 = 0;
+      while (nodeToFaces(nf1+1, nn) != -1)
+      {
+        int nf2 = nf1+1;
+        int face1 = nodeToFaces(nf1, nn);
+        while (nodeToFaces(nf2, nn) != -1)
+        {
+          int face2 = nodeToFaces(nf2, nn);
+          int hash = hashPairOrdered(face1, face2);
+          if (mapFaces.count(hash) == 0)
+          {
+            faceToFaces(nbNghbFaces(face1), face1) = face2;
+            faceToFaces(nbNghbFaces(face2), face2) = face1;
+            ++nbNghbFaces(face1);
+            ++nbNghbFaces(face2);
+          }
+          else
+            mapFaces[hash] = hash;
+          ++nf2;
+        }
+        ++nf1;
       }
     }
+
   } //getLocNghbFacesEdgeType
 

XIOS/trunk/src/transformation/algo_types.hpp

@@ -25 +25 @@
 #include "domain_algorithm_compute_connectivity.hpp"
 #include "compute_connectivity_domain.hpp"
+#include "expand_domain.hpp"
+#include "domain_algorithm_expand.hpp"
 
 #endif // __XIOS_ALGORITHM_TRANSFORMATION_TYPES_HPP__

XIOS/trunk/src/transformation/domain_algorithm_expand.cpp

@@ -3 +3 @@
    \author Ha NGUYEN
    \since 08 Aug 2016
-   \date 08 Aug 2016
-
-   \brief Algorithm for expand on an domain.
+   \date 19 Sep 2016
+
+   \brief Algorithm for expanding an domain.
  */
 #include "domain_algorithm_expand.hpp"
@@ -13 +13 @@
 #include "grid.hpp"
 #include "grid_transformation_factory_impl.hpp"
+#include "context.hpp"
+#include "context_client.hpp"
 
 namespace xios {
@@ -45 +47 @@
 : CDomainAlgorithmTransformation(domainDestination, domainSource)
 {
+  if (domainDestination == domainSource)
+  {
+    ERROR("CDomainAlgorithmExpand::CDomainAlgorithmExpand(CDomain* domainDestination,CDomain* domainSource, CExpandDomain* expandDomain)",
+           << "Domain source and domain destination are the same. Please make sure domain destination refers to domain source" << std::endl
+           << "Domain source " <<domainSource->getId() << std::endl
+           << "Domain destination " <<domainDestination->getId() << std::endl);
+  }
+
+  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);
-}
+
+  switch (expandDomain->type)
+  {
+    case CExpandDomain::type_attr::node :
+      expandDomainNodeConnectivity(domainDestination,
+                                   domainSource);
+      break;
+    case CExpandDomain::type_attr::edge :
+      expandDomainEdgeConnectivity(domainDestination,
+                                   domainSource);
+      break;
+    default:
+      break;
+  }
+}
+
+/*!
+ *  Expand domain with edge-type neighbor
+ *  \param[in/out] domainDestination domain destination and will be modified
+ *  \param[in] domainSource domain source
+ */
+void CDomainAlgorithmExpand::expandDomainEdgeConnectivity(CDomain* domainDestination,
+                                                          CDomain* domainSource)
+{
+  CContext* context = CContext::getCurrent();
+  CContextClient* client=context->client;
+
+  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);
+}
+
+/*!
+ *  Expand domain with node-type neighbor
+ *  \param[in/out] domainDestination domain destination and will be modified
+ *  \param[in] domainSource domain source
+ */
+void CDomainAlgorithmExpand::expandDomainNodeConnectivity(CDomain* domainDestination,
+                                                          CDomain* domainSource)
+{
+  CContext* context = CContext::getCurrent();
+  CContextClient* client=context->client;
+
+  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);
+}
+
+/*!
+ *  Extend 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 domain extended part
+ */
+void CDomainAlgorithmExpand::updateDomainAttributes(CDomain* domainDestination,
+                                                    CDomain* domainSource,
+                                                    CArray<int,2>& neighborsDomainSrc)
+{
+  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()
+                                                                      : "";
+  domainDestination->domain_ref.setValue(domainSource->getId());
+  domainDestination->solveRefInheritance(true);
+  if (domainDstRef.empty()) domainDestination->domain_ref.reset();
+  else domainDestination->domain_ref.setValue(domainDstRef);
+
+  // Now extend domain destination
+  int niGlob = domainSource->ni_glo;
+  CArray<bool,1>& mask_1d_src = domainSource->mask_1d;
+  CArray<int,1>& i_index_src = domainSource->i_index;
+  CArray<double,1>& lon_src = domainSource->lonvalue_1d;
+  CArray<double,1>& lat_src = domainSource->latvalue_1d;
+  CArray<double,2>& bounds_lon_src = domainSource->bounds_lon_1d;
+  CArray<double,2>& bounds_lat_src = domainSource->bounds_lat_1d;
+  CArray<int,1>& data_i_index_src = domainSource->data_i_index;
+
+  int oldNbLocal = i_index_src.numElements(), index, globalIndex;
+  // Uncompress data_i_index
+  CArray<int,1> data_i_index_src_full(oldNbLocal);
+  int nbUnMaskedPointOnLocalDomain = 0;
+  data_i_index_src_full = -1; // Suppose all values are masked
+  for (int idx = 0; idx < data_i_index_src.numElements(); ++idx)
+  {
+    int dataIdx = data_i_index_src(idx);
+    if ((0 <= dataIdx) && (dataIdx < oldNbLocal))
+    {
+      data_i_index_src_full(nbUnMaskedPointOnLocalDomain) = dataIdx;
+      ++nbUnMaskedPointOnLocalDomain;
+    }
+  }
+
+  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<double,1>& lon_dst = domainDestination->lonvalue_1d;
+  CArray<double,1>& lat_dst = domainDestination->latvalue_1d;
+  CArray<double,2>& bounds_lon_dst = domainDestination->bounds_lon_1d;
+  CArray<double,2>& bounds_lat_dst = domainDestination->bounds_lat_1d;
+  CArray<int,1>& data_i_index_dst = domainDestination->data_i_index;
+  CArray<int,1>& data_j_index_dst = domainDestination->data_j_index;
+
+  // Resize all array-like attributes of domain destination
+  int nbNeighbor    = neighborsDomainSrc.shape()[1];
+  int newNbLocalDst = nbNeighbor + oldNbLocal;
+  int nVertex       = bounds_lon_dst.shape()[0];
+
+  mask_1d_dst.resizeAndPreserve(newNbLocalDst);
+  i_index_dst.resizeAndPreserve(newNbLocalDst);
+  j_index_dst.resizeAndPreserve(newNbLocalDst);
+  lon_dst.resizeAndPreserve(newNbLocalDst);
+  lat_dst.resizeAndPreserve(newNbLocalDst);
+  bounds_lon_dst.resizeAndPreserve(nVertex, newNbLocalDst);
+  bounds_lat_dst.resizeAndPreserve(nVertex, newNbLocalDst);
+  CArray<int,1> data_i_index_dst_full(newNbLocalDst);
+  data_i_index_dst_full(Range(0,oldNbLocal-1)) = data_i_index_src_full;
+  data_i_index_dst_full(Range(oldNbLocal,newNbLocalDst-1)) = -1;
+
+  // 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()));
+  // First, index mapping for local domain
+  for (int idx = 0; idx < oldNbLocal; ++idx)
+  {
+    index = i_index_dst(idx);
+    transMap[index].push_back(index);
+    transWeight[index].push_back(1.0);
+  }
+  // Then, index mapping for extended part
+  for (int idx = 0; idx < nbNeighbor; ++idx)
+  {
+    index = idx + oldNbLocal;
+    globalIndex = neighborsDomainSrc(0,idx);
+    i_index_dst(index) = globalIndex;
+    j_index_dst(index) = 0;
+    transMap[globalIndex].push_back(globalIndex);
+    transWeight[globalIndex].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; // data_i_index
+  // Initialize database
+  for (int idx = 0; idx < oldNbLocal; ++idx)
+  {
+    index = i_index_src(idx);
+    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; ++dataIdx;
+    data[dataIdx] = data_i_index_src_full(idx);
+  }
+
+  CClientClientDHTDouble dhtData(localData,client->intraComm);
+  CArray<size_t,1> neighborInd(nbNeighbor);
+  for (int idx = 0; idx < nbNeighbor; ++idx)
+    neighborInd(idx) = neighborsDomainSrc(0,idx);
+
+  // Compute local data on other domains
+  dhtData.computeIndexInfoMapping(neighborInd);
+  CClientClientDHTDouble::Index2VectorInfoTypeMap& neighborData = dhtData.getInfoIndexMap();
+  CClientClientDHTDouble::Index2VectorInfoTypeMap::iterator ite = neighborData.end(), it;
+  // Ok get neighbor data
+  size_t nIdx;
+  int nbUnMaskedPointOnExtendedPart = 0;
+  for (int idx = 0; idx < nbNeighbor; ++idx)
+  {
+    nIdx  = neighborInd(idx);
+    it = neighborData.find(nIdx);
+    if (ite != it)
+    {
+      index = idx + oldNbLocal;
+      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;
+      }
+      mask_1d_dst(index) = (1.0 == data[dataIdx]) ? true : false; ++dataIdx;
+      data_i_index_dst_full(index) = (int)(data[dataIdx]);
+      if (0 <= data_i_index_dst_full(index))
+      {
+        data_i_index_dst_full(index) = index;
+        ++nbUnMaskedPointOnExtendedPart;
+      }
+    }
+  }
+
+
+  // Finally, update data_i_index
+  int nbUnMaskedPointOnNewDstDomain = (nbUnMaskedPointOnExtendedPart + nbUnMaskedPointOnLocalDomain);
+  int count = 0, dataIdx;
+  for (int idx = 0; idx < newNbLocalDst; ++idx)
+  {
+    dataIdx = data_i_index_dst_full(idx);
+    if ((0 <= dataIdx) && (dataIdx < newNbLocalDst))
+    {
+      ++count;
+    }
+  }
+
+  data_i_index_dst.resize(count);
+  data_j_index_dst.resize(count);
+  data_j_index_dst = 0;
+
+  count = 0;
+  for (int idx = 0; idx < newNbLocalDst; ++idx)
+  {
+    dataIdx = data_i_index_dst_full(idx);
+    if ((0 <= dataIdx) && (dataIdx < newNbLocalDst))
+    {
+      data_i_index_dst(count) = dataIdx;
+      ++count;
+    }
+  }
+
+  // Update ni
+  domainDestination->ni.setValue(newNbLocalDst);
+
+}
+
 
 /*!

XIOS/trunk/src/transformation/domain_algorithm_expand.hpp

@@ -28 +28 @@
 
   static bool registerTrans();
+
+protected:
+  void expandDomainEdgeConnectivity(CDomain* domainDestination, CDomain* domainSource);
+  void expandDomainNodeConnectivity(CDomain* domainDestination, CDomain* domainSource);
+  void updateDomainAttributes(CDomain* domainDestination,
+                              CDomain* domainSource,
+                              CArray<int,2>& neighborsDomainSrc);
+
 protected:
   void computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs);

XIOS/trunk/src/transformation/grid_transformation.cpp

@@ -340 +340 @@
       algo = algoTransformation_[std::distance(itb, it)];
 
-    if ((0 != algo) && (CGenericAlgorithmTransformation::ELEMENT_NO_MODIFICATION_WITH_DATA == algo->type())) // Only registered transformation can be executed
+    if ((0 != algo) &&
+        ((CGenericAlgorithmTransformation::ELEMENT_NO_MODIFICATION_WITH_DATA == algo->type()) ||
+        (CGenericAlgorithmTransformation::ELEMENT_MODIFICATION_WITH_DATA == algo->type()))) // Only registered transformation can be executed
     {
       algo->computeIndexSourceMapping(dataAuxInputs);

XIOS/trunk/src/transformation/grid_transformation_selector.cpp

@@ -325 +325 @@
   CDomainAlgorithmInterpolate::registerTrans();
   CDomainAlgorithmZoom::registerTrans();
-
-}
-
-}
+  CDomainAlgorithmExpand::registerTrans();
+
+}
+
+}