Changeset 783

Timestamp:

11/12/15 16:33:08 (7 years ago)

Author:

mhnguyen

Message:

Adding attribute reading of unstructured grid

+) Fix a minor bug relating to unstructured grid detection
+) Add attribute reading for unstructured domain

Test
+) On Curie
+) test_remap passes

Location:

XIOS/trunk

Files:

• inputs/REMAP/iodef.xml (modified) (6 diffs)
• src/io/inetcdf4.cpp (modified) (1 diff)
• src/io/nc4_data_input.cpp (modified) (10 diffs)
• src/io/nc4_data_input.hpp (modified) (1 diff)
• src/node/domain.cpp (modified) (2 diffs)
• src/test/test_remap.f90 (modified) (3 diffs)

XIOS/trunk/inputs/REMAP/iodef.xml

@@ -8 +8 @@
    <field_definition level="1" >
      <field id="src_field" operation="instant" domain_ref="src_domain"/>
-<!--     <field id="dst_field" operation="instant" field_ref="src_field" domain_ref="dst_domain"/>-->
+     <field id="dst_field" operation="instant" field_ref="src_field" domain_ref="dst_domain"/>
      <field id="tmp_field" operation="instant" domain_ref="src_domain_regular_tmp"/>
      <field id="tmp_field_1" operation="instant" grid_ref="src_grid_curvilinear" />
+     <field id="tmp_field_2" operation="instant" domain_ref="domain_src_unstructured" />
      <field id="dst_field_regular" operation="instant" field_ref="tmp_field" domain_ref="dst_domain_regular" read_access="true"/>
-     <field id="dst_field_regular_pole" operation="instant" field_ref="src_field" domain_ref="dst_domain_regular_pole" read_access="true"/>
+     <field id="dst_field_regular_pole" operation="instant" field_ref="src_field" domain_ref="dst_domain_regular_pole"/>
    </field_definition>
 
@@ -32 +33 @@
 <!--        <field field_ref="tmp_field_1" operation="instant"/>-->
      </file>
+     <file id="output_dst_unstructured" name="output_dst_unstructured" type="one_file">
+        <field field_ref="tmp_field_2" operation="instant"/>
+     </file>
 
      <file id="output_src_regular" name="output_src_regular" mode="read" type="multiple_file">
@@ -41 +45 @@
 
      <file id="output_src_curvilinear" name="output_src_curvilinear" mode="read" type="one_file">
-<!--        <field id="src_field_curvilinear" name="field_A" domain_ref="src_domain_curvilinear" axis_ref="src_axis_curvilinear" operation="instant"/>-->
         <field id="src_field_curvilinear" name="field_A" grid_ref="src_grid_curvilinear" operation="instant"/>
      </file>
-     <file id="output_surface_curvi" name="output_surface_curvi" mode="read" type="one_file">
-        <field id="field_surface_curvi" name="field_A" domain_ref="domain_surface_curvi" axis_ref="axis_surface_curvi" operation="instant"/>
+     <file id="output_unstructured" name="output_unstructured" mode="read" type="one_file">
+        <field id="field_src_unstructred" name="field" domain_ref="domain_src_unstructured" operation="instant"/>
      </file>
    </file_definition>
@@ -59 +62 @@
    <domain_definition>
      <domain id="src_domain" />
-     <domain id="dst_domain" domain_src="src_domain_regular">
+     <domain id="dst_domain" domain_src="src_domain">
        <interpolate_domain/>
      </domain>
@@ -66 +69 @@
        <interpolate_domain/>
      </domain>
-<!--     <domain id="dst_domain_regular" domain_src="src_domain_regular" ni_glo="90" nj_glo="45" type="rectilinear">-->
      <domain id="dst_domain_regular" domain_src="src_domain_regular_tmp" ni_glo="90" nj_glo="45" type="rectilinear">
        <generate_rectilinear_domain />
@@ -80 +82 @@
        <generate_rectilinear_domain />
      </domain>
-     <domain id="domain_surface_curvi" type="curvilinear">
+     <domain id="domain_src_unstructured" type="unstructured">
        <generate_rectilinear_domain />
      </domain>

XIOS/trunk/src/io/inetcdf4.cpp

@@ -551 +551 @@
           return false;
       }
-      else return false;
     }
 

XIOS/trunk/src/io/nc4_data_input.cpp

@@ -143 +143 @@
     // Verify the compatibility of dimension of declared grid and real grid in file
     int realGridDim = 1;
+    bool isUnstructuredGrid = SuperClassWriter::isUnstructured(fieldId);
     std::map<StdString, StdSize> dimSizeMap = SuperClassWriter::getDimensions(&fieldId);
     realGridDim = SuperClassWriter::isTemporal(fieldId) ? dimSizeMap.size() - 1 : dimSizeMap.size();
+    if (isUnstructuredGrid) ++realGridDim;
 
     if (gridDim != realGridDim)
@@ -154 +156 @@
     if (SuperClassWriter::isTemporal(fieldId)) dimSizeMap.erase(SuperClassWriter::getUnlimitedDimensionName());
     int mapSize = dimSizeMap.size() - 1;
+    std::list<std::pair<StdString, StdSize> > listDimSize;
+    for (std::map<StdString, StdSize>::const_iterator itMap = dimSizeMap.begin(); itMap != dimSizeMap.end(); ++itMap)
+      listDimSize.push_front(*itMap);
 
     // Now process domain and axis
@@ -169 +174 @@
         {
            it = readValueDomains_.insert(domainP[idxDomain]->getId());
-           if (it.second) readDomainAttributeValueFromFile(domainP[idxDomain], dimSizeMap, mapSize - 1 - elementPosition, fieldId);
+           if (it.second) readDomainAttributeValueFromFile(domainP[idxDomain], listDimSize, elementPosition, fieldId);
         }
         else
         {
           it = readMetaDataDomains_.insert(domainP[idxDomain]->getId());
-          if (it.second) readDomainAttributesFromFile(domainP[idxDomain], dimSizeMap, mapSize - 1 - elementPosition, fieldId);
+          if (it.second) readDomainAttributesFromFile(domainP[idxDomain], listDimSize, elementPosition, fieldId);
         }
         ++idxDomain;
-        elementPosition += 2;
+        if (isUnstructuredGrid) ++elementPosition;
+        else elementPosition += 2;
       }
       else
@@ -184 +190 @@
         {
           it = readValueAxis_.insert(axisP[idxAxis]->getId());
-          if (it.second) readAxisAttributeValueFromFile(axisP[idxAxis], dimSizeMap, mapSize - elementPosition, fieldId);
+          if (it.second) readAxisAttributeValueFromFile(axisP[idxAxis], listDimSize, elementPosition, fieldId);
         }
         else
         {
           it = readMetaDataAxis_.insert(axisP[idxAxis]->getId());
-          if (it.second) readAxisAttributesFromFile(axisP[idxAxis], dimSizeMap, mapSize - elementPosition, fieldId);
+          if (it.second) readAxisAttributesFromFile(axisP[idxAxis], listDimSize, elementPosition, fieldId);
         }
         ++idxAxis;
@@ -204 +210 @@
     \param [in] fieldId id (or name) associated with the grid
   */
-  void CNc4DataInput::readDomainAttributeValueFromFile(CDomain* domain, std::map<StdString, StdSize>& dimSizeMap,
+  void CNc4DataInput::readDomainAttributeValueFromFile(CDomain* domain, std::list<std::pair<StdString, StdSize> >& dimSizeMap,
                                                        int elementPosition, const StdString& fieldId)
   {
     // There are some optional attributes of a domain to retrieve from file    // + lon lat?
-    std::map<StdString, StdSize>::const_iterator itMapNj = dimSizeMap.begin(), itMapNi,
-                                                 iteMap  = dimSizeMap.end();
-
-    for (int i = 0; i < elementPosition; ++i, ++itMapNj) {}
-    itMapNi = itMapNj; ++itMapNi;
+    std::list<std::pair<StdString, StdSize> >::const_iterator itMapNi = dimSizeMap.begin(), itMapNj,
+                                                              iteMap  = dimSizeMap.end();
+
+    for (int i = 0; i < elementPosition; ++i, ++itMapNi) {}
+    itMapNj = itMapNi; ++itMapNj;
 
     if (this->isRectilinear(fieldId))
@@ -225 +231 @@
       readFieldVariableValue(domain->lonvalue_rectilinear_read_from_file, itMapNi->first, nBeginLon, nSizeLon, true);
       domain->fillInRectilinearLonLat();
+    }
+    else if (this->isUnstructured(fieldId))
+    {
+      int ni = domain->ni;
+      std::vector<StdSize> nBeginLatLon(1,0), nSizeLatLon(1,0);
+      nBeginLatLon[0] = domain->ibegin.getValue();
+      nSizeLatLon[0]  = ni;
+
+      StdString latName = this->getLatCoordName(fieldId);
+      domain->latvalue_1d.resize(ni);
+      readFieldVariableValue(domain->latvalue_1d, latName, nBeginLatLon, nSizeLatLon);
+      StdString lonName = this->getLonCoordName(fieldId);
+      domain->lonvalue_1d.resize(ni);
+      readFieldVariableValue(domain->lonvalue_1d, lonName, nBeginLatLon, nSizeLatLon);
+
+      StdString boundsLatName = this->getBoundsId(latName);
+      if (0 == boundsLatName.compare(""))
+         ERROR("CNc4DataInput::readDomainAttributeValueFromFile(...)",
+              << "Field '" << fieldId << std::endl
+              << "Trying to read attributes from unstructured grid."
+              << "Latitude variable " << latName << " does not have bounds.");
+      StdString boundsLonName = this->getBoundsId(lonName);
+      if (0 == boundsLonName.compare(""))
+         ERROR("CNc4DataInput::readDomainAttributeValueFromFile(...)",
+              << "Field '" << fieldId << std::endl
+              << "Trying to read attributes from unstructured grid."
+              << "Longitude variable " << lonName << " does not have bounds.");
+
+      int nbVertex = this->getNbVertex(fieldId);
+      domain->nvertex.setValue(nbVertex);
+      std::vector<StdSize> nBeginBndsLatLon(2), nSizeBndsLatLon(2);
+      nBeginBndsLatLon[0] = domain->ibegin.getValue(); nSizeBndsLatLon[0] = ni;
+      nBeginBndsLatLon[1] = 0; nSizeBndsLatLon[1] = nbVertex;
+
+      domain->bounds_lat_1d.resize(nbVertex,ni);
+      readFieldVariableValue(domain->bounds_lat_1d, boundsLatName, nBeginBndsLatLon, nSizeBndsLatLon);
+      domain->bounds_lon_1d.resize(nbVertex,ni);
+      readFieldVariableValue(domain->bounds_lon_1d, boundsLonName, nBeginBndsLatLon, nSizeBndsLatLon);
     }
     else if (this->isCurvilinear(fieldId))
@@ -254 +298 @@
               << "Longitude variable " << lonName << " does not have bounds.");
 
-      domain->nvertex.setValue(4);
+      int nbVertex = this->getNbVertex(fieldId);
+      domain->nvertex.setValue(nbVertex);
       std::vector<StdSize> nBeginBndsLatLon(3), nSizeBndsLatLon(3);
       nBeginBndsLatLon[0] = domain->jbegin.getValue(); nSizeBndsLatLon[0] = nj;
       nBeginBndsLatLon[1] = domain->ibegin.getValue(); nSizeBndsLatLon[1] = ni;
-      nBeginBndsLatLon[2] = 0; nSizeBndsLatLon[2] = 4;
-
-      domain->bounds_lat_2d.resize(4,ni,nj);
+      nBeginBndsLatLon[2] = 0; nSizeBndsLatLon[2] = nbVertex;
+
+      domain->bounds_lat_2d.resize(nbVertex,ni,nj);
       readFieldVariableValue(domain->bounds_lat_2d, boundsLatName, nBeginBndsLatLon, nSizeBndsLatLon);
-      domain->bounds_lon_2d.resize(4,ni,nj);
+      domain->bounds_lon_2d.resize(nbVertex,ni,nj);
       readFieldVariableValue(domain->bounds_lon_2d, boundsLonName, nBeginBndsLatLon, nSizeBndsLatLon);
-    }
-    else if (this->isUnstructured(fieldId))
-    {
-
     }
   }
@@ -278 +319 @@
     \param [in] fieldId id (or name) associated with the grid
   */
-  void CNc4DataInput::readDomainAttributesFromFile(CDomain* domain, std::map<StdString, StdSize>& dimSizeMap,
+  void CNc4DataInput::readDomainAttributesFromFile(CDomain* domain, std::list<std::pair<StdString, StdSize> >& dimSizeMap,
                                                    int elementPosition, const StdString& fieldId)
   {
     // There are some mandatory attributes of a domain to retrieve from file
     // + ni_glo, nj_glo
-    std::map<StdString, StdSize>::const_iterator itMapNj = dimSizeMap.begin(), itMapNi,
-                                                 iteMap  = dimSizeMap.end();
-    for (int i = 0; i < elementPosition; ++i, ++itMapNj) {}
-    itMapNi = itMapNj; ++itMapNi;
+    std::list<std::pair<StdString, StdSize> >::const_iterator itMapNi = dimSizeMap.begin(), itMapNj,
+                                                              iteMap  = dimSizeMap.end();
+    for (int i = 0; i < elementPosition; ++i, ++itMapNi) {}
+    itMapNj = itMapNi; ++itMapNj;
 
     if (this->isRectilinear(fieldId) || this->isCurvilinear(fieldId))
     {
       domain->nj_glo.setValue(itMapNj->second);
-      domain->ni_glo.setValue((itMapNi)->second);
+      domain->ni_glo.setValue(itMapNi->second);
     }
     else if (this->isUnstructured(fieldId))
     {
-
+      domain->nj_glo.setValue(1);
+      domain->ni_glo.setValue(itMapNi->second);
     }
   }
@@ -306 +348 @@
     \param [in] fieldId id (or name) associated with the grid
   */
-  void CNc4DataInput::readAxisAttributesFromFile(CAxis* axis, std::map<StdString, StdSize>& dimSizeMap,
+  void CNc4DataInput::readAxisAttributesFromFile(CAxis* axis, std::list<std::pair<StdString, StdSize> >& dimSizeMap,
                                                  int elementPosition, const StdString& fieldId)
   {
-    std::map<StdString, StdSize>::const_iterator itMapN = dimSizeMap.begin(),
-                                                 iteMap  = dimSizeMap.end();
+    std::list<std::pair<StdString, StdSize> >::const_iterator itMapN = dimSizeMap.begin(),
+                                                              iteMap = dimSizeMap.end();
     for (int i = 0; i < elementPosition; ++i, ++itMapN) {}
     axis->n_glo.setValue(itMapN->second);
@@ -322 +364 @@
     \param [in] fieldId id (or name) associated with the grid
   */
-  void CNc4DataInput::readAxisAttributeValueFromFile(CAxis* axis, std::map<StdString, StdSize>& dimSizeMap,
+  void CNc4DataInput::readAxisAttributeValueFromFile(CAxis* axis, std::list<std::pair<StdString, StdSize> >& dimSizeMap,
                                                     int elementPosition, const StdString& fieldId)
   {
-    std::map<StdString, StdSize>::const_iterator itMapN = dimSizeMap.begin(),
-                                                 iteMap  = dimSizeMap.end();
+    std::list<std::pair<StdString, StdSize> >::const_iterator itMapN = dimSizeMap.begin(),
+                                                              iteMap = dimSizeMap.end();
     for (int i = 0; i < elementPosition; ++i, ++itMapN) {}
 

XIOS/trunk/src/io/nc4_data_input.hpp

@@ -40 +40 @@
 
   private:
-    void readDomainAttributesFromFile(CDomain* domain, std::map<StdString, StdSize>& dimSizeMap,
+    void readDomainAttributesFromFile(CDomain* domain, std::list<std::pair<StdString, StdSize> >& dimSizeMap,
                                       int elementPosition, const StdString& fieldId);
-    void readDomainAttributeValueFromFile(CDomain* domain, std::map<StdString, StdSize>& dimSizeMap,
+    void readDomainAttributeValueFromFile(CDomain* domain, std::list<std::pair<StdString, StdSize> >& dimSizeMap,
                                           int elementPosition, const StdString& fieldId);
 
-    void readAxisAttributesFromFile(CAxis* axis, std::map<StdString, StdSize>& dimSizeMap,
+    void readAxisAttributesFromFile(CAxis* axis, std::list<std::pair<StdString, StdSize> >& dimSizeMap,
                                     int elementPosition, const StdString& fieldId);
-    void readAxisAttributeValueFromFile(CAxis* axis, std::map<StdString, StdSize>& dimSizeMap,
+    void readAxisAttributeValueFromFile(CAxis* axis, std::list<std::pair<StdString, StdSize> >& dimSizeMap,
                                         int elementPosition, const StdString& fieldId);
 

XIOS/trunk/src/node/domain.cpp

@@ -228 +228 @@
    {
      if (this->isRedistributed_) return;
+
+     this->isRedistributed_ = true;
+     CContext* context = CContext::getCurrent();
+     CContextClient* client = context->client;
+     int rankClient = client->clientRank;
+     int rankOnDomain = rankClient%nbLocalDomain;
+
+     if (ni_glo.isEmpty() || ni_glo <= 0 )
+     {
+        ERROR("CDomain::redistribute(int nbLocalDomain)",
+           << "[ Id = " << this->getId() << " ] "
+           << "The global domain is badly defined,"
+           << " check the \'ni_glo\'  value !")
+     }
+
+     if (nj_glo.isEmpty() || nj_glo <= 0 )
+     {
+        ERROR("CDomain::redistribute(int nbLocalDomain)",
+           << "[ Id = " << this->getId() << " ] "
+           << "The global domain is badly defined,"
+           << " check the \'nj_glo\'  value !")
+     }
+
      if ((type_attr::rectilinear == type)  || (type_attr::curvilinear == type))
      {
-        this->isRedistributed_ = true;
-
-        CContext* context = CContext::getCurrent();
-        CContextClient* client = context->client;
-        int rankClient = client->clientRank;
-        int rankOnDomain = rankClient%nbLocalDomain;
-
-        if (ni_glo.isEmpty() || ni_glo <= 0 )
-        {
-           ERROR("CDomain::redistribute(int nbLocalDomain)",
-              << "[ Id = " << this->getId() << " ] "
-              << "The global domain is badly defined,"
-              << " check the \'ni_glo\'  value !")
-        }
-
-        if (nj_glo.isEmpty() || nj_glo <= 0 )
-        {
-           ERROR("CDomain::redistribute(int nbLocalDomain)",
-              << "[ Id = " << this->getId() << " ] "
-              << "The global domain is badly defined,"
-              << " check the \'nj_glo\'  value !")
-        }
-
         int globalDomainSize = ni_glo * nj_glo;
         if (globalDomainSize <= nbLocalDomain)
@@ -329 +329 @@
         // Now fill other attributes
         if (type_attr::rectilinear == type) fillInRectilinearLonLat();
+     }
+     else  // unstructured domain
+     {
+        int globalDomainSize = ni_glo * nj_glo;
+        if (globalDomainSize <= nbLocalDomain)
+        {
+          for (int idx = 0; idx < nbLocalDomain; ++idx)
+          {
+            if (rankOnDomain < globalDomainSize)
+            {
+              int iIdx = rankOnDomain % ni_glo;
+              int jIdx = rankOnDomain / ni_glo;
+              ibegin.setValue(iIdx); jbegin.setValue(jIdx);
+              ni.setValue(1); nj.setValue(1);
+            }
+            else
+            {
+              ibegin.setValue(0); jbegin.setValue(0);
+              ni.setValue(0); nj.setValue(0);
+            }
+          }
+        }
+        else
+        {
+          float njGlo = nj_glo.getValue();
+          float niGlo = ni_glo.getValue();
+          std::vector<int> ibeginVec(nbLocalDomain,0);
+          std::vector<int> niVec(nbLocalDomain);
+          for (int i = 1; i < nbLocalDomain; ++i)
+          {
+            int range = ni_glo / nbLocalDomain;
+            if (i < (ni_glo%nbLocalDomain)) ++range;
+            niVec[i-1] = range;
+            ibeginVec[i] = ibeginVec[i-1] + niVec[i-1];
+          }
+          niVec[nbLocalDomain-1] = ni_glo - ibeginVec[nbLocalDomain-1];
+
+          int iIdx = rankOnDomain % nbLocalDomain;
+          ibegin.setValue(ibeginVec[iIdx]);
+          jbegin.setValue(0);
+          ni.setValue(niVec[iIdx]);
+          nj.setValue(1);
+        }
      }
    }

XIOS/trunk/src/test/test_remap.f90

@@ -20 +20 @@
   DOUBLE PRECISION,ALLOCATABLE :: src_boundslon(:,:), dst_boundslon(:,:)
   DOUBLE PRECISION,ALLOCATABLE :: src_boundslat(:,:), dst_boundslat(:,:)
-  DOUBLE PRECISION,ALLOCATABLE :: src_field(:), tmp_field(:), tmp_field_1(:)
+  DOUBLE PRECISION,ALLOCATABLE :: src_field(:), tmp_field(:), tmp_field_1(:), tmp_field_2(:)
   INTEGER :: src_ni_glo, dst_ni_glo;
   INTEGER :: src_nvertex, dst_nvertex;
@@ -132 +132 @@
   ALLOCATE(tmp_field_1(src_tmp_ni*src_tmp_nj*src_tmp_n))
 
+  CALL xios_get_domain_attr("domain_src_unstructured", ni=src_tmp_ni, nj=src_tmp_nj)
+  ALLOCATE(tmp_field_2(src_tmp_ni*src_tmp_nj))
+
   DO ts=1,1
     CALL xios_recv_field("src_field_regular_tmp", tmp_field)
     CALL xios_recv_field("src_field_curvilinear", tmp_field_1)
+    CALL xios_recv_field("field_src_unstructred", tmp_field_2)
     CALL xios_update_calendar(ts)
     CALL xios_send_field("src_field",src_field)
     CALL xios_send_field("tmp_field",tmp_field)
     CALL xios_send_field("tmp_field_1",tmp_field_1)
+    CALL xios_send_field("tmp_field_2",tmp_field_2)
     CALL wait_us(5000) ;
   ENDDO
@@ -146 +151 @@
   DEALLOCATE(src_lon, src_lat, src_boundslon,src_boundslat, src_field)
   DEALLOCATE(dst_lon, dst_lat, dst_boundslon,dst_boundslat)
+  DEALLOCATE(tmp_field, tmp_field_1, tmp_field_2)
 
   CALL MPI_COMM_FREE(comm, ierr)