source: XIOS/dev/branch_openmp/extern/remap/src/mapper.cpp @ 1328

Last change on this file since 1328 was 1328, checked in by yushan, 6 years ago

dev_omp

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