source: XIOS/trunk/extern/remap/src/mapper.cpp @ 1614

Last change on this file since 1614 was 1614, checked in by ymipsl, 3 years ago

Interpolation : enhancement : you have now the possibility to give computed area to the remaper for better global conservation.

YM

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