Changeset 688


Ignore:
Timestamp:
09/15/15 17:31:06 (9 years ago)
Author:
mhnguyen
Message:

Integrating remap library into XIOS

+) Change name of some files of remap library to be compatible with XIOS
+) Implement function to fill in automatically boundary longitude and latitude

Test
+) On Curie
+) test_remap correct

Location:
XIOS/trunk
Files:
56 added
10 edited

Legend:

Unmodified
Added
Removed

  • XIOS/trunk/bld.cfg

    r660 r688  
    2828src::blitz $PWD/extern/blitz/src 
    2929src::netcdf $PWD/extern/netcdf4 
     30src::remap $PWD/extern/remap/src 
    3031bld::lib xios 
    31 bld::target libxios.a 
     32bld::lib::remap mapper 
     33bld::target libxios.a libmapper.a 
    3234#bld::target generate_fortran_interface.exe  
    3335bld::target xios_server.exe test_remap.exe 
    34 bld::target test_new_features.exe test_unstruct_complete.exe  
    35 bld::target test_client.exe test_complete.exe 
     36#bld::target test_new_features.exe test_unstruct_complete.exe  
     37#bld::target test_client.exe test_complete.exe 
    3638bld::exe_dep 
    3739 

  • XIOS/trunk/inputs/REMAP/iodef.xml

    r666 r688  
    2727   <domain_definition> 
    2828     <domain id="src_domain" /> 
    29      <domain id="dst_domain" domain_ref="src_domain"> 
     29     <domain id="dst_domain" domain_src="src_domain"> 
     30       <generate_rectilinear_domain /> 
    3031       <interpolate_from_file_domain file="weight.nc" /> 
    3132     </domain> 

  • XIOS/trunk/src/filter/spatial_transform_filter.cpp

    r653 r688  
    8585    double* sendBuff; 
    8686    if (0 != sendBuffSize) sendBuff = new double[sendBuffSize]; 
     87    std::vector<MPI_Request> sendRequest; 
    8788    for (itSend = itbSend; itSend != iteSend; ++itSend) 
    8889    { 
    … …  
    9495        sendBuff[idx] = dataSrc(localIndex_p(idx)); 
    9596      } 
    96       MPI_Send(sendBuff, countSize, MPI_DOUBLE, destRank, 12, client->intraComm); 
     97      sendRequest.push_back(MPI_Request()); 
     98      MPI_Isend(sendBuff, countSize, MPI_DOUBLE, destRank, 12, client->intraComm, &sendRequest.back()); 
    9799    } 
    98100 
    … …  
    123125    } 
    124126 
     127    std::vector<MPI_Status> requestStatus(sendRequest.size()); 
     128    MPI_Wait(&sendRequest[0], &requestStatus[0]); 
    125129    if (0 != sendBuffSize) delete [] sendBuff; 
    126130    if (0 != recvBuffSize) delete [] recvBuff; 

  • XIOS/trunk/src/node/domain.cpp

    r687 r688  
    215215        else 
    216216        { 
     217          float njGlo = nj_glo.getValue(); 
     218          float niGlo = ni_glo.getValue(); 
     219          int nbProcOnX, nbProcOnY, range; 
     220 
    217221          // Compute (approximately) number of segment on x and y axis 
    218           float yOverXRatio = (nj_glo.getValue())/(ni_glo.getValue()); 
    219           int nbProcOnX, nbProcOnY, range; 
     222          float yOverXRatio = njGlo/niGlo; 
     223 
    220224          nbProcOnX = std::ceil(std::sqrt(nbLocalDomain/yOverXRatio)); 
    221225          nbProcOnY = std::ceil(((float)nbLocalDomain)/nbProcOnX); 
    … …  
    266270        // Now fill other attributes 
    267271        fillInRectilinearLonLat(); 
     272//        fillInRectilinearBoundLonLat(); 
    268273     } 
    269274   } 
    … …  
    280285     lonvalue_1d.resize(ni); 
    281286     latvalue_1d.resize(nj); 
     287 
    282288     double lonStep = double(360/ni_glo.getValue()); 
    283289     double latStep = double(180/nj_glo.getValue()); 
    … …  
    286292     for (int i = 0; i < ni; ++i) 
    287293     { 
    288        lonvalue_1d(i) = static_cast<double>(ibegin + i) * lonStep; 
     294       lonvalue_1d(i) = static_cast<double>(ibegin + i) * lonStep -180 + lonStep/2; 
    289295     } 
    290296 
    291297     for (int j = 0; j < nj; ++j) 
    292298     { 
    293        latvalue_1d(j) = static_cast<double>(jbegin + j) * latStep - 90; 
    294      } 
     299       latvalue_1d(j) = static_cast<double>(jbegin + j) * latStep - 90 + latStep/2; 
     300     } 
     301   } 
     302 
     303   void CDomain::fillInRectilinearBoundLonLat(CArray<double,2>& boundsLon, CArray<double,2>& boundsLat) 
     304   { 
     305     int i,j,k; 
     306     const int nvertexValue = 4; 
     307 
     308     boundsLon.resize(nvertexValue,ni*nj); 
     309     boundsLat.resize(nvertexValue,nj*ni); 
     310 
     311     double lonStep = double(360/ni_glo.getValue()); 
     312     double latStep = double(180/nj_glo.getValue()); 
     313 
     314     for(j=0;j<nj;++j) 
     315       for(i=0;i<ni;++i) 
     316       { 
     317         k=j*ni+i; 
     318         boundsLon(0,k) = boundsLon(1,k) = (0 != (ibegin + i)) ? (ibegin + i) * lonStep -180 
     319                                                               : -180; 
     320         boundsLon(2,k) = boundsLon(3,k) = ((ibegin + i + 1) != ni_glo) ? (ibegin + i +1) * lonStep -180 
     321                                                                        : 180; 
     322 
     323         boundsLat(1,k) = boundsLat(2,k) = (0 != (jbegin + j)) ? (jbegin + j) * latStep - 90 
     324                                                               : -90; 
     325         boundsLat(0,k) = boundsLat(3,k) = ((jbegin + j +1) != nj_glo) ? (jbegin + j +1) * latStep - 90 
     326                                                               : +90; 
     327       } 
    295328   } 
    296329 
    … …  
    625658            } 
    626659          } 
     660        lonvalue_1d.free(); 
    627661     } 
    628662 
    … …  
    666700         } 
    667701       } 
     702       lonvalue_2d.free(); 
    668703     } 
    669704 

  • XIOS/trunk/src/node/domain.hpp

    r687 r688  
    124124         void sendLonLatArea(void); 
    125125         void computeConnectedServer(void) ; 
     126         void fillInRectilinearBoundLonLat(CArray<double,2>& boundsLon, CArray<double,2>& boundsLat); 
    126127 
    127128         static bool dispatchEvent(CEventServer& event); 

  • XIOS/trunk/src/test/test_remap.f90

    r673 r688  
    115115                            bounds_lon_1D=src_boundslon, bounds_lat_1D=src_boundslat, nvertex=src_nvertex) 
    116116 
    117   CALL xios_set_domain_attr("dst_domain", ni_glo=dst_ni_glo, ibegin=dst_ibegin, ni=dst_ni, type="unstructured") 
     117!  CALL xios_set_domain_attr("dst_domain", ni_glo=dst_ni_glo, ibegin=dst_ibegin, ni=dst_ni, type="unstructured") 
    118118  CALL xios_set_domain_attr("dst_domain", lonvalue_1D=dst_lon, latvalue_1D=dst_lat, & 
    119119                            bounds_lon_1D=dst_boundslon, bounds_lat_1D=dst_boundslat, nvertex=dst_nvertex) 
     120 
    120121 
    121122  dtime%second = 3600 

  • XIOS/trunk/src/timer.hpp

    r652 r688  
    1 #ifndef __TIMER_HPP__ 
    2 #define __TIMER_HPP__ 
     1#ifndef __XIOS_TIMER_HPP__ 
     2#define __XIOS_TIMER_HPP__ 
    33 
    44#include <string> 
    … …  
    1414      bool suspended; 
    1515      std::string name; 
    16        
     16 
    1717      CTimer(const std::string& name); 
    1818      void suspend(void); 

  • XIOS/trunk/src/transformation/domain_algorithm_interpolate_from_file.cpp

    r663 r688  
    33   \author Ha NGUYEN 
    44   \since 09 Jul 2015 
    5    \date 17 Jul 2015 
     5   \date 11 Sep 2015 
    66 
    77   \brief Algorithm for interpolation on a domain. 
    … …  
    1414#include "client_server_mapping_distributed.hpp" 
    1515#include "netcdf.hpp" 
     16#include "mapper.hpp" 
    1617 
    1718namespace xios { 
    … …  
    2324  computeIndexSourceMapping(); 
    2425} 
     26 
     27void CDomainAlgorithmInterpolateFromFile::computeRemap() 
     28{ 
     29  using namespace sphereRemap; 
     30 
     31  CContext* context = CContext::getCurrent(); 
     32  CContextClient* client=context->client; 
     33  int clientRank = client->clientRank; 
     34  int i, j, k, idx; 
     35        double srcPole[] = {0, 0, 0}; 
     36        double dstPole[] = {0, 0, 1}; 
     37        int orderInterp = 2; 
     38 
     39  int constNVertex = 4; // Value by default number of vertex for rectangular domain 
     40  int nVertexSrc, nVertexDest; 
     41  nVertexSrc = nVertexDest = constNVertex; 
     42 
     43  // First of all, try to retrieve the boundary values of domain source and domain destination 
     44  int localDomainSrcSize = domainSrc_->i_index.numElements(); 
     45  int niSrc = domainSrc_->ni.getValue(), njSrc = domainSrc_->nj.getValue(); 
     46  bool hasBoundSrc = domainSrc_->hasBounds; 
     47  if (hasBoundSrc) nVertexSrc = domainSrc_->nvertex.getValue(); 
     48  CArray<double,2> boundsLonSrc(nVertexSrc,localDomainSrcSize); 
     49  CArray<double,2> boundsLatSrc(nVertexSrc,localDomainSrcSize); 
     50 
     51  if (hasBoundSrc)  // Suppose that domain source is curvilinear or unstructured 
     52  { 
     53    if (!domainSrc_->bounds_lon_2d.isEmpty()) 
     54    { 
     55      for (j = 0; j < njSrc; ++j) 
     56        for (i = 0; i < niSrc; ++i) 
     57        { 
     58          k=j*niSrc+i; 
     59          for(int n=0;n<nVertexSrc;++n) 
     60          { 
     61            boundsLonSrc(n,k) = domainSrc_->bounds_lon_2d(n,i,j); 
     62            boundsLatSrc(n,k) = domainSrc_->bounds_lat_2d(n,i,j); 
     63          } 
     64        } 
     65    } 
     66    else 
     67    { 
     68      boundsLonSrc = domainSrc_->bounds_lon_1d; 
     69      boundsLatSrc = domainSrc_->bounds_lat_1d; 
     70    } 
     71  } 
     72  else // if domain source is rectilinear, not do anything now 
     73  { 
     74    nVertexSrc = constNVertex; 
     75  } 
     76 
     77  int localDomainDestSize = domainDest_->i_index.numElements(); 
     78  int niDest = domainDest_->ni.getValue(), njDest = domainDest_->nj.getValue(); 
     79  bool hasBoundDest = domainDest_->hasBounds; 
     80  if (hasBoundDest) nVertexDest = domainDest_->nvertex.getValue(); 
     81  CArray<double,2> boundsLonDest(nVertexDest,localDomainDestSize); 
     82  CArray<double,2> boundsLatDest(nVertexDest,localDomainDestSize); 
     83 
     84  if (hasBoundDest) 
     85  { 
     86    if (!domainDest_->bounds_lon_2d.isEmpty()) 
     87    { 
     88      for (j = 0; j < njDest; ++j) 
     89        for (i = 0; i < niDest; ++i) 
     90        { 
     91          k=j*niDest+i; 
     92          for(int n=0;n<nVertexDest;++n) 
     93          { 
     94            boundsLonDest(n,k) = domainDest_->bounds_lon_2d(n,i,j); 
     95            boundsLatDest(n,k) = domainDest_->bounds_lat_2d(n,i,j); 
     96          } 
     97        } 
     98    } 
     99    else 
     100    { 
     101      boundsLonDest = domainDest_->bounds_lon_1d; 
     102      boundsLatDest = domainDest_->bounds_lat_1d; 
     103    } 
     104  } 
     105  else 
     106  { 
     107    // Ok, fill in boundary values for rectangular domain 
     108    domainDest_->fillInRectilinearBoundLonLat(boundsLonDest, boundsLatDest); 
     109    nVertexDest = 4; 
     110  } 
     111 
     112  // Ok, now use mapper to calculate 
     113  int nSrcLocal = domainSrc_->i_index.numElements(); 
     114  int nDstLocal = domainDest_->i_index.numElements(); 
     115  Mapper mapper(client->intraComm); 
     116  mapper.setVerbosity(PROGRESS) ; 
     117  mapper.setSourceMesh(boundsLonSrc.dataFirst(), boundsLatSrc.dataFirst(), nVertexSrc, nSrcLocal, srcPole); 
     118  mapper.setTargetMesh(boundsLonDest.dataFirst(), boundsLatDest.dataFirst(), nVertexDest, nDstLocal, dstPole); 
     119  std::vector<double> timings = mapper.computeWeights(orderInterp); 
     120 
     121  for (int idx = 0;  idx < mapper.nWeights; ++idx) 
     122  { 
     123    transformationMapping_[mapper.targetWeightId[idx]].push_back(mapper.sourceWeightId[idx]); 
     124    transformationWeight_[mapper.targetWeightId[idx]].push_back(mapper.remapMatrix[idx]); 
     125  } 
     126} 
     127 
     128void CDomainAlgorithmInterpolateFromFile::computeIndexSourceMapping() 
     129{ 
     130  computeRemap(); 
     131} 
     132 
    25133 
    26134/*! 
    27135  Compute the index mapping between domain on grid source and one on grid destination 
    28136*/ 
    29 void CDomainAlgorithmInterpolateFromFile::computeIndexSourceMapping() 
     137void CDomainAlgorithmInterpolateFromFile::readRemapInfo() 
    30138{ 
    31139  CContext* context = CContext::getCurrent(); 
    … …  
    37145  readInterpolationInfo(filename, interpMapValue); 
    38146 
    39   //randomizeInterpolationInfo(interpMapValue); 
    40147  boost::unordered_map<size_t,int> globalIndexOfDomainDest; 
    41148  int ni = domainDest_->ni.getValue(); 
    … …  
    199306} 
    200307 
    201 void CDomainAlgorithmInterpolateFromFile::randomizeInterpolationInfo(std::map<int,std::vector<std::pair<int,double> > >& interpMapValue) 
    202 { 
    203   int iDestBegin = 2; 
    204   int jDestBegin = 2; 
    205  
    206   int ni_glo_src = domainSrc_->ni_glo.getValue(); 
    207   int ni_glo = domainDest_->ni_glo.getValue(); 
    208   size_t globalIndex; 
    209   int nIndexSize = domainDest_->i_index.numElements(), i_ind, j_ind; 
    210   for (int idx = 0; idx < nIndexSize; ++idx) 
    211   { 
    212     i_ind=domainDest_->i_index(idx) ; 
    213     j_ind=domainDest_->j_index(idx) ; 
    214  
    215     globalIndex = i_ind + j_ind * ni_glo; 
    216  
    217     interpMapValue[globalIndex].push_back(make_pair((i_ind+iDestBegin) + (j_ind+jDestBegin-1)*ni_glo_src, 0.25)); 
    218     interpMapValue[globalIndex].push_back(make_pair((i_ind+iDestBegin+1) + (j_ind+jDestBegin)*ni_glo_src, 0.25)); 
    219     interpMapValue[globalIndex].push_back(make_pair((i_ind+iDestBegin) + (j_ind+jDestBegin+1)*ni_glo_src, 0.25)); 
    220     interpMapValue[globalIndex].push_back(make_pair((i_ind+iDestBegin-1) + (j_ind+jDestBegin)*ni_glo_src, 0.25)); 
    221  
    222  
    223   } 
    224 } 
    225  
    226308/*! 
    227309  Read interpolation information from a file 

  • XIOS/trunk/src/transformation/domain_algorithm_interpolate_from_file.hpp

    r657 r688  
    33   \author Ha NGUYEN 
    44   \since 09 July 2015 
    5    \date 09 July 2015 
     5   \date 09 Sep 2015 
    66 
    77   \brief Algorithm for interpolation on a domain. 
    … …  
    3030private: 
    3131  void readInterpolationInfo(std::string& filename, std::map<int,std::vector<std::pair<int,double> > >& interpMapValue); 
    32   void randomizeInterpolationInfo(std::map<int,std::vector<std::pair<int,double> > >& interpMapValue); 
     32  void computeRemap(); 
     33  void readRemapInfo(); 
     34 
    3335private: 
    3436  CInterpolateFromFileDomain* interpDomain_; 

  • XIOS/trunk/src/transformation/grid_transformation.cpp

    r687 r688  
    5454    CDomain* domain = CDomain::createDomain(); 
    5555    domain->setAttributes(domainSrcTmp[idx]); 
     56    domain->checkAttributesOnClient(); 
    5657    domainSrc.push_back(domain); 
    5758  } 
Note: See TracChangeset for help on using the changeset viewer.