Changeset 2538 for XIOS3/trunk/extern/remap/src/mapper.cpp
- Timestamp:
- 07/24/23 11:55:12 (11 months ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
XIOS3/trunk/extern/remap/src/mapper.cpp
r2269 r2538 154 154 nIntersections++; 155 155 } 156 156 157 /* overallocate for NMAX neighbours for each elements */ 158 /* 157 159 remapMatrix = new double[nIntersections*NMAX]; 158 160 srcAddress = new int[nIntersections*NMAX]; … … 161 163 sourceWeightId =new long[nIntersections*NMAX]; 162 164 targetWeightId =new long[nIntersections*NMAX]; 163 165 */ 164 166 165 167 if (mpiRank == 0 && verbose) cout << "Remapping..." << endl; … … 195 197 196 198 int *nbSendElement = new int[mpiSize]; 199 int *nbSendVertex = new int[mpiSize]; 197 200 int **sendElement = new int*[mpiSize]; /* indices of elements required from other rank */ 198 201 double **recvValue = new double*[mpiSize]; 199 202 double **recvArea = new double*[mpiSize]; 200 203 double **recvGivenArea = new double*[mpiSize]; 204 int **recvVertexOffset = new int*[mpiSize]; 201 205 Coord **recvGrad = new Coord*[mpiSize]; 202 206 GloId **recvNeighIds = new GloId*[mpiSize]; /* ids of the of the source neighbours which also contribute through gradient */ … … 219 223 { 220 224 sendElement[rank] = new int[nbSendElement[rank]]; 221 recvValue[rank] = new double[nbSendElement[rank]];222 recvArea[rank] = new double[nbSendElement[rank]];223 recvGivenArea[rank] = new double[nbSendElement[rank]];224 if (order == 2)225 {226 recvNeighIds[rank] = new GloId[nbSendElement[rank]*(NMAX+1)];227 recvGrad[rank] = new Coord[nbSendElement[rank]*(NMAX+1)];228 }229 else230 recvNeighIds[rank] = new GloId[nbSendElement[rank]];231 232 225 last = -1; 233 226 index = -1; … … 244 237 /* communicate sizes of source elements to be sent (index lists and later values and gradients) */ 245 238 int *nbRecvElement = new int[mpiSize]; 239 int *nbRecvVertex = new int[mpiSize]; 246 240 MPI_Alltoall(nbSendElement, 1, MPI_INT, nbRecvElement, 1, MPI_INT, communicator); 247 241 … … 253 247 double **sendArea = new double*[mpiSize]; 254 248 double **sendGivenArea = new double*[mpiSize]; 249 int **sendVertexOffset = new int*[mpiSize]; 255 250 Coord **sendGrad = new Coord*[mpiSize]; 256 251 GloId **sendNeighIds = new GloId*[mpiSize]; 257 MPI_Request *sendRequest = new MPI_Request[ 5*mpiSize];258 MPI_Request *recvRequest = new MPI_Request[ 5*mpiSize];252 MPI_Request *sendRequest = new MPI_Request[6*mpiSize]; 253 MPI_Request *recvRequest = new MPI_Request[6*mpiSize]; 259 254 for (int rank = 0; rank < mpiSize; rank++) 260 255 { … … 268 263 { 269 264 recvElement[rank] = new int[nbRecvElement[rank]]; 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[6*mpiSize]; 270 271 MPI_Waitall(nbSendRequest, sendRequest, status); 272 MPI_Waitall(nbRecvRequest, recvRequest, status); 273 274 nbSendRequest = 0; 275 nbRecvRequest = 0; 276 277 for (int rank = 0; rank < mpiSize; rank++) 278 { 279 if (nbRecvElement[rank] > 0) 280 { 281 nbRecvVertex[rank]=0 ; 282 for (int j = 0; j < nbRecvElement[rank]; j++) nbRecvVertex[rank]+=sstree.localElements[recvElement[rank][j]].n+1; 283 MPI_Issend(&nbRecvVertex[rank], 1, MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 284 nbSendRequest++; 285 } 286 287 if (nbSendElement[rank] > 0) 288 { 289 MPI_Irecv(&nbSendVertex[rank], 1, MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 290 nbRecvRequest++; 291 } 292 } 293 294 MPI_Waitall(nbSendRequest, sendRequest, status); 295 MPI_Waitall(nbRecvRequest, recvRequest, status); 296 297 298 /* for all indices that have been received from requesting ranks: pack values and gradients, then send */ 299 nbSendRequest = 0; 300 nbRecvRequest = 0; 301 for (int rank = 0; rank < mpiSize; rank++) 302 { 303 if (nbRecvElement[rank] > 0) 304 { 270 305 sendValue[rank] = new double[nbRecvElement[rank]]; 271 306 sendArea[rank] = new double[nbRecvElement[rank]]; 272 307 sendGivenArea[rank] = new double[nbRecvElement[rank]]; 308 sendVertexOffset[rank] = new int[nbRecvElement[rank]]; 309 273 310 if (order == 2) 274 311 { 275 sendNeighIds[rank] = new GloId[nbRecvElement[rank]*(NMAX+1)]; 276 sendGrad[rank] = new Coord[nbRecvElement[rank]*(NMAX+1)]; 277 } 278 else 279 { 280 sendNeighIds[rank] = new GloId[nbRecvElement[rank]]; 281 } 282 MPI_Irecv(recvElement[rank], nbRecvElement[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 283 nbRecvRequest++; 284 } 285 } 286 MPI_Status *status = new MPI_Status[5*mpiSize]; 287 288 MPI_Waitall(nbSendRequest, sendRequest, status); 289 MPI_Waitall(nbRecvRequest, recvRequest, status); 290 291 /* for all indices that have been received from requesting ranks: pack values and gradients, then send */ 292 nbSendRequest = 0; 293 nbRecvRequest = 0; 294 for (int rank = 0; rank < mpiSize; rank++) 295 { 296 if (nbRecvElement[rank] > 0) 297 { 312 sendNeighIds[rank] = new GloId[nbRecvVertex[rank]]; 313 sendGrad[rank] = new Coord[nbRecvVertex[rank]]; 314 } 315 else sendNeighIds[rank] = new GloId[nbRecvElement[rank]]; 316 298 317 int jj = 0; // jj == j if no weight writing 318 sendVertexOffset[rank][0]=0 ; 299 319 for (int j = 0; j < nbRecvElement[rank]; j++) 300 320 { … … 302 322 sendArea[rank][j] = sstree.localElements[recvElement[rank][j]].area; 303 323 sendGivenArea[rank][j] = sstree.localElements[recvElement[rank][j]].given_area; 324 325 if (j==0) 326 { 327 if (order==2) sendVertexOffset[rank][j]=sstree.localElements[recvElement[rank][j]].n+1 ; 328 else sendVertexOffset[rank][j]=1 ; 329 } 330 else 331 { 332 if (order == 2) sendVertexOffset[rank][j] = sendVertexOffset[rank][j-1]+sstree.localElements[recvElement[rank][j]].n+1 ; 333 else sendVertexOffset[rank][j] = sendVertexOffset[rank][j-1]+1 ; 334 } 335 304 336 if (order == 2) 305 337 { … … 308 340 sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].src_id; 309 341 jj++; 310 for (int i = 0; i < NMAX; i++)342 for (int i = 0; i < sstree.localElements[recvElement[rank][j]].n ; i++) 311 343 { 312 344 sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].gradNeigh[i]; … … 325 357 MPI_Issend(sendGivenArea[rank], nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 326 358 nbSendRequest++; 359 MPI_Issend(sendVertexOffset[rank], nbRecvElement[rank], MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 360 nbSendRequest++; 361 327 362 if (order == 2) 328 363 { 329 MPI_Issend(sendGrad[rank], 3*nbRecv Element[rank]*(NMAX+1), MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]);364 MPI_Issend(sendGrad[rank], 3*nbRecvVertex[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 330 365 nbSendRequest++; 331 MPI_Issend(sendNeighIds[rank], 4*nbRecv Element[rank]*(NMAX+1), MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]);366 MPI_Issend(sendNeighIds[rank], 4*nbRecvVertex[rank], MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 332 367 //ym --> attention taille GloId 333 368 nbSendRequest++; … … 342 377 if (nbSendElement[rank] > 0) 343 378 { 379 recvValue[rank] = new double[nbSendElement[rank]]; 380 recvArea[rank] = new double[nbSendElement[rank]]; 381 recvGivenArea[rank] = new double[nbSendElement[rank]]; 382 recvVertexOffset[rank] = new int[nbSendElement[rank]]; 383 384 if (order == 2) 385 { 386 recvNeighIds[rank] = new GloId[nbSendVertex[rank]]; 387 recvGrad[rank] = new Coord[nbSendVertex[rank]]; 388 } 389 else recvNeighIds[rank] = new GloId[nbSendElement[rank]]; 390 391 344 392 MPI_Irecv(recvValue[rank], nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 345 393 nbRecvRequest++; … … 348 396 MPI_Irecv(recvGivenArea[rank], nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 349 397 nbRecvRequest++; 398 MPI_Irecv(recvVertexOffset[rank], nbSendElement[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 399 nbRecvRequest++; 400 350 401 if (order == 2) 351 402 { 352 MPI_Irecv(recvGrad[rank], 3*nbSendElement[rank]*(NMAX+1), 353 MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 403 MPI_Irecv(recvGrad[rank], 3*nbSendVertex[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 354 404 nbRecvRequest++; 355 MPI_Irecv(recvNeighIds[rank], 4*nbSend Element[rank]*(NMAX+1), MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]);405 MPI_Irecv(recvNeighIds[rank], 4*nbSendVertex[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 356 406 //ym --> attention taille GloId 357 407 nbRecvRequest++; … … 366 416 } 367 417 368 418 MPI_Waitall(nbSendRequest, sendRequest, status); 369 419 MPI_Waitall(nbRecvRequest, recvRequest, status); 370 420 … … 397 447 398 448 /* first order: src value times weight (weight = supermesh area), later divide by target area */ 399 int kk = (order == 2) ? n1 * (NMAX + 1) : n1; 400 GloId neighID = recvNeighIds[rank][kk]; 449 // int kk = (order == 2) ? n1 * (NMAX + 1) : n1; 450 int offset ; 451 int nvertex ; 452 if (n1==0) 453 { 454 offset=0 ; 455 nvertex=recvVertexOffset[rank][0] ; 456 } 457 else 458 { 459 offset = recvVertexOffset[rank][n1-1]; 460 nvertex = recvVertexOffset[rank][n1]-recvVertexOffset[rank][n1-1]; 461 } 462 463 GloId neighID = recvNeighIds[rank][offset]; 401 464 wgt_map[neighID] += w; 402 465 403 466 if (order == 2) 404 467 { 405 for (int k = 0; k < NMAX+1; k++)468 for (int k = 0; k < nvertex; k++) 406 469 { 407 int kk = n1 * (NMAX + 1)+ k;470 int kk = offset + k; 408 471 GloId neighID = recvNeighIds[rank][kk]; 409 472 if (neighID.ind != -1) wgt_map[neighID] += w * scalarprod(recvGrad[rank][kk], (*it)->x); … … 423 486 for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) 424 487 { 425 if (quantity) this->remapMatrix [i] = (it->second ) / renorm;426 else this->remapMatrix [i] = (it->second / e.area) / renorm;427 this->srcAddress [i] = it->first.ind;428 this->srcRank [i] = it->first.rank;429 this->dstAddress [i] = j;430 this->sourceWeightId [i]= it->first.globalId;431 this->targetWeightId [i]= targetGlobalId[j];488 if (quantity) this->remapMatrix.push_back((it->second ) / renorm); 489 else this->remapMatrix.push_back((it->second / e.area) / renorm); 490 this->srcAddress.push_back(it->first.ind); 491 this->srcRank.push_back(it->first.rank); 492 this->dstAddress.push_back(j); 493 this->sourceWeightId.push_back(it->first.globalId) ; 494 this->targetWeightId.push_back(targetGlobalId[j]) ; 432 495 i++; 433 496 } … … 443 506 delete[] recvArea[rank]; 444 507 delete[] recvGivenArea[rank]; 508 delete[] recvVertexOffset[rank]; 445 509 if (order == 2) 446 510 { … … 455 519 delete[] sendArea[rank]; 456 520 delete[] sendGivenArea[rank]; 521 delete[] sendVertexOffset[rank]; 457 522 if (order == 2) 458 523 delete[] sendGrad[rank]; … … 466 531 delete[] nbSendElement; 467 532 delete[] nbRecvElement; 533 delete[] nbSendVertex; 534 delete[] nbRecvVertex; 468 535 delete[] sendElement; 469 536 delete[] recvElement; … … 956 1023 Mapper::~Mapper() 957 1024 { 958 delete [] remapMatrix;959 delete [] srcAddress;960 delete [] srcRank;961 delete [] dstAddress;962 delete [] sourceWeightId ;963 delete [] targetWeightId ;964 if (neighbourElements) delete [] neighbourElements;1025 //delete [] remapMatrix; 1026 //delete [] srcAddress; 1027 //delete [] srcRank; 1028 //delete [] dstAddress; 1029 //delete [] sourceWeightId ; 1030 //delete [] targetWeightId ; 1031 //if (neighbourElements) delete [] neighbourElements; 965 1032 CTimer::release() ; 966 1033 }
Note: See TracChangeset
for help on using the changeset viewer.