Changeset 944


Ignore:
Timestamp:
09/22/16 10:59:32 (4 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:
5 edited

Legend:

Unmodified
Added
Removed
  • XIOS/trunk/inputs/Unstruct/iodef.xml

    r934 r944  
    77   <field_definition level="1" enabled=".TRUE." default_value="1000"> 
    88     <field id="field_A_srf"  operation="average" freq_op="3600s" grid_ref="grid_A"/> 
     9     <field id="field_A_expand"  operation="average" grid_ref="grid_dst" field_ref="field_A_srf"/> 
    910   </field_definition> 
    1011 
    11    <file_definition type="one_file" par_access="collective" output_freq="6h" output_level="10" enabled=".TRUE." > 
     12   <file_definition type="one_file" par_access="collective" output_freq="1h" output_level="10" enabled=".TRUE." > 
    1213     <file id="output" name="output"> 
    13         <field field_ref="field_A_srf" name="field_A_srf"/> 
     14        <field field_ref="field_A_expand" name="field"/> 
    1415     </file> 
    1516   </file_definition> 
    1617 
    1718   <axis_definition> 
    18      <axis id="axis_srf"/> 
     19     <axis id="axis_srf" positive="up"/> 
    1920   </axis_definition> 
    2021 
    2122   <domain_definition> 
    22      <domain id="domain_srf" > 
    23        <compute_connectivity_domain id="compute" type="edge" /> 
     23     <domain id="domain_srf" />      
     24     <domain id="domain_dst" domain_ref="domain_srf" > 
     25       <expand_domain/> 
     26       <compute_connectivity_domain id="compute" type="node"/> 
    2427     </domain> 
    2528   </domain_definition> 
     
    2831     <grid id="grid_A"> 
    2932       <domain domain_ref="domain_srf" /> 
     33       <axis axis_ref="axis_srf" /> 
     34     </grid> 
     35     <grid id="grid_dst"> 
     36       <domain domain_ref="domain_dst" /> 
    3037       <axis axis_ref="axis_srf" /> 
    3138     </grid> 
  • XIOS/trunk/src/node/compute_connectivity_domain.cpp

    r934 r944  
    3838  void CComputeConnectivityDomain::checkValid(CDomain* domainDst) 
    3939  { 
    40     if (type.isEmpty()) 
     40    if (CDomain::type_attr::unstructured != domainDst->type) 
     41    { 
    4142      ERROR("CComputeConnectivityDomain::checkValid(CDomain* domainDst)", 
    42        << "Connectivity type is not defined. Chose 'node' or 'edge' for the type." 
    43        << "Domain " <<domainDst->getId() << std::endl 
    44        << "Connectivity object " << this->getId()); 
     43            << "Domain connectivity computation is only supported for unstructured" << std::endl 
     44            << "Check type of domain destination, id = " << domainDst->getId()); 
     45    } 
    4546 
     47    if (type.isEmpty()) type.setValue(CComputeConnectivityDomain::type_attr::edge); 
    4648    if (n_neighbor_max.isEmpty()) n_neighbor_max.setValue(0); 
    4749    if (n_neighbor.isEmpty()) n_neighbor.resize(domainDst->i_index.numElements()); 
  • XIOS/trunk/src/test/test_unstruct_complete.f90

    r934 r944  
    1717  INTEGER,PARAMETER :: ni_glo=100 
    1818  INTEGER,PARAMETER :: nj_glo=100 
    19   INTEGER,PARAMETER :: llm=5 
     19  INTEGER,PARAMETER :: llm=1 
    2020  DOUBLE PRECISION  :: lval(llm)=1 
    2121  TYPE(xios_field) :: field_hdl 
    2222  TYPE(xios_fieldgroup) :: fieldgroup_hdl 
    2323  TYPE(xios_file) :: file_hdl 
     24  TYPE(xios_domain) :: domain_hdl 
    2425  LOGICAL :: ok 
    2526 
     
    179180      IF (MOD(ind,8)>=0 .AND. MOD(ind,8)<2) THEN 
    180181        mask(ncell)=.FALSE. 
     182!        data_n_index=data_n_index+1 
    181183      ELSE 
    182184        mask(ncell)=.TRUE. 
     
    215217  CALL xios_close_context_definition() 
    216218 
     219  CALL xios_field_get_domain("field_A_expand",domain_hdl) 
     220  CALL xios_get_attr(domain_hdl,ni=ni) 
     221  PRINT *, "Domain destination ni = ",ni, "Domain source n= ",ncell 
    217222  CALL xios_get_compute_connectivity_domain_attr("compute", n_neighbor_max=nbMax) 
    218223  ALLOCATE(local_neighbor(nbMax,ncell)) 
    219224  CALL xios_get_compute_connectivity_domain_attr("compute", n_neighbor=n_local, local_neighbor=local_neighbor) 
    220  
    221    DO ts=1,24*10 
     225  PRINT *, "Connectivity max = ",nbMax 
     226 
     227   DO ts=1,1 
    222228     CALL xios_update_calendar(ts) 
    223229     CALL xios_send_field("field_A_srf",field_A_compressed) 
  • XIOS/trunk/src/transformation/domain_algorithm_compute_connectivity.cpp

    r934 r944  
    4949  CArray<int,1>& nbNeighbor = compute_connectivityDomain->n_neighbor; 
    5050  CArray<int,2>& localNeighbors = compute_connectivityDomain->local_neighbor; 
     51  int type = 1; // Edge type 
    5152  switch (compute_connectivityDomain->type) 
    5253  { 
    5354    case CComputeConnectivityDomain::type_attr::node : 
    54       computeNodeConnectivity(domainDestination, 
    55                               nbNeighborMax, 
    56                               nbNeighbor, 
    57                               localNeighbors); 
     55      type = 0; 
    5856      break; 
    5957    case CComputeConnectivityDomain::type_attr::edge : 
    60       computeEdgeConnectivity(domainDestination, 
    61                               nbNeighborMax, 
    62                               nbNeighbor, 
    63                               localNeighbors); 
     58      type = 1; 
    6459      break; 
    6560    default: 
    6661      break; 
    6762  } 
     63 
     64  computeLocalConnectivity(type, domainDestination, nbNeighborMax, nbNeighbor, localNeighbors); 
    6865} 
    6966 
    70 void CDomainAlgorithmComputeConnectivity::computeNodeConnectivity(CDomain* domain, 
     67/*! 
     68 *  Compute local connectivity of a domain 
     69 *  \param[in] type type of connectivity (node or edge) 
     70 *  \param[in] domain domain on which we calculate local connectivity 
     71 *  \param[in/out] nbConnectivityMax maximum number of neighbor a cell of domain has 
     72 *  \param[in/out] nbConnectivity number of neighbor a cell has 
     73 *  \param[in/out] localConnectivity localConnectivity local index of neighbor of a cell 
     74 */ 
     75void CDomainAlgorithmComputeConnectivity::computeLocalConnectivity(int type, 
     76                                                                  CDomain* domain, 
    7177                                                                  int& nbConnectivityMax, 
    7278                                                                  CArray<int,1>& nbConnectivity, 
    7379                                                                  CArray<int,2>& localConnectivity) 
    7480{ 
     81 
    7582  CMesh mesh; 
    76   CArray<double,1>& lon = domain->lonvalue_1d; 
    77   CArray<double,1>& lat = domain->latvalue_1d; 
     83 
    7884  CArray<double,2>& bounds_lon = domain->bounds_lon_1d; 
    7985  CArray<double,2>& bounds_lat = domain->bounds_lat_1d; 
    80   mesh.createMesh(lon, lat, bounds_lon, bounds_lat); 
    81   CArray<int, 2>& face_nodes = mesh.face_nodes; 
    82   int nvertex = mesh.nvertex; 
    83   int nFaces  = mesh.nbFaces; 
    84   int nNodes  = mesh.nbNodes; 
    85   std::vector<std::vector<size_t> > nodeFaceConnectivity(nNodes, std::vector<size_t>()); 
    86   for (int nf = 0; nf < nFaces; ++nf) 
    87   { 
    88     for (int nv = 0; nv < nvertex; ++nv) 
    89     { 
    90       nodeFaceConnectivity[face_nodes(nv,nf)].push_back(nf); 
    91     } 
    92   } 
     86  int ncell = bounds_lon.shape()[1]; 
     87  CArray<int,1> localIndex(ncell); 
     88  for (int idx = 0; idx <ncell; ++idx) localIndex(idx) = idx; 
    9389 
    94   if (nFaces != nbConnectivity.numElements()) nbConnectivity.resize(nFaces); 
    95   nbConnectivityMax = 1; 
    96   for (int nf = 0; nf < nFaces; ++nf) 
    97   { 
    98     std::set<size_t> neighFaces; 
    99     for (int nv = 0; nv < nvertex; ++nv) 
    100     { 
    101       std::vector<size_t>& tmpFaces = nodeFaceConnectivity[face_nodes(nv,nf)]; 
    102       for (int nFn = 0; nFn < tmpFaces.size(); ++nFn) 
    103       { 
    104         neighFaces.insert(tmpFaces[nFn]); 
    105       } 
    106     } 
    107     if (nbConnectivityMax < (neighFaces.size()-1)) nbConnectivityMax = neighFaces.size()-1; 
    108   } 
    109   if ((nbConnectivityMax != localConnectivity.shape()[0]) || (nFaces != localConnectivity.shape()[1])) 
    110     localConnectivity.resize(nbConnectivityMax, nFaces); 
    111  
    112   for (int nf = 0; nf < nFaces; ++nf) 
    113   { 
    114     std::set<size_t> neighFaces; 
    115     for (int nv = 0; nv < nvertex; ++nv) 
    116     { 
    117       std::vector<size_t>& tmpFaces = nodeFaceConnectivity[face_nodes(nv,nf)]; 
    118       for (int nFn = 0; nFn < tmpFaces.size(); ++nFn) 
    119       { 
    120         neighFaces.insert(tmpFaces[nFn]); 
    121       } 
    122     } 
    123  
    124     neighFaces.erase(nf); 
    125     nbConnectivity(nf) = neighFaces.size(); 
    126     std::set<size_t>::iterator it = neighFaces.begin(), ite = neighFaces.end(); 
    127     for (int idx = 0; it != ite; ++it, ++idx) 
    128     { 
    129       localConnectivity(idx, nf) = *it; 
    130     } 
    131   } 
     90  mesh.getLocalNghbFaces(type, localIndex, bounds_lon, bounds_lat, localConnectivity, nbConnectivity); 
     91  nbConnectivityMax = 0; 
     92  for (int idx =0; idx < nbConnectivity.numElements(); ++idx) 
     93    if (nbConnectivityMax < nbConnectivity(idx)) nbConnectivityMax = nbConnectivity(idx); 
    13294} 
    13395 
    134 void CDomainAlgorithmComputeConnectivity::computeEdgeConnectivity(CDomain* domain, 
    135                                                                   int& nbConnectivityMax, 
    136                                                                   CArray<int,1>& nbConnectivity, 
    137                                                                   CArray<int,2>& localConnectivity) 
    138 { 
    139   CMesh mesh; 
    140   CArray<double,1>& lon = domain->lonvalue_1d; 
    141   CArray<double,1>& lat = domain->latvalue_1d; 
    142   CArray<double,2>& bounds_lon = domain->bounds_lon_1d; 
    143   CArray<double,2>& bounds_lat = domain->bounds_lat_1d; 
    144   mesh.createMesh(lon, lat, bounds_lon, bounds_lat); 
    14596 
    146   CArray<int, 2>& face_edges = mesh.face_edges; 
    147   int nvertex = mesh.nvertex; 
    148   int nFaces  = mesh.nbFaces; 
    149   int nEdges  = mesh.nbEdges; 
    150   std::vector<std::vector<size_t> > edgeFaceConnectivity(nEdges, std::vector<size_t>()); 
    151   for (int nf = 0; nf < nFaces; ++nf) 
    152   { 
    153     for (int nv = 0; nv < nvertex; ++nv) 
    154     { 
    155       edgeFaceConnectivity[face_edges(nv,nf)].push_back(nf); 
    156     } 
    157   } 
    158  
    159   if (nFaces != nbConnectivity.numElements()) nbConnectivity.resize(nFaces); 
    160   nbConnectivityMax = 1; 
    161   for (int nf = 0; nf < nFaces; ++nf) 
    162   { 
    163     std::set<size_t> neighFaces; 
    164     for (int nv = 0; nv < nvertex; ++nv) 
    165     { 
    166       std::vector<size_t>& tmpFaces = edgeFaceConnectivity[face_edges(nv,nf)]; 
    167       for (int nFn = 0; nFn < tmpFaces.size(); ++nFn) 
    168       { 
    169         neighFaces.insert(tmpFaces[nFn]); 
    170       } 
    171     } 
    172     if (nbConnectivityMax < (neighFaces.size()-1)) nbConnectivityMax = neighFaces.size() - 1; 
    173   } 
    174   if ((nbConnectivityMax != localConnectivity.shape()[0]) || (nFaces != localConnectivity.shape()[1])) 
    175     localConnectivity.resize(nbConnectivityMax, nFaces); 
    176  
    177   for (int nf = 0; nf < nFaces; ++nf) 
    178   { 
    179     std::set<size_t> neighFaces; 
    180     for (int nv = 0; nv < nvertex; ++nv) 
    181     { 
    182       std::vector<size_t>& tmpFaces = edgeFaceConnectivity[face_edges(nv,nf)]; 
    183       for (int nFn = 0; nFn < tmpFaces.size(); ++nFn) 
    184       { 
    185         neighFaces.insert(tmpFaces[nFn]); 
    186       } 
    187     } 
    188  
    189     neighFaces.erase(nf); 
    190     nbConnectivity(nf) = neighFaces.size(); 
    191     std::set<size_t>::iterator it = neighFaces.begin(), ite = neighFaces.end(); 
    192     for (int idx = 0; it != ite; ++it, ++idx) 
    193     { 
    194       localConnectivity(idx, nf) = *it; 
    195     } 
    196   } 
    197  
    198 } 
    19997 
    20098/*! 
     
    203101void CDomainAlgorithmComputeConnectivity::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs) 
    204102{ 
    205  
    206103} 
    207104 
  • XIOS/trunk/src/transformation/domain_algorithm_compute_connectivity.hpp

    r934 r944  
    3232 
    3333protected: 
    34   void computeNodeConnectivity(CDomain* domain, int& nbConnectivityMax, 
    35                                CArray<int,1>& nbConnectivity, CArray<int,2>& localConnectivity); 
    36  
    37   void computeEdgeConnectivity(CDomain* domain, int& nbConnectivityMax, 
    38                                CArray<int,1>& nbConnectivity, CArray<int,2>& localConnectivity); 
     34  void computeLocalConnectivity(int type, 
     35                                CDomain* domain, 
     36                                int& nbConnectivityMax, 
     37                                CArray<int,1>& nbConnectivity, 
     38                                CArray<int,2>& localConnectivity); 
    3939 
    4040private: 
Note: See TracChangeset for help on using the changeset viewer.