Changeset 944

Timestamp:

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

Author:

mhnguyen

Message:

Update compute_connectivity_domain with new functions of class CMesh

+) Make use of the function computing local neighbor of class CMesh

Test
+) On Curie
+) Work

Location:

XIOS/trunk

Files:

• inputs/Unstruct/iodef.xml (modified) (2 diffs)
• src/node/compute_connectivity_domain.cpp (modified) (1 diff)
• src/test/test_unstruct_complete.f90 (modified) (3 diffs)
• src/transformation/domain_algorithm_compute_connectivity.cpp (modified) (2 diffs)
• src/transformation/domain_algorithm_compute_connectivity.hpp (modified) (1 diff)

XIOS/trunk/inputs/Unstruct/iodef.xml

@@ -7 +7 @@
    <field_definition level="1" enabled=".TRUE." default_value="1000">
      <field id="field_A_srf"  operation="average" freq_op="3600s" grid_ref="grid_A"/>
+     <field id="field_A_expand"  operation="average" grid_ref="grid_dst" field_ref="field_A_srf"/>
    </field_definition>
 
-   <file_definition type="one_file" par_access="collective" output_freq="6h" output_level="10" enabled=".TRUE." >
+   <file_definition type="one_file" par_access="collective" output_freq="1h" output_level="10" enabled=".TRUE." >
      <file id="output" name="output">
-        <field field_ref="field_A_srf" name="field_A_srf"/>
+        <field field_ref="field_A_expand" name="field"/>
      </file>
    </file_definition>
 
    <axis_definition>
-     <axis id="axis_srf"/>
+     <axis id="axis_srf" positive="up"/>
    </axis_definition>
 
    <domain_definition>
-     <domain id="domain_srf" >
-       <compute_connectivity_domain id="compute" type="edge" />
+     <domain id="domain_srf" />     
+     <domain id="domain_dst" domain_ref="domain_srf" >
+       <expand_domain/>
+       <compute_connectivity_domain id="compute" type="node"/>
      </domain>
    </domain_definition>
@@ -28 +31 @@
      <grid id="grid_A">
        <domain domain_ref="domain_srf" />
+       <axis axis_ref="axis_srf" />
+     </grid>
+     <grid id="grid_dst">
+       <domain domain_ref="domain_dst" />
        <axis axis_ref="axis_srf" />
      </grid>

XIOS/trunk/src/node/compute_connectivity_domain.cpp

@@ -38 +38 @@
   void CComputeConnectivityDomain::checkValid(CDomain* domainDst)
   {
-    if (type.isEmpty())
+    if (CDomain::type_attr::unstructured != domainDst->type)
+    {
       ERROR("CComputeConnectivityDomain::checkValid(CDomain* domainDst)",
-       << "Connectivity type is not defined. Chose 'node' or 'edge' for the type."
-       << "Domain " <<domainDst->getId() << std::endl
-       << "Connectivity object " << this->getId());
+            << "Domain connectivity computation is only supported for unstructured" << std::endl
+            << "Check type of domain destination, id = " << domainDst->getId());
+    }
 
+    if (type.isEmpty()) type.setValue(CComputeConnectivityDomain::type_attr::edge);
     if (n_neighbor_max.isEmpty()) n_neighbor_max.setValue(0);
     if (n_neighbor.isEmpty()) n_neighbor.resize(domainDst->i_index.numElements());

XIOS/trunk/src/test/test_unstruct_complete.f90

@@ -17 +17 @@
   INTEGER,PARAMETER :: ni_glo=100
   INTEGER,PARAMETER :: nj_glo=100
-  INTEGER,PARAMETER :: llm=5
+  INTEGER,PARAMETER :: llm=1
   DOUBLE PRECISION  :: lval(llm)=1
   TYPE(xios_field) :: field_hdl
   TYPE(xios_fieldgroup) :: fieldgroup_hdl
   TYPE(xios_file) :: file_hdl
+  TYPE(xios_domain) :: domain_hdl
   LOGICAL :: ok
 
@@ -179 +180 @@
       IF (MOD(ind,8)>=0 .AND. MOD(ind,8)<2) THEN
         mask(ncell)=.FALSE.
+!        data_n_index=data_n_index+1
       ELSE
         mask(ncell)=.TRUE.
@@ -215 +217 @@
   CALL xios_close_context_definition()
 
+  CALL xios_field_get_domain("field_A_expand",domain_hdl)
+  CALL xios_get_attr(domain_hdl,ni=ni)
+  PRINT *, "Domain destination ni = ",ni, "Domain source n= ",ncell
   CALL xios_get_compute_connectivity_domain_attr("compute", n_neighbor_max=nbMax)
   ALLOCATE(local_neighbor(nbMax,ncell))
   CALL xios_get_compute_connectivity_domain_attr("compute", n_neighbor=n_local, local_neighbor=local_neighbor)
-
-   DO ts=1,24*10
+  PRINT *, "Connectivity max = ",nbMax
+
+   DO ts=1,1
      CALL xios_update_calendar(ts)
      CALL xios_send_field("field_A_srf",field_A_compressed)

XIOS/trunk/src/transformation/domain_algorithm_compute_connectivity.cpp

@@ -49 +49 @@
   CArray<int,1>& nbNeighbor = compute_connectivityDomain->n_neighbor;
   CArray<int,2>& localNeighbors = compute_connectivityDomain->local_neighbor;
+  int type = 1; // Edge type
   switch (compute_connectivityDomain->type)
   {
     case CComputeConnectivityDomain::type_attr::node :
-      computeNodeConnectivity(domainDestination,
-                              nbNeighborMax,
-                              nbNeighbor,
-                              localNeighbors);
+      type = 0;
       break;
     case CComputeConnectivityDomain::type_attr::edge :
-      computeEdgeConnectivity(domainDestination,
-                              nbNeighborMax,
-                              nbNeighbor,
-                              localNeighbors);
+      type = 1;
       break;
     default:
       break;
   }
+
+  computeLocalConnectivity(type, domainDestination, nbNeighborMax, nbNeighbor, localNeighbors);
 }
 
-void CDomainAlgorithmComputeConnectivity::computeNodeConnectivity(CDomain* domain,
+/*!
+ *  Compute local connectivity of a domain
+ *  \param[in] type type of connectivity (node or edge)
+ *  \param[in] domain domain on which we calculate local connectivity
+ *  \param[in/out] nbConnectivityMax maximum number of neighbor a cell of domain has
+ *  \param[in/out] nbConnectivity number of neighbor a cell has
+ *  \param[in/out] localConnectivity localConnectivity local index of neighbor of a cell
+ */
+void CDomainAlgorithmComputeConnectivity::computeLocalConnectivity(int type,
+                                                                  CDomain* domain,
                                                                   int& nbConnectivityMax,
                                                                   CArray<int,1>& nbConnectivity,
                                                                   CArray<int,2>& localConnectivity)
 {
+
   CMesh mesh;
-  CArray<double,1>& lon = domain->lonvalue_1d;
-  CArray<double,1>& lat = domain->latvalue_1d;
+
   CArray<double,2>& bounds_lon = domain->bounds_lon_1d;
   CArray<double,2>& bounds_lat = domain->bounds_lat_1d;
-  mesh.createMesh(lon, lat, bounds_lon, bounds_lat);
-  CArray<int, 2>& face_nodes = mesh.face_nodes;
-  int nvertex = mesh.nvertex;
-  int nFaces  = mesh.nbFaces;
-  int nNodes  = mesh.nbNodes;
-  std::vector<std::vector<size_t> > nodeFaceConnectivity(nNodes, std::vector<size_t>());
-  for (int nf = 0; nf < nFaces; ++nf)
-  {
-    for (int nv = 0; nv < nvertex; ++nv)
-    {
-      nodeFaceConnectivity[face_nodes(nv,nf)].push_back(nf);
-    }
-  }
+  int ncell = bounds_lon.shape()[1];
+  CArray<int,1> localIndex(ncell);
+  for (int idx = 0; idx <ncell; ++idx) localIndex(idx) = idx;
 
-  if (nFaces != nbConnectivity.numElements()) nbConnectivity.resize(nFaces);
-  nbConnectivityMax = 1;
-  for (int nf = 0; nf < nFaces; ++nf)
-  {
-    std::set<size_t> neighFaces;
-    for (int nv = 0; nv < nvertex; ++nv)
-    {
-      std::vector<size_t>& tmpFaces = nodeFaceConnectivity[face_nodes(nv,nf)];
-      for (int nFn = 0; nFn < tmpFaces.size(); ++nFn)
-      {
-        neighFaces.insert(tmpFaces[nFn]);
-      }
-    }
-    if (nbConnectivityMax < (neighFaces.size()-1)) nbConnectivityMax = neighFaces.size()-1;
-  }
-  if ((nbConnectivityMax != localConnectivity.shape()[0]) || (nFaces != localConnectivity.shape()[1]))
-    localConnectivity.resize(nbConnectivityMax, nFaces);
-
-  for (int nf = 0; nf < nFaces; ++nf)
-  {
-    std::set<size_t> neighFaces;
-    for (int nv = 0; nv < nvertex; ++nv)
-    {
-      std::vector<size_t>& tmpFaces = nodeFaceConnectivity[face_nodes(nv,nf)];
-      for (int nFn = 0; nFn < tmpFaces.size(); ++nFn)
-      {
-        neighFaces.insert(tmpFaces[nFn]);
-      }
-    }
-
-    neighFaces.erase(nf);
-    nbConnectivity(nf) = neighFaces.size();
-    std::set<size_t>::iterator it = neighFaces.begin(), ite = neighFaces.end();
-    for (int idx = 0; it != ite; ++it, ++idx)
-    {
-      localConnectivity(idx, nf) = *it;
-    }
-  }
+  mesh.getLocalNghbFaces(type, localIndex, bounds_lon, bounds_lat, localConnectivity, nbConnectivity);
+  nbConnectivityMax = 0;
+  for (int idx =0; idx < nbConnectivity.numElements(); ++idx)
+    if (nbConnectivityMax < nbConnectivity(idx)) nbConnectivityMax = nbConnectivity(idx);
 }
 
-void CDomainAlgorithmComputeConnectivity::computeEdgeConnectivity(CDomain* domain,
-                                                                  int& nbConnectivityMax,
-                                                                  CArray<int,1>& nbConnectivity,
-                                                                  CArray<int,2>& localConnectivity)
-{
-  CMesh mesh;
-  CArray<double,1>& lon = domain->lonvalue_1d;
-  CArray<double,1>& lat = domain->latvalue_1d;
-  CArray<double,2>& bounds_lon = domain->bounds_lon_1d;
-  CArray<double,2>& bounds_lat = domain->bounds_lat_1d;
-  mesh.createMesh(lon, lat, bounds_lon, bounds_lat);
 
-  CArray<int, 2>& face_edges = mesh.face_edges;
-  int nvertex = mesh.nvertex;
-  int nFaces  = mesh.nbFaces;
-  int nEdges  = mesh.nbEdges;
-  std::vector<std::vector<size_t> > edgeFaceConnectivity(nEdges, std::vector<size_t>());
-  for (int nf = 0; nf < nFaces; ++nf)
-  {
-    for (int nv = 0; nv < nvertex; ++nv)
-    {
-      edgeFaceConnectivity[face_edges(nv,nf)].push_back(nf);
-    }
-  }
-
-  if (nFaces != nbConnectivity.numElements()) nbConnectivity.resize(nFaces);
-  nbConnectivityMax = 1;
-  for (int nf = 0; nf < nFaces; ++nf)
-  {
-    std::set<size_t> neighFaces;
-    for (int nv = 0; nv < nvertex; ++nv)
-    {
-      std::vector<size_t>& tmpFaces = edgeFaceConnectivity[face_edges(nv,nf)];
-      for (int nFn = 0; nFn < tmpFaces.size(); ++nFn)
-      {
-        neighFaces.insert(tmpFaces[nFn]);
-      }
-    }
-    if (nbConnectivityMax < (neighFaces.size()-1)) nbConnectivityMax = neighFaces.size() - 1;
-  }
-  if ((nbConnectivityMax != localConnectivity.shape()[0]) || (nFaces != localConnectivity.shape()[1]))
-    localConnectivity.resize(nbConnectivityMax, nFaces);
-
-  for (int nf = 0; nf < nFaces; ++nf)
-  {
-    std::set<size_t> neighFaces;
-    for (int nv = 0; nv < nvertex; ++nv)
-    {
-      std::vector<size_t>& tmpFaces = edgeFaceConnectivity[face_edges(nv,nf)];
-      for (int nFn = 0; nFn < tmpFaces.size(); ++nFn)
-      {
-        neighFaces.insert(tmpFaces[nFn]);
-      }
-    }
-
-    neighFaces.erase(nf);
-    nbConnectivity(nf) = neighFaces.size();
-    std::set<size_t>::iterator it = neighFaces.begin(), ite = neighFaces.end();
-    for (int idx = 0; it != ite; ++it, ++idx)
-    {
-      localConnectivity(idx, nf) = *it;
-    }
-  }
-
-}
 
 /*!
@@ -203 +101 @@
 void CDomainAlgorithmComputeConnectivity::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
 {
-
 }
 

XIOS/trunk/src/transformation/domain_algorithm_compute_connectivity.hpp

@@ -32 +32 @@
 
 protected:
-  void computeNodeConnectivity(CDomain* domain, int& nbConnectivityMax,
-                               CArray<int,1>& nbConnectivity, CArray<int,2>& localConnectivity);
-
-  void computeEdgeConnectivity(CDomain* domain, int& nbConnectivityMax,
-                               CArray<int,1>& nbConnectivity, CArray<int,2>& localConnectivity);
+  void computeLocalConnectivity(int type,
+                                CDomain* domain,
+                                int& nbConnectivityMax,
+                                CArray<int,1>& nbConnectivity,
+                                CArray<int,2>& localConnectivity);
 
 private: