Changeset 809


Ignore:
Timestamp:
12/14/15 09:21:50 (8 years ago)
Author:
ymipsl
Message:

Fix problems for interpolation onto regular domain.

YM

Location:
XIOS/trunk/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r808 r809  
    457457   } 
    458458 
    459    void CDomain::fillInRectilinearBoundLonLat(CArray<double,2>& boundsLon, CArray<double,2>& boundsLat, 
    460                                               bool isNorthPole, bool isSouthPole) 
     459 
     460 
     461   void CDomain::AllgatherRectilinearLonLat(CArray<double,1>& lon, CArray<double,1>& lat, CArray<double,1>& lon_g, CArray<double,1>& lat_g) 
     462   { 
     463          CContext* context = CContext::getCurrent(); 
     464      CContextClient* client = context->client; 
     465          lon_g.resize(ni_glo) ; 
     466          lat_g.resize(nj_glo) ; 
     467           
     468           
     469          int* ibegin_g = new int[client->clientSize] ; 
     470          int* jbegin_g = new int[client->clientSize] ; 
     471          int* ni_g = new int[client->clientSize] ; 
     472          int* nj_g = new int[client->clientSize] ; 
     473          int v ; 
     474          v=ibegin ; 
     475          MPI_Allgather(&v,1,MPI_INT,ibegin_g,1,MPI_INT,client->intraComm) ; 
     476          v=jbegin ; 
     477          MPI_Allgather(&v,1,MPI_INT,jbegin_g,1,MPI_INT,client->intraComm) ; 
     478          v=ni ; 
     479          MPI_Allgather(&v,1,MPI_INT,ni_g,1,MPI_INT,client->intraComm) ; 
     480          v=nj ; 
     481          MPI_Allgather(&v,1,MPI_INT,nj_g,1,MPI_INT,client->intraComm) ; 
     482           
     483          MPI_Allgatherv(lon.dataFirst(),ni,MPI_DOUBLE,lon_g.dataFirst(),ni_g, ibegin_g,MPI_DOUBLE,client->intraComm) ; 
     484          MPI_Allgatherv(lat.dataFirst(),nj,MPI_DOUBLE,lat_g.dataFirst(),nj_g, jbegin_g,MPI_DOUBLE,client->intraComm) ; 
     485    
     486      delete[] ibegin_g ; 
     487      delete[] jbegin_g ; 
     488      delete[] ni_g ; 
     489      delete[] nj_g ; 
     490   } 
     491 
     492   void CDomain::fillInRectilinearBoundLonLat(CArray<double,1>& lon, CArray<double,1>& lat, 
     493                                              CArray<double,2>& boundsLon, CArray<double,2>& boundsLat) 
    461494   { 
    462495     int i,j,k; 
     496      
    463497     const int nvertexValue = 4; 
    464498     boundsLon.resize(nvertexValue,ni*nj); 
    465499 
    466      if (!lonvalue_rectilinear_read_from_file.isEmpty()) 
    467      { 
    468        double lonStepStart = lonvalue_rectilinear_read_from_file(1)-lonvalue_rectilinear_read_from_file(0); 
    469        bounds_lon_start.setValue(lonvalue_rectilinear_read_from_file(0) - lonStepStart/2); 
    470        double lonStepEnd = (lonvalue_rectilinear_read_from_file(ni_glo-1)-lonvalue_rectilinear_read_from_file(ni_glo-2)); 
    471        bounds_lon_end.setValue(lonvalue_rectilinear_read_from_file(ni_glo-1) + lonStepEnd/2); 
    472        double errorBoundsLon = std::abs(360-std::abs(bounds_lon_end-bounds_lon_start)); 
    473        if (errorBoundsLon > NumTraits<double>::epsilon()) bounds_lon_end.setValue(bounds_lon_start+360); 
    474        for(j=0;j<nj;++j) 
    475          for(i=0;i<ni;++i) 
    476          { 
    477            k=j*ni+i; 
    478            boundsLon(0,k) = boundsLon(1,k) = (0 == (ibegin + i)) ? bounds_lon_start 
    479                                                                  : (lonvalue_rectilinear_read_from_file(ibegin + i)+lonvalue_rectilinear_read_from_file(ibegin + i-1))/2; 
    480            boundsLon(2,k) = boundsLon(3,k) = ((ibegin + i + 1) == ni_glo) ? bounds_lon_end 
    481                                                                           : (lonvalue_rectilinear_read_from_file(ibegin + i + 1)+lonvalue_rectilinear_read_from_file(ibegin + i))/2; 
    482          } 
    483      } 
    484      else 
    485      { 
    486        double boundsLonRange = bounds_lon_end - bounds_lon_start; 
    487        double lonStep = boundsLonRange/double(ni_glo.getValue()); 
    488        for(j=0;j<nj;++j) 
    489          for(i=0;i<ni;++i) 
    490          { 
    491            k=j*ni+i; 
    492            boundsLon(0,k) = boundsLon(1,k) = (0 != (ibegin + i)) ? (ibegin + i) * lonStep + bounds_lon_start 
    493                                                                  : bounds_lon_start; 
    494            boundsLon(2,k) = boundsLon(3,k) = ((ibegin + i + 1) != ni_glo) ? (ibegin + i +1) * lonStep + bounds_lon_start 
    495                                                                           : bounds_lon_end; 
    496          } 
    497      } 
    498  
    499      boundsLat.resize(nvertexValue,nj*ni); 
    500      if (!latvalue_rectilinear_read_from_file.isEmpty()) 
    501      { 
    502        double latStepStart = latvalue_rectilinear_read_from_file(1)-latvalue_rectilinear_read_from_file(0); 
    503        if (isNorthPole) bounds_lat_start.setValue(latvalue_rectilinear_read_from_file(0) ); 
    504        else bounds_lat_start.setValue(latvalue_rectilinear_read_from_file(0)-latStepStart/2 ); 
    505  
    506  
    507        double latStepEnd = (latvalue_rectilinear_read_from_file(nj_glo-1)-latvalue_rectilinear_read_from_file(nj_glo-2)); 
    508        if (isSouthPole) bounds_lat_end.setValue(latvalue_rectilinear_read_from_file(nj_glo-1)); 
    509        else bounds_lat_end.setValue(latvalue_rectilinear_read_from_file(nj_glo-1)+latStepEnd/2); 
     500     double lonStepStart = lon(1)-lon(0); 
     501     bounds_lon_start=lon(0) - lonStepStart/2; 
     502     double lonStepEnd = lon(ni_glo-1)-lon(ni_glo-2); 
     503     bounds_lon_end=lon(ni_glo-1) + lonStepEnd/2; 
     504     double errorBoundsLon = std::abs(360-std::abs(bounds_lon_end-bounds_lon_start)); 
     505     if (errorBoundsLon > NumTraits<double>::epsilon()) bounds_lon_end=bounds_lon_start+360; 
     506     for(j=0;j<nj;++j) 
     507       for(i=0;i<ni;++i) 
     508       { 
     509         k=j*ni+i; 
     510         boundsLon(0,k) = boundsLon(1,k) = (0 == (ibegin + i)) ? bounds_lon_start 
     511                                                               : (lon(ibegin + i)+lon(ibegin + i-1))/2; 
     512         boundsLon(2,k) = boundsLon(3,k) = ((ibegin + i + 1) == ni_glo) ? bounds_lon_end 
     513                                                                        : (lon(ibegin + i + 1)+lon(ibegin + i))/2; 
     514       } 
     515     
     516 
     517    boundsLat.resize(nvertexValue,nj*ni); 
     518    bool isNorthPole, isSouthPole ; 
     519    if (std::abs(90 - std::abs(lat(0))) < NumTraits<double>::epsilon()) isNorthPole = true; 
     520    if (std::abs(90 - std::abs(lat(nj_glo-1))) < NumTraits<double>::epsilon()) isSouthPole = true; 
     521 
     522    double latStepStart = lat(1)-lat(0); 
     523    if (isNorthPole) bounds_lat_start=lat(0); 
     524    else bounds_lat_start=lat(0)-latStepStart/2; 
     525 
     526    double latStepEnd = lat(nj_glo-1)-lat(nj_glo-2); 
     527    if (isSouthPole) bounds_lat_end=lat(nj_glo-1); 
     528    else bounds_lat_end=lat(nj_glo-1)+latStepEnd/2; 
    510529        
    511        if (bounds_lat_start > 90.-1e-3) bounds_lat_start=90 ; 
    512        if (bounds_lat_start < -90.+1e-3) bounds_lat_start=-90 ; 
    513        if (bounds_lat_end > 90.-1e-3) bounds_lat_end=90 ; 
    514        if (bounds_lat_end < -90.+1e-3) bounds_lat_end=-90 ; 
     530// Work arround for small value close to pole, not too good for remapping     
     531    if (bounds_lat_start > 90.-1e-3) bounds_lat_start=90 ; 
     532    if (bounds_lat_start < -90.+1e-3) bounds_lat_start=-90 ; 
     533    if (bounds_lat_end > 90.-1e-3) bounds_lat_end=90 ; 
     534    if (bounds_lat_end < -90.+1e-3) bounds_lat_end=-90 ; 
    515535       
    516        for(j=0;j<nj;++j) 
    517          for(i=0;i<ni;++i) 
    518          { 
    519            k=j*ni+i; 
    520            boundsLat(1,k) = boundsLat(2,k) = (0 == (jbegin + j)) ? bounds_lat_start 
    521                                                                  : (latvalue_rectilinear_read_from_file(jbegin + j)+latvalue_rectilinear_read_from_file(jbegin + j-1))/2; 
    522            boundsLat(0,k) = boundsLat(3,k) = ((jbegin + j +1) == nj_glo) ? bounds_lat_end 
    523                                                                  : (latvalue_rectilinear_read_from_file(jbegin + j + 1)+latvalue_rectilinear_read_from_file(jbegin + j))/2; 
    524          } 
    525      } 
    526      else 
    527      { 
    528        double boundsLatRange = bounds_lat_end - bounds_lat_start; 
    529        double latStep = boundsLatRange/double(nj_glo.getValue()); 
    530        double bounds_lat_start_pole = bounds_lat_start; 
    531        double bounds_lat_end_pole   = bounds_lat_end; 
    532        if (isNorthPole) bounds_lat_start_pole = lat_start; 
    533        if (isSouthPole) bounds_lat_end_pole   = lat_end; 
    534  
    535        for(j=0;j<nj;++j) 
    536          for(i=0;i<ni;++i) 
    537          { 
    538            k=j*ni+i; 
    539            boundsLat(1,k) = boundsLat(2,k) = (0 != (jbegin + j)) ? (jbegin + j) * latStep + bounds_lat_start 
    540                                                                  : bounds_lat_start_pole; 
    541            boundsLat(0,k) = boundsLat(3,k) = ((jbegin + j +1) != nj_glo) ? (jbegin + j +1) * latStep + bounds_lat_start 
    542                                                                  : bounds_lat_end_pole; 
    543          } 
    544      } 
    545  
     536    for(j=0;j<nj;++j) 
     537      for(i=0;i<ni;++i) 
     538      { 
     539        k=j*ni+i; 
     540        boundsLat(1,k) = boundsLat(2,k) = (0 == (jbegin + j)) ? bounds_lat_start 
     541                                                              : (lat(jbegin + j)+lat(jbegin + j-1))/2; 
     542        boundsLat(0,k) = boundsLat(3,k) = ((jbegin + j +1) == nj_glo) ? bounds_lat_end 
     543                                                                      : (lat(jbegin + j + 1)+lat(jbegin + j))/2; 
     544      } 
    546545   } 
    547546 

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

    r795 r809  
    124124         void sendLonLatArea(void); 
    125125         void computeConnectedServer(void) ; 
    126          void fillInRectilinearBoundLonLat(CArray<double,2>& boundsLon, CArray<double,2>& boundsLat, 
    127                                            bool isNorthPole = false, bool isSouthPole = false); 
     126          
     127         void AllgatherRectilinearLonLat(CArray<double,1>& lon, CArray<double,1>& lat,  
     128                                         CArray<double,1>& lon_g, CArray<double,1>& lat_g); 
     129          
     130         void fillInRectilinearBoundLonLat(CArray<double,1>& lon, CArray<double,1>& lat,  
     131                                           CArray<double,2>& boundsLon, CArray<double,2>& boundsLat); 
    128132         void fillInRectilinearLonLat(); 
    129133 

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

    r808 r809  
    7878    bool isNorthPole = false; 
    7979    bool isSouthPole = false; 
    80  
    81     if (domainSrc_->latvalue_rectilinear_read_from_file.isEmpty()) 
    82     { 
    83       if (std::abs(poleValue - std::abs(domainSrc_->lat_start)) < NumTraits<double>::epsilon()) isNorthPole = true; 
    84       if (std::abs(poleValue - std::abs(domainSrc_->lat_end)) < NumTraits<double>::epsilon()) isSouthPole = true; 
    85     } 
    86     else 
    87     { 
    88       if (std::abs(poleValue - std::abs(domainSrc_->latvalue_rectilinear_read_from_file(0))) < NumTraits<double>::epsilon()) isNorthPole = true; 
    89       if (std::abs(poleValue - std::abs(domainSrc_->latvalue_rectilinear_read_from_file(domainSrc_->nj_glo-1))) < NumTraits<double>::epsilon()) isSouthPole = true; 
    90     } 
     80    CArray<double,1> lon_g ; 
     81    CArray<double,1> lat_g ; 
     82 
     83    if (!domainSrc_->lonvalue_1d.isEmpty() && !domainSrc_->latvalue_1d.isEmpty()) 
     84    { 
     85                domainSrc_->AllgatherRectilinearLonLat(domainSrc_->lonvalue_1d,domainSrc_->latvalue_1d, lon_g,lat_g) ; 
     86        } 
     87        else if (! domainSrc_->latvalue_rectilinear_read_from_file.isEmpty() && ! domainSrc_->lonvalue_rectilinear_read_from_file.isEmpty() ) 
     88    { 
     89                lat_g=domainSrc_->latvalue_rectilinear_read_from_file ; 
     90                lon_g=domainSrc_->lonvalue_rectilinear_read_from_file ; 
     91        } 
     92        else if (!domainSrc_->lon_start.isEmpty() && !domainSrc_->lon_end.isEmpty() && 
     93                 !domainSrc_->lat_start.isEmpty() && !domainSrc_->lat_end.isEmpty()) 
     94        { 
     95          double step=(domainSrc_->lon_end-domainSrc_->lon_start)/domainSrc_->ni_glo ; 
     96          for(int i=0; i<domainSrc_->ni_glo; ++i) lon_g(i)=domainSrc_->lon_start+i*step ; 
     97          step=(domainSrc_->lat_end-domainSrc_->lat_start)/domainSrc_->nj_glo ;  
     98          for(int i=0; i<domainSrc_->ni_glo; ++i) lat_g(i)=domainSrc_->lat_start+i*step ; 
     99        } 
     100        else ERROR("void CDomainAlgorithmInterpolate::computeRemap()",<<"Cannot compute bounds for rectilinear domain") ; 
    91101 
    92102    nVertexSrc = constNVertex; 
    93     domainSrc_->fillInRectilinearBoundLonLat(boundsLonSrc, boundsLatSrc, isNorthPole, isSouthPole); 
     103    domainSrc_->fillInRectilinearBoundLonLat(lon_g,lat_g, boundsLonSrc, boundsLatSrc); 
    94104  } 
    95105 
    … …  
    132142    if (std::abs(poleValue - std::abs(domainDest_->lat_start)) < NumTraits<double>::epsilon()) isNorthPole = true; 
    133143    if (std::abs(poleValue - std::abs(domainDest_->lat_end)) < NumTraits<double>::epsilon()) isSouthPole = true; 
     144 
     145    CArray<double,1> lon_g ; 
     146    CArray<double,1> lat_g ; 
     147 
     148    if (!domainDest_->lonvalue_1d.isEmpty() && !domainDest_->latvalue_1d.isEmpty()) 
     149    { 
     150                domainDest_->AllgatherRectilinearLonLat(domainDest_->lonvalue_1d,domainDest_->latvalue_1d, lon_g,lat_g) ; 
     151        } 
     152        else if (! domainDest_->latvalue_rectilinear_read_from_file.isEmpty() && ! domainDest_->lonvalue_rectilinear_read_from_file.isEmpty() ) 
     153    { 
     154                lat_g=domainDest_->latvalue_rectilinear_read_from_file ; 
     155                lon_g=domainDest_->lonvalue_rectilinear_read_from_file ; 
     156        } 
     157        else if (!domainDest_->lon_start.isEmpty() && !domainDest_->lon_end.isEmpty() && 
     158                 !domainDest_->lat_start.isEmpty() && !domainDest_->lat_end.isEmpty()) 
     159        { 
     160          double step=(domainDest_->lon_end-domainDest_->lon_start)/domainDest_->ni_glo ; 
     161          for(int i=0; i<domainDest_->ni_glo; ++i) lon_g(i)=domainDest_->lon_start+i*step ; 
     162          step=(domainDest_->lat_end-domainDest_->lat_start)/domainDest_->nj_glo ;  
     163          for(int i=0; i<domainDest_->ni_glo; ++i) lat_g(i)=domainDest_->lat_start+i*step ; 
     164        } 
     165        else ERROR("void CDomainAlgorithmInterpolate::computeRemap()",<<"Cannot compute bounds for rectilinear domain") ; 
     166    if (std::abs(poleValue - std::abs(lat_g(0))) < NumTraits<double>::epsilon()) isNorthPole = true; 
     167    if (std::abs(poleValue - std::abs(lat_g(domainDest_->nj_glo-1))) < NumTraits<double>::epsilon()) isSouthPole = true; 
     168 
     169     
     170     
     171     
    134172    if (isNorthPole && (0 == domainDest_->jbegin.getValue())) 
    135173    { 
    … …  
    155193    // Ok, fill in boundary values for rectangular domain 
    156194    nVertexDest = constNVertex; 
    157     domainDest_->fillInRectilinearBoundLonLat(boundsLonDest, boundsLatDest, isNorthPole, isSouthPole); 
     195    domainDest_->fillInRectilinearBoundLonLat(lon_g,lat_g, boundsLonDest, boundsLatDest); 
    158196  } 
    159197 
Note: See TracChangeset for help on using the changeset viewer.