source: XIOS/trunk/extern/remap/src/mapper.cpp @ 923

Last change on this file since 923 was 923, checked in by ymipsl, 5 years ago

Improve remapper :

  • increase stability for smal grid remapping
  • increase stability for large number of processes

YM

File size: 28.2 KB
Line 
1#include "mpi.hpp"
2#include <map>
3#include "cputime.hpp"       /* time */
4#include "meshutil.hpp"
5#include "polyg.hpp"
6#include "circle.hpp"
7#include "intersect.hpp"
8#include "intersection_ym.hpp"
9#include "errhandle.hpp"
10#include "mpi_routing.hpp"
11#include "gridRemap.hpp"
12
13#include "mapper.hpp"
14
15namespace sphereRemap {
16
17/* A subdivition of an array into N sub-arays
18   can be represented by the length of the N arrays
19   or by the offsets when each of the N arrays starts.
20   This function convertes from the former to the latter. */
21void cptOffsetsFromLengths(const int *lengths, int *offsets, int sz)
22{
23        offsets[0] = 0;
24        for (int i = 1; i < sz; i++)
25                offsets[i] = offsets[i-1] + lengths[i-1];
26}
27
28
29void Mapper::setSourceMesh(const double* boundsLon, const double* boundsLat, int nVertex, int nbCells, const double* pole, const long int* globalId)
30{
31  srcGrid.pole = Coord(pole[0], pole[1], pole[2]);
32
33        int mpiRank, mpiSize;
34        MPI_Comm_rank(communicator, &mpiRank);
35        MPI_Comm_size(communicator, &mpiSize);
36
37        sourceElements.reserve(nbCells);
38        sourceMesh.reserve(nbCells);
39  sourceGlobalId.resize(nbCells) ;
40
41  if (globalId==NULL)
42  {
43    long int offset ;
44    long int nb=nbCells ;
45    MPI_Scan(&nb,&offset,1,MPI_LONG,MPI_SUM,communicator) ;
46    offset=offset-nb ;
47    for(int i=0;i<nbCells;i++) sourceGlobalId[i]=offset+i ;
48  }
49  else sourceGlobalId.assign(globalId,globalId+nbCells);
50
51        for (int i = 0; i < nbCells; i++)
52        {
53                int offs = i*nVertex;
54                Elt elt(boundsLon + offs, boundsLat + offs, nVertex);
55                elt.src_id.rank = mpiRank;
56                elt.src_id.ind = i;
57    elt.src_id.globalId = sourceGlobalId[i];
58                sourceElements.push_back(elt);
59                sourceMesh.push_back(Node(elt.x, cptRadius(elt), &sourceElements.back()));
60                cptEltGeom(sourceElements[i], Coord(pole[0], pole[1], pole[2]));
61        }
62
63}
64
65void Mapper::setTargetMesh(const double* boundsLon, const double* boundsLat, int nVertex, int nbCells, const double* pole, const long int* globalId)
66{
67  tgtGrid.pole = Coord(pole[0], pole[1], pole[2]);
68
69        int mpiRank, mpiSize;
70        MPI_Comm_rank(communicator, &mpiRank);
71        MPI_Comm_size(communicator, &mpiSize);
72
73        targetElements.reserve(nbCells);
74        targetMesh.reserve(nbCells);
75
76  targetGlobalId.resize(nbCells) ;
77  if (globalId==NULL)
78  {
79    long int offset ;
80    long int nb=nbCells ;
81    MPI_Scan(&nb,&offset,1,MPI_LONG,MPI_SUM,communicator) ;
82    offset=offset-nb ;
83    for(int i=0;i<nbCells;i++) targetGlobalId[i]=offset+i ;
84  }
85  else targetGlobalId.assign(globalId,globalId+nbCells);
86
87        for (int i = 0; i < nbCells; i++)
88        {
89                int offs = i*nVertex;
90                Elt elt(boundsLon + offs, boundsLat + offs, nVertex);
91                targetElements.push_back(elt);
92                targetMesh.push_back(Node(elt.x, cptRadius(elt), &sourceElements.back()));
93                cptEltGeom(targetElements[i], Coord(pole[0], pole[1], pole[2]));
94        }
95
96
97}
98
99void Mapper::setSourceValue(const double* val)
100{
101  int size=sourceElements.size() ;
102  for(int i=0;i<size;++i) sourceElements[i].val=val[i] ;
103}
104
105void Mapper::getTargetValue(double* val)
106{
107  int size=targetElements.size() ;
108  for(int i=0;i<size;++i) val[i]=targetElements[i].val ;
109}
110
111vector<double> Mapper::computeWeights(int interpOrder, bool renormalize)
112{
113        vector<double> timings;
114        int mpiSize, mpiRank;
115        MPI_Comm_size(communicator, &mpiSize);
116        MPI_Comm_rank(communicator, &mpiRank);
117
118  this->buildSSTree(sourceMesh, targetMesh);
119
120        if (mpiRank == 0 && verbose) cout << "Computing intersections ..." << endl;
121        double tic = cputime();
122        computeIntersection(&targetElements[0], targetElements.size());
123        timings.push_back(cputime() - tic);
124
125        tic = cputime();
126        if (interpOrder == 2) {
127                if (mpiRank == 0 && verbose) cout << "Computing grads ..." << endl;
128                buildMeshTopology();
129                computeGrads();
130        }
131        timings.push_back(cputime() - tic);
132
133        /* Prepare computation of weights */
134        /* compute number of intersections which for the first order case
135           corresponds to the number of edges in the remap matrix */
136        int nIntersections = 0;
137        for (int j = 0; j < targetElements.size(); j++)
138        {
139                Elt &elt = targetElements[j];
140                for (list<Polyg*>::iterator it = elt.is.begin(); it != elt.is.end(); it++)
141                        nIntersections++;
142        }
143        /* overallocate for NMAX neighbours for each elements */
144        remapMatrix = new double[nIntersections*NMAX];
145        srcAddress = new int[nIntersections*NMAX];
146        srcRank = new int[nIntersections*NMAX];
147        dstAddress = new int[nIntersections*NMAX];
148  sourceWeightId =new long[nIntersections*NMAX];
149  targetWeightId =new long[nIntersections*NMAX];
150
151
152        if (mpiRank == 0 && verbose) cout << "Remapping..." << endl;
153        tic = cputime();
154        nWeights = remap(&targetElements[0], targetElements.size(), interpOrder, renormalize);
155        timings.push_back(cputime() - tic);
156
157  for (int i = 0; i < targetElements.size(); i++) targetElements[i].delete_intersections();
158
159        return timings;
160}
161
162/**
163   @param elements are cells of the target grid that are distributed over CPUs
164          indepentently of the distribution of the SS-tree.
165   @param nbElements is the size of the elements array.
166   @param order is the order of interpolaton (must be 1 or 2).
167*/
168int Mapper::remap(Elt *elements, int nbElements, int order, bool renormalize)
169{
170        int mpiSize, mpiRank;
171        MPI_Comm_size(communicator, &mpiSize);
172        MPI_Comm_rank(communicator, &mpiRank);
173
174        /* create list of intersections (super mesh elements) for each rank */
175        multimap<int, Polyg *> *elementList = new multimap<int, Polyg *>[mpiSize];
176        for (int j = 0; j < nbElements; j++)
177        {
178                Elt& e = elements[j];
179                for (list<Polyg *>::iterator it = e.is.begin(); it != e.is.end(); it++)
180                        elementList[(*it)->id.rank].insert(pair<int, Polyg *>((*it)->id.ind, *it));
181        }
182
183        int *nbSendElement = new int[mpiSize];
184        int **sendElement = new int*[mpiSize]; /* indices of elements required from other rank */
185        double **recvValue = new double*[mpiSize];
186        Coord **recvGrad = new Coord*[mpiSize];
187        GloId **recvNeighIds = new GloId*[mpiSize]; /* ids of the of the source neighbours which also contribute through gradient */
188        for (int rank = 0; rank < mpiSize; rank++)
189        {
190                /* get size for allocation */
191                int last = -1; /* compares unequal to any index */
192                int index = -1; /* increased to starting index 0 in first iteration */
193                for (multimap<int, Polyg *>::iterator it = elementList[rank].begin(); it != elementList[rank].end(); ++it)
194                {
195                        if (last != it->first)
196                                index++;
197                        (it->second)->id.ind = index;
198                        last = it->first;
199                }
200                nbSendElement[rank] = index + 1;
201
202                /* if size is non-zero allocate and collect indices of elements on other ranks that we intersect */
203                if (nbSendElement[rank] > 0)
204                {
205                        sendElement[rank] = new int[nbSendElement[rank]];
206                        recvValue[rank]   = new double[nbSendElement[rank]];
207                        if (order == 2)
208                        {
209                                recvNeighIds[rank] = new GloId[nbSendElement[rank]*(NMAX+1)];
210                                recvGrad[rank]    = new Coord[nbSendElement[rank]*(NMAX+1)];
211                        }
212                        else
213                                recvNeighIds[rank] = new GloId[nbSendElement[rank]];
214
215                        last = -1;
216                        index = -1;
217                        for (multimap<int, Polyg *>::iterator it = elementList[rank].begin(); it != elementList[rank].end(); ++it)
218                        {
219                                if (last != it->first)
220                                        index++;
221                                sendElement[rank][index] = it->first;
222                                last = it->first;
223                        }
224                }
225        }
226
227        /* communicate sizes of source elements to be sent (index lists and later values and gradients) */
228        int *nbRecvElement = new int[mpiSize];
229        MPI_Alltoall(nbSendElement, 1, MPI_INT, nbRecvElement, 1, MPI_INT, communicator);
230
231        /* communicate indices of source elements on other ranks whoes value and gradient we need (since intersection) */
232        int nbSendRequest = 0;
233        int nbRecvRequest = 0;
234        int **recvElement = new int*[mpiSize];
235        double **sendValue = new double*[mpiSize];
236        Coord **sendGrad = new Coord*[mpiSize];
237        GloId **sendNeighIds = new GloId*[mpiSize];
238        MPI_Request *sendRequest = new MPI_Request[3*mpiSize];
239        MPI_Request *recvRequest = new MPI_Request[3*mpiSize];
240        for (int rank = 0; rank < mpiSize; rank++)
241        {
242                if (nbSendElement[rank] > 0)
243                {
244                        MPI_Issend(sendElement[rank], nbSendElement[rank], MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]);
245                        nbSendRequest++;
246                }
247
248                if (nbRecvElement[rank] > 0)
249                {
250                        recvElement[rank] = new int[nbRecvElement[rank]];
251                        sendValue[rank]   = new double[nbRecvElement[rank]];
252                        if (order == 2)
253                        {
254                                sendNeighIds[rank] = new GloId[nbRecvElement[rank]*(NMAX+1)];
255                                sendGrad[rank]    = new Coord[nbRecvElement[rank]*(NMAX+1)];
256                        }
257                        else
258                        {
259                                sendNeighIds[rank] = new GloId[nbRecvElement[rank]];
260                        }
261                        MPI_Irecv(recvElement[rank], nbRecvElement[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]);
262                        nbRecvRequest++;
263                }
264        }
265        MPI_Status *status = new MPI_Status[3*mpiSize];
266        MPI_Waitall(nbRecvRequest, recvRequest, status);
267        MPI_Waitall(nbSendRequest, sendRequest, status);
268
269        /* for all indices that have been received from requesting ranks: pack values and gradients, then send */
270        nbSendRequest = 0;
271        nbRecvRequest = 0;
272        for (int rank = 0; rank < mpiSize; rank++)
273        {
274                if (nbRecvElement[rank] > 0)
275                {
276                        int jj = 0; // jj == j if no weight writing
277                        for (int j = 0; j < nbRecvElement[rank]; j++)
278                        {
279                                sendValue[rank][j] = sstree.localElements[recvElement[rank][j]].val;
280                                if (order == 2)
281                                {
282                                        sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].grad;
283//          cout<<"grad  "<<jj<<"  "<<recvElement[rank][j]<<"  "<<sendGrad[rank][jj]<<" "<<sstree.localElements[recvElement[rank][j]].grad<<endl ;
284                                        sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].src_id;
285                                        jj++;
286                                        for (int i = 0; i < NMAX; i++)
287                                        {
288                                                sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].gradNeigh[i];
289//            cout<<"grad  "<<jj<<"  "<<sendGrad[rank][jj]<<" "<<sstree.localElements[recvElement[rank][j]].grad<<endl ;
290            sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].neighId[i];
291                                                jj++;
292                                        }
293                                }
294                                else
295                                        sendNeighIds[rank][j] = sstree.localElements[recvElement[rank][j]].src_id;
296                        }
297                        MPI_Issend(sendValue[rank],  nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]);
298                        nbSendRequest++;
299                        if (order == 2)
300                        {
301                                MPI_Issend(sendGrad[rank], 3*nbRecvElement[rank]*(NMAX+1),
302                                                                MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]);
303                                nbSendRequest++;
304                                MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank]*(NMAX+1), MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]);
305//ym  --> attention taille GloId
306                                nbSendRequest++;
307                        }
308                        else
309                        {
310                                MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank], MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]);
311//ym  --> attention taille GloId
312                                nbSendRequest++;
313                        }
314                }
315                if (nbSendElement[rank] > 0)
316                {
317                        MPI_Irecv(recvValue[rank],  nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]);
318                        nbRecvRequest++;
319                        if (order == 2)
320                        {
321                                MPI_Irecv(recvGrad[rank], 3*nbSendElement[rank]*(NMAX+1),
322                                                MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]);
323                                nbRecvRequest++;
324                                MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank]*(NMAX+1), MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]);
325//ym  --> attention taille GloId
326                                nbRecvRequest++;
327                        }
328                        else
329                        {
330                                MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]);
331//ym  --> attention taille GloId
332                                nbRecvRequest++;
333                        }
334                }
335        }
336        MPI_Waitall(nbRecvRequest, recvRequest, status);
337        MPI_Waitall(nbSendRequest, sendRequest, status);
338
339        /* now that all values and gradients are available use them to computed interpolated values on target
340           and also to compute weights */
341        int i = 0;
342        for (int j = 0; j < nbElements; j++)
343        {
344                Elt& e = elements[j];
345
346                /* since for the 2nd order case source grid elements can contribute to a destination grid element over several "paths"
347                   (step1: gradient is computed using neighbours on same grid, step2: intersection uses several elements on other grid)
348                   accumulate them so that there is only one final weight between two elements */
349                map<GloId,double> wgt_map;
350
351                /* for destination element `e` loop over all intersetions/the corresponding source elements */
352                for (list<Polyg *>::iterator it = e.is.begin(); it != e.is.end(); it++)
353                {
354                        /* it is the intersection element, so it->x and it->area are barycentre and area of intersection element (super mesh)
355                        but it->id is id of the source element that it intersects */
356                        int n1 = (*it)->id.ind;
357                        int rank = (*it)->id.rank;
358                        double fk = recvValue[rank][n1];
359                        double w = (*it)->area;
360
361                        /* first order: src value times weight (weight = supermesh area), later divide by target area */
362                        int kk = (order == 2) ? n1 * (NMAX + 1) : n1;
363                        GloId neighID = recvNeighIds[rank][kk];
364                        wgt_map[neighID] += (*it)->area;
365
366                        if (order == 2)
367                        {
368                                for (int k = 0; k < NMAX+1; k++)
369                                {
370                                        int kk = n1 * (NMAX + 1) + k;
371                                        GloId neighID = recvNeighIds[rank][kk];
372                                        if (neighID.ind != -1)  wgt_map[neighID] += w * scalarprod(recvGrad[rank][kk], (*it)->x);
373                                }
374
375                        }
376                }
377
378    double renorm=0;
379    if (renormalize) 
380      for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) renorm+=it->second / e.area;
381    else renorm=1. ;
382
383    for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++)
384                {
385                        this->remapMatrix[i] = (it->second / e.area) / renorm;
386                        this->srcAddress[i] = it->first.ind;
387                        this->srcRank[i] = it->first.rank;
388                        this->dstAddress[i] = j;
389      this->sourceWeightId[i]= it->first.globalId ;
390      this->targetWeightId[i]= targetGlobalId[j] ;
391                        i++;
392                }
393        }
394
395        /* free all memory allocated in this function */
396        for (int rank = 0; rank < mpiSize; rank++)
397        {
398                if (nbSendElement[rank] > 0)
399                {
400                        delete[] sendElement[rank];
401                        delete[] recvValue[rank];
402                        if (order == 2)
403                        {
404                                delete[] recvGrad[rank];
405                        }
406                        delete[] recvNeighIds[rank];
407                }
408                if (nbRecvElement[rank] > 0)
409                {
410                        delete[] recvElement[rank];
411                        delete[] sendValue[rank];
412                        if (order == 2)
413                                delete[] sendGrad[rank];
414                        delete[] sendNeighIds[rank];
415                }
416        }
417        delete[] status;
418        delete[] sendRequest;
419        delete[] recvRequest;
420        delete[] elementList;
421        delete[] nbSendElement;
422        delete[] nbRecvElement;
423        delete[] sendElement;
424        delete[] recvElement;
425        delete[] sendValue;
426        delete[] recvValue;
427        delete[] sendGrad;
428        delete[] recvGrad;
429        delete[] sendNeighIds;
430        delete[] recvNeighIds;
431        return i;
432}
433
434void Mapper::computeGrads()
435{
436        /* array of pointers to collect local elements and elements received from other cpu */
437        vector<Elt*> globalElements(sstree.nbLocalElements + nbNeighbourElements);
438        int index = 0;
439        for (int i = 0; i < sstree.nbLocalElements; i++, index++)
440                globalElements[index] = &(sstree.localElements[i]);
441        for (int i = 0; i < nbNeighbourElements; i++, index++)
442                globalElements[index] = &neighbourElements[i];
443
444        update_baryc(sstree.localElements, sstree.nbLocalElements);
445        computeGradients(&globalElements[0], sstree.nbLocalElements);
446}
447
448/** for each element of the source grid, finds all the neighbouring elements that share an edge
449    (filling array neighbourElements). This is used later to compute gradients */
450void Mapper::buildMeshTopology()
451{
452        int mpiSize, mpiRank;
453        MPI_Comm_size(communicator, &mpiSize);
454        MPI_Comm_rank(communicator, &mpiRank);
455
456        vector<Node> *routingList = new vector<Node>[mpiSize];
457        vector<vector<int> > routes(sstree.localTree.leafs.size());
458
459        sstree.routeIntersections(routes, sstree.localTree.leafs);
460
461        for (int i = 0; i < routes.size(); ++i)
462                for (int k = 0; k < routes[i].size(); ++k)
463                        routingList[routes[i][k]].push_back(sstree.localTree.leafs[i]);
464        routingList[mpiRank].clear();
465
466
467        CMPIRouting mpiRoute(communicator);
468        mpiRoute.init(routes);
469        int nRecv = mpiRoute.getTotalSourceElement();
470// cout << mpiRank << " NRECV " << nRecv << "(" << routes.size() << ")"<< endl;
471
472        int *nbSendNode = new int[mpiSize];
473        int *nbRecvNode = new int[mpiSize];
474        int *sendMessageSize = new int[mpiSize];
475        int *recvMessageSize = new int[mpiSize];
476
477        for (int rank = 0; rank < mpiSize; rank++)
478        {
479                nbSendNode[rank] = routingList[rank].size();
480                sendMessageSize[rank] = 0;
481                for (size_t j = 0; j < routingList[rank].size(); j++)
482                {
483                        Elt *elt = (Elt *) (routingList[rank][j].data);
484                        sendMessageSize[rank] += packedPolygonSize(*elt);
485                }
486        }
487
488        MPI_Alltoall(nbSendNode, 1, MPI_INT, nbRecvNode, 1, MPI_INT, communicator);
489        MPI_Alltoall(sendMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator);
490
491        char **sendBuffer = new char*[mpiSize];
492        char **recvBuffer = new char*[mpiSize];
493        int *pos = new int[mpiSize];
494
495        for (int rank = 0; rank < mpiSize; rank++)
496        {
497                if (nbSendNode[rank] > 0) sendBuffer[rank] = new char[sendMessageSize[rank]];
498                if (nbRecvNode[rank] > 0) recvBuffer[rank] = new char[recvMessageSize[rank]];
499        }
500
501        for (int rank = 0; rank < mpiSize; rank++)
502        {
503                pos[rank] = 0;
504                for (size_t j = 0; j < routingList[rank].size(); j++)
505                {
506                        Elt *elt = (Elt *) (routingList[rank][j].data);
507                        packPolygon(*elt, sendBuffer[rank], pos[rank]);
508                }
509        }
510        delete [] routingList;
511
512
513        int nbSendRequest = 0;
514        int nbRecvRequest = 0;
515        MPI_Request *sendRequest = new MPI_Request[mpiSize];
516        MPI_Request *recvRequest = new MPI_Request[mpiSize];
517        MPI_Status  *status      = new MPI_Status[mpiSize];
518
519        for (int rank = 0; rank < mpiSize; rank++)
520        {
521                if (nbSendNode[rank] > 0)
522                {
523                        MPI_Issend(sendBuffer[rank], sendMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]);
524                        nbSendRequest++;
525                }
526                if (nbRecvNode[rank] > 0)
527                {
528                        MPI_Irecv(recvBuffer[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]);
529                        nbRecvRequest++;
530                }
531        }
532
533        MPI_Waitall(nbRecvRequest, recvRequest, status);
534        MPI_Waitall(nbSendRequest, sendRequest, status);
535
536        for (int rank = 0; rank < mpiSize; rank++)
537                if (nbSendNode[rank] > 0) delete [] sendBuffer[rank];
538        delete [] sendBuffer;
539
540        char **sendBuffer2 = new char*[mpiSize];
541        char **recvBuffer2 = new char*[mpiSize];
542
543        for (int rank = 0; rank < mpiSize; rank++)
544        {
545                nbSendNode[rank] = 0;
546                sendMessageSize[rank] = 0;
547
548                if (nbRecvNode[rank] > 0)
549                {
550                        set<NodePtr> neighbourList;
551                        pos[rank] = 0;
552                        for (int j = 0; j < nbRecvNode[rank]; j++)
553                        {
554                                Elt elt;
555                                unpackPolygon(elt, recvBuffer[rank], pos[rank]);
556                                Node node(elt.x, cptRadius(elt), &elt);
557                                findNeighbour(sstree.localTree.root, &node, neighbourList);
558                        }
559                        nbSendNode[rank] = neighbourList.size();
560                        for (set<NodePtr>::iterator it = neighbourList.begin(); it != neighbourList.end(); it++)
561                        {
562                                Elt *elt = (Elt *) ((*it)->data);
563                                sendMessageSize[rank] += packedPolygonSize(*elt);
564                        }
565
566                        sendBuffer2[rank] = new char[sendMessageSize[rank]];
567                        pos[rank] = 0;
568
569                        for (set<NodePtr>::iterator it = neighbourList.begin(); it != neighbourList.end(); it++)
570                        {
571                                Elt *elt = (Elt *) ((*it)->data);
572                                packPolygon(*elt, sendBuffer2[rank], pos[rank]);
573                        }
574                }
575        }
576        for (int rank = 0; rank < mpiSize; rank++)
577                if (nbRecvNode[rank] > 0) delete [] recvBuffer[rank];
578        delete [] recvBuffer;
579
580
581        MPI_Barrier(communicator);
582        MPI_Alltoall(nbSendNode, 1, MPI_INT, nbRecvNode, 1, MPI_INT, communicator);
583        MPI_Alltoall(sendMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator);
584
585        for (int rank = 0; rank < mpiSize; rank++)
586                if (nbRecvNode[rank] > 0) recvBuffer2[rank] = new char[recvMessageSize[rank]];
587
588        nbSendRequest = 0;
589        nbRecvRequest = 0;
590
591        for (int rank = 0; rank < mpiSize; rank++)
592        {
593                if (nbSendNode[rank] > 0)
594                {
595                        MPI_Issend(sendBuffer2[rank], sendMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]);
596                        nbSendRequest++;
597                }
598                if (nbRecvNode[rank] > 0)
599                {
600                        MPI_Irecv(recvBuffer2[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]);
601                        nbRecvRequest++;
602                }
603        }
604
605        MPI_Waitall(nbRecvRequest, recvRequest, status);
606        MPI_Waitall(nbSendRequest, sendRequest, status);
607
608        int nbNeighbourNodes = 0;
609        for (int rank = 0; rank < mpiSize; rank++)
610                nbNeighbourNodes += nbRecvNode[rank];
611
612        neighbourElements = new Elt[nbNeighbourNodes];
613        nbNeighbourElements = nbNeighbourNodes;
614
615        int index = 0;
616        for (int rank = 0; rank < mpiSize; rank++)
617        {
618                pos[rank] = 0;
619                for (int j = 0; j < nbRecvNode[rank]; j++)
620                {
621                        unpackPolygon(neighbourElements[index], recvBuffer2[rank], pos[rank]);
622                        neighbourElements[index].id.ind = sstree.localTree.leafs.size() + index;
623                        index++;
624                }
625        }
626        for (int rank = 0; rank < mpiSize; rank++)
627        {
628                if (nbRecvNode[rank] > 0) delete [] recvBuffer2[rank];
629                if (nbSendNode[rank] > 0) delete [] sendBuffer2[rank];
630        }
631        delete [] recvBuffer2;
632        delete [] sendBuffer2;
633        delete [] sendMessageSize;
634        delete [] recvMessageSize;
635        delete [] nbSendNode;
636        delete [] nbRecvNode;
637        delete [] sendRequest;
638        delete [] recvRequest;
639        delete [] status;
640        delete [] pos;
641
642        /* re-compute on received elements to avoid having to send this information */
643        neighbourNodes.resize(nbNeighbourNodes);
644        setCirclesAndLinks(neighbourElements, neighbourNodes);
645        cptAllEltsGeom(neighbourElements, nbNeighbourNodes, srcGrid.pole);
646
647        /* the local SS tree must include nodes from other cpus if they are potential
648           intersector of a local node */
649        sstree.localTree.insertNodes(neighbourNodes);
650
651        /* for every local element,
652           use the SS-tree to find all elements (including neighbourElements)
653           who are potential neighbours because their circles intersect,
654           then check all canditates for common edges to build up connectivity information
655        */
656        for (int j = 0; j < sstree.localTree.leafs.size(); j++)
657        {
658                Node& node = sstree.localTree.leafs[j];
659
660                /* find all leafs whoes circles that intersect node's circle and save into node->intersectors */
661                node.search(sstree.localTree.root);
662
663                Elt *elt = (Elt *)(node.data);
664
665                for (int i = 0; i < elt->n; i++) elt->neighbour[i] = NOT_FOUND;
666
667                /* for element `elt` loop through all nodes in the SS-tree
668                   whoes circles intersect with the circle around `elt` (the SS intersectors)
669                   and check if they are neighbours in the sense that the two elements share an edge.
670                   If they do, save this information for elt */
671                for (list<NodePtr>::iterator it = (node.intersectors).begin(); it != (node.intersectors).end(); ++it)
672                {
673                        Elt *elt2 = (Elt *)((*it)->data);
674                        set_neighbour(*elt, *elt2);
675                }
676
677/*
678                for (int i = 0; i < elt->n; i++)
679                {
680                        if (elt->neighbour[i] == NOT_FOUND)
681                                error_exit("neighbour not found");
682                }
683*/
684        }
685}
686
687/** @param elements are the target grid elements */
688void Mapper::computeIntersection(Elt *elements, int nbElements)
689{
690        int mpiSize, mpiRank;
691        MPI_Comm_size(communicator, &mpiSize);
692        MPI_Comm_rank(communicator, &mpiRank);
693
694        MPI_Barrier(communicator);
695
696        vector<Node> *routingList = new vector<Node>[mpiSize];
697
698        vector<Node> routeNodes;  routeNodes.reserve(nbElements);
699        for (int j = 0; j < nbElements; j++)
700        {
701                elements[j].id.ind = j;
702                elements[j].id.rank = mpiRank;
703                routeNodes.push_back(Node(elements[j].x, cptRadius(elements[j]), &elements[j]));
704        }
705
706        vector<vector<int> > routes(routeNodes.size());
707        sstree.routeIntersections(routes, routeNodes);
708        for (int i = 0; i < routes.size(); ++i)
709                for (int k = 0; k < routes[i].size(); ++k)
710                        routingList[routes[i][k]].push_back(routeNodes[i]);
711
712        if (verbose >= 2)
713        {
714                cout << " --> rank  " << mpiRank << " nbElements " << nbElements << " : ";
715                for (int rank = 0; rank < mpiSize; rank++)
716                        cout << routingList[rank].size() << "   ";
717                cout << endl;
718        }
719        MPI_Barrier(communicator);
720
721        int *nbSendNode = new int[mpiSize];
722        int *nbRecvNode = new int[mpiSize];
723        int *sentMessageSize = new int[mpiSize];
724        int *recvMessageSize = new int[mpiSize];
725
726        for (int rank = 0; rank < mpiSize; rank++)
727        {
728                nbSendNode[rank] = routingList[rank].size();
729                sentMessageSize[rank] = 0;
730                for (size_t j = 0; j < routingList[rank].size(); j++)
731                {
732                        Elt *elt = (Elt *) (routingList[rank][j].data);
733                        sentMessageSize[rank] += packedPolygonSize(*elt);
734                }
735        }
736
737        MPI_Alltoall(nbSendNode, 1, MPI_INT, nbRecvNode, 1, MPI_INT, communicator);
738        MPI_Alltoall(sentMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator);
739
740        int total = 0;
741
742        for (int rank = 0; rank < mpiSize; rank++)
743        {
744                total = total + nbRecvNode[rank];
745        }
746
747        if (verbose >= 2) cout << "---> rank " << mpiRank << " : compute intersection : total received nodes  " << total << endl;
748        char **sendBuffer = new char*[mpiSize];
749        char **recvBuffer = new char*[mpiSize];
750        int *pos = new int[mpiSize];
751
752        for (int rank = 0; rank < mpiSize; rank++)
753        {
754                if (nbSendNode[rank] > 0) sendBuffer[rank] = new char[sentMessageSize[rank]];
755                if (nbRecvNode[rank] > 0) recvBuffer[rank] = new char[recvMessageSize[rank]];
756        }
757
758        for (int rank = 0; rank < mpiSize; rank++)
759        {
760                pos[rank] = 0;
761                for (size_t j = 0; j < routingList[rank].size(); j++)
762                {
763                        Elt* elt = (Elt *) (routingList[rank][j].data);
764                        packPolygon(*elt, sendBuffer[rank], pos[rank]);
765                }
766        }
767        delete [] routingList;
768
769        int nbSendRequest = 0;
770        int nbRecvRequest = 0;
771        MPI_Request *sendRequest = new MPI_Request[mpiSize];
772        MPI_Request *recvRequest = new MPI_Request[mpiSize];
773        MPI_Status   *status = new MPI_Status[mpiSize];
774
775        for (int rank = 0; rank < mpiSize; rank++)
776        {
777                if (nbSendNode[rank] > 0)
778                {
779                        MPI_Issend(sendBuffer[rank], sentMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]);
780                        nbSendRequest++;
781                }
782                if (nbRecvNode[rank] > 0)
783                {
784                        MPI_Irecv(recvBuffer[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]);
785                        nbRecvRequest++;
786                }
787        }
788
789        MPI_Waitall(nbRecvRequest, recvRequest, status);
790        MPI_Waitall(nbSendRequest, sendRequest, status);
791        char **sendBuffer2 = new char*[mpiSize];
792        char **recvBuffer2 = new char*[mpiSize];
793
794        double tic = cputime();
795        for (int rank = 0; rank < mpiSize; rank++)
796        {
797                sentMessageSize[rank] = 0;
798
799                if (nbRecvNode[rank] > 0)
800                {
801                        Elt *recvElt = new Elt[nbRecvNode[rank]];
802                        pos[rank] = 0;
803                        for (int j = 0; j < nbRecvNode[rank]; j++)
804                        {
805                                unpackPolygon(recvElt[j], recvBuffer[rank], pos[rank]);
806                                cptEltGeom(recvElt[j], tgtGrid.pole);
807                                Node recvNode(recvElt[j].x, cptRadius(recvElt[j]), &recvElt[j]);
808                                recvNode.search(sstree.localTree.root);
809                                /* for a node holding an element of the target, loop throught candidates for intersecting source */
810                                for (list<NodePtr>::iterator it = (recvNode.intersectors).begin(); it != (recvNode.intersectors).end(); ++it)
811                                {
812                                        Elt *elt2 = (Elt *) ((*it)->data);
813                                        /* recvElt is target, elt2 is source */
814//                                      intersect(&recvElt[j], elt2);
815                                        intersect_ym(&recvElt[j], elt2);
816                                }
817
818                                if (recvElt[j].is.size() > 0) sentMessageSize[rank] += packIntersectionSize(recvElt[j]);
819
820                                // here recvNode goes out of scope
821                        }
822
823                        if (sentMessageSize[rank] > 0)
824                        {
825                                sentMessageSize[rank] += sizeof(int);
826                                sendBuffer2[rank] = new char[sentMessageSize[rank]];
827                                *((int *) sendBuffer2[rank]) = 0;
828                                pos[rank] = sizeof(int);
829                                for (int j = 0; j < nbRecvNode[rank]; j++)
830                                {
831                                        packIntersection(recvElt[j], sendBuffer2[rank], pos[rank]);
832                                        //FIXME should be deleted: recvElt[j].delete_intersections(); // intersection areas have been packed to buffer and won't be used any more
833                                }
834                        }
835                        delete [] recvElt;
836
837                }
838        }
839        delete [] pos;
840
841        for (int rank = 0; rank < mpiSize; rank++)
842        {
843                if (nbSendNode[rank] > 0) delete [] sendBuffer[rank];
844                if (nbRecvNode[rank] > 0) delete [] recvBuffer[rank];
845                nbSendNode[rank] = 0;
846        }
847
848        if (verbose >= 2) cout << "Rank " << mpiRank << "  Compute (internal) intersection " << cputime() - tic << " s" << endl;
849        MPI_Alltoall(sentMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator);
850
851        for (int rank = 0; rank < mpiSize; rank++)
852                if (recvMessageSize[rank] > 0)
853                        recvBuffer2[rank] = new char[recvMessageSize[rank]];
854
855        nbSendRequest = 0;
856        nbRecvRequest = 0;
857
858        for (int rank = 0; rank < mpiSize; rank++)
859        {
860                if (sentMessageSize[rank] > 0)
861                {
862                        MPI_Issend(sendBuffer2[rank], sentMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]);
863                        nbSendRequest++;
864                }
865                if (recvMessageSize[rank] > 0)
866                {
867                        MPI_Irecv(recvBuffer2[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]);
868                        nbRecvRequest++;
869                }
870        }
871
872        MPI_Waitall(nbRecvRequest, recvRequest, status);
873        MPI_Waitall(nbSendRequest, sendRequest, status);
874
875        delete [] sendRequest;
876        delete [] recvRequest;
877        delete [] status;
878        for (int rank = 0; rank < mpiSize; rank++)
879        {
880                if (nbRecvNode[rank] > 0)
881                {
882                        if (sentMessageSize[rank] > 0)
883                                delete [] sendBuffer2[rank];
884                }
885
886                if (recvMessageSize[rank] > 0)
887                {
888                        unpackIntersection(elements, recvBuffer2[rank]);
889                        delete [] recvBuffer2[rank];
890                }
891        }
892        delete [] sendBuffer2;
893        delete [] recvBuffer2;
894        delete [] sendBuffer;
895        delete [] recvBuffer;
896
897        delete [] nbSendNode;
898        delete [] nbRecvNode;
899        delete [] sentMessageSize;
900        delete [] recvMessageSize;
901}
902
903Mapper::~Mapper()
904{
905        delete [] remapMatrix;
906        delete [] srcAddress;
907        delete [] srcRank;
908        delete [] dstAddress;
909        if (neighbourElements) delete [] neighbourElements;
910}
911
912}
Note: See TracBrowser for help on using the repository browser.