Ignore:
Timestamp:
03/22/18 10:43:20 (6 years ago)
Author:
yushan
Message:

branch_openmp merged with XIOS_DEV_CMIP6@1459

Location:
XIOS/dev/branch_openmp/extern/remap/src
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • XIOS/dev/branch_openmp/extern/remap/src/gridRemap.hpp

    r1335 r1460  
    1414Coord readPole(std::istream&); 
    1515 
    16  
     16//extern CRemapGrid srcGrid; 
     17//extern CRemapGrid tgtGrid; 
    1718 
    1819} 

  • XIOS/dev/branch_openmp/extern/remap/src/mapper.cpp

    r1341 r1460  
    2929void cptOffsetsFromLengths(const int *lengths, int *offsets, int sz) 
    3030{ 
    31   offsets[0] = 0; 
    32   for (int i = 1; i < sz; i++) 
    33     offsets[i] = offsets[i-1] + lengths[i-1]; 
     31        offsets[0] = 0; 
     32        for (int i = 1; i < sz; i++) 
     33                offsets[i] = offsets[i-1] + lengths[i-1]; 
    3434} 
    3535 
    … …  
    3939  srcGrid.pole = Coord(pole[0], pole[1], pole[2]); 
    4040 
    41   int mpiRank, mpiSize; 
    42   MPI_Comm_rank(communicator, &mpiRank); 
    43   MPI_Comm_size(communicator, &mpiSize); 
    44  
    45   sourceElements.reserve(nbCells); 
    46   sourceMesh.reserve(nbCells); 
     41        int mpiRank, mpiSize; 
     42        MPI_Comm_rank(communicator, &mpiRank); 
     43        MPI_Comm_size(communicator, &mpiSize); 
     44 
     45        sourceElements.reserve(nbCells); 
     46        sourceMesh.reserve(nbCells); 
    4747  sourceGlobalId.resize(nbCells) ; 
    4848 
    … …  
    5757  else sourceGlobalId.assign(globalId,globalId+nbCells); 
    5858 
    59   for (int i = 0; i < nbCells; i++) 
    60   { 
    61     int offs = i*nVertex; 
    62     Elt elt(boundsLon + offs, boundsLat + offs, nVertex); 
    63     elt.src_id.rank = mpiRank; 
    64     elt.src_id.ind = i; 
     59        for (int i = 0; i < nbCells; i++) 
     60        { 
     61                int offs = i*nVertex; 
     62                Elt elt(boundsLon + offs, boundsLat + offs, nVertex); 
     63                elt.src_id.rank = mpiRank; 
     64                elt.src_id.ind = i; 
    6565    elt.src_id.globalId = sourceGlobalId[i]; 
    66     sourceElements.push_back(elt); 
    67     sourceMesh.push_back(Node(elt.x, cptRadius(elt), &sourceElements.back())); 
    68     cptEltGeom(sourceElements[i], Coord(pole[0], pole[1], pole[2])); 
    69   } 
     66                sourceElements.push_back(elt); 
     67                sourceMesh.push_back(Node(elt.x, cptRadius(elt), &sourceElements.back())); 
     68                cptEltGeom(sourceElements[i], Coord(pole[0], pole[1], pole[2])); 
     69        } 
    7070 
    7171} 
    … …  
    7575  tgtGrid.pole = Coord(pole[0], pole[1], pole[2]); 
    7676 
    77   int mpiRank, mpiSize; 
    78   MPI_Comm_rank(communicator, &mpiRank); 
    79   MPI_Comm_size(communicator, &mpiSize); 
    80  
    81   targetElements.reserve(nbCells); 
    82   targetMesh.reserve(nbCells); 
     77        int mpiRank, mpiSize; 
     78        MPI_Comm_rank(communicator, &mpiRank); 
     79        MPI_Comm_size(communicator, &mpiSize); 
     80 
     81        targetElements.reserve(nbCells); 
     82        targetMesh.reserve(nbCells); 
    8383 
    8484  targetGlobalId.resize(nbCells) ; 
    … …  
    9393  else targetGlobalId.assign(globalId,globalId+nbCells); 
    9494 
    95   for (int i = 0; i < nbCells; i++) 
    96   { 
    97     int offs = i*nVertex; 
    98     Elt elt(boundsLon + offs, boundsLat + offs, nVertex); 
    99     targetElements.push_back(elt); 
    100     targetMesh.push_back(Node(elt.x, cptRadius(elt), &sourceElements.back())); 
    101     cptEltGeom(targetElements[i], Coord(pole[0], pole[1], pole[2])); 
    102   } 
     95        for (int i = 0; i < nbCells; i++) 
     96        { 
     97                int offs = i*nVertex; 
     98                Elt elt(boundsLon + offs, boundsLat + offs, nVertex); 
     99                targetElements.push_back(elt); 
     100                targetMesh.push_back(Node(elt.x, cptRadius(elt), &sourceElements.back())); 
     101                cptEltGeom(targetElements[i], Coord(pole[0], pole[1], pole[2])); 
     102        } 
    103103 
    104104 
    … …  
    119119vector<double> Mapper::computeWeights(int interpOrder, bool renormalize, bool quantity) 
    120120{ 
    121   vector<double> timings; 
    122   int mpiSize, mpiRank; 
    123   MPI_Comm_size(communicator, &mpiSize); 
    124   MPI_Comm_rank(communicator, &mpiRank); 
    125  
     121        vector<double> timings; 
     122        int mpiSize, mpiRank; 
     123        MPI_Comm_size(communicator, &mpiSize); 
     124        MPI_Comm_rank(communicator, &mpiRank); 
    126125 
    127126  this->buildSSTree(sourceMesh, targetMesh); 
    128127 
    129   if (mpiRank == 0 && verbose) cout << "Computing intersections ..." << endl; 
    130   double tic = cputime(); 
    131   computeIntersection(&targetElements[0], targetElements.size()); 
    132   timings.push_back(cputime() - tic); 
    133  
    134   tic = cputime(); 
    135   if (interpOrder == 2)  
    136   { 
    137     if (mpiRank == 0 && verbose) cout << "Computing grads ..." << endl; 
    138     buildMeshTopology(); 
    139     computeGrads(); 
    140   } 
    141   timings.push_back(cputime() - tic); 
    142  
    143   /* Prepare computation of weights */ 
    144   /* compute number of intersections which for the first order case 
    145      corresponds to the number of edges in the remap matrix */ 
    146   int nIntersections = 0; 
    147   for (int j = 0; j < targetElements.size(); j++) 
    148   { 
    149     Elt &elt = targetElements[j]; 
    150     for (list<Polyg*>::iterator it = elt.is.begin(); it != elt.is.end(); it++) 
    151       nIntersections++; 
    152   } 
    153   /* overallocate for NMAX neighbours for each elements */ 
    154   remapMatrix = new double[nIntersections*NMAX]; 
    155   srcAddress = new int[nIntersections*NMAX]; 
    156   srcRank = new int[nIntersections*NMAX]; 
    157   dstAddress = new int[nIntersections*NMAX]; 
     128        if (mpiRank == 0 && verbose) cout << "Computing intersections ..." << endl; 
     129        double tic = cputime(); 
     130        computeIntersection(&targetElements[0], targetElements.size()); 
     131        timings.push_back(cputime() - tic); 
     132 
     133        tic = cputime(); 
     134        if (interpOrder == 2) { 
     135                if (mpiRank == 0 && verbose) cout << "Computing grads ..." << endl; 
     136                buildMeshTopology(); 
     137                computeGrads(); 
     138        } 
     139        timings.push_back(cputime() - tic); 
     140 
     141        /* Prepare computation of weights */ 
     142        /* compute number of intersections which for the first order case 
     143           corresponds to the number of edges in the remap matrix */ 
     144        int nIntersections = 0; 
     145        for (int j = 0; j < targetElements.size(); j++) 
     146        { 
     147                Elt &elt = targetElements[j]; 
     148                for (list<Polyg*>::iterator it = elt.is.begin(); it != elt.is.end(); it++) 
     149                        nIntersections++; 
     150        } 
     151        /* overallocate for NMAX neighbours for each elements */ 
     152        remapMatrix = new double[nIntersections*NMAX]; 
     153        srcAddress = new int[nIntersections*NMAX]; 
     154        srcRank = new int[nIntersections*NMAX]; 
     155        dstAddress = new int[nIntersections*NMAX]; 
    158156  sourceWeightId =new long[nIntersections*NMAX]; 
    159157  targetWeightId =new long[nIntersections*NMAX]; 
    160158 
    161159 
    162   if (mpiRank == 0 && verbose) cout << "Remapping..." << endl; 
    163   tic = cputime(); 
    164   nWeights = remap(&targetElements[0], targetElements.size(), interpOrder, renormalize, quantity); 
    165   timings.push_back(cputime() - tic); 
     160        if (mpiRank == 0 && verbose) cout << "Remapping..." << endl; 
     161        tic = cputime(); 
     162        nWeights = remap(&targetElements[0], targetElements.size(), interpOrder, renormalize, quantity); 
     163        timings.push_back(cputime() - tic); 
    166164 
    167165  for (int i = 0; i < targetElements.size(); i++) targetElements[i].delete_intersections(); 
    168166 
    169   return timings; 
     167        return timings; 
    170168} 
    171169 
    … …  
    178176int Mapper::remap(Elt *elements, int nbElements, int order, bool renormalize, bool quantity) 
    179177{ 
    180   int mpiSize, mpiRank; 
    181   MPI_Comm_size(communicator, &mpiSize); 
    182   MPI_Comm_rank(communicator, &mpiRank); 
    183  
    184   /* create list of intersections (super mesh elements) for each rank */ 
    185   multimap<int, Polyg *> *elementList = new multimap<int, Polyg *>[mpiSize]; 
    186   for (int j = 0; j < nbElements; j++) 
    187   { 
    188     Elt& e = elements[j]; 
    189     for (list<Polyg *>::iterator it = e.is.begin(); it != e.is.end(); it++) 
    190     { 
    191       elementList[(*it)->id.rank].insert(pair<int, Polyg *>((*it)->id.ind, *it)); 
    192       //std::cout<<"elementList["<<(*it)->id.rank<<"].size = "<< elementList[(*it)->id.rank].size()<<std::endl; 
    193     } 
    194   } 
    195  
    196   int *nbSendElement = new int[mpiSize]; 
    197   int **sendElement = new int*[mpiSize]; /* indices of elements required from other rank */ 
    198   double **recvValue = new double*[mpiSize]; 
    199   double **recvArea = new double*[mpiSize]; 
    200   Coord **recvGrad = new Coord*[mpiSize]; 
    201   GloId **recvNeighIds = new GloId*[mpiSize]; /* ids of the of the source neighbours which also contribute through gradient */ 
    202   for (int rank = 0; rank < mpiSize; rank++) 
    203   { 
    204     /* get size for allocation */ 
    205     int last = -1; /* compares unequal to any index */ 
    206     int index = -1; /* increased to starting index 0 in first iteration */ 
    207     for (multimap<int, Polyg *>::iterator it = elementList[rank].begin(); it != elementList[rank].end(); ++it) 
    208     { 
    209       if (last != it->first) 
    210         index++; 
    211       (it->second)->id.ind = index; 
    212       last = it->first; 
    213     } 
    214     nbSendElement[rank] = index + 1; 
    215  
    216     /* if size is non-zero allocate and collect indices of elements on other ranks that we intersect */ 
    217     if (nbSendElement[rank] > 0) 
    218     { 
    219       sendElement[rank] = new int[nbSendElement[rank]]; 
    220       recvValue[rank]   = new double[nbSendElement[rank]]; 
    221       recvArea[rank]    = new double[nbSendElement[rank]]; 
    222       if (order == 2) 
    223       { 
    224         recvNeighIds[rank] = new GloId[nbSendElement[rank]*(NMAX+1)]; 
    225         recvGrad[rank]    = new Coord[nbSendElement[rank]*(NMAX+1)]; 
    226       } 
    227       else 
    228         recvNeighIds[rank] = new GloId[nbSendElement[rank]]; 
    229  
    230       last = -1; 
    231       index = -1; 
    232       for (multimap<int, Polyg *>::iterator it = elementList[rank].begin(); it != elementList[rank].end(); ++it) 
    233       { 
    234         if (last != it->first) 
    235           index++; 
    236         sendElement[rank][index] = it->first; 
    237         last = it->first; 
    238       } 
    239     } 
    240   } 
    241  
    242   /* communicate sizes of source elements to be sent (index lists and later values and gradients) */ 
    243   int *nbRecvElement = new int[mpiSize]; 
    244   MPI_Alltoall(nbSendElement, 1, MPI_INT, nbRecvElement, 1, MPI_INT, communicator); 
    245    
    246  
    247   /* communicate indices of source elements on other ranks whoes value and gradient we need (since intersection) */ 
    248   int nbSendRequest = 0; 
    249   int nbRecvRequest = 0; 
    250   int **recvElement = new int*[mpiSize]; 
    251   double **sendValue = new double*[mpiSize]; 
    252   double **sendArea = new double*[mpiSize]; 
    253   Coord **sendGrad = new Coord*[mpiSize]; 
    254   GloId **sendNeighIds = new GloId*[mpiSize]; 
    255   MPI_Request *sendRequest = new MPI_Request[4*mpiSize]; 
    256   MPI_Request *recvRequest = new MPI_Request[4*mpiSize]; 
    257   for (int rank = 0; rank < mpiSize; rank++) 
    258   { 
    259     if (nbSendElement[rank] > 0) 
    260     { 
    261       MPI_Issend(sendElement[rank], nbSendElement[rank], MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 
    262       nbSendRequest++; 
    263     } 
    264  
    265     if (nbRecvElement[rank] > 0) 
    266     { 
    267       recvElement[rank] = new int[nbRecvElement[rank]]; 
    268       sendValue[rank]   = new double[nbRecvElement[rank]]; 
    269       sendArea[rank]   = new double[nbRecvElement[rank]]; 
    270       if (order == 2) 
    271       { 
    272         sendNeighIds[rank] = new GloId[nbRecvElement[rank]*(NMAX+1)]; 
    273         sendGrad[rank]    = new Coord[nbRecvElement[rank]*(NMAX+1)]; 
    274       } 
    275       else 
    276       { 
    277         sendNeighIds[rank] = new GloId[nbRecvElement[rank]]; 
    278       } 
    279       MPI_Irecv(recvElement[rank], nbRecvElement[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
    280       nbRecvRequest++; 
    281     } 
    282   } 
    283  
    284   MPI_Status *status = new MPI_Status[4*mpiSize]; 
    285  
    286   MPI_Waitall(nbSendRequest, sendRequest, status); 
    287   MPI_Waitall(nbRecvRequest, recvRequest, status); 
    288  
    289   /* for all indices that have been received from requesting ranks: pack values and gradients, then send */ 
    290   nbSendRequest = 0; 
    291   nbRecvRequest = 0; 
    292   for (int rank = 0; rank < mpiSize; rank++) 
    293   { 
    294     if (nbRecvElement[rank] > 0) 
    295     { 
    296       int jj = 0; // jj == j if no weight writing 
    297       for (int j = 0; j < nbRecvElement[rank]; j++) 
    298       { 
    299         sendValue[rank][j] = sstree.localElements[recvElement[rank][j]].val; 
    300         sendArea[rank][j] = sstree.localElements[recvElement[rank][j]].area; 
    301         if (order == 2) 
    302         { 
    303           sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].grad; 
    304           sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].src_id; 
    305           jj++; 
    306           for (int i = 0; i < NMAX; i++) 
    307           { 
    308             sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].gradNeigh[i]; 
     178        int mpiSize, mpiRank; 
     179        MPI_Comm_size(communicator, &mpiSize); 
     180        MPI_Comm_rank(communicator, &mpiRank); 
     181 
     182        /* create list of intersections (super mesh elements) for each rank */ 
     183        multimap<int, Polyg *> *elementList = new multimap<int, Polyg *>[mpiSize]; 
     184        for (int j = 0; j < nbElements; j++) 
     185        { 
     186                Elt& e = elements[j]; 
     187                for (list<Polyg *>::iterator it = e.is.begin(); it != e.is.end(); it++) 
     188                        elementList[(*it)->id.rank].insert(pair<int, Polyg *>((*it)->id.ind, *it)); 
     189        } 
     190 
     191        int *nbSendElement = new int[mpiSize]; 
     192        int **sendElement = new int*[mpiSize]; /* indices of elements required from other rank */ 
     193        double **recvValue = new double*[mpiSize]; 
     194        double **recvArea = new double*[mpiSize]; 
     195        Coord **recvGrad = new Coord*[mpiSize]; 
     196        GloId **recvNeighIds = new GloId*[mpiSize]; /* ids of the of the source neighbours which also contribute through gradient */ 
     197        for (int rank = 0; rank < mpiSize; rank++) 
     198        { 
     199                /* get size for allocation */ 
     200                int last = -1; /* compares unequal to any index */ 
     201                int index = -1; /* increased to starting index 0 in first iteration */ 
     202                for (multimap<int, Polyg *>::iterator it = elementList[rank].begin(); it != elementList[rank].end(); ++it) 
     203                { 
     204                        if (last != it->first) 
     205                                index++; 
     206                        (it->second)->id.ind = index; 
     207                        last = it->first; 
     208                } 
     209                nbSendElement[rank] = index + 1; 
     210 
     211                /* if size is non-zero allocate and collect indices of elements on other ranks that we intersect */ 
     212                if (nbSendElement[rank] > 0) 
     213                { 
     214                        sendElement[rank] = new int[nbSendElement[rank]]; 
     215                        recvValue[rank]   = new double[nbSendElement[rank]]; 
     216                        recvArea[rank]    = new double[nbSendElement[rank]]; 
     217                        if (order == 2) 
     218                        { 
     219                                recvNeighIds[rank] = new GloId[nbSendElement[rank]*(NMAX+1)]; 
     220                                recvGrad[rank]    = new Coord[nbSendElement[rank]*(NMAX+1)]; 
     221                        } 
     222                        else 
     223                                recvNeighIds[rank] = new GloId[nbSendElement[rank]]; 
     224 
     225                        last = -1; 
     226                        index = -1; 
     227                        for (multimap<int, Polyg *>::iterator it = elementList[rank].begin(); it != elementList[rank].end(); ++it) 
     228                        { 
     229                                if (last != it->first) 
     230                                        index++; 
     231                                sendElement[rank][index] = it->first; 
     232                                last = it->first; 
     233                        } 
     234                } 
     235        } 
     236 
     237        /* communicate sizes of source elements to be sent (index lists and later values and gradients) */ 
     238        int *nbRecvElement = new int[mpiSize]; 
     239        MPI_Alltoall(nbSendElement, 1, MPI_INT, nbRecvElement, 1, MPI_INT, communicator); 
     240 
     241        /* communicate indices of source elements on other ranks whoes value and gradient we need (since intersection) */ 
     242        int nbSendRequest = 0; 
     243        int nbRecvRequest = 0; 
     244        int **recvElement = new int*[mpiSize]; 
     245        double **sendValue = new double*[mpiSize]; 
     246        double **sendArea = new double*[mpiSize]; 
     247        Coord **sendGrad = new Coord*[mpiSize]; 
     248        GloId **sendNeighIds = new GloId*[mpiSize]; 
     249        MPI_Request *sendRequest = new MPI_Request[4*mpiSize]; 
     250        MPI_Request *recvRequest = new MPI_Request[4*mpiSize]; 
     251        for (int rank = 0; rank < mpiSize; rank++) 
     252        { 
     253                if (nbSendElement[rank] > 0) 
     254                { 
     255                        MPI_Issend(sendElement[rank], nbSendElement[rank], MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 
     256                        nbSendRequest++; 
     257                } 
     258 
     259                if (nbRecvElement[rank] > 0) 
     260                { 
     261                        recvElement[rank] = new int[nbRecvElement[rank]]; 
     262                        sendValue[rank]   = new double[nbRecvElement[rank]]; 
     263                        sendArea[rank]   = new double[nbRecvElement[rank]]; 
     264                        if (order == 2) 
     265                        { 
     266                                sendNeighIds[rank] = new GloId[nbRecvElement[rank]*(NMAX+1)]; 
     267                                sendGrad[rank]    = new Coord[nbRecvElement[rank]*(NMAX+1)]; 
     268                        } 
     269                        else 
     270                        { 
     271                                sendNeighIds[rank] = new GloId[nbRecvElement[rank]]; 
     272                        } 
     273                        MPI_Irecv(recvElement[rank], nbRecvElement[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
     274                        nbRecvRequest++; 
     275                } 
     276        } 
     277        MPI_Status *status = new MPI_Status[4*mpiSize]; 
     278         
     279        MPI_Waitall(nbSendRequest, sendRequest, status); 
     280        MPI_Waitall(nbRecvRequest, recvRequest, status); 
     281 
     282        /* for all indices that have been received from requesting ranks: pack values and gradients, then send */ 
     283        nbSendRequest = 0; 
     284        nbRecvRequest = 0; 
     285        for (int rank = 0; rank < mpiSize; rank++) 
     286        { 
     287                if (nbRecvElement[rank] > 0) 
     288                { 
     289                        int jj = 0; // jj == j if no weight writing 
     290                        for (int j = 0; j < nbRecvElement[rank]; j++) 
     291                        { 
     292                                sendValue[rank][j] = sstree.localElements[recvElement[rank][j]].val; 
     293                                sendArea[rank][j] = sstree.localElements[recvElement[rank][j]].area; 
     294                                if (order == 2) 
     295                                { 
     296                                        sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].grad; 
     297//          cout<<"grad  "<<jj<<"  "<<recvElement[rank][j]<<"  "<<sendGrad[rank][jj]<<" "<<sstree.localElements[recvElement[rank][j]].grad<<endl ; 
     298                                        sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].src_id; 
     299                                        jj++; 
     300                                        for (int i = 0; i < NMAX; i++) 
     301                                        { 
     302                                                sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].gradNeigh[i]; 
     303//            cout<<"grad  "<<jj<<"  "<<sendGrad[rank][jj]<<" "<<sstree.localElements[recvElement[rank][j]].grad<<endl ; 
    309304            sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].neighId[i]; 
    310             jj++; 
    311           } 
    312         } 
    313         else 
    314           sendNeighIds[rank][j] = sstree.localElements[recvElement[rank][j]].src_id; 
    315       } 
    316       MPI_Issend(sendValue[rank],  nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 
    317       nbSendRequest++; 
    318       MPI_Issend(sendArea[rank],  nbRecvElement[rank], MPI_DOUBLE, rank, 1, communicator, &sendRequest[nbSendRequest]); 
    319       nbSendRequest++; 
    320       if (order == 2) 
    321       { 
    322         MPI_Issend(sendGrad[rank], 3*nbRecvElement[rank]*(NMAX+1), MPI_DOUBLE, rank, 2, communicator, &sendRequest[nbSendRequest]); 
    323         nbSendRequest++; 
    324         MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank]*(NMAX+1), MPI_INT, rank, 3, communicator, &sendRequest[nbSendRequest]); 
     305                                                jj++; 
     306                                        } 
     307                                } 
     308                                else 
     309                                        sendNeighIds[rank][j] = sstree.localElements[recvElement[rank][j]].src_id; 
     310                        } 
     311                        MPI_Issend(sendValue[rank],  nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 
     312                        nbSendRequest++; 
     313                        MPI_Issend(sendArea[rank],  nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 
     314                        nbSendRequest++; 
     315                        if (order == 2) 
     316                        { 
     317                                MPI_Issend(sendGrad[rank], 3*nbRecvElement[rank]*(NMAX+1), 
     318                                                                MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 
     319                                nbSendRequest++; 
     320                                MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank]*(NMAX+1), MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 
    325321//ym  --> attention taille GloId 
    326         nbSendRequest++; 
    327       } 
    328       else 
    329       { 
    330         MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank], MPI_INT, rank, 4, communicator, &sendRequest[nbSendRequest]); 
     322                                nbSendRequest++; 
     323                        } 
     324                        else 
     325                        { 
     326                                MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank], MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 
    331327//ym  --> attention taille GloId 
    332         nbSendRequest++; 
    333       } 
    334     } 
    335     if (nbSendElement[rank] > 0) 
    336     { 
    337       MPI_Irecv(recvValue[rank],  nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
    338       nbRecvRequest++; 
    339       MPI_Irecv(recvArea[rank],  nbSendElement[rank], MPI_DOUBLE, rank, 1, communicator, &recvRequest[nbRecvRequest]); 
    340       nbRecvRequest++; 
    341       if (order == 2) 
    342       { 
    343         MPI_Irecv(recvGrad[rank], 3*nbSendElement[rank]*(NMAX+1), MPI_DOUBLE, rank, 2, communicator, &recvRequest[nbRecvRequest]); 
    344         nbRecvRequest++; 
    345         MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank]*(NMAX+1), MPI_INT, rank, 3, communicator, &recvRequest[nbRecvRequest]); 
     328                                nbSendRequest++; 
     329                        } 
     330                } 
     331                if (nbSendElement[rank] > 0) 
     332                { 
     333                        MPI_Irecv(recvValue[rank],  nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
     334                        nbRecvRequest++; 
     335                        MPI_Irecv(recvArea[rank],  nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
     336                        nbRecvRequest++; 
     337                        if (order == 2) 
     338                        { 
     339                                MPI_Irecv(recvGrad[rank], 3*nbSendElement[rank]*(NMAX+1), 
     340                                                MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
     341                                nbRecvRequest++; 
     342                                MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank]*(NMAX+1), MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
    346343//ym  --> attention taille GloId 
    347         nbRecvRequest++; 
    348       } 
    349       else 
    350       { 
    351         MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank], MPI_INT, rank, 4, communicator, &recvRequest[nbRecvRequest]); 
     344                                nbRecvRequest++; 
     345                        } 
     346                        else 
     347                        { 
     348                                MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 
    352349//ym  --> attention taille GloId 
    353         nbRecvRequest++; 
    354       } 
    355     } 
    356   } 
     350                                nbRecvRequest++; 
     351                        } 
     352                } 
     353        } 
    357354         
    358   MPI_Waitall(nbSendRequest, sendRequest, status); 
    359   MPI_Waitall(nbRecvRequest, recvRequest, status); 
     355        MPI_Waitall(nbSendRequest, sendRequest, status); 
     356        MPI_Waitall(nbRecvRequest, recvRequest, status); 
    360357         
    361358 
    362   /* now that all values and gradients are available use them to computed interpolated values on target 
    363     and also to compute weights */ 
    364   int i = 0; 
    365   for (int j = 0; j < nbElements; j++) 
    366   { 
    367     Elt& e = elements[j]; 
    368  
    369     /* since for the 2nd order case source grid elements can contribute to a destination grid element over several "paths" 
    370     (step1: gradient is computed using neighbours on same grid, step2: intersection uses several elements on other grid) 
    371     accumulate them so that there is only one final weight between two elements */ 
    372     map<GloId,double> wgt_map; 
    373  
    374     /* for destination element `e` loop over all intersetions/the corresponding source elements */ 
    375     for (list<Polyg *>::iterator it = e.is.begin(); it != e.is.end(); it++) 
    376     { 
    377       /* it is the intersection element, so it->x and it->area are barycentre and area of intersection element (super mesh) 
    378       but it->id is id of the source element that it intersects */ 
    379       int n1 = (*it)->id.ind; 
    380       int rank = (*it)->id.rank; 
    381       double fk = recvValue[rank][n1]; 
    382       double srcArea = recvArea[rank][n1]; 
    383       double w = (*it)->area; 
     359        /* now that all values and gradients are available use them to computed interpolated values on target 
     360           and also to compute weights */ 
     361        int i = 0; 
     362        for (int j = 0; j < nbElements; j++) 
     363        { 
     364                Elt& e = elements[j]; 
     365 
     366                /* since for the 2nd order case source grid elements can contribute to a destination grid element over several "paths" 
     367                   (step1: gradient is computed using neighbours on same grid, step2: intersection uses several elements on other grid) 
     368                   accumulate them so that there is only one final weight between two elements */ 
     369                map<GloId,double> wgt_map; 
     370 
     371                /* for destination element `e` loop over all intersetions/the corresponding source elements */ 
     372                for (list<Polyg *>::iterator it = e.is.begin(); it != e.is.end(); it++) 
     373                { 
     374                        /* it is the intersection element, so it->x and it->area are barycentre and area of intersection element (super mesh) 
     375                        but it->id is id of the source element that it intersects */ 
     376                        int n1 = (*it)->id.ind; 
     377                        int rank = (*it)->id.rank; 
     378                        double fk = recvValue[rank][n1]; 
     379                        double srcArea = recvArea[rank][n1]; 
     380                        double w = (*it)->area; 
    384381      if (quantity) w/=srcArea ; 
    385382 
    386       /* first order: src value times weight (weight = supermesh area), later divide by target area */ 
    387       int kk = (order == 2) ? n1 * (NMAX + 1) : n1; 
    388       GloId neighID = recvNeighIds[rank][kk]; 
    389       wgt_map[neighID] += w; 
    390  
    391       if (order == 2) 
    392       { 
    393         for (int k = 0; k < NMAX+1; k++) 
    394         { 
    395           int kk = n1 * (NMAX + 1) + k; 
    396           GloId neighID = recvNeighIds[rank][kk]; 
    397           if (neighID.ind != -1)  wgt_map[neighID] += w * scalarprod(recvGrad[rank][kk], (*it)->x); 
    398         } 
    399  
    400       } 
    401     } 
     383                        /* first order: src value times weight (weight = supermesh area), later divide by target area */ 
     384                        int kk = (order == 2) ? n1 * (NMAX + 1) : n1; 
     385                        GloId neighID = recvNeighIds[rank][kk]; 
     386                        wgt_map[neighID] += w; 
     387 
     388                        if (order == 2) 
     389                        { 
     390                                for (int k = 0; k < NMAX+1; k++) 
     391                                { 
     392                                        int kk = n1 * (NMAX + 1) + k; 
     393                                        GloId neighID = recvNeighIds[rank][kk]; 
     394                                        if (neighID.ind != -1)  wgt_map[neighID] += w * scalarprod(recvGrad[rank][kk], (*it)->x); 
     395                                } 
     396 
     397                        } 
     398                } 
    402399 
    403400    double renorm=0; 
    … …  
    407404 
    408405    for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) 
    409     { 
     406                { 
    410407      if (quantity)  this->remapMatrix[i] = (it->second ) / renorm; 
    411       else this->remapMatrix[i] = (it->second / e.area) / renorm; 
    412       this->srcAddress[i] = it->first.ind; 
    413       this->srcRank[i] = it->first.rank; 
    414       this->dstAddress[i] = j; 
     408                        else this->remapMatrix[i] = (it->second / e.area) / renorm; 
     409                        this->srcAddress[i] = it->first.ind; 
     410                        this->srcRank[i] = it->first.rank; 
     411                        this->dstAddress[i] = j; 
    415412      this->sourceWeightId[i]= it->first.globalId ; 
    416413      this->targetWeightId[i]= targetGlobalId[j] ; 
    417       i++; 
    418     } 
    419   } 
    420  
    421   /* free all memory allocated in this function */ 
    422   for (int rank = 0; rank < mpiSize; rank++) 
    423   { 
    424     if (nbSendElement[rank] > 0) 
    425     { 
    426       delete[] sendElement[rank]; 
    427       delete[] recvValue[rank]; 
    428       delete[] recvArea[rank]; 
    429       if (order == 2) 
    430       { 
    431         delete[] recvGrad[rank]; 
    432       } 
    433       delete[] recvNeighIds[rank]; 
    434     } 
    435     if (nbRecvElement[rank] > 0) 
    436     { 
    437       delete[] recvElement[rank]; 
    438       delete[] sendValue[rank]; 
    439       delete[] sendArea[rank]; 
    440       if (order == 2) 
    441         delete[] sendGrad[rank]; 
    442       delete[] sendNeighIds[rank]; 
    443     } 
    444   } 
    445   delete[] status; 
    446   delete[] sendRequest; 
    447   delete[] recvRequest; 
    448   delete[] elementList; 
    449   delete[] nbSendElement; 
    450   delete[] nbRecvElement; 
    451   delete[] sendElement; 
    452   delete[] recvElement; 
    453   delete[] sendValue; 
    454   delete[] recvValue; 
    455   delete[] sendGrad; 
    456   delete[] recvGrad; 
    457   delete[] sendNeighIds; 
    458   delete[] recvNeighIds; 
    459   return i; 
     414                        i++; 
     415                } 
     416        } 
     417 
     418        /* free all memory allocated in this function */ 
     419        for (int rank = 0; rank < mpiSize; rank++) 
     420        { 
     421                if (nbSendElement[rank] > 0) 
     422                { 
     423                        delete[] sendElement[rank]; 
     424                        delete[] recvValue[rank]; 
     425                        delete[] recvArea[rank]; 
     426                        if (order == 2) 
     427                        { 
     428                                delete[] recvGrad[rank]; 
     429                        } 
     430                        delete[] recvNeighIds[rank]; 
     431                } 
     432                if (nbRecvElement[rank] > 0) 
     433                { 
     434                        delete[] recvElement[rank]; 
     435                        delete[] sendValue[rank]; 
     436                        delete[] sendArea[rank]; 
     437                        if (order == 2) 
     438                                delete[] sendGrad[rank]; 
     439                        delete[] sendNeighIds[rank]; 
     440                } 
     441        } 
     442        delete[] status; 
     443        delete[] sendRequest; 
     444        delete[] recvRequest; 
     445        delete[] elementList; 
     446        delete[] nbSendElement; 
     447        delete[] nbRecvElement; 
     448        delete[] sendElement; 
     449        delete[] recvElement; 
     450        delete[] sendValue; 
     451        delete[] recvValue; 
     452        delete[] sendGrad; 
     453        delete[] recvGrad; 
     454        delete[] sendNeighIds; 
     455        delete[] recvNeighIds; 
     456        return i; 
    460457} 
    461458 
    … …  
    496493        mpiRoute.init(routes); 
    497494        int nRecv = mpiRoute.getTotalSourceElement(); 
     495// cout << mpiRank << " NRECV " << nRecv << "(" << routes.size() << ")"<< endl; 
    498496 
    499497        int *nbSendNode = new int[mpiSize]; 
    … …  
    558556        } 
    559557 
    560   MPI_Waitall(nbRecvRequest, recvRequest, status); 
    561   MPI_Waitall(nbSendRequest, sendRequest, status); 
     558        MPI_Waitall(nbRecvRequest, recvRequest, status); 
     559        MPI_Waitall(nbSendRequest, sendRequest, status); 
    562560 
    563561        for (int rank = 0; rank < mpiSize; rank++) 
    … …  
    702700                } 
    703701 
     702/* 
     703                for (int i = 0; i < elt->n; i++) 
     704                { 
     705                        if (elt->neighbour[i] == NOT_FOUND) 
     706                                error_exit("neighbour not found"); 
     707                } 
     708*/ 
    704709        } 
    705710} 
    … …  
    809814        MPI_Waitall(nbRecvRequest, recvRequest, status); 
    810815        MPI_Waitall(nbSendRequest, sendRequest, status); 
    811  
    812816        char **sendBuffer2 = new char*[mpiSize]; 
    813817        char **recvBuffer2 = new char*[mpiSize]; 
    … …  
    836840                                        intersect_ym(&recvElt[j], elt2); 
    837841                                } 
     842 
    838843                                if (recvElt[j].is.size() > 0) sentMessageSize[rank] += packIntersectionSize(recvElt[j]); 
    839844 
    … …  
    854859                        } 
    855860                        delete [] recvElt; 
     861 
    856862                } 
    857863        } 
    … …  
    889895        } 
    890896 
    891   MPI_Waitall(nbRecvRequest, recvRequest, status); 
    892   MPI_Waitall(nbSendRequest, sendRequest, status); 
     897        MPI_Waitall(nbRecvRequest, recvRequest, status); 
     898        MPI_Waitall(nbSendRequest, sendRequest, status); 
    893899 
    894900        delete [] sendRequest; 

  • XIOS/dev/branch_openmp/extern/remap/src/node.cpp

    r1328 r1460  
    281281        } 
    282282 
     283 
    283284  return q; 
    284285} 

  • XIOS/dev/branch_openmp/extern/remap/src/parallel_tree.cpp

    r1355 r1460  
    197197        delete[] displs; 
    198198 
     199        /* unpack */ 
     200/* 
     201        randomArray.resize(nrecv); 
     202        randomizeArray(randomArray); 
     203        tree.leafs.resize(nrecv); 
     204        index = 0; 
     205 
     206 
     207  for (int i = 0; i < nrecv; i++) 
     208        { 
     209                Coord x = *(Coord *)(&recvBuffer[index]); 
     210                index += sizeof(Coord)/sizeof(*recvBuffer); 
     211                double radius = recvBuffer[index++]; 
     212                tree.leafs[randomArray[i]].centre = x; 
     213                tree.leafs[randomArray[i]].radius = radius; 
     214 
     215        } 
     216*/ 
    199217 
    200218  randomArray.resize(blocSize); 
    … …  
    231249    cerr << comm.rank << ": PROBLEM: (node assign)" << tree.levelSize[assignLevel] << " != " << comm.group_size << " (keepNodes)"  
    232250         << "   node size : "<<node.size()<<"   bloc size : "<<blocSize<<"  total number of leaf : "<<tree.leafs.size()<<endl ; 
    233  
     251/* 
     252        MPI_Allreduce(&ok, &allok, 1, MPI_INT, MPI_PROD, communicator); 
     253        if (!allok) { 
     254                MPI_Finalize(); 
     255                exit(1); 
     256        } 
     257*/ 
    234258    MPI_Abort(MPI_COMM_WORLD,-1) ; 
    235259  } 
    … …  
    302326  nb=nb1+nb2 ; 
    303327  MPI_Allreduce(&nb, &nbTot, 1, MPI_LONG, MPI_SUM, communicator) ; 
    304  
    305  
    306328  int commSize ; 
    307329  MPI_Comm_size(communicator,&commSize) ; 
    … …  
    326348        randomizeArray(randomArray2); 
    327349 
     350/*       
     351        int s1,s2 ; 
     352         
     353        if (node.size()< nbSampleNodes/2) 
     354        {  
     355          s1 = node.size() ; 
     356          s2 = nbSampleNodes-s1 ; 
     357        } 
     358        else if (node2.size()< nbSampleNodes/2) 
     359        { 
     360          s2 = node.size() ; 
     361          s1 = nbSampleNodes-s2 ; 
     362        } 
     363        else 
     364        { 
     365          s1=nbSampleNodes/2 ; 
     366          s2=nbSampleNodes/2 ; 
     367        } 
     368*/ 
    328369        for (int i = 0; i <nbSampleNodes1; i++) sampleNodes.push_back(Node(node[randomArray1[i%nb1]].centre,  node[randomArray1[i%nb1]].radius, NULL)); 
    329370        for (int i = 0; i <nbSampleNodes2; i++) sampleNodes.push_back(Node(node2[randomArray2[i%nb2]].centre, node2[randomArray2[i%nb2]].radius, NULL)); 
    330371 
     372/*           
     373        for (int i = 0; i < nbSampleNodes/2; i++) 
     374        { 
     375          sampleNodes.push_back(Node(node[randomArray1[i]].centre,  node[randomArray1[i]].radius, NULL)); 
     376          sampleNodes.push_back(Node(node2[randomArray2[i]].centre, node2[randomArray2[i]].radius, NULL)); 
     377        } 
     378*/ 
    331379        CTimer::get("buildParallelSampleTree").resume(); 
    332380        //sampleTree.buildParallelSampleTree(sampleNodes, cascade); 

  • XIOS/dev/branch_openmp/extern/remap/src/parallel_tree.hpp

    r1328 r1460  
    66#include "mpi_cascade.hpp" 
    77#include "mpi.hpp" 
     8 
    89namespace sphereRemap { 
    910 
Note: See TracChangeset for help on using the changeset viewer.