source: XIOS/dev/dev_trunk_omp/extern/remap/src/mapper.cpp @ 1619

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

branch_openmp merged and NOT tested with trunk r1618

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