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

File:
1 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]; 
Note: See TracChangeset for help on using the changeset viewer.