Changeset 743


Ignore:
Timestamp:
10/19/15 17:41:35 (6 years ago)
Author:
mhnguyen
Message:

Processing pole of rectangular grid with latitude = 90

+) Use average value of all points with latitude 90 (-90)

Test
+) On Curie
+) test_remap pass

Location:
XIOS/trunk
Files:
5 edited

Legend:

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

    r742 r743  
    1212     <field id="dst_field_regular_grid"  operation="instant" field_ref="tmp_field" grid_ref="dst_grid_regular"  read_access="true"/> 
    1313     <field id="dst_field_regular"  operation="instant" field_ref="tmp_field" domain_ref="dst_domain_regular"  read_access="true"/> 
     14     <field id="dst_field_regular_pole"  operation="instant" field_ref="src_field" domain_ref="dst_domain_regular_pole"  read_access="true"/> 
    1415 
    1516   </field_definition> 
     
    2223     <file id="output_dst" name="output_dst" type="one_file"> 
    2324        <field field_ref="dst_field" name="field" /> 
     25     </file> 
     26     <file id="out_dst_regular_pole" name="out_dst_regular_pole" type="one_file"> 
     27        <field field_ref="dst_field_regular_pole" name="field" /> 
    2428     </file> 
    2529     <file id="output_dst_regular" name="output_dst_regular" type="one_file"> 
     
    4246     <domain id="src_domain" /> 
    4347     <domain id="dst_domain" domain_src="src_domain"> 
     48       <interpolate_domain/> 
     49     </domain> 
     50     <domain id="dst_domain_regular_pole" domain_src="src_domain" ni_glo="180" nj_glo="90" type="rectilinear"> 
     51       <generate_rectilinear_domain bounds_lat_start="-90" bounds_lat_end="90" lon_start="2" lon_end="360"/> 
    4452       <interpolate_domain/> 
    4553     </domain> 
  • XIOS/trunk/src/node/domain.cpp

    r734 r743  
    375375   } 
    376376 
    377    void CDomain::fillInRectilinearBoundLonLat(CArray<double,2>& boundsLon, CArray<double,2>& boundsLat) 
     377   void CDomain::fillInRectilinearBoundLonLat(CArray<double,2>& boundsLon, CArray<double,2>& boundsLat, 
     378                                              bool isNorthPole, bool isSouthPole) 
    378379   { 
    379380     int i,j,k; 
     
    397398         boundsLon(2,k) = boundsLon(3,k) = ((ibegin + i + 1) != ni_glo) ? (ibegin + i +1) * lonStep + bounds_lon_start 
    398399                                                                        : bounds_lon_end; 
    399  
     400       } 
     401 
     402     double bounds_lat_start_pole = bounds_lat_start; 
     403     double bounds_lat_end_pole   = bounds_lat_end; 
     404     if (isNorthPole) bounds_lat_start_pole = lat_start; 
     405     if (isSouthPole) bounds_lat_end_pole   = lat_end; 
     406 
     407     for(j=0;j<nj;++j) 
     408       for(i=0;i<ni;++i) 
     409       { 
     410         k=j*ni+i; 
    400411         boundsLat(1,k) = boundsLat(2,k) = (0 != (jbegin + j)) ? (jbegin + j) * latStep + bounds_lat_start 
    401                                                                : bounds_lat_start; 
     412                                                               : bounds_lat_start_pole; 
    402413         boundsLat(0,k) = boundsLat(3,k) = ((jbegin + j +1) != nj_glo) ? (jbegin + j +1) * latStep + bounds_lat_start 
    403                                                                : bounds_lat_end; 
     414                                                               : bounds_lat_end_pole; 
    404415       } 
    405416   } 
  • XIOS/trunk/src/node/domain.hpp

    r731 r743  
    131131         void sendLonLatArea(void); 
    132132         void computeConnectedServer(void) ; 
    133          void fillInRectilinearBoundLonLat(CArray<double,2>& boundsLon, CArray<double,2>& boundsLat); 
     133         void fillInRectilinearBoundLonLat(CArray<double,2>& boundsLon, CArray<double,2>& boundsLat, 
     134                                           bool isNorthPole = false, bool isSouthPole = false); 
    134135 
    135136         static bool dispatchEvent(CEventServer& event); 
  • XIOS/trunk/src/transformation/domain_algorithm_interpolate.cpp

    r709 r743  
    3939        int orderInterp = interpDomain_->order.getValue(); 
    4040 
    41   int constNVertex = 4; // Value by default number of vertex for rectangular domain 
     41  const double poleValue = 90.0; 
     42  const int constNVertex = 4; // Value by default number of vertex for rectangular domain 
    4243  int nVertexSrc, nVertexDest; 
    4344  nVertexSrc = nVertexDest = constNVertex; 
     
    7980  } 
    8081 
     82  bool isNorthPole = false; 
     83  bool isSouthPole = false; 
     84  std::map<int,std::vector<std::pair<int,double> > > interpMapValueNorthPole; 
     85  std::map<int,std::vector<std::pair<int,double> > > interpMapValueSouthPole; 
     86 
    8187  int localDomainDestSize = domainDest_->i_index.numElements(); 
    8288  int niDest = domainDest_->ni.getValue(), njDest = domainDest_->nj.getValue(); 
     
    110116  else 
    111117  { 
     118    if (poleValue == std::abs(domainDest_->lat_start)) isNorthPole = true; 
     119    if (poleValue == std::abs(domainDest_->lat_end)) isSouthPole = true; 
     120    if (isNorthPole && (0 == domainDest_->jbegin.getValue())) 
     121    { 
     122      int ibegin = domainDest_->ibegin.getValue(); 
     123      for (i = 0; i < niDest; ++i) 
     124      { 
     125        interpMapValueNorthPole[i+ibegin]; 
     126      } 
     127    } 
     128 
     129    if (isSouthPole && (domainDest_->nj_glo.getValue() == (domainDest_->jbegin.getValue() + njDest))) 
     130    { 
     131      int ibegin = domainDest_->ibegin.getValue(); 
     132      int njGlo = domainDest_->nj_glo.getValue(); 
     133      int niGlo = domainDest_->ni_glo.getValue(); 
     134      for (i = 0; i < niDest; ++i) 
     135      { 
     136        k = (njGlo - 1)*niGlo + i + ibegin; 
     137        interpMapValueSouthPole[k]; 
     138      } 
     139    } 
     140 
    112141    // Ok, fill in boundary values for rectangular domain 
    113     domainDest_->fillInRectilinearBoundLonLat(boundsLonDest, boundsLatDest); 
    114142    nVertexDest = constNVertex; 
     143    domainDest_->fillInRectilinearBoundLonLat(boundsLonDest, boundsLatDest, isNorthPole, isSouthPole); 
    115144  } 
    116145 
     
    152181 
    153182  std::map<int,std::vector<std::pair<int,double> > > interpMapValue; 
     183  std::map<int,std::vector<std::pair<int,double> > >::const_iterator iteNorthPole = interpMapValueNorthPole.end(), 
     184                                                                     iteSouthPole = interpMapValueSouthPole.end(); 
    154185  for (int idx = 0;  idx < mapper.nWeights; ++idx) 
    155186  { 
    156187    interpMapValue[mapper.targetWeightId[idx]].push_back(make_pair(mapper.sourceWeightId[idx],mapper.remapMatrix[idx])); 
    157   } 
     188    if (iteNorthPole != interpMapValueNorthPole.find(mapper.targetWeightId[idx])) 
     189    { 
     190      interpMapValueNorthPole[mapper.targetWeightId[idx]].push_back(make_pair(mapper.sourceWeightId[idx],mapper.remapMatrix[idx])); 
     191    } 
     192 
     193    if (iteSouthPole != interpMapValueSouthPole.find(mapper.targetWeightId[idx])) 
     194    { 
     195      interpMapValueSouthPole[mapper.targetWeightId[idx]].push_back(make_pair(mapper.sourceWeightId[idx],mapper.remapMatrix[idx])); 
     196    } 
     197  } 
     198  int niGloDst = domainDest_->ni_glo.getValue(); 
     199  processPole(interpMapValueNorthPole, niGloDst); 
     200  processPole(interpMapValueSouthPole, niGloDst); 
     201 
     202  if (!interpMapValueNorthPole.empty()) 
     203  { 
     204     std::map<int,std::vector<std::pair<int,double> > >::iterator itNorthPole = interpMapValueNorthPole.begin(); 
     205     for (; itNorthPole != iteNorthPole; ++itNorthPole) 
     206     { 
     207       if (!(itNorthPole->second.empty())) 
     208        itNorthPole->second.swap(interpMapValue[itNorthPole->first]); 
     209     } 
     210  } 
     211 
     212  if (!interpMapValueSouthPole.empty()) 
     213  { 
     214     std::map<int,std::vector<std::pair<int,double> > >::iterator itSouthPole = interpMapValueSouthPole.begin(); 
     215     for (; itSouthPole != iteSouthPole; ++itSouthPole) 
     216     { 
     217       if (!(itSouthPole->second.empty())) 
     218        itSouthPole->second.swap(interpMapValue[itSouthPole->first]); 
     219     } 
     220  } 
     221 
    158222  exchangeRemapInfo(interpMapValue); 
    159223 
    160224  delete [] globalSrc; 
    161225  delete [] globalDst; 
     226} 
     227 
     228void CDomainAlgorithmInterpolate::processPole(std::map<int,std::vector<std::pair<int,double> > >& interMapValuePole, 
     229                                              int nbGlobalPointOnPole) 
     230{ 
     231  CContext* context = CContext::getCurrent(); 
     232  CContextClient* client=context->client; 
     233 
     234  MPI_Comm poleComme(MPI_COMM_NULL); 
     235  MPI_Comm_split(client->intraComm, interMapValuePole.empty() ? MPI_UNDEFINED : 1, 0, &poleComme); 
     236  if (MPI_COMM_NULL != poleComme) 
     237  { 
     238    int nbClientPole; 
     239    MPI_Comm_size(poleComme, &nbClientPole); 
     240 
     241    std::map<int,std::vector<std::pair<int,double> > >::iterator itePole = interMapValuePole.end(), itPole, 
     242                                                                 itbPole = interMapValuePole.begin(); 
     243 
     244    int nbWeight = 0; 
     245    for (itPole = itbPole; itPole != itePole; ++itPole) 
     246       nbWeight += itPole->second.size(); 
     247 
     248    std::vector<int> recvCount(nbClientPole,0); 
     249    std::vector<int> displ(nbClientPole,0); 
     250    MPI_Allgather(&nbWeight,1,MPI_INT,&recvCount[0],1,MPI_INT,poleComme) ; 
     251 
     252    displ[0]=0; 
     253    for(int n=1;n<nbClientPole;++n) displ[n]=displ[n-1]+recvCount[n-1] ; 
     254    int recvSize=displ[nbClientPole-1]+recvCount[nbClientPole-1] ; 
     255 
     256    std::vector<int> sendSourceIndexBuff(nbWeight); 
     257    std::vector<double> sendSourceWeightBuff(nbWeight); 
     258    int k = 0; 
     259    for (itPole = itbPole; itPole != itePole; ++itPole) 
     260    { 
     261      for (int idx = 0; idx < itPole->second.size(); ++idx) 
     262      { 
     263        sendSourceIndexBuff[k] = (itPole->second)[idx].first; 
     264        sendSourceWeightBuff[k] = (itPole->second)[idx].second; 
     265        ++k; 
     266      } 
     267    } 
     268 
     269    std::vector<int> recvSourceIndexBuff(recvSize); 
     270    std::vector<double> recvSourceWeightBuff(recvSize); 
     271 
     272    // Gather all index and weight for pole 
     273    MPI_Allgatherv(&sendSourceIndexBuff[0],nbWeight,MPI_INT,&recvSourceIndexBuff[0],&recvCount[0],&displ[0],MPI_INT,poleComme); 
     274    MPI_Allgatherv(&sendSourceWeightBuff[0],nbWeight,MPI_DOUBLE,&recvSourceWeightBuff[0],&recvCount[0],&displ[0],MPI_DOUBLE,poleComme); 
     275 
     276    std::map<int,double> recvTemp; 
     277    for (int idx = 0; idx < recvSize; ++idx) 
     278    { 
     279      if (recvTemp.end() != recvTemp.find(recvSourceIndexBuff[idx])) 
     280        recvTemp[recvSourceIndexBuff[idx]] += recvSourceWeightBuff[idx]/nbGlobalPointOnPole; 
     281      else 
     282        recvTemp[recvSourceIndexBuff[idx]] = 0.0; 
     283    } 
     284 
     285    std::map<int,double>::const_iterator itRecvTemp, itbRecvTemp = recvTemp.begin(), iteRecvTemp = recvTemp.end(); 
     286 
     287    for (itPole = itbPole; itPole != itePole; ++itPole) 
     288    { 
     289      itPole->second.clear(); 
     290      for (itRecvTemp = itbRecvTemp; itRecvTemp != iteRecvTemp; ++itRecvTemp) 
     291          itPole->second.push_back(make_pair(itRecvTemp->first, itRecvTemp->second)); 
     292    } 
     293  } 
     294 
    162295} 
    163296 
  • XIOS/trunk/src/transformation/domain_algorithm_interpolate.hpp

    r709 r743  
    3030private: 
    3131  void readInterpolationInfo(std::string& filename, std::map<int,std::vector<std::pair<int,double> > >& interpMapValue); 
     32  void processPole(std::map<int,std::vector<std::pair<int,double> > >& interMapValuePole, 
     33                   int nbGlobalPointOnPole); 
    3234  void computeRemap(); 
    3335  void readRemapInfo(); 
Note: See TracChangeset for help on using the changeset viewer.