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

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

dev_omp OK

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