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

Last change on this file since 1646 was 1646, checked in by yushan, 5 years ago

branch merged with trunk @1645. arch file (ep&mpi) added for ADA

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