Changeset 2200

Timestamp:

08/19/21 13:28:53 (3 years ago)

Author:

ymipsl

Message:

Bugs fix for UGRID convention output

• global indexation was not taking into account
• coherence problem in connectivity for node, edge and face mesh
• add UGRID testcase based on sphere cube

YM

Location:

XIOS/trunk/src

Files:

• io/nc4_data_output.cpp (modified) (5 diffs)
• node/file.cpp (modified) (3 diffs)
• node/mesh.cpp (modified) (5 diffs)
• test/test_ugrid.f90 (added)

XIOS/trunk/src/io/nc4_data_output.cpp

@@ -556 +556 @@
       // The first domain for the same mesh that will be written is that with the highest value of nvertex.
       // Thus the entire mesh connectivity will be generated at once.
+
       if (isWrittenDomain(domid)) return ;
       else setWrittenDomain(domid);
@@ -570 +571 @@
       std::vector<StdString> dim0;
       StdString domainName = domain->name;
-      domain->assignMesh(domainName, domain->nvertex);
-      domain->mesh->createMeshEpsilon(server->intraComm, domain->lonvalue, domain->latvalue, domain->bounds_lonvalue, domain->bounds_latvalue);
 
       StdString node_x = domainName + "_node_x";
@@ -600 +599 @@
       SuperClassWriter::addAttribute("topology_dimension", 2, &domainName);
       SuperClassWriter::addAttribute("node_coordinates", node_x + " " + node_y, &domainName);
+
 
       try
@@ -810 +810 @@
                  startNodes[0] = domain->mesh->node_start;
                  countNodes[0] = domain->mesh->node_count;
+                 if (countNodes[0]==0) startNodes[0]=0 ; // for netcdf error
                  startEdgeNodes[0] = domain->ibegin;
                  startEdgeNodes[1] = 0;
@@ -846 +847 @@
                  startNodes[0] = domain->mesh->node_start;
                  countNodes[0] = domain->mesh->node_count;
+                 if (countNodes[0]==0) startNodes[0]=0;
+                 
                  startEdges[0] = domain->mesh->edge_start;
                  countEdges[0] = domain->mesh->edge_count;
+                 if (countEdges[0]==0) startEdges[0]=0 ; // for netcdf error
+                 
                  startEdgeNodes[0] = domain->mesh->edge_start;
+                 countEdgeNodes[0] = domain->mesh->edge_count;
+                 if (countEdgeNodes[0]==0) startEdgeNodes[0]=0; // for netcdf error
                  startEdgeNodes[1] = 0;
-                 countEdgeNodes[0] = domain->mesh->edge_count;
                  countEdgeNodes[1]= 2;
+
                  startEdgeFaces[0] = domain->mesh->edge_start;
+                 countEdgeFaces[0] = domain->mesh->edge_count;
+                 if (countEdgeFaces[0]==0) startEdgeFaces[0]=0 ; // for netcdf error
                  startEdgeFaces[1]= 0;
-                 countEdgeFaces[0] = domain->mesh->edge_count;
                  countEdgeFaces[1]= 2;
+                
                  startFaceConctv[0] = domain->ibegin;
+                 countFaceConctv[0] = domain->ni;
                  startFaceConctv[1] = 0;
-                 countFaceConctv[0] = domain->ni;
                  countFaceConctv[1] = domain->nvertex;
                }

XIOS/trunk/src/node/file.cpp

@@ -867 +867 @@
    \brief Sorting domains with the same name (= describing the same mesh) in the decreasing order of nvertex for UGRID files.
    This insures that the domain with the highest nvertex is written first and thus all known mesh connectivity is generated at once by this domain.
+   Secondly, we associate the UGRID mesh to the domain in incresing order, i.e. nvertex=1 first, nvertex=2 and then nvertex>2.
+   In this case the connectivity of each component of the mesh (node, edge and face) are fully coherent.
    */
    void CFile::sortEnabledFieldsForUgrid()
@@ -874 +876 @@
      std::vector<int> domainNvertices;
      std::vector<StdString> domainNames;
+     std::map<string, tuple<CDomain*,CDomain*,CDomain*>> registeredDomains ;
 
      for (int i = 0; i < size; ++i)
@@ -916 +919 @@
          domainNvertices.push_back(nvertex);
        }
+
+       if (nvertex==1)  std::get<0>(registeredDomains[domainName])=domain[0] ;
+       else if (nvertex==2) std::get<1>(registeredDomains[domainName])=domain[0] ;
+       else  std::get<2>(registeredDomains[domainName])=domain[0] ;
+     }
+
+     for(auto& it:registeredDomains)
+     {
+       list<CDomain*> domains ;
+       string domainName=it.first ;
+
+       if (std::get<0>(it.second)!=nullptr) domains.push_back(std::get<0>(it.second)) ;
+       if (std::get<1>(it.second)!=nullptr) domains.push_back(std::get<1>(it.second)) ; 
+       if (std::get<2>(it.second)!=nullptr) domains.push_back(std::get<2>(it.second)) ; 
+       
+       // for each component of a given mesh (i.e. domains with same name but different number of vertices)
+       // associate the UGRID mesh in increasing order
+       for(auto& domain : domains )
+       {
+         domain-> computeWrittenIndex();
+         CArray<int, 1>& indexToWrite = domain->localIndexToWriteOnServer;
+         int nbWritten = indexToWrite.numElements();
+         CArray<double,1> writtenLat, writtenLon;
+         CArray<double,2> writtenBndsLat, writtenBndsLon;
+
+         writtenLat.resize(nbWritten);
+         writtenLon.resize(nbWritten);
+         for (int idx = 0; idx < nbWritten; ++idx)
+         {
+           writtenLat(idx) = domain->latvalue(indexToWrite(idx));
+           writtenLon(idx) = domain->lonvalue(indexToWrite(idx));
+         }
+    
+         int nvertex = domain->nvertex, idx;
+        if (nvertex>1)
+         {
+           writtenBndsLat.resize(nvertex, nbWritten);
+           writtenBndsLon.resize(nvertex, nbWritten);
+           CArray<double,2>& boundslat = domain->bounds_latvalue;
+           CArray<double,2>& boundslon = domain->bounds_lonvalue;
+           for (idx = 0; idx < nbWritten; ++idx)
+             for (int nv = 0; nv < nvertex; ++nv)
+             {
+                writtenBndsLat(nv, idx) = boundslat(nv, int(indexToWrite(idx)));
+                writtenBndsLon(nv, idx) = boundslon(nv, int(indexToWrite(idx)));
+             }
+         }
+         domain->assignMesh(domainName, domain->nvertex);
+         CContextServer* server=CContext::getCurrent()->server ;
+         domain->mesh->createMeshEpsilon(server->intraComm, writtenLon, writtenLat, writtenBndsLon, writtenBndsLat);
+       }
+
      }
    }

XIOS/trunk/src/node/mesh.cpp

@@ -433 +433 @@
       node_lat = latvalue;
 
+      unsigned long nodeCount = nbNodes_;
+      unsigned long nodeStart, nbNodes;
+      MPI_Scan(&nodeCount, &nodeStart, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
+      int nNodes = nodeStart;
+      MPI_Bcast(&nNodes, 1, MPI_UNSIGNED_LONG, mpiSize-1, comm);
+      nbNodesGlo = nNodes;
+
+      nodeStart -= nodeCount;
+      node_start = nodeStart;
+      node_count = nodeCount;
+
       // Global node indexes
       vector<size_t> hashValues(4);
@@ -441 +452 @@
         for (size_t nh = 0; nh < 4; ++nh)
         {
-          nodeHash2IdxGlo[hashValues[nh]].push_back(mpiRank*nbNodes_ + nn);
+          nodeHash2IdxGlo[hashValues[nh]].push_back(nodeStart + nn); 
         }
       }
       pNodeGlobalIndex = new CClientClientDHTSizet (nodeHash2IdxGlo, comm);
+
       nodesAreWritten = true;
     }
@@ -460 +472 @@
       unsigned long nbEdgesOnProc = nbEdges_;
       unsigned long nbEdgesAccum;
+      unsigned long nbEdgesGlo;
       MPI_Scan(&nbEdgesOnProc, &nbEdgesAccum, 1, MPI_UNSIGNED_LONG, MPI_SUM, comm);
+      nbEdgesGlo = nbEdgesAccum ;
+      MPI_Bcast(&nbEdgesGlo, 1, MPI_UNSIGNED_LONG, mpiSize-1, comm);
       nbEdgesAccum -= nbEdges_;
+      edge_start = nbEdgesAccum ;
+      edge_count = nbEdgesOnProc ;
 
       CClientClientDHTSizet::Index2VectorInfoTypeMap edgeHash2IdxGlo;
@@ -670 +687 @@
     } //nvertex = 2
 
-    else
+    else // nvertex > 2
     {
       nbFaces_ = bounds_lon.shape()[1];
@@ -1601 +1618 @@
   } // createMeshEpsilon
 
+
   ///----------------------------------------------------------------
   /*!