Changeset 1114


Ignore:
Timestamp:
05/03/17 15:07:53 (4 years ago)
Author:
ymipsl
Message:

Add conservative remapping option for quantity, ie, the integrated value over a cell area.

YM

Location:
XIOS/trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • XIOS/trunk/extern/remap/src/mapper.cpp

    r923 r1114  
    109109} 
    110110 
    111 vector<double> Mapper::computeWeights(int interpOrder, bool renormalize) 
     111vector<double> Mapper::computeWeights(int interpOrder, bool renormalize, bool quantity) 
    112112{ 
    113113        vector<double> timings; 
     
    152152        if (mpiRank == 0 && verbose) cout << "Remapping..." << endl; 
    153153        tic = cputime(); 
    154         nWeights = remap(&targetElements[0], targetElements.size(), interpOrder, renormalize); 
     154        nWeights = remap(&targetElements[0], targetElements.size(), interpOrder, renormalize, quantity); 
    155155        timings.push_back(cputime() - tic); 
    156156 
     
    166166   @param order is the order of interpolaton (must be 1 or 2). 
    167167*/ 
    168 int Mapper::remap(Elt *elements, int nbElements, int order, bool renormalize) 
     168int Mapper::remap(Elt *elements, int nbElements, int order, bool renormalize, bool quantity) 
    169169{ 
    170170        int mpiSize, mpiRank; 
     
    184184        int **sendElement = new int*[mpiSize]; /* indices of elements required from other rank */ 
    185185        double **recvValue = new double*[mpiSize]; 
     186        double **recvArea = new double*[mpiSize]; 
    186187        Coord **recvGrad = new Coord*[mpiSize]; 
    187188        GloId **recvNeighIds = new GloId*[mpiSize]; /* ids of the of the source neighbours which also contribute through gradient */ 
     
    205206                        sendElement[rank] = new int[nbSendElement[rank]]; 
    206207                        recvValue[rank]   = new double[nbSendElement[rank]]; 
     208                        recvArea[rank]    = new double[nbSendElement[rank]]; 
    207209                        if (order == 2) 
    208210                        { 
     
    234236        int **recvElement = new int*[mpiSize]; 
    235237        double **sendValue = new double*[mpiSize]; 
     238        double **sendArea = new double*[mpiSize]; 
    236239        Coord **sendGrad = new Coord*[mpiSize]; 
    237240        GloId **sendNeighIds = new GloId*[mpiSize]; 
     
    250253                        recvElement[rank] = new int[nbRecvElement[rank]]; 
    251254                        sendValue[rank]   = new double[nbRecvElement[rank]]; 
     255                        sendArea[rank]   = new double[nbRecvElement[rank]]; 
    252256                        if (order == 2) 
    253257                        { 
     
    278282                        { 
    279283                                sendValue[rank][j] = sstree.localElements[recvElement[rank][j]].val; 
     284                                sendArea[rank][j] = sstree.localElements[recvElement[rank][j]].area; 
    280285                                if (order == 2) 
    281286                                { 
     
    297302                        MPI_Issend(sendValue[rank],  nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 
    298303                        nbSendRequest++; 
     304                        MPI_Issend(sendArea[rank],  nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 
     305                        nbSendRequest++; 
    299306                        if (order == 2) 
    300307                        { 
     
    317324                        MPI_Irecv(recvValue[rank],  nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
    318325                        nbRecvRequest++; 
     326                        MPI_Irecv(recvArea[rank],  nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
     327                        nbRecvRequest++; 
    319328                        if (order == 2) 
    320329                        { 
     
    357366                        int rank = (*it)->id.rank; 
    358367                        double fk = recvValue[rank][n1]; 
     368                        double srcArea = recvArea[rank][n1]; 
    359369                        double w = (*it)->area; 
     370      if (quantity) w/=srcArea ; 
    360371 
    361372                        /* first order: src value times weight (weight = supermesh area), later divide by target area */ 
    362373                        int kk = (order == 2) ? n1 * (NMAX + 1) : n1; 
    363374                        GloId neighID = recvNeighIds[rank][kk]; 
    364                         wgt_map[neighID] += (*it)->area; 
     375                        wgt_map[neighID] += w; 
    365376 
    366377                        if (order == 2) 
     
    383394    for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) 
    384395                { 
    385                         this->remapMatrix[i] = (it->second / e.area) / renorm; 
     396      if (quantity)  this->remapMatrix[i] = (it->second ) / renorm; 
     397                        else this->remapMatrix[i] = (it->second / e.area) / renorm; 
    386398                        this->srcAddress[i] = it->first.ind; 
    387399                        this->srcRank[i] = it->first.rank; 
     
    400412                        delete[] sendElement[rank]; 
    401413                        delete[] recvValue[rank]; 
     414                        delete[] recvArea[rank]; 
    402415                        if (order == 2) 
    403416                        { 
     
    410423                        delete[] recvElement[rank]; 
    411424                        delete[] sendValue[rank]; 
     425                        delete[] sendArea[rank]; 
    412426                        if (order == 2) 
    413427                                delete[] sendGrad[rank]; 
  • XIOS/trunk/extern/remap/src/mapper.hpp

    r844 r1114  
    3434       /** @param trgElts are the elements of the unstructured target grid 
    3535           Returns the timings for substeps: */ 
    36        vector<double> computeWeights(int interpOrder, bool renormalize=false); 
     36       vector<double> computeWeights(int interpOrder, bool renormalize=false, bool quantity=false); 
    3737       int getNbWeights(void) { return nWeights ; } 
    3838/* 
     
    5151private: 
    5252       /** @return number of weights (local to cpu) */ 
    53        int remap(Elt* elements, int nbElements, int order, bool renormalize=false); 
     53       int remap(Elt* elements, int nbElements, int order, bool renormalize=false, bool quantity=false); 
    5454 
    5555       void buildMeshTopology(); 
  • XIOS/trunk/src/config/interpolate_domain_attribute.conf

    r1014 r1114  
    22DECLARE_ATTRIBUTE(int, order) 
    33DECLARE_ATTRIBUTE(bool, renormalize) 
     4DECLARE_ATTRIBUTE(bool, quantity) 
    45 
    56/* Write interpolation weights into file */ 
  • XIOS/trunk/src/transformation/domain_algorithm_interpolate.cpp

    r1046 r1114  
    100100  int orderInterp = interpDomain_->order.getValue(); 
    101101  bool renormalize ; 
     102  bool quantity ; 
    102103 
    103104  if (interpDomain_->renormalize.isEmpty()) renormalize=true; 
    104105  else renormalize = interpDomain_->renormalize; 
     106 
     107  if (interpDomain_->quantity.isEmpty()) quantity=false; 
     108  else quantity = interpDomain_->quantity; 
    105109 
    106110  const double poleValue = 90.0; 
     
    342346  mapper.setTargetMesh(boundsLonDestUnmasked.dataFirst(), boundsLatDestUnmasked.dataFirst(), nVertexDest, nDstLocalUnmasked, &dstPole[0], globalDstUnmasked); 
    343347 
    344   std::vector<double> timings = mapper.computeWeights(orderInterp,renormalize); 
     348  std::vector<double> timings = mapper.computeWeights(orderInterp,renormalize,quantity); 
    345349 
    346350  std::map<int,std::vector<std::pair<int,double> > > interpMapValue; 
Note: See TracChangeset for help on using the changeset viewer.