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

Last change on this file since 1147 was 1147, checked in by yushan, 7 years ago

bug fixed in MPI_Gather(v)

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