source: XIOS/dev/branch_yushan_merged/extern/remap/src/mapper.cpp @ 1134

Last change on this file since 1134 was 1134, checked in by yushan, 4 years ago

branch merged with trunk r1130

File size: 29.0 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, bool quantity)
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, quantity);
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, bool quantity)
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        double **recvArea = new double*[mpiSize];
187        Coord **recvGrad = new Coord*[mpiSize];
188        GloId **recvNeighIds = new GloId*[mpiSize]; /* ids of the of the source neighbours which also contribute through gradient */
189        for (int rank = 0; rank < mpiSize; rank++)
190        {
191                /* get size for allocation */
192                int last = -1; /* compares unequal to any index */
193                int index = -1; /* increased to starting index 0 in first iteration */
194                for (multimap<int, Polyg *>::iterator it = elementList[rank].begin(); it != elementList[rank].end(); ++it)
195                {
196                        if (last != it->first)
197                                index++;
198                        (it->second)->id.ind = index;
199                        last = it->first;
200                }
201                nbSendElement[rank] = index + 1;
202
203                /* if size is non-zero allocate and collect indices of elements on other ranks that we intersect */
204                if (nbSendElement[rank] > 0)
205                {
206                        sendElement[rank] = new int[nbSendElement[rank]];
207                        recvValue[rank]   = new double[nbSendElement[rank]];
208                        recvArea[rank]    = new double[nbSendElement[rank]];
209                        if (order == 2)
210                        {
211                                recvNeighIds[rank] = new GloId[nbSendElement[rank]*(NMAX+1)];
212                                recvGrad[rank]    = new Coord[nbSendElement[rank]*(NMAX+1)];
213                        }
214                        else
215                                recvNeighIds[rank] = new GloId[nbSendElement[rank]];
216
217                        last = -1;
218                        index = -1;
219                        for (multimap<int, Polyg *>::iterator it = elementList[rank].begin(); it != elementList[rank].end(); ++it)
220                        {
221                                if (last != it->first)
222                                        index++;
223                                sendElement[rank][index] = it->first;
224                                last = it->first;
225                        }
226                }
227        }
228
229        /* communicate sizes of source elements to be sent (index lists and later values and gradients) */
230        int *nbRecvElement = new int[mpiSize];
231        MPI_Alltoall(nbSendElement, 1, MPI_INT, nbRecvElement, 1, MPI_INT, communicator);
232
233        /* communicate indices of source elements on other ranks whoes value and gradient we need (since intersection) */
234        int nbSendRequest = 0;
235        int nbRecvRequest = 0;
236        int **recvElement = new int*[mpiSize];
237        double **sendValue = new double*[mpiSize];
238        double **sendArea = new double*[mpiSize];
239        Coord **sendGrad = new Coord*[mpiSize];
240        GloId **sendNeighIds = new GloId*[mpiSize];
241        MPI_Request *sendRequest = new MPI_Request[3*mpiSize];
242        MPI_Request *recvRequest = new MPI_Request[3*mpiSize];
243        for (int rank = 0; rank < mpiSize; rank++)
244        {
245                if (nbSendElement[rank] > 0)
246                {
247                        MPI_Issend(sendElement[rank], nbSendElement[rank], MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]);
248                        nbSendRequest++;
249                }
250
251                if (nbRecvElement[rank] > 0)
252                {
253                        recvElement[rank] = new int[nbRecvElement[rank]];
254                        sendValue[rank]   = new double[nbRecvElement[rank]];
255                        sendArea[rank]   = new double[nbRecvElement[rank]];
256                        if (order == 2)
257                        {
258                                sendNeighIds[rank] = new GloId[nbRecvElement[rank]*(NMAX+1)];
259                                sendGrad[rank]    = new Coord[nbRecvElement[rank]*(NMAX+1)];
260                        }
261                        else
262                        {
263                                sendNeighIds[rank] = new GloId[nbRecvElement[rank]];
264                        }
265                        MPI_Irecv(recvElement[rank], nbRecvElement[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]);
266                        nbRecvRequest++;
267                }
268        }
269        MPI_Status *status = new MPI_Status[3*mpiSize];
270        MPI_Waitall(nbRecvRequest, recvRequest, status);
271        MPI_Waitall(nbSendRequest, sendRequest, status);
272
273        /* for all indices that have been received from requesting ranks: pack values and gradients, then send */
274        nbSendRequest = 0;
275        nbRecvRequest = 0;
276        for (int rank = 0; rank < mpiSize; rank++)
277        {
278                if (nbRecvElement[rank] > 0)
279                {
280                        int jj = 0; // jj == j if no weight writing
281                        for (int j = 0; j < nbRecvElement[rank]; j++)
282                        {
283                                sendValue[rank][j] = sstree.localElements[recvElement[rank][j]].val;
284                                sendArea[rank][j] = sstree.localElements[recvElement[rank][j]].area;
285                                if (order == 2)
286                                {
287                                        sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].grad;
288//          cout<<"grad  "<<jj<<"  "<<recvElement[rank][j]<<"  "<<sendGrad[rank][jj]<<" "<<sstree.localElements[recvElement[rank][j]].grad<<endl ;
289                                        sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].src_id;
290                                        jj++;
291                                        for (int i = 0; i < NMAX; i++)
292                                        {
293                                                sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].gradNeigh[i];
294//            cout<<"grad  "<<jj<<"  "<<sendGrad[rank][jj]<<" "<<sstree.localElements[recvElement[rank][j]].grad<<endl ;
295            sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].neighId[i];
296                                                jj++;
297                                        }
298                                }
299                                else
300                                        sendNeighIds[rank][j] = sstree.localElements[recvElement[rank][j]].src_id;
301                        }
302                        MPI_Issend(sendValue[rank],  nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]);
303                        nbSendRequest++;
304                        MPI_Issend(sendArea[rank],  nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]);
305                        nbSendRequest++;
306                        if (order == 2)
307                        {
308                                MPI_Issend(sendGrad[rank], 3*nbRecvElement[rank]*(NMAX+1),
309                                                                MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]);
310                                nbSendRequest++;
311                                MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank]*(NMAX+1), MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]);
312//ym  --> attention taille GloId
313                                nbSendRequest++;
314                        }
315                        else
316                        {
317                                MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank], MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]);
318//ym  --> attention taille GloId
319                                nbSendRequest++;
320                        }
321                }
322                if (nbSendElement[rank] > 0)
323                {
324                        MPI_Irecv(recvValue[rank],  nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]);
325                        nbRecvRequest++;
326                        MPI_Irecv(recvArea[rank],  nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]);
327                        nbRecvRequest++;
328                        if (order == 2)
329                        {
330                                MPI_Irecv(recvGrad[rank], 3*nbSendElement[rank]*(NMAX+1),
331                                                MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]);
332                                nbRecvRequest++;
333                                MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank]*(NMAX+1), MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]);
334//ym  --> attention taille GloId
335                                nbRecvRequest++;
336                        }
337                        else
338                        {
339                                MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]);
340//ym  --> attention taille GloId
341                                nbRecvRequest++;
342                        }
343                }
344        }
345        MPI_Waitall(nbRecvRequest, recvRequest, status);
346        MPI_Waitall(nbSendRequest, sendRequest, status);
347
348        /* now that all values and gradients are available use them to computed interpolated values on target
349           and also to compute weights */
350        int i = 0;
351        for (int j = 0; j < nbElements; j++)
352        {
353                Elt& e = elements[j];
354
355                /* since for the 2nd order case source grid elements can contribute to a destination grid element over several "paths"
356                   (step1: gradient is computed using neighbours on same grid, step2: intersection uses several elements on other grid)
357                   accumulate them so that there is only one final weight between two elements */
358                map<GloId,double> wgt_map;
359
360                /* for destination element `e` loop over all intersetions/the corresponding source elements */
361                for (list<Polyg *>::iterator it = e.is.begin(); it != e.is.end(); it++)
362                {
363                        /* it is the intersection element, so it->x and it->area are barycentre and area of intersection element (super mesh)
364                        but it->id is id of the source element that it intersects */
365                        int n1 = (*it)->id.ind;
366                        int rank = (*it)->id.rank;
367                        double fk = recvValue[rank][n1];
368                        double srcArea = recvArea[rank][n1];
369                        double w = (*it)->area;
370                        if (quantity) w/=srcArea ;
371
372                        /* first order: src value times weight (weight = supermesh area), later divide by target area */
373                        int kk = (order == 2) ? n1 * (NMAX + 1) : n1;
374                        GloId neighID = recvNeighIds[rank][kk];
375                        wgt_map[neighID] += w;
376
377                        if (order == 2)
378                        {
379                                for (int k = 0; k < NMAX+1; k++)
380                                {
381                                        int kk = n1 * (NMAX + 1) + k;
382                                        GloId neighID = recvNeighIds[rank][kk];
383                                        if (neighID.ind != -1)  wgt_map[neighID] += w * scalarprod(recvGrad[rank][kk], (*it)->x);
384                                }
385
386                        }
387                }
388
389    double renorm=0;
390    if (renormalize) 
391      for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) renorm+=it->second / e.area;
392    else renorm=1. ;
393
394    for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++)
395                {
396      if (quantity)  this->remapMatrix[i] = (it->second ) / renorm;
397                        else this->remapMatrix[i] = (it->second / e.area) / renorm;
398                        this->srcAddress[i] = it->first.ind;
399                        this->srcRank[i] = it->first.rank;
400                        this->dstAddress[i] = j;
401      this->sourceWeightId[i]= it->first.globalId ;
402      this->targetWeightId[i]= targetGlobalId[j] ;
403                        i++;
404                }
405        }
406
407        /* free all memory allocated in this function */
408        for (int rank = 0; rank < mpiSize; rank++)
409        {
410                if (nbSendElement[rank] > 0)
411                {
412                        delete[] sendElement[rank];
413                        delete[] recvValue[rank];
414                        delete[] recvArea[rank];
415                        if (order == 2)
416                        {
417                                delete[] recvGrad[rank];
418                        }
419                        delete[] recvNeighIds[rank];
420                }
421                if (nbRecvElement[rank] > 0)
422                {
423                        delete[] recvElement[rank];
424                        delete[] sendValue[rank];
425                        delete[] sendArea[rank];
426                        if (order == 2)
427                                delete[] sendGrad[rank];
428                        delete[] sendNeighIds[rank];
429                }
430        }
431        delete[] status;
432        delete[] sendRequest;
433        delete[] recvRequest;
434        delete[] elementList;
435        delete[] nbSendElement;
436        delete[] nbRecvElement;
437        delete[] sendElement;
438        delete[] recvElement;
439        delete[] sendValue;
440        delete[] recvValue;
441        delete[] sendGrad;
442        delete[] recvGrad;
443        delete[] sendNeighIds;
444        delete[] recvNeighIds;
445        return i;
446}
447
448void Mapper::computeGrads()
449{
450        /* array of pointers to collect local elements and elements received from other cpu */
451        vector<Elt*> globalElements(sstree.nbLocalElements + nbNeighbourElements);
452        int index = 0;
453        for (int i = 0; i < sstree.nbLocalElements; i++, index++)
454                globalElements[index] = &(sstree.localElements[i]);
455        for (int i = 0; i < nbNeighbourElements; i++, index++)
456                globalElements[index] = &neighbourElements[i];
457
458        update_baryc(sstree.localElements, sstree.nbLocalElements);
459        computeGradients(&globalElements[0], sstree.nbLocalElements);
460}
461
462/** for each element of the source grid, finds all the neighbouring elements that share an edge
463    (filling array neighbourElements). This is used later to compute gradients */
464void Mapper::buildMeshTopology()
465{
466        int mpiSize, mpiRank;
467        MPI_Comm_size(communicator, &mpiSize);
468        MPI_Comm_rank(communicator, &mpiRank);
469
470        vector<Node> *routingList = new vector<Node>[mpiSize];
471        vector<vector<int> > routes(sstree.localTree.leafs.size());
472
473        sstree.routeIntersections(routes, sstree.localTree.leafs);
474
475        for (int i = 0; i < routes.size(); ++i)
476                for (int k = 0; k < routes[i].size(); ++k)
477                        routingList[routes[i][k]].push_back(sstree.localTree.leafs[i]);
478        routingList[mpiRank].clear();
479
480
481        CMPIRouting mpiRoute(communicator);
482        mpiRoute.init(routes);
483        int nRecv = mpiRoute.getTotalSourceElement();
484// cout << mpiRank << " NRECV " << nRecv << "(" << routes.size() << ")"<< endl;
485
486        int *nbSendNode = new int[mpiSize];
487        int *nbRecvNode = new int[mpiSize];
488        int *sendMessageSize = new int[mpiSize];
489        int *recvMessageSize = new int[mpiSize];
490
491        for (int rank = 0; rank < mpiSize; rank++)
492        {
493                nbSendNode[rank] = routingList[rank].size();
494                sendMessageSize[rank] = 0;
495                for (size_t j = 0; j < routingList[rank].size(); j++)
496                {
497                        Elt *elt = (Elt *) (routingList[rank][j].data);
498                        sendMessageSize[rank] += packedPolygonSize(*elt);
499                }
500        }
501
502        MPI_Alltoall(nbSendNode, 1, MPI_INT, nbRecvNode, 1, MPI_INT, communicator);
503        MPI_Alltoall(sendMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator);
504
505        char **sendBuffer = new char*[mpiSize];
506        char **recvBuffer = new char*[mpiSize];
507        int *pos = new int[mpiSize];
508
509        for (int rank = 0; rank < mpiSize; rank++)
510        {
511                if (nbSendNode[rank] > 0) sendBuffer[rank] = new char[sendMessageSize[rank]];
512                if (nbRecvNode[rank] > 0) recvBuffer[rank] = new char[recvMessageSize[rank]];
513        }
514
515        for (int rank = 0; rank < mpiSize; rank++)
516        {
517                pos[rank] = 0;
518                for (size_t j = 0; j < routingList[rank].size(); j++)
519                {
520                        Elt *elt = (Elt *) (routingList[rank][j].data);
521                        packPolygon(*elt, sendBuffer[rank], pos[rank]);
522                }
523        }
524        delete [] routingList;
525
526
527        int nbSendRequest = 0;
528        int nbRecvRequest = 0;
529        MPI_Request *sendRequest = new MPI_Request[mpiSize];
530        MPI_Request *recvRequest = new MPI_Request[mpiSize];
531        MPI_Status  *status      = new MPI_Status[mpiSize];
532
533        for (int rank = 0; rank < mpiSize; rank++)
534        {
535                if (nbSendNode[rank] > 0)
536                {
537                        MPI_Issend(sendBuffer[rank], sendMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]);
538                        nbSendRequest++;
539                }
540                if (nbRecvNode[rank] > 0)
541                {
542                        MPI_Irecv(recvBuffer[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]);
543                        nbRecvRequest++;
544                }
545        }
546
547        MPI_Waitall(nbRecvRequest, recvRequest, status);
548        MPI_Waitall(nbSendRequest, sendRequest, status);
549
550        for (int rank = 0; rank < mpiSize; rank++)
551                if (nbSendNode[rank] > 0) delete [] sendBuffer[rank];
552        delete [] sendBuffer;
553
554        char **sendBuffer2 = new char*[mpiSize];
555        char **recvBuffer2 = new char*[mpiSize];
556
557        for (int rank = 0; rank < mpiSize; rank++)
558        {
559                nbSendNode[rank] = 0;
560                sendMessageSize[rank] = 0;
561
562                if (nbRecvNode[rank] > 0)
563                {
564                        set<NodePtr> neighbourList;
565                        pos[rank] = 0;
566                        for (int j = 0; j < nbRecvNode[rank]; j++)
567                        {
568                                Elt elt;
569                                unpackPolygon(elt, recvBuffer[rank], pos[rank]);
570                                Node node(elt.x, cptRadius(elt), &elt);
571                                findNeighbour(sstree.localTree.root, &node, neighbourList);
572                        }
573                        nbSendNode[rank] = neighbourList.size();
574                        for (set<NodePtr>::iterator it = neighbourList.begin(); it != neighbourList.end(); it++)
575                        {
576                                Elt *elt = (Elt *) ((*it)->data);
577                                sendMessageSize[rank] += packedPolygonSize(*elt);
578                        }
579
580                        sendBuffer2[rank] = new char[sendMessageSize[rank]];
581                        pos[rank] = 0;
582
583                        for (set<NodePtr>::iterator it = neighbourList.begin(); it != neighbourList.end(); it++)
584                        {
585                                Elt *elt = (Elt *) ((*it)->data);
586                                packPolygon(*elt, sendBuffer2[rank], pos[rank]);
587                        }
588                }
589        }
590        for (int rank = 0; rank < mpiSize; rank++)
591                if (nbRecvNode[rank] > 0) delete [] recvBuffer[rank];
592        delete [] recvBuffer;
593
594
595        MPI_Barrier(communicator);
596        MPI_Alltoall(nbSendNode, 1, MPI_INT, nbRecvNode, 1, MPI_INT, communicator);
597        MPI_Alltoall(sendMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator);
598
599        for (int rank = 0; rank < mpiSize; rank++)
600                if (nbRecvNode[rank] > 0) recvBuffer2[rank] = new char[recvMessageSize[rank]];
601
602        nbSendRequest = 0;
603        nbRecvRequest = 0;
604
605        for (int rank = 0; rank < mpiSize; rank++)
606        {
607                if (nbSendNode[rank] > 0)
608                {
609                        MPI_Issend(sendBuffer2[rank], sendMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]);
610                        nbSendRequest++;
611                }
612                if (nbRecvNode[rank] > 0)
613                {
614                        MPI_Irecv(recvBuffer2[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]);
615                        nbRecvRequest++;
616                }
617        }
618
619        MPI_Waitall(nbRecvRequest, recvRequest, status);
620        MPI_Waitall(nbSendRequest, sendRequest, status);
621
622        int nbNeighbourNodes = 0;
623        for (int rank = 0; rank < mpiSize; rank++)
624                nbNeighbourNodes += nbRecvNode[rank];
625
626        neighbourElements = new Elt[nbNeighbourNodes];
627        nbNeighbourElements = nbNeighbourNodes;
628
629        int index = 0;
630        for (int rank = 0; rank < mpiSize; rank++)
631        {
632                pos[rank] = 0;
633                for (int j = 0; j < nbRecvNode[rank]; j++)
634                {
635                        unpackPolygon(neighbourElements[index], recvBuffer2[rank], pos[rank]);
636                        neighbourElements[index].id.ind = sstree.localTree.leafs.size() + index;
637                        index++;
638                }
639        }
640        for (int rank = 0; rank < mpiSize; rank++)
641        {
642                if (nbRecvNode[rank] > 0) delete [] recvBuffer2[rank];
643                if (nbSendNode[rank] > 0) delete [] sendBuffer2[rank];
644        }
645        delete [] recvBuffer2;
646        delete [] sendBuffer2;
647        delete [] sendMessageSize;
648        delete [] recvMessageSize;
649        delete [] nbSendNode;
650        delete [] nbRecvNode;
651        delete [] sendRequest;
652        delete [] recvRequest;
653        delete [] status;
654        delete [] pos;
655
656        /* re-compute on received elements to avoid having to send this information */
657        neighbourNodes.resize(nbNeighbourNodes);
658        setCirclesAndLinks(neighbourElements, neighbourNodes);
659        cptAllEltsGeom(neighbourElements, nbNeighbourNodes, srcGrid.pole);
660
661        /* the local SS tree must include nodes from other cpus if they are potential
662           intersector of a local node */
663        sstree.localTree.insertNodes(neighbourNodes);
664
665        /* for every local element,
666           use the SS-tree to find all elements (including neighbourElements)
667           who are potential neighbours because their circles intersect,
668           then check all canditates for common edges to build up connectivity information
669        */
670        for (int j = 0; j < sstree.localTree.leafs.size(); j++)
671        {
672                Node& node = sstree.localTree.leafs[j];
673
674                /* find all leafs whoes circles that intersect node's circle and save into node->intersectors */
675                node.search(sstree.localTree.root);
676
677                Elt *elt = (Elt *)(node.data);
678
679                for (int i = 0; i < elt->n; i++) elt->neighbour[i] = NOT_FOUND;
680
681                /* for element `elt` loop through all nodes in the SS-tree
682                   whoes circles intersect with the circle around `elt` (the SS intersectors)
683                   and check if they are neighbours in the sense that the two elements share an edge.
684                   If they do, save this information for elt */
685                for (list<NodePtr>::iterator it = (node.intersectors).begin(); it != (node.intersectors).end(); ++it)
686                {
687                        Elt *elt2 = (Elt *)((*it)->data);
688                        set_neighbour(*elt, *elt2);
689                }
690
691/*
692                for (int i = 0; i < elt->n; i++)
693                {
694                        if (elt->neighbour[i] == NOT_FOUND)
695                                error_exit("neighbour not found");
696                }
697*/
698        }
699}
700
701/** @param elements are the target grid elements */
702void Mapper::computeIntersection(Elt *elements, int nbElements)
703{
704        int mpiSize, mpiRank;
705        MPI_Comm_size(communicator, &mpiSize);
706        MPI_Comm_rank(communicator, &mpiRank);
707
708        MPI_Barrier(communicator);
709
710        vector<Node> *routingList = new vector<Node>[mpiSize];
711
712        vector<Node> routeNodes;  routeNodes.reserve(nbElements);
713        for (int j = 0; j < nbElements; j++)
714        {
715                elements[j].id.ind = j;
716                elements[j].id.rank = mpiRank;
717                routeNodes.push_back(Node(elements[j].x, cptRadius(elements[j]), &elements[j]));
718        }
719
720        vector<vector<int> > routes(routeNodes.size());
721        sstree.routeIntersections(routes, routeNodes);
722        for (int i = 0; i < routes.size(); ++i)
723                for (int k = 0; k < routes[i].size(); ++k)
724                        routingList[routes[i][k]].push_back(routeNodes[i]);
725
726        if (verbose >= 2)
727        {
728                cout << " --> rank  " << mpiRank << " nbElements " << nbElements << " : ";
729                for (int rank = 0; rank < mpiSize; rank++)
730                        cout << routingList[rank].size() << "   ";
731                cout << endl;
732        }
733        MPI_Barrier(communicator);
734
735        int *nbSendNode = new int[mpiSize];
736        int *nbRecvNode = new int[mpiSize];
737        int *sentMessageSize = new int[mpiSize];
738        int *recvMessageSize = new int[mpiSize];
739
740        for (int rank = 0; rank < mpiSize; rank++)
741        {
742                nbSendNode[rank] = routingList[rank].size();
743                sentMessageSize[rank] = 0;
744                for (size_t j = 0; j < routingList[rank].size(); j++)
745                {
746                        Elt *elt = (Elt *) (routingList[rank][j].data);
747                        sentMessageSize[rank] += packedPolygonSize(*elt);
748                }
749        }
750
751        MPI_Alltoall(nbSendNode, 1, MPI_INT, nbRecvNode, 1, MPI_INT, communicator);
752        MPI_Alltoall(sentMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator);
753
754        int total = 0;
755
756        for (int rank = 0; rank < mpiSize; rank++)
757        {
758                total = total + nbRecvNode[rank];
759        }
760
761        if (verbose >= 2) cout << "---> rank " << mpiRank << " : compute intersection : total received nodes  " << total << endl;
762        char **sendBuffer = new char*[mpiSize];
763        char **recvBuffer = new char*[mpiSize];
764        int *pos = new int[mpiSize];
765
766        for (int rank = 0; rank < mpiSize; rank++)
767        {
768                if (nbSendNode[rank] > 0) sendBuffer[rank] = new char[sentMessageSize[rank]];
769                if (nbRecvNode[rank] > 0) recvBuffer[rank] = new char[recvMessageSize[rank]];
770        }
771
772        for (int rank = 0; rank < mpiSize; rank++)
773        {
774                pos[rank] = 0;
775                for (size_t j = 0; j < routingList[rank].size(); j++)
776                {
777                        Elt* elt = (Elt *) (routingList[rank][j].data);
778                        packPolygon(*elt, sendBuffer[rank], pos[rank]);
779                }
780        }
781        delete [] routingList;
782
783        int nbSendRequest = 0;
784        int nbRecvRequest = 0;
785        MPI_Request *sendRequest = new MPI_Request[mpiSize];
786        MPI_Request *recvRequest = new MPI_Request[mpiSize];
787        MPI_Status   *status = new MPI_Status[mpiSize];
788
789        for (int rank = 0; rank < mpiSize; rank++)
790        {
791                if (nbSendNode[rank] > 0)
792                {
793                        MPI_Issend(sendBuffer[rank], sentMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]);
794                        nbSendRequest++;
795                }
796                if (nbRecvNode[rank] > 0)
797                {
798                        MPI_Irecv(recvBuffer[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]);
799                        nbRecvRequest++;
800                }
801        }
802
803        MPI_Waitall(nbRecvRequest, recvRequest, status);
804        MPI_Waitall(nbSendRequest, sendRequest, status);
805        char **sendBuffer2 = new char*[mpiSize];
806        char **recvBuffer2 = new char*[mpiSize];
807
808        double tic = cputime();
809        for (int rank = 0; rank < mpiSize; rank++)
810        {
811                sentMessageSize[rank] = 0;
812
813                if (nbRecvNode[rank] > 0)
814                {
815                        Elt *recvElt = new Elt[nbRecvNode[rank]];
816                        pos[rank] = 0;
817                        for (int j = 0; j < nbRecvNode[rank]; j++)
818                        {
819                                unpackPolygon(recvElt[j], recvBuffer[rank], pos[rank]);
820                                cptEltGeom(recvElt[j], tgtGrid.pole);
821                                Node recvNode(recvElt[j].x, cptRadius(recvElt[j]), &recvElt[j]);
822                                recvNode.search(sstree.localTree.root);
823                                /* for a node holding an element of the target, loop throught candidates for intersecting source */
824                                for (list<NodePtr>::iterator it = (recvNode.intersectors).begin(); it != (recvNode.intersectors).end(); ++it)
825                                {
826                                        Elt *elt2 = (Elt *) ((*it)->data);
827                                        /* recvElt is target, elt2 is source */
828//                                      intersect(&recvElt[j], elt2);
829                                        intersect_ym(&recvElt[j], elt2);
830                                }
831
832                                if (recvElt[j].is.size() > 0) sentMessageSize[rank] += packIntersectionSize(recvElt[j]);
833
834                                // here recvNode goes out of scope
835                        }
836
837                        if (sentMessageSize[rank] > 0)
838                        {
839                                sentMessageSize[rank] += sizeof(int);
840                                sendBuffer2[rank] = new char[sentMessageSize[rank]];
841                                *((int *) sendBuffer2[rank]) = 0;
842                                pos[rank] = sizeof(int);
843                                for (int j = 0; j < nbRecvNode[rank]; j++)
844                                {
845                                        packIntersection(recvElt[j], sendBuffer2[rank], pos[rank]);
846                                        //FIXME should be deleted: recvElt[j].delete_intersections(); // intersection areas have been packed to buffer and won't be used any more
847                                }
848                        }
849                        delete [] recvElt;
850
851                }
852        }
853        delete [] pos;
854
855        for (int rank = 0; rank < mpiSize; rank++)
856        {
857                if (nbSendNode[rank] > 0) delete [] sendBuffer[rank];
858                if (nbRecvNode[rank] > 0) delete [] recvBuffer[rank];
859                nbSendNode[rank] = 0;
860        }
861
862        if (verbose >= 2) cout << "Rank " << mpiRank << "  Compute (internal) intersection " << cputime() - tic << " s" << endl;
863        MPI_Alltoall(sentMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator);
864
865        for (int rank = 0; rank < mpiSize; rank++)
866                if (recvMessageSize[rank] > 0)
867                        recvBuffer2[rank] = new char[recvMessageSize[rank]];
868
869        nbSendRequest = 0;
870        nbRecvRequest = 0;
871
872        for (int rank = 0; rank < mpiSize; rank++)
873        {
874                if (sentMessageSize[rank] > 0)
875                {
876                        MPI_Issend(sendBuffer2[rank], sentMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]);
877                        nbSendRequest++;
878                }
879                if (recvMessageSize[rank] > 0)
880                {
881                        MPI_Irecv(recvBuffer2[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]);
882                        nbRecvRequest++;
883                }
884        }
885
886        MPI_Waitall(nbRecvRequest, recvRequest, status);
887        MPI_Waitall(nbSendRequest, sendRequest, status);
888
889        delete [] sendRequest;
890        delete [] recvRequest;
891        delete [] status;
892        for (int rank = 0; rank < mpiSize; rank++)
893        {
894                if (nbRecvNode[rank] > 0)
895                {
896                        if (sentMessageSize[rank] > 0)
897                                delete [] sendBuffer2[rank];
898                }
899
900                if (recvMessageSize[rank] > 0)
901                {
902                        unpackIntersection(elements, recvBuffer2[rank]);
903                        delete [] recvBuffer2[rank];
904                }
905        }
906        delete [] sendBuffer2;
907        delete [] recvBuffer2;
908        delete [] sendBuffer;
909        delete [] recvBuffer;
910
911        delete [] nbSendNode;
912        delete [] nbRecvNode;
913        delete [] sentMessageSize;
914        delete [] recvMessageSize;
915}
916
917Mapper::~Mapper()
918{
919        delete [] remapMatrix;
920        delete [] srcAddress;
921        delete [] srcRank;
922        delete [] dstAddress;
923        if (neighbourElements) delete [] neighbourElements;
924}
925
926}
Note: See TracBrowser for help on using the repository browser.