Changeset 1338 for XIOS/dev/branch_openmp/extern/remap
 Timestamp:
 11/21/17 10:47:57 (4 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

XIOS/dev/branch_openmp/extern/remap/src/mapper.cpp
r1335 r1338 132 132 timings.push_back(cputime()  tic); 133 133 134 tic = cputime(); 135 if (interpOrder == 2) { 136 if (mpiRank == 0 && verbose) cout << "Computing grads ..." << endl; 137 buildMeshTopology(); 138 computeGrads(); 139 } 140 timings.push_back(cputime()  tic); 141 142 /* Prepare computation of weights */ 143 /* compute number of intersections which for the first order case 144 corresponds to the number of edges in the remap matrix */ 145 int nIntersections = 0; 146 for (int j = 0; j < targetElements.size(); j++) 147 { 148 Elt &elt = targetElements[j]; 149 for (list<Polyg*>::iterator it = elt.is.begin(); it != elt.is.end(); it++) 150 nIntersections++; 151 } 152 /* overallocate for NMAX neighbours for each elements */ 153 remapMatrix = new double[nIntersections*NMAX]; 154 srcAddress = new int[nIntersections*NMAX]; 155 srcRank = new int[nIntersections*NMAX]; 156 dstAddress = new int[nIntersections*NMAX]; 134 tic = cputime(); 135 if (interpOrder == 2) 136 { 137 if (mpiRank == 0 && verbose) cout << "Computing grads ..." << endl; 138 buildMeshTopology(); 139 computeGrads(); 140 } 141 timings.push_back(cputime()  tic); 142 143 /* Prepare computation of weights */ 144 /* compute number of intersections which for the first order case 145 corresponds to the number of edges in the remap matrix */ 146 int nIntersections = 0; 147 for (int j = 0; j < targetElements.size(); j++) 148 { 149 Elt &elt = targetElements[j]; 150 for (list<Polyg*>::iterator it = elt.is.begin(); it != elt.is.end(); it++) 151 nIntersections++; 152 } 153 /* overallocate for NMAX neighbours for each elements */ 154 remapMatrix = new double[nIntersections*NMAX]; 155 srcAddress = new int[nIntersections*NMAX]; 156 srcRank = new int[nIntersections*NMAX]; 157 dstAddress = new int[nIntersections*NMAX]; 157 158 sourceWeightId =new long[nIntersections*NMAX]; 158 159 targetWeightId =new long[nIntersections*NMAX]; 159 160 160 161 161 162 163 164 162 if (mpiRank == 0 && verbose) cout << "Remapping..." << endl; 163 tic = cputime(); 164 nWeights = remap(&targetElements[0], targetElements.size(), interpOrder, renormalize, quantity); 165 timings.push_back(cputime()  tic); 165 166 166 167 for (int i = 0; i < targetElements.size(); i++) targetElements[i].delete_intersections(); 167 168 168 169 return timings; 169 170 } 170 171 … … 177 178 int Mapper::remap(Elt *elements, int nbElements, int order, bool renormalize, bool quantity) 178 179 { 179 int mpiSize, mpiRank; 180 MPI_Comm_size(communicator, &mpiSize); 181 MPI_Comm_rank(communicator, &mpiRank); 182 183 /* create list of intersections (super mesh elements) for each rank */ 184 multimap<int, Polyg *> *elementList = new multimap<int, Polyg *>[mpiSize]; 185 for (int j = 0; j < nbElements; j++) 186 { 187 Elt& e = elements[j]; 188 for (list<Polyg *>::iterator it = e.is.begin(); it != e.is.end(); it++) 189 elementList[(*it)>id.rank].insert(pair<int, Polyg *>((*it)>id.ind, *it)); 190 } 191 192 int *nbSendElement = new int[mpiSize]; 193 int **sendElement = new int*[mpiSize]; /* indices of elements required from other rank */ 194 double **recvValue = new double*[mpiSize]; 195 double **recvArea = new double*[mpiSize]; 196 Coord **recvGrad = new Coord*[mpiSize]; 197 GloId **recvNeighIds = new GloId*[mpiSize]; /* ids of the of the source neighbours which also contribute through gradient */ 198 for (int rank = 0; rank < mpiSize; rank++) 199 { 200 /* get size for allocation */ 201 int last = 1; /* compares unequal to any index */ 202 int index = 1; /* increased to starting index 0 in first iteration */ 203 for (multimap<int, Polyg *>::iterator it = elementList[rank].begin(); it != elementList[rank].end(); ++it) 204 { 205 if (last != it>first) 206 index++; 207 (it>second)>id.ind = index; 208 last = it>first; 209 } 210 nbSendElement[rank] = index + 1; 211 212 /* if size is nonzero allocate and collect indices of elements on other ranks that we intersect */ 213 if (nbSendElement[rank] > 0) 214 { 215 sendElement[rank] = new int[nbSendElement[rank]]; 216 recvValue[rank] = new double[nbSendElement[rank]]; 217 recvArea[rank] = new double[nbSendElement[rank]]; 218 if (order == 2) 219 { 220 recvNeighIds[rank] = new GloId[nbSendElement[rank]*(NMAX+1)]; 221 recvGrad[rank] = new Coord[nbSendElement[rank]*(NMAX+1)]; 222 } 223 else 224 recvNeighIds[rank] = new GloId[nbSendElement[rank]]; 225 226 last = 1; 227 index = 1; 228 for (multimap<int, Polyg *>::iterator it = elementList[rank].begin(); it != elementList[rank].end(); ++it) 229 { 230 if (last != it>first) 231 index++; 232 sendElement[rank][index] = it>first; 233 last = it>first; 234 } 235 } 236 } 237 238 /* communicate sizes of source elements to be sent (index lists and later values and gradients) */ 239 int *nbRecvElement = new int[mpiSize]; 240 MPI_Alltoall(nbSendElement, 1, MPI_INT, nbRecvElement, 1, MPI_INT, communicator); 241 242 /* communicate indices of source elements on other ranks whoes value and gradient we need (since intersection) */ 243 int nbSendRequest = 0; 244 int nbRecvRequest = 0; 245 int **recvElement = new int*[mpiSize]; 246 double **sendValue = new double*[mpiSize]; 247 double **sendArea = new double*[mpiSize]; 248 Coord **sendGrad = new Coord*[mpiSize]; 249 GloId **sendNeighIds = new GloId*[mpiSize]; 250 MPI_Request *sendRequest = new MPI_Request[4*mpiSize]; 251 MPI_Request *recvRequest = new MPI_Request[4*mpiSize]; 252 for (int rank = 0; rank < mpiSize; rank++) 253 { 254 if (nbSendElement[rank] > 0) 255 { 256 MPI_Issend(sendElement[rank], nbSendElement[rank], MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 257 nbSendRequest++; 258 } 259 260 if (nbRecvElement[rank] > 0) 261 { 262 recvElement[rank] = new int[nbRecvElement[rank]]; 263 sendValue[rank] = new double[nbRecvElement[rank]]; 264 sendArea[rank] = new double[nbRecvElement[rank]]; 265 if (order == 2) 266 { 267 sendNeighIds[rank] = new GloId[nbRecvElement[rank]*(NMAX+1)]; 268 sendGrad[rank] = new Coord[nbRecvElement[rank]*(NMAX+1)]; 269 } 270 else 271 { 272 sendNeighIds[rank] = new GloId[nbRecvElement[rank]]; 273 } 274 MPI_Irecv(recvElement[rank], nbRecvElement[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 275 nbRecvRequest++; 276 } 277 } 180 int mpiSize, mpiRank; 181 MPI_Comm_size(communicator, &mpiSize); 182 MPI_Comm_rank(communicator, &mpiRank); 183 184 /* create list of intersections (super mesh elements) for each rank */ 185 multimap<int, Polyg *> *elementList = new multimap<int, Polyg *>[mpiSize]; 186 for (int j = 0; j < nbElements; j++) 187 { 188 Elt& e = elements[j]; 189 for (list<Polyg *>::iterator it = e.is.begin(); it != e.is.end(); it++) 190 { 191 elementList[(*it)>id.rank].insert(pair<int, Polyg *>((*it)>id.ind, *it)); 192 //std::cout<<"elementList["<<(*it)>id.rank<<"].size = "<< elementList[(*it)>id.rank].size()<<std::endl; 193 } 194 } 195 196 int *nbSendElement = new int[mpiSize]; 197 int **sendElement = new int*[mpiSize]; /* indices of elements required from other rank */ 198 double **recvValue = new double*[mpiSize]; 199 double **recvArea = new double*[mpiSize]; 200 Coord **recvGrad = new Coord*[mpiSize]; 201 GloId **recvNeighIds = new GloId*[mpiSize]; /* ids of the of the source neighbours which also contribute through gradient */ 202 for (int rank = 0; rank < mpiSize; rank++) 203 { 204 /* get size for allocation */ 205 int last = 1; /* compares unequal to any index */ 206 int index = 1; /* increased to starting index 0 in first iteration */ 207 for (multimap<int, Polyg *>::iterator it = elementList[rank].begin(); it != elementList[rank].end(); ++it) 208 { 209 if (last != it>first) 210 index++; 211 (it>second)>id.ind = index; 212 last = it>first; 213 } 214 nbSendElement[rank] = index + 1; 215 216 /* if size is nonzero allocate and collect indices of elements on other ranks that we intersect */ 217 if (nbSendElement[rank] > 0) 218 { 219 sendElement[rank] = new int[nbSendElement[rank]]; 220 recvValue[rank] = new double[nbSendElement[rank]]; 221 recvArea[rank] = new double[nbSendElement[rank]]; 222 if (order == 2) 223 { 224 recvNeighIds[rank] = new GloId[nbSendElement[rank]*(NMAX+1)]; 225 recvGrad[rank] = new Coord[nbSendElement[rank]*(NMAX+1)]; 226 } 227 else 228 recvNeighIds[rank] = new GloId[nbSendElement[rank]]; 229 230 last = 1; 231 index = 1; 232 for (multimap<int, Polyg *>::iterator it = elementList[rank].begin(); it != elementList[rank].end(); ++it) 233 { 234 if (last != it>first) 235 index++; 236 sendElement[rank][index] = it>first; 237 last = it>first; 238 } 239 } 240 } 241 242 /* communicate sizes of source elements to be sent (index lists and later values and gradients) */ 243 int *nbRecvElement = new int[mpiSize]; 244 MPI_Alltoall(nbSendElement, 1, MPI_INT, nbRecvElement, 1, MPI_INT, communicator); 245 246 /* communicate indices of source elements on other ranks whoes value and gradient we need (since intersection) */ 247 int nbSendRequest = 0; 248 int nbRecvRequest = 0; 249 int **recvElement = new int*[mpiSize]; 250 double **sendValue = new double*[mpiSize]; 251 double **sendArea = new double*[mpiSize]; 252 Coord **sendGrad = new Coord*[mpiSize]; 253 GloId **sendNeighIds = new GloId*[mpiSize]; 254 MPI_Request *sendRequest = new MPI_Request[4*mpiSize]; 255 MPI_Request *recvRequest = new MPI_Request[4*mpiSize]; 256 for (int rank = 0; rank < mpiSize; rank++) 257 { 258 if (nbSendElement[rank] > 0) 259 { 260 MPI_Issend(sendElement[rank], nbSendElement[rank], MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 261 nbSendRequest++; 262 } 263 264 if (nbRecvElement[rank] > 0) 265 { 266 recvElement[rank] = new int[nbRecvElement[rank]]; 267 sendValue[rank] = new double[nbRecvElement[rank]]; 268 sendArea[rank] = new double[nbRecvElement[rank]]; 269 if (order == 2) 270 { 271 sendNeighIds[rank] = new GloId[nbRecvElement[rank]*(NMAX+1)]; 272 sendGrad[rank] = new Coord[nbRecvElement[rank]*(NMAX+1)]; 273 } 274 else 275 { 276 sendNeighIds[rank] = new GloId[nbRecvElement[rank]]; 277 } 278 MPI_Irecv(recvElement[rank], nbRecvElement[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 279 nbRecvRequest++; 280 } 281 } 278 282 279 283 MPI_Status *status = new MPI_Status[4*mpiSize]; … … 282 286 MPI_Waitall(nbRecvRequest, recvRequest, status); 283 287 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 288 /* for all indices that have been received from requesting ranks: pack values and gradients, then send */ 289 nbSendRequest = 0; 290 nbRecvRequest = 0; 291 for (int rank = 0; rank < mpiSize; rank++) 292 { 293 if (nbRecvElement[rank] > 0) 294 { 295 int jj = 0; // jj == j if no weight writing 296 for (int j = 0; j < nbRecvElement[rank]; j++) 297 { 298 sendValue[rank][j] = sstree.localElements[recvElement[rank][j]].val; 299 sendArea[rank][j] = sstree.localElements[recvElement[rank][j]].area; 300 if (order == 2) 301 { 302 sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].grad; 303 sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].src_id; 304 jj++; 305 for (int i = 0; i < NMAX; i++) 306 { 307 sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].gradNeigh[i]; 304 308 sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].neighId[i]; 305 jj++; 306 } 307 } 308 else 309 sendNeighIds[rank][j] = sstree.localElements[recvElement[rank][j]].src_id; 310 } 311 MPI_Issend(sendValue[rank], nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 312 nbSendRequest++; 313 MPI_Issend(sendArea[rank], nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 314 nbSendRequest++; 315 if (order == 2) 316 { 317 MPI_Issend(sendGrad[rank], 3*nbRecvElement[rank]*(NMAX+1), 318 MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 319 nbSendRequest++; 320 MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank]*(NMAX+1), MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 309 jj++; 310 } 311 } 312 else 313 sendNeighIds[rank][j] = sstree.localElements[recvElement[rank][j]].src_id; 314 } 315 MPI_Issend(sendValue[rank], nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 316 nbSendRequest++; 317 MPI_Issend(sendArea[rank], nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 318 nbSendRequest++; 319 if (order == 2) 320 { 321 MPI_Issend(sendGrad[rank], 3*nbRecvElement[rank]*(NMAX+1), MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 322 nbSendRequest++; 323 MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank]*(NMAX+1), MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 321 324 //ym > attention taille GloId 322 323 324 325 326 325 nbSendRequest++; 326 } 327 else 328 { 329 MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank], MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 327 330 //ym > attention taille GloId 328 nbSendRequest++; 329 } 330 } 331 if (nbSendElement[rank] > 0) 332 { 333 MPI_Irecv(recvValue[rank], nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 334 nbRecvRequest++; 335 MPI_Irecv(recvArea[rank], nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 336 nbRecvRequest++; 337 if (order == 2) 338 { 339 MPI_Irecv(recvGrad[rank], 3*nbSendElement[rank]*(NMAX+1), 340 MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 341 nbRecvRequest++; 342 MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank]*(NMAX+1), MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 331 nbSendRequest++; 332 } 333 } 334 if (nbSendElement[rank] > 0) 335 { 336 MPI_Irecv(recvValue[rank], nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 337 nbRecvRequest++; 338 MPI_Irecv(recvArea[rank], nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 339 nbRecvRequest++; 340 if (order == 2) 341 { 342 MPI_Irecv(recvGrad[rank], 3*nbSendElement[rank]*(NMAX+1), MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 343 nbRecvRequest++; 344 MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank]*(NMAX+1), MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 343 345 //ym > attention taille GloId 344 345 346 347 348 346 nbRecvRequest++; 347 } 348 else 349 { 350 MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 349 351 //ym > attention taille GloId 350 351 352 353 352 nbRecvRequest++; 353 } 354 } 355 } 354 356 355 356 357 MPI_Waitall(nbSendRequest, sendRequest, status); 358 MPI_Waitall(nbRecvRequest, recvRequest, status); 357 359 358 360 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 361 /* now that all values and gradients are available use them to computed interpolated values on target 362 and also to compute weights */ 363 int i = 0; 364 for (int j = 0; j < nbElements; j++) 365 { 366 Elt& e = elements[j]; 367 368 /* since for the 2nd order case source grid elements can contribute to a destination grid element over several "paths" 369 (step1: gradient is computed using neighbours on same grid, step2: intersection uses several elements on other grid) 370 accumulate them so that there is only one final weight between two elements */ 371 map<GloId,double> wgt_map; 372 373 /* for destination element `e` loop over all intersetions/the corresponding source elements */ 374 for (list<Polyg *>::iterator it = e.is.begin(); it != e.is.end(); it++) 375 { 376 /* it is the intersection element, so it>x and it>area are barycentre and area of intersection element (super mesh) 377 but it>id is id of the source element that it intersects */ 378 int n1 = (*it)>id.ind; 379 int rank = (*it)>id.rank; 380 double fk = recvValue[rank][n1]; 381 double srcArea = recvArea[rank][n1]; 382 double w = (*it)>area; 381 383 if (quantity) w/=srcArea ; 382 384 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 385 /* first order: src value times weight (weight = supermesh area), later divide by target area */ 386 int kk = (order == 2) ? n1 * (NMAX + 1) : n1; 387 GloId neighID = recvNeighIds[rank][kk]; 388 wgt_map[neighID] += w; 389 390 if (order == 2) 391 { 392 for (int k = 0; k < NMAX+1; k++) 393 { 394 int kk = n1 * (NMAX + 1) + k; 395 GloId neighID = recvNeighIds[rank][kk]; 396 if (neighID.ind != 1) wgt_map[neighID] += w * scalarprod(recvGrad[rank][kk], (*it)>x); 397 } 398 399 } 400 } 399 401 400 402 double renorm=0; … … 404 406 405 407 for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) 406 408 { 407 409 if (quantity) this>remapMatrix[i] = (it>second ) / renorm; 408 409 410 411 410 else this>remapMatrix[i] = (it>second / e.area) / renorm; 411 this>srcAddress[i] = it>first.ind; 412 this>srcRank[i] = it>first.rank; 413 this>dstAddress[i] = j; 412 414 this>sourceWeightId[i]= it>first.globalId ; 413 415 this>targetWeightId[i]= targetGlobalId[j] ; 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 delete[] recvNeighIds; 456 416 i++; 417 } 418 } 419 420 /* free all memory allocated in this function */ 421 /* for (int rank = 0; rank < mpiSize; rank++) 422 { 423 if (nbSendElement[rank] > 0) 424 { 425 delete[] sendElement[rank]; 426 delete[] recvValue[rank]; 427 delete[] recvArea[rank]; 428 if (order == 2) 429 { 430 delete[] recvGrad[rank]; 431 } 432 delete[] recvNeighIds[rank]; 433 } 434 if (nbRecvElement[rank] > 0) 435 { 436 delete[] recvElement[rank]; 437 delete[] sendValue[rank]; 438 delete[] sendArea[rank]; 439 if (order == 2) 440 delete[] sendGrad[rank]; 441 delete[] sendNeighIds[rank]; 442 } 443 } 444 delete[] status; 445 delete[] sendRequest; 446 delete[] recvRequest; 447 delete[] elementList; 448 delete[] nbSendElement; 449 delete[] nbRecvElement; 450 delete[] sendElement; 451 delete[] recvElement; 452 delete[] sendValue; 453 delete[] recvValue; 454 delete[] sendGrad; 455 delete[] recvGrad; 456 delete[] sendNeighIds; 457 delete[] recvNeighIds;*/ 458 return i; 457 459 } 458 460
Note: See TracChangeset
for help on using the changeset viewer.