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