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

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

save modif

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