Changeset 1615


Ignore:
Timestamp:
11/26/18 10:35:26 (9 months ago)
Author:
ymipsl
Message:

Interpolation enhancement :
Domain area can be used to improved global conservation when there is a discrepency between model area and xios computed area.

New domain attribute :

radius : radius af the spherical domain, used to compute area ond the sphere with a normalized radius of 1 (for remapper).

New domain_interpolate attribute :
use_area : remapper will take model computed area in order to perform a global conservation for flux integrated over the cell (mass for example).
In this case the domain attributes "area" and "radius" must be defined on the source or target domain to be taking into account.

YM

Location:
XIOS/trunk/src
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • XIOS/trunk/src/config/domain_attribute.conf

    r1553 r1615  
    5858 
    5959DECLARE_ARRAY(double, 2, area) 
     60DECLARE_ATTRIBUTE(double, radius) 
    6061 
    6162DECLARE_ENUM4(type,rectilinear,curvilinear,unstructured, gaussian) 
  • XIOS/trunk/src/config/interpolate_domain_attribute.conf

    r1269 r1615  
    1010DECLARE_ATTRIBUTE(bool, write_weight) 
    1111DECLARE_ENUM2(read_write_convention, c, fortran) 
     12DECLARE_ATTRIBUTE(bool, use_area) 
  • XIOS/trunk/src/context_client.cpp

    r1475 r1615  
    9595    { 
    9696      list<int> ranks = event.getRanks(); 
    97  
     97      info(100)<<"Event "<<timeLine<<" of context "<<context->getId()<<endl ; 
    9898      if (CXios::checkEventSync) 
    9999      { 
     
    124124        { 
    125125          event.send(timeLine, sizes, buffList); 
     126          info(100)<<"Event "<<timeLine<<" of context "<<context->getId()<<"  sent"<<endl ; 
    126127 
    127128          checkBuffers(ranks); 
     
    142143          info(100)<<"DEBUG : temporaly event created : timeline "<<timeLine<<endl ; 
    143144          event.send(timeLine, tmpBufferedEvent.sizes, tmpBufferedEvent.buffers); 
     145          info(100)<<"Event "<<timeLine<<" of context "<<context->getId()<<"  sent"<<endl ; 
    144146        } 
    145147      } 
  • XIOS/trunk/src/node/interpolate_domain.cpp

    r1269 r1615  
    6565    if (this->read_write_convention.isEmpty()) this->read_write_convention.setValue(read_write_convention_attr::fortran); 
    6666 
     67 
    6768  } 
    6869 
  • XIOS/trunk/src/transformation/domain_algorithm_interpolate.cpp

    r1542 r1615  
    304304  CArray<double,2> boundsLonSrcUnmasked(nVertexSrc,nSrcLocalUnmasked); 
    305305  CArray<double,2> boundsLatSrcUnmasked(nVertexSrc,nSrcLocalUnmasked); 
     306  CArray<double,1> areaSrcUnmasked(nSrcLocalUnmasked); 
     307   
    306308  long int * globalSrcUnmasked = new long int [nSrcLocalUnmasked]; 
    307309 
    308310  nSrcLocalUnmasked=0 ; 
     311  bool hasSrcArea=!domainSrc_->area.isEmpty() && !domainSrc_->radius.isEmpty() && !interpDomain_->use_area.isEmpty() && interpDomain_->use_area==true  ; 
     312  double srcAreaFactor ; 
     313  if (hasSrcArea) srcAreaFactor=1./(domainSrc_->radius*domainSrc_->radius) ; 
     314   
    309315  for (int idx=0 ; idx < nSrcLocal; idx++) 
    310316  { 
     
    316322        boundsLatSrcUnmasked(n,nSrcLocalUnmasked) = boundsLatSrc(n,idx) ; 
    317323      } 
     324      if (hasSrcArea) areaSrcUnmasked(nSrcLocalUnmasked) = domainSrc_->area(idx)*srcAreaFactor ; 
    318325      globalSrcUnmasked[nSrcLocalUnmasked]=globalSrc[idx] ; 
    319326      ++nSrcLocalUnmasked ; 
    320327    } 
    321328  } 
    322  
     329  
    323330 
    324331  int nDstLocalUnmasked = 0 ; 
     
    327334  CArray<double,2> boundsLonDestUnmasked(nVertexDest,nDstLocalUnmasked); 
    328335  CArray<double,2> boundsLatDestUnmasked(nVertexDest,nDstLocalUnmasked); 
     336  CArray<double,1>   areaDstUnmasked(nSrcLocalUnmasked); 
     337 
    329338  long int * globalDstUnmasked = new long int [nDstLocalUnmasked]; 
    330339 
    331340  nDstLocalUnmasked=0 ; 
     341  bool hasDstArea=!domainDest_->area.isEmpty() && !domainDest_->radius.isEmpty() && !interpDomain_->use_area.isEmpty() && interpDomain_->use_area==true ; 
     342  double dstAreaFactor ; 
     343  if (hasDstArea) dstAreaFactor=1./(domainDest_->radius*domainDest_->radius) ; 
    332344  for (int idx=0 ; idx < nDstLocal; idx++) 
    333345  { 
     
    339351        boundsLatDestUnmasked(n,nDstLocalUnmasked) = boundsLatDest(n,idx) ; 
    340352      } 
     353      if (hasDstArea) areaDstUnmasked(nDstLocalUnmasked) = domainDest_->area(idx)*dstAreaFactor ; 
    341354      globalDstUnmasked[nDstLocalUnmasked]=globalDst[idx] ; 
    342355      ++nDstLocalUnmasked ; 
     
    344357  } 
    345358 
    346   mapper.setSourceMesh(boundsLonSrcUnmasked.dataFirst(), boundsLatSrcUnmasked.dataFirst(), nVertexSrc, nSrcLocalUnmasked, &srcPole[0], globalSrcUnmasked); 
    347   mapper.setTargetMesh(boundsLonDestUnmasked.dataFirst(), boundsLatDestUnmasked.dataFirst(), nVertexDest, nDstLocalUnmasked, &dstPole[0], globalDstUnmasked); 
     359  double* ptAreaSrcUnmasked = NULL ; 
     360  if (hasSrcArea) ptAreaSrcUnmasked=areaSrcUnmasked.dataFirst() ; 
     361 
     362  double* ptAreaDstUnmasked = NULL ; 
     363  if (hasDstArea) ptAreaDstUnmasked=areaDstUnmasked.dataFirst() ; 
     364 
     365  mapper.setSourceMesh(boundsLonSrcUnmasked.dataFirst(), boundsLatSrcUnmasked.dataFirst(), ptAreaSrcUnmasked, nVertexSrc, nSrcLocalUnmasked, &srcPole[0], globalSrcUnmasked); 
     366  mapper.setTargetMesh(boundsLonDestUnmasked.dataFirst(), boundsLatDestUnmasked.dataFirst(), ptAreaDstUnmasked, nVertexDest, nDstLocalUnmasked, &dstPole[0], globalDstUnmasked); 
    348367 
    349368  std::vector<double> timings = mapper.computeWeights(orderInterp,renormalize,quantity); 
     
    949968    size_t nbIndex=dataOut.numElements() ;  
    950969 
    951     for (int idx = 0; idx < nbIndex; ++idx) 
    952     { 
    953       if (allMissing(idx)) dataOut(idx) = defaultValue; // If all data source are nan then data destination must be nan 
     970    if (allMissing.numElements()>0) 
     971    { 
     972      for (int idx = 0; idx < nbIndex; ++idx) 
     973      { 
     974        if (allMissing(idx)) dataOut(idx) = defaultValue; // If all data source are nan then data destination must be nan 
     975      } 
    954976    } 
    955977     
Note: See TracChangeset for help on using the changeset viewer.