Changeset 1460 for XIOS/dev/branch_openmp/extern
- Timestamp:
- 03/22/18 10:43:20 (6 years ago)
- Location:
- XIOS/dev/branch_openmp/extern
- Files:
-
- 8 edited
Legend:
- Unmodified
- Added
- Removed
-
XIOS/dev/branch_openmp/extern/remap/src/gridRemap.hpp
r1335 r1460 14 14 Coord readPole(std::istream&); 15 15 16 16 //extern CRemapGrid srcGrid; 17 //extern CRemapGrid tgtGrid; 17 18 18 19 } -
XIOS/dev/branch_openmp/extern/remap/src/mapper.cpp
r1341 r1460 29 29 void cptOffsetsFromLengths(const int *lengths, int *offsets, int sz) 30 30 { 31 32 33 31 offsets[0] = 0; 32 for (int i = 1; i < sz; i++) 33 offsets[i] = offsets[i-1] + lengths[i-1]; 34 34 } 35 35 … … 39 39 srcGrid.pole = Coord(pole[0], pole[1], pole[2]); 40 40 41 42 43 44 45 46 41 int mpiRank, mpiSize; 42 MPI_Comm_rank(communicator, &mpiRank); 43 MPI_Comm_size(communicator, &mpiSize); 44 45 sourceElements.reserve(nbCells); 46 sourceMesh.reserve(nbCells); 47 47 sourceGlobalId.resize(nbCells) ; 48 48 … … 57 57 else sourceGlobalId.assign(globalId,globalId+nbCells); 58 58 59 60 61 62 63 64 59 for (int i = 0; i < nbCells; i++) 60 { 61 int offs = i*nVertex; 62 Elt elt(boundsLon + offs, boundsLat + offs, nVertex); 63 elt.src_id.rank = mpiRank; 64 elt.src_id.ind = i; 65 65 elt.src_id.globalId = sourceGlobalId[i]; 66 67 68 69 66 sourceElements.push_back(elt); 67 sourceMesh.push_back(Node(elt.x, cptRadius(elt), &sourceElements.back())); 68 cptEltGeom(sourceElements[i], Coord(pole[0], pole[1], pole[2])); 69 } 70 70 71 71 } … … 75 75 tgtGrid.pole = Coord(pole[0], pole[1], pole[2]); 76 76 77 78 79 80 81 82 77 int mpiRank, mpiSize; 78 MPI_Comm_rank(communicator, &mpiRank); 79 MPI_Comm_size(communicator, &mpiSize); 80 81 targetElements.reserve(nbCells); 82 targetMesh.reserve(nbCells); 83 83 84 84 targetGlobalId.resize(nbCells) ; … … 93 93 else targetGlobalId.assign(globalId,globalId+nbCells); 94 94 95 96 97 98 99 100 101 102 95 for (int i = 0; i < nbCells; i++) 96 { 97 int offs = i*nVertex; 98 Elt elt(boundsLon + offs, boundsLat + offs, nVertex); 99 targetElements.push_back(elt); 100 targetMesh.push_back(Node(elt.x, cptRadius(elt), &sourceElements.back())); 101 cptEltGeom(targetElements[i], Coord(pole[0], pole[1], pole[2])); 102 } 103 103 104 104 … … 119 119 vector<double> Mapper::computeWeights(int interpOrder, bool renormalize, bool quantity) 120 120 { 121 vector<double> timings; 122 int mpiSize, mpiRank; 123 MPI_Comm_size(communicator, &mpiSize); 124 MPI_Comm_rank(communicator, &mpiRank); 125 121 vector<double> timings; 122 int mpiSize, mpiRank; 123 MPI_Comm_size(communicator, &mpiSize); 124 MPI_Comm_rank(communicator, &mpiRank); 126 125 127 126 this->buildSSTree(sourceMesh, targetMesh); 128 127 129 if (mpiRank == 0 && verbose) cout << "Computing intersections ..." << endl; 130 double tic = cputime(); 131 computeIntersection(&targetElements[0], targetElements.size()); 132 timings.push_back(cputime() - tic); 133 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]; 128 if (mpiRank == 0 && verbose) cout << "Computing intersections ..." << endl; 129 double tic = cputime(); 130 computeIntersection(&targetElements[0], targetElements.size()); 131 timings.push_back(cputime() - tic); 132 133 tic = cputime(); 134 if (interpOrder == 2) { 135 if (mpiRank == 0 && verbose) cout << "Computing grads ..." << endl; 136 buildMeshTopology(); 137 computeGrads(); 138 } 139 timings.push_back(cputime() - tic); 140 141 /* Prepare computation of weights */ 142 /* compute number of intersections which for the first order case 143 corresponds to the number of edges in the remap matrix */ 144 int nIntersections = 0; 145 for (int j = 0; j < targetElements.size(); j++) 146 { 147 Elt &elt = targetElements[j]; 148 for (list<Polyg*>::iterator it = elt.is.begin(); it != elt.is.end(); it++) 149 nIntersections++; 150 } 151 /* overallocate for NMAX neighbours for each elements */ 152 remapMatrix = new double[nIntersections*NMAX]; 153 srcAddress = new int[nIntersections*NMAX]; 154 srcRank = new int[nIntersections*NMAX]; 155 dstAddress = new int[nIntersections*NMAX]; 158 156 sourceWeightId =new long[nIntersections*NMAX]; 159 157 targetWeightId =new long[nIntersections*NMAX]; 160 158 161 159 162 163 164 165 160 if (mpiRank == 0 && verbose) cout << "Remapping..." << endl; 161 tic = cputime(); 162 nWeights = remap(&targetElements[0], targetElements.size(), interpOrder, renormalize, quantity); 163 timings.push_back(cputime() - tic); 166 164 167 165 for (int i = 0; i < targetElements.size(); i++) targetElements[i].delete_intersections(); 168 166 169 167 return timings; 170 168 } 171 169 … … 178 176 int Mapper::remap(Elt *elements, int nbElements, int order, bool renormalize, bool quantity) 179 177 { 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 non-zero 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 247 /* communicate indices of source elements on other ranks whoes value and gradient we need (since intersection) */ 248 int nbSendRequest = 0; 249 int nbRecvRequest = 0; 250 int **recvElement = new int*[mpiSize]; 251 double **sendValue = new double*[mpiSize]; 252 double **sendArea = new double*[mpiSize]; 253 Coord **sendGrad = new Coord*[mpiSize]; 254 GloId **sendNeighIds = new GloId*[mpiSize]; 255 MPI_Request *sendRequest = new MPI_Request[4*mpiSize]; 256 MPI_Request *recvRequest = new MPI_Request[4*mpiSize]; 257 for (int rank = 0; rank < mpiSize; rank++) 258 { 259 if (nbSendElement[rank] > 0) 260 { 261 MPI_Issend(sendElement[rank], nbSendElement[rank], MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 262 nbSendRequest++; 263 } 264 265 if (nbRecvElement[rank] > 0) 266 { 267 recvElement[rank] = new int[nbRecvElement[rank]]; 268 sendValue[rank] = new double[nbRecvElement[rank]]; 269 sendArea[rank] = new double[nbRecvElement[rank]]; 270 if (order == 2) 271 { 272 sendNeighIds[rank] = new GloId[nbRecvElement[rank]*(NMAX+1)]; 273 sendGrad[rank] = new Coord[nbRecvElement[rank]*(NMAX+1)]; 274 } 275 else 276 { 277 sendNeighIds[rank] = new GloId[nbRecvElement[rank]]; 278 } 279 MPI_Irecv(recvElement[rank], nbRecvElement[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 280 nbRecvRequest++; 281 } 282 } 283 284 MPI_Status *status = new MPI_Status[4*mpiSize]; 285 286 MPI_Waitall(nbSendRequest, sendRequest, status); 287 MPI_Waitall(nbRecvRequest, recvRequest, status); 288 289 /* for all indices that have been received from requesting ranks: pack values and gradients, then send */ 290 nbSendRequest = 0; 291 nbRecvRequest = 0; 292 for (int rank = 0; rank < mpiSize; rank++) 293 { 294 if (nbRecvElement[rank] > 0) 295 { 296 int jj = 0; // jj == j if no weight writing 297 for (int j = 0; j < nbRecvElement[rank]; j++) 298 { 299 sendValue[rank][j] = sstree.localElements[recvElement[rank][j]].val; 300 sendArea[rank][j] = sstree.localElements[recvElement[rank][j]].area; 301 if (order == 2) 302 { 303 sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].grad; 304 sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].src_id; 305 jj++; 306 for (int i = 0; i < NMAX; i++) 307 { 308 sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].gradNeigh[i]; 178 int mpiSize, mpiRank; 179 MPI_Comm_size(communicator, &mpiSize); 180 MPI_Comm_rank(communicator, &mpiRank); 181 182 /* create list of intersections (super mesh elements) for each rank */ 183 multimap<int, Polyg *> *elementList = new multimap<int, Polyg *>[mpiSize]; 184 for (int j = 0; j < nbElements; j++) 185 { 186 Elt& e = elements[j]; 187 for (list<Polyg *>::iterator it = e.is.begin(); it != e.is.end(); it++) 188 elementList[(*it)->id.rank].insert(pair<int, Polyg *>((*it)->id.ind, *it)); 189 } 190 191 int *nbSendElement = new int[mpiSize]; 192 int **sendElement = new int*[mpiSize]; /* indices of elements required from other rank */ 193 double **recvValue = new double*[mpiSize]; 194 double **recvArea = new double*[mpiSize]; 195 Coord **recvGrad = new Coord*[mpiSize]; 196 GloId **recvNeighIds = new GloId*[mpiSize]; /* ids of the of the source neighbours which also contribute through gradient */ 197 for (int rank = 0; rank < mpiSize; rank++) 198 { 199 /* get size for allocation */ 200 int last = -1; /* compares unequal to any index */ 201 int index = -1; /* increased to starting index 0 in first iteration */ 202 for (multimap<int, Polyg *>::iterator it = elementList[rank].begin(); it != elementList[rank].end(); ++it) 203 { 204 if (last != it->first) 205 index++; 206 (it->second)->id.ind = index; 207 last = it->first; 208 } 209 nbSendElement[rank] = index + 1; 210 211 /* if size is non-zero allocate and collect indices of elements on other ranks that we intersect */ 212 if (nbSendElement[rank] > 0) 213 { 214 sendElement[rank] = new int[nbSendElement[rank]]; 215 recvValue[rank] = new double[nbSendElement[rank]]; 216 recvArea[rank] = new double[nbSendElement[rank]]; 217 if (order == 2) 218 { 219 recvNeighIds[rank] = new GloId[nbSendElement[rank]*(NMAX+1)]; 220 recvGrad[rank] = new Coord[nbSendElement[rank]*(NMAX+1)]; 221 } 222 else 223 recvNeighIds[rank] = new GloId[nbSendElement[rank]]; 224 225 last = -1; 226 index = -1; 227 for (multimap<int, Polyg *>::iterator it = elementList[rank].begin(); it != elementList[rank].end(); ++it) 228 { 229 if (last != it->first) 230 index++; 231 sendElement[rank][index] = it->first; 232 last = it->first; 233 } 234 } 235 } 236 237 /* communicate sizes of source elements to be sent (index lists and later values and gradients) */ 238 int *nbRecvElement = new int[mpiSize]; 239 MPI_Alltoall(nbSendElement, 1, MPI_INT, nbRecvElement, 1, MPI_INT, communicator); 240 241 /* communicate indices of source elements on other ranks whoes value and gradient we need (since intersection) */ 242 int nbSendRequest = 0; 243 int nbRecvRequest = 0; 244 int **recvElement = new int*[mpiSize]; 245 double **sendValue = new double*[mpiSize]; 246 double **sendArea = new double*[mpiSize]; 247 Coord **sendGrad = new Coord*[mpiSize]; 248 GloId **sendNeighIds = new GloId*[mpiSize]; 249 MPI_Request *sendRequest = new MPI_Request[4*mpiSize]; 250 MPI_Request *recvRequest = new MPI_Request[4*mpiSize]; 251 for (int rank = 0; rank < mpiSize; rank++) 252 { 253 if (nbSendElement[rank] > 0) 254 { 255 MPI_Issend(sendElement[rank], nbSendElement[rank], MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 256 nbSendRequest++; 257 } 258 259 if (nbRecvElement[rank] > 0) 260 { 261 recvElement[rank] = new int[nbRecvElement[rank]]; 262 sendValue[rank] = new double[nbRecvElement[rank]]; 263 sendArea[rank] = new double[nbRecvElement[rank]]; 264 if (order == 2) 265 { 266 sendNeighIds[rank] = new GloId[nbRecvElement[rank]*(NMAX+1)]; 267 sendGrad[rank] = new Coord[nbRecvElement[rank]*(NMAX+1)]; 268 } 269 else 270 { 271 sendNeighIds[rank] = new GloId[nbRecvElement[rank]]; 272 } 273 MPI_Irecv(recvElement[rank], nbRecvElement[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 274 nbRecvRequest++; 275 } 276 } 277 MPI_Status *status = new MPI_Status[4*mpiSize]; 278 279 MPI_Waitall(nbSendRequest, sendRequest, status); 280 MPI_Waitall(nbRecvRequest, recvRequest, status); 281 282 /* for all indices that have been received from requesting ranks: pack values and gradients, then send */ 283 nbSendRequest = 0; 284 nbRecvRequest = 0; 285 for (int rank = 0; rank < mpiSize; rank++) 286 { 287 if (nbRecvElement[rank] > 0) 288 { 289 int jj = 0; // jj == j if no weight writing 290 for (int j = 0; j < nbRecvElement[rank]; j++) 291 { 292 sendValue[rank][j] = sstree.localElements[recvElement[rank][j]].val; 293 sendArea[rank][j] = sstree.localElements[recvElement[rank][j]].area; 294 if (order == 2) 295 { 296 sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].grad; 297 // cout<<"grad "<<jj<<" "<<recvElement[rank][j]<<" "<<sendGrad[rank][jj]<<" "<<sstree.localElements[recvElement[rank][j]].grad<<endl ; 298 sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].src_id; 299 jj++; 300 for (int i = 0; i < NMAX; i++) 301 { 302 sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].gradNeigh[i]; 303 // cout<<"grad "<<jj<<" "<<sendGrad[rank][jj]<<" "<<sstree.localElements[recvElement[rank][j]].grad<<endl ; 309 304 sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].neighId[i]; 310 jj++; 311 } 312 } 313 else 314 sendNeighIds[rank][j] = sstree.localElements[recvElement[rank][j]].src_id; 315 } 316 MPI_Issend(sendValue[rank], nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 317 nbSendRequest++; 318 MPI_Issend(sendArea[rank], nbRecvElement[rank], MPI_DOUBLE, rank, 1, communicator, &sendRequest[nbSendRequest]); 319 nbSendRequest++; 320 if (order == 2) 321 { 322 MPI_Issend(sendGrad[rank], 3*nbRecvElement[rank]*(NMAX+1), MPI_DOUBLE, rank, 2, communicator, &sendRequest[nbSendRequest]); 323 nbSendRequest++; 324 MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank]*(NMAX+1), MPI_INT, rank, 3, communicator, &sendRequest[nbSendRequest]); 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]); 325 321 //ym --> attention taille GloId 326 327 328 329 330 MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank], MPI_INT, rank, 4, communicator, &sendRequest[nbSendRequest]);322 nbSendRequest++; 323 } 324 else 325 { 326 MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank], MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 331 327 //ym --> attention taille GloId 332 nbSendRequest++; 333 } 334 } 335 if (nbSendElement[rank] > 0) 336 { 337 MPI_Irecv(recvValue[rank], nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 338 nbRecvRequest++; 339 MPI_Irecv(recvArea[rank], nbSendElement[rank], MPI_DOUBLE, rank, 1, communicator, &recvRequest[nbRecvRequest]); 340 nbRecvRequest++; 341 if (order == 2) 342 { 343 MPI_Irecv(recvGrad[rank], 3*nbSendElement[rank]*(NMAX+1), MPI_DOUBLE, rank, 2, communicator, &recvRequest[nbRecvRequest]); 344 nbRecvRequest++; 345 MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank]*(NMAX+1), MPI_INT, rank, 3, communicator, &recvRequest[nbRecvRequest]); 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]); 346 343 //ym --> attention taille GloId 347 348 349 350 351 MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank], MPI_INT, rank, 4, communicator, &recvRequest[nbRecvRequest]);344 nbRecvRequest++; 345 } 346 else 347 { 348 MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 352 349 //ym --> attention taille GloId 353 354 355 356 350 nbRecvRequest++; 351 } 352 } 353 } 357 354 358 MPI_Waitall(nbSendRequest, sendRequest, status);359 355 MPI_Waitall(nbSendRequest, sendRequest, status); 356 MPI_Waitall(nbRecvRequest, recvRequest, status); 360 357 361 358 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 359 /* now that all values and gradients are available use them to computed interpolated values on target 360 and also to compute weights */ 361 int i = 0; 362 for (int j = 0; j < nbElements; j++) 363 { 364 Elt& e = elements[j]; 365 366 /* since for the 2nd order case source grid elements can contribute to a destination grid element over several "paths" 367 (step1: gradient is computed using neighbours on same grid, step2: intersection uses several elements on other grid) 368 accumulate them so that there is only one final weight between two elements */ 369 map<GloId,double> wgt_map; 370 371 /* for destination element `e` loop over all intersetions/the corresponding source elements */ 372 for (list<Polyg *>::iterator it = e.is.begin(); it != e.is.end(); it++) 373 { 374 /* it is the intersection element, so it->x and it->area are barycentre and area of intersection element (super mesh) 375 but it->id is id of the source element that it intersects */ 376 int n1 = (*it)->id.ind; 377 int rank = (*it)->id.rank; 378 double fk = recvValue[rank][n1]; 379 double srcArea = recvArea[rank][n1]; 380 double w = (*it)->area; 384 381 if (quantity) w/=srcArea ; 385 382 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 383 /* first order: src value times weight (weight = supermesh area), later divide by target area */ 384 int kk = (order == 2) ? n1 * (NMAX + 1) : n1; 385 GloId neighID = recvNeighIds[rank][kk]; 386 wgt_map[neighID] += w; 387 388 if (order == 2) 389 { 390 for (int k = 0; k < NMAX+1; k++) 391 { 392 int kk = n1 * (NMAX + 1) + k; 393 GloId neighID = recvNeighIds[rank][kk]; 394 if (neighID.ind != -1) wgt_map[neighID] += w * scalarprod(recvGrad[rank][kk], (*it)->x); 395 } 396 397 } 398 } 402 399 403 400 double renorm=0; … … 407 404 408 405 for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) 409 406 { 410 407 if (quantity) this->remapMatrix[i] = (it->second ) / renorm; 411 412 413 414 408 else this->remapMatrix[i] = (it->second / e.area) / renorm; 409 this->srcAddress[i] = it->first.ind; 410 this->srcRank[i] = it->first.rank; 411 this->dstAddress[i] = j; 415 412 this->sourceWeightId[i]= it->first.globalId ; 416 413 this->targetWeightId[i]= targetGlobalId[j] ; 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 456 457 458 459 414 i++; 415 } 416 } 417 418 /* free all memory allocated in this function */ 419 for (int rank = 0; rank < mpiSize; rank++) 420 { 421 if (nbSendElement[rank] > 0) 422 { 423 delete[] sendElement[rank]; 424 delete[] recvValue[rank]; 425 delete[] recvArea[rank]; 426 if (order == 2) 427 { 428 delete[] recvGrad[rank]; 429 } 430 delete[] recvNeighIds[rank]; 431 } 432 if (nbRecvElement[rank] > 0) 433 { 434 delete[] recvElement[rank]; 435 delete[] sendValue[rank]; 436 delete[] sendArea[rank]; 437 if (order == 2) 438 delete[] sendGrad[rank]; 439 delete[] sendNeighIds[rank]; 440 } 441 } 442 delete[] status; 443 delete[] sendRequest; 444 delete[] recvRequest; 445 delete[] elementList; 446 delete[] nbSendElement; 447 delete[] nbRecvElement; 448 delete[] sendElement; 449 delete[] recvElement; 450 delete[] sendValue; 451 delete[] recvValue; 452 delete[] sendGrad; 453 delete[] recvGrad; 454 delete[] sendNeighIds; 455 delete[] recvNeighIds; 456 return i; 460 457 } 461 458 … … 496 493 mpiRoute.init(routes); 497 494 int nRecv = mpiRoute.getTotalSourceElement(); 495 // cout << mpiRank << " NRECV " << nRecv << "(" << routes.size() << ")"<< endl; 498 496 499 497 int *nbSendNode = new int[mpiSize]; … … 558 556 } 559 557 560 561 558 MPI_Waitall(nbRecvRequest, recvRequest, status); 559 MPI_Waitall(nbSendRequest, sendRequest, status); 562 560 563 561 for (int rank = 0; rank < mpiSize; rank++) … … 702 700 } 703 701 702 /* 703 for (int i = 0; i < elt->n; i++) 704 { 705 if (elt->neighbour[i] == NOT_FOUND) 706 error_exit("neighbour not found"); 707 } 708 */ 704 709 } 705 710 } … … 809 814 MPI_Waitall(nbRecvRequest, recvRequest, status); 810 815 MPI_Waitall(nbSendRequest, sendRequest, status); 811 812 816 char **sendBuffer2 = new char*[mpiSize]; 813 817 char **recvBuffer2 = new char*[mpiSize]; … … 836 840 intersect_ym(&recvElt[j], elt2); 837 841 } 842 838 843 if (recvElt[j].is.size() > 0) sentMessageSize[rank] += packIntersectionSize(recvElt[j]); 839 844 … … 854 859 } 855 860 delete [] recvElt; 861 856 862 } 857 863 } … … 889 895 } 890 896 891 892 897 MPI_Waitall(nbRecvRequest, recvRequest, status); 898 MPI_Waitall(nbSendRequest, sendRequest, status); 893 899 894 900 delete [] sendRequest; -
XIOS/dev/branch_openmp/extern/remap/src/node.cpp
r1328 r1460 281 281 } 282 282 283 283 284 return q; 284 285 } -
XIOS/dev/branch_openmp/extern/remap/src/parallel_tree.cpp
r1355 r1460 197 197 delete[] displs; 198 198 199 /* unpack */ 200 /* 201 randomArray.resize(nrecv); 202 randomizeArray(randomArray); 203 tree.leafs.resize(nrecv); 204 index = 0; 205 206 207 for (int i = 0; i < nrecv; i++) 208 { 209 Coord x = *(Coord *)(&recvBuffer[index]); 210 index += sizeof(Coord)/sizeof(*recvBuffer); 211 double radius = recvBuffer[index++]; 212 tree.leafs[randomArray[i]].centre = x; 213 tree.leafs[randomArray[i]].radius = radius; 214 215 } 216 */ 199 217 200 218 randomArray.resize(blocSize); … … 231 249 cerr << comm.rank << ": PROBLEM: (node assign)" << tree.levelSize[assignLevel] << " != " << comm.group_size << " (keepNodes)" 232 250 << " node size : "<<node.size()<<" bloc size : "<<blocSize<<" total number of leaf : "<<tree.leafs.size()<<endl ; 233 251 /* 252 MPI_Allreduce(&ok, &allok, 1, MPI_INT, MPI_PROD, communicator); 253 if (!allok) { 254 MPI_Finalize(); 255 exit(1); 256 } 257 */ 234 258 MPI_Abort(MPI_COMM_WORLD,-1) ; 235 259 } … … 302 326 nb=nb1+nb2 ; 303 327 MPI_Allreduce(&nb, &nbTot, 1, MPI_LONG, MPI_SUM, communicator) ; 304 305 306 328 int commSize ; 307 329 MPI_Comm_size(communicator,&commSize) ; … … 326 348 randomizeArray(randomArray2); 327 349 350 /* 351 int s1,s2 ; 352 353 if (node.size()< nbSampleNodes/2) 354 { 355 s1 = node.size() ; 356 s2 = nbSampleNodes-s1 ; 357 } 358 else if (node2.size()< nbSampleNodes/2) 359 { 360 s2 = node.size() ; 361 s1 = nbSampleNodes-s2 ; 362 } 363 else 364 { 365 s1=nbSampleNodes/2 ; 366 s2=nbSampleNodes/2 ; 367 } 368 */ 328 369 for (int i = 0; i <nbSampleNodes1; i++) sampleNodes.push_back(Node(node[randomArray1[i%nb1]].centre, node[randomArray1[i%nb1]].radius, NULL)); 329 370 for (int i = 0; i <nbSampleNodes2; i++) sampleNodes.push_back(Node(node2[randomArray2[i%nb2]].centre, node2[randomArray2[i%nb2]].radius, NULL)); 330 371 372 /* 373 for (int i = 0; i < nbSampleNodes/2; i++) 374 { 375 sampleNodes.push_back(Node(node[randomArray1[i]].centre, node[randomArray1[i]].radius, NULL)); 376 sampleNodes.push_back(Node(node2[randomArray2[i]].centre, node2[randomArray2[i]].radius, NULL)); 377 } 378 */ 331 379 CTimer::get("buildParallelSampleTree").resume(); 332 380 //sampleTree.buildParallelSampleTree(sampleNodes, cascade); -
XIOS/dev/branch_openmp/extern/remap/src/parallel_tree.hpp
r1328 r1460 6 6 #include "mpi_cascade.hpp" 7 7 #include "mpi.hpp" 8 8 9 namespace sphereRemap { 9 10 -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_declaration.cpp
r1373 r1460 24 24 ::MPI_Datatype MPI_UNSIGNED_LONG_STD = MPI_UNSIGNED_LONG; 25 25 ::MPI_Datatype MPI_UNSIGNED_CHAR_STD = MPI_UNSIGNED_CHAR; 26 ::MPI_Datatype MPI_UINT64_T_STD = MPI_UINT64_T; 26 27 27 28 #undef MPI_INT … … 32 33 #undef MPI_UNSIGNED_LONG 33 34 #undef MPI_UNSIGNED_CHAR 35 #undef MPI_UINT64_T 34 36 35 37 … … 37 39 ::MPI_Op MPI_MAX_STD = MPI_MAX; 38 40 ::MPI_Op MPI_MIN_STD = MPI_MIN; 41 ::MPI_Op MPI_LOR_STD = MPI_LOR; 39 42 40 43 #undef MPI_SUM 41 44 #undef MPI_MAX 42 45 #undef MPI_MIN 43 44 45 /*#undef MPI_INT 46 #undef MPI_FLOAT 47 #undef MPI_DOUBLE 48 #undef MPI_CHAR 49 #undef MPI_LONG 50 #undef MPI_UNSIGNED_LONG 51 #undef MPI_UNSIGNED_CHAR 52 53 #undef MPI_SUM 54 #undef MPI_MAX 55 #undef MPI_MIN 56 57 #undef MPI_COMM_WORLD 58 #undef MPI_COMM_NULL 59 60 #undef MPI_STATUS_IGNORE 61 #undef MPI_REQUEST_NULL 62 #undef MPI_INFO_NULL 63 */ 46 #undef MPI_LOR 64 47 65 48 … … 73 56 extern ::MPI_Datatype MPI_UNSIGNED_LONG_STD; 74 57 extern ::MPI_Datatype MPI_UNSIGNED_CHAR_STD; 58 extern ::MPI_Datatype MPI_UINT64_T_STD; 59 75 60 76 61 extern ::MPI_Op MPI_SUM_STD; 77 62 extern ::MPI_Op MPI_MAX_STD; 78 63 extern ::MPI_Op MPI_MIN_STD; 64 extern ::MPI_Op MPI_LOR_STD; 79 65 80 66 extern ::MPI_Comm MPI_COMM_WORLD_STD; … … 92 78 ep_lib::MPI_Datatype MPI_UNSIGNED_LONG = &MPI_UNSIGNED_LONG_STD; 93 79 ep_lib::MPI_Datatype MPI_UNSIGNED_CHAR = &MPI_UNSIGNED_CHAR_STD; 80 ep_lib::MPI_Datatype MPI_UINT64_T = &MPI_UINT64_T_STD; 81 94 82 95 83 ep_lib::MPI_Op MPI_SUM = &MPI_SUM_STD; 96 84 ep_lib::MPI_Op MPI_MAX = &MPI_MAX_STD; 97 85 ep_lib::MPI_Op MPI_MIN = &MPI_MIN_STD; 86 ep_lib::MPI_Op MPI_LOR = &MPI_LOR_STD; 98 87 99 88 ep_lib::MPI_Comm MPI_COMM_WORLD(&MPI_COMM_WORLD_STD); -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_declaration.hpp
r1369 r1460 9 9 #undef MPI_UNSIGNED_LONG 10 10 #undef MPI_UNSIGNED_CHAR 11 #undef MPI_UINT64_T 11 12 12 13 #undef MPI_SUM 13 14 #undef MPI_MAX 14 15 #undef MPI_MIN 16 #undef MPI_LOR 15 17 16 18 #undef MPI_COMM_WORLD … … 28 30 extern ep_lib::MPI_Datatype MPI_UNSIGNED_LONG; 29 31 extern ep_lib::MPI_Datatype MPI_UNSIGNED_CHAR; 32 extern ep_lib::MPI_Datatype MPI_UINT64_T; 30 33 31 34 extern ep_lib::MPI_Op MPI_SUM; 32 35 extern ep_lib::MPI_Op MPI_MAX; 33 36 extern ep_lib::MPI_Op MPI_MIN; 37 extern ep_lib::MPI_Op MPI_LOR; 34 38 35 39 extern ep_lib::MPI_Comm MPI_COMM_WORLD; -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_reduce.cpp
r1365 r1460 29 29 30 30 template<typename T> 31 T lor_op(T a, T b) 32 { 33 return a||b; 34 } 35 36 template<typename T> 31 37 void reduce_max(const T * buffer, T* recvbuf, int count) 32 38 { … … 44 50 { 45 51 transform(buffer, buffer+count, recvbuf, recvbuf, std::plus<T>()); 52 } 53 54 template<typename T> 55 void reduce_lor(const T * buffer, T* recvbuf, int count) 56 { 57 transform(buffer, buffer+count, recvbuf, recvbuf, lor_op<T>); 46 58 } 47 59 … … 112 124 } 113 125 126 else if(datatype == MPI_UINT64_T) 127 { 128 assert(datasize == sizeof(uint64_t)); 129 for(int i=1; i<num_ep; i++) 130 reduce_max<uint64_t>(static_cast<uint64_t*>(comm.my_buffer->void_buffer[i]), static_cast<uint64_t*>(recvbuf), count); 131 } 132 114 133 else printf("datatype Error\n"); 115 134 … … 160 179 } 161 180 181 else if(datatype == MPI_UINT64_T) 182 { 183 assert(datasize == sizeof(uint64_t)); 184 for(int i=1; i<num_ep; i++) 185 reduce_min<uint64_t>(static_cast<uint64_t*>(comm.my_buffer->void_buffer[i]), static_cast<uint64_t*>(recvbuf), count); 186 } 187 162 188 else printf("datatype Error\n"); 163 189 … … 209 235 } 210 236 237 else if(datatype ==MPI_UINT64_T) 238 { 239 assert(datasize == sizeof(uint64_t)); 240 for(int i=1; i<num_ep; i++) 241 reduce_sum<uint64_t>(static_cast<uint64_t*>(comm.my_buffer->void_buffer[i]), static_cast<uint64_t*>(recvbuf), count); 242 } 243 211 244 else printf("datatype Error\n"); 212 245 246 } 247 248 if(op == MPI_SUM) 249 { 250 if(datatype != MPI_INT) 251 printf("datatype Error, must be MPI_INT\n"); 252 else 253 { 254 assert(datasize == sizeof(int)); 255 for(int i=1; i<num_ep; i++) 256 reduce_lor<int>(static_cast<int*>(comm.my_buffer->void_buffer[i]), static_cast<int*>(recvbuf), count); 257 } 213 258 } 214 259 }
Note: See TracChangeset
for help on using the changeset viewer.