Changeset 1614
 Timestamp:
 11/26/18 10:16:06 (5 months ago)
 Location:
 XIOS/trunk/extern/remap/src
 Files:

 5 edited
Legend:
 Unmodified
 Added
 Removed

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