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