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

Last change on this file since 1460 was 1460, checked in by yushan, 3 years ago

branch_openmp merged with XIOS_DEV_CMIP6@1459

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