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

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

tests (client, complete, toy, rempa) work fine with 2-level server mode. Tested on ADA (2*2+2)

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