Changeset 1328 for XIOS/dev/branch_openmp/extern
- Timestamp:
- 11/15/17 12:14:34 (7 years ago)
- Location:
- XIOS/dev/branch_openmp/extern
- Files:
-
- 28 edited
Legend:
- Unmodified
- Added
- Removed
-
XIOS/dev/branch_openmp/extern/remap/src/clipper.cpp
r1205 r1328 4274 4274 { 4275 4275 //The equation of a line in general form (Ax + By + C = 0) 4276 //given 2 points (x ,y) & (x,y) is ...4277 //(y - y)x + (x - x)y + (y - y)x - (x - x)y= 04278 //A = (y - y); B = (x - x); C = (y - y)x - (x - x)y4279 //perpendicular distance of point (x ,y) = (Ax + By + C)/Sqrt(A + B)4276 //given 2 points (x¹,y¹) & (x²,y²) is ... 4277 //(y¹ - y²)x + (x² - x¹)y + (y² - y¹)x¹ - (x² - x¹)y¹ = 0 4278 //A = (y¹ - y²); B = (x² - x¹); C = (y² - y¹)x¹ - (x² - x¹)y¹ 4279 //perpendicular distance of point (x³,y³) = (Ax³ + By³ + C)/Sqrt(A² + B²) 4280 4280 //see http://en.wikipedia.org/wiki/Perpendicular_distance 4281 4281 double A = double(ln1.Y - ln2.Y); -
XIOS/dev/branch_openmp/extern/remap/src/cputime.cpp
r694 r1328 1 1 #include "mpi.hpp" 2 using namespace ep_lib; 2 3 3 4 namespace sphereRemap { -
XIOS/dev/branch_openmp/extern/remap/src/elt.hpp
r1172 r1328 53 53 struct Elt : Polyg 54 54 { 55 Elt() {} 56 Elt(const double *bounds_lon, const double *bounds_lat, int max_num_vert) 57 { 58 int k = 0; 59 vertex[k++] = xyz(bounds_lon[0], bounds_lat[0]); 60 for (int i = 1; i < max_num_vert; i++) 61 { 62 vertex[k] = xyz(bounds_lon[i], bounds_lat[i]); 63 /* netCDF convention: if first vertex repeats element is finished (at least three vertices == triagle) */ 64 if (k >= 3 && squaredist(vertex[k], vertex[0]) < EPS*EPS) 65 break; 66 /* eliminate zero edges: move to next vertex only if it is different */ 67 if (squaredist(vertex[k], vertex[k-1]) > EPS*EPS) 68 k++; 69 //else cout << "Removed edge " << k << " due to zero length (coinciding endpoints)." << endl ; 70 } 71 n = k; 72 x = barycentre(vertex, n); 73 } 55 Elt() {} 56 Elt(const double *bounds_lon, const double *bounds_lat, int max_num_vert) 57 { 58 int k = 0; 59 vertex[k++] = xyz(bounds_lon[0], bounds_lat[0]); 60 for (int i = 1; i < max_num_vert; i++) 61 { 62 vertex[k] = xyz(bounds_lon[i], bounds_lat[i]); 63 /* netCDF convention: if first vertex repeats element is finished (at least three vertices == triagle) */ 64 if (k >= 3 && squaredist(vertex[k], vertex[0]) < EPS*EPS) 65 break; 66 /* eliminate zero edges: move to next vertex only if it is different */ 67 if (squaredist(vertex[k], vertex[k-1]) > EPS*EPS) 68 k++; 69 else 70 /* cout << "Removed edge " << k << " due to zero length (coinciding endpoints)." << endl */ ; 71 } 72 n = k; 73 x = barycentre(vertex, n); 74 } 74 75 75 76 Elt& operator=(const Elt& rhs) … … 95 96 } 96 97 97 98 99 100 101 102 103 104 98 void delete_intersections() 99 { 100 for (list<Polyg*>::iterator it = this->is.begin(); it != this->is.end(); it++) 101 { 102 Polyg* poly = *it; 103 delete poly; 104 } 105 } 105 106 106 107 void insert_vertex(int i, const Coord& v) -
XIOS/dev/branch_openmp/extern/remap/src/gridRemap.cpp
r1155 r1328 11 11 12 12 CRemapGrid srcGrid; 13 #pragma omp threadprivate(srcGrid)14 15 13 CRemapGrid tgtGrid; 16 #pragma omp threadprivate(tgtGrid)17 14 18 15 } -
XIOS/dev/branch_openmp/extern/remap/src/gridRemap.hpp
r1220 r1328 14 14 Coord readPole(std::istream&); 15 15 16 16 extern CRemapGrid srcGrid; 17 extern CRemapGrid tgtGrid; 17 18 18 19 } -
XIOS/dev/branch_openmp/extern/remap/src/intersect.cpp
r1220 r1328 15 15 namespace sphereRemap { 16 16 17 extern CRemapGrid srcGrid;18 #pragma omp threadprivate(srcGrid)19 20 extern CRemapGrid tgtGrid;21 #pragma omp threadprivate(tgtGrid)22 23 17 using namespace std; 24 18 … … 27 21 int neighbour_idx(const Elt& a, const Elt& b) 28 22 { 29 30 31 32 33 34 35 36 37 38 39 40 41 42 23 for (int i = 0; i < a.n; i++) 24 { 25 for (int j = 0; j < b.n; j++) 26 { 27 assert(squaredist(a.vertex[ i ], b.vertex[ j ]) > EPS*EPS || 28 squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+1)%b.n]) > EPS*EPS); 29 if ( squaredist(a.vertex[ i ], b.vertex[ j ]) < 1e-13*1e-13 && 30 squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+b.n-1)%b.n]) < 1e-13*1e-13) 31 { 32 return i; 33 } 34 } 35 } 36 return NOT_FOUND; 43 37 } 44 38 … … 211 205 bool isNeighbour(Elt& a, const Elt& b) 212 206 { 213 207 // return neighbour_idx(a, b) != NOT_FOUND; 214 208 return insertNeighbour(a,b,false) ; 215 209 } … … 219 213 void intersect(Elt *a, Elt *b) 220 214 { 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 // 215 int na = a->n; /* vertices of a */ 216 int nb = b->n; /* vertices of b */ 217 Coord *c = new Coord[na+nb]; 218 Coord *c2 = new Coord[na+nb]; 219 Coord *xc = new Coord[na+nb]; 220 Coord *xc2 = new Coord[na+nb]; 221 Coord gc, gc2; 222 double *d = new double[na+nb]; 223 double *d2 = new double[na+nb]; 224 double are, are2; 225 Ipt ipt[NMAX*NMAX]; 226 Ipt ipt2[NMAX*NMAX]; 227 ptsec(a, b, ipt); 228 /* make ipt2 transpose of ipt */ 229 for (int ii = 0; ii < na; ii++) 230 for (int jj = 0; jj < nb; jj++) 231 ipt2[jj*na+ii] = ipt[ii*nb+jj]; 232 list<Sgm> iscot; 233 recense(a, b, ipt, iscot, 0); 234 recense(b, a, ipt2, iscot, 1); 235 236 int nseg = iscot.size(); 237 int nc = 0; 238 int nc2 = 0; 239 while (iscot.size() && nc < 2) 240 nc = assemble(iscot, c, d, xc); 241 while (iscot.size() && nc2 < 2) 242 nc2 = assemble(iscot, c2, d2, xc2); 243 // assert(nseg == nc + nc2 || nseg == 1); // unused segment 250 244 251 245 if (!(nseg == nc + nc2 || nseg == 1)) … … 266 260 267 261 // intersect_ym(a,b) ; 268 269 270 271 272 273 274 275 276 277 278 279 262 if (nc == 1) nc = 0; 263 if (nc2 == 1) nc2 = 0; 264 gc = barycentre(xc, nc); 265 gc2 = barycentre(xc2, nc2); 266 orient(nc, xc, c, d, gc); 267 268 Coord pole = srcGrid.pole; 269 if (pole == ORIGIN) pole = tgtGrid.pole; 270 const double MINBASE = 1e-11; 271 if (nc == 2) /* nc is the number of vertices of super mesh element */ 272 { 273 double base = arcdist(xc[0], xc[1]); 280 274 cerr << "DID ARRIVE " << base << xc[0] << xc[1] << endl; 281 282 283 284 285 286 287 288 289 290 291 275 gc = midpoint(gc, midpointSC(xc[0], xc[1])); 276 /* intersection area `are` must be zero here unless we have one great and one small circle */ 277 are = alun(base, fabs(scalarprod(xc[0], pole))); 278 } 279 else 280 { 281 are = airbar(nc, xc, c, d, pole, gc); 282 } 283 if (nc2 == 2) 284 { 285 double base = arcdist(xc2[0], xc2[1]); 292 286 cerr << "DID ARRIVE " << base << xc2[0] << xc2[1] << endl; 293 294 295 296 297 298 299 300 287 assert(base > MINBASE); 288 gc2 = midpoint(gc2, midpointSC(xc2[0], xc2[1])); 289 are2 = alun(base, fabs(scalarprod(xc2[0], pole))); // 0 290 } 291 else 292 { 293 are2 = airbar(nc2, xc2, c2, d2, pole, gc2); 294 } 301 295 302 296 // double ym_area=intersect_ym(a,b) ; 303 297 304 298 if (nc > 1) 305 306 307 308 309 310 311 312 313 314 299 { 300 /* create one super mesh polygon that src and dest point to */ 301 Polyg *is = new Polyg; 302 is->x = gc; 303 is->area = are; 304 is->id = b->id; 305 is->src_id = b->src_id; 306 is->n = nc; 307 (a->is).push_back(is); 308 (b->is).push_back(is); 315 309 /* 316 if (2*fabs(are-ym_area)/(are+ym_area) > 1.1 && ym_area>1e-8)310 if ( 2*fabs(are-ym_area)/(are+ym_area) > 1.1 && ym_area>1e-8) 317 311 { 318 312 cout<<"Big area difference : "<<are<<" "<<ym_area<<endl ; … … 320 314 } 321 315 */ 322 // 323 324 325 326 327 328 329 330 331 332 333 316 // cout<<"intersection : "<<are<<" "<< ym_area<<" diff : "<<fabs(are-ym_area)<<" ratio : "<<fabs(are-ym_area)/(0.5*(are+ym_area))<<endl ; 317 } 318 if (nc2 > 1) 319 { 320 Polyg *is = new Polyg; 321 is->x = gc2; 322 is->area = are2; 323 is->id = b->id; /* intersection holds id of corresponding source element (see Elt class definition for details about id) */ 324 is->src_id = b->src_id; 325 is->n = nc2; 326 (a->is).push_back(is); 327 (b->is).push_back(is); 334 328 /* 335 if ( 329 if ( 2*fabs(are-ym_area)/(are+ym_area) > 1.1 && ym_area>1e-8 ) 336 330 { 337 331 cout<<"Big area difference : "<<are<<" "<<ym_area<<endl ; … … 339 333 } 340 334 */ 341 // 342 335 // cout<<"intersection : "<<are2<<" "<< ym_area<<" diff : "<<fabs(are-ym_area)<<" ratio : "<<fabs(are-ym_area)/(0.5*(are+ym_area))<<endl ; 336 } 343 337 /* 344 338 if (nc<=1 && nc2<=1) … … 350 344 } 351 345 */ 352 353 354 355 356 357 358 } 359 360 } 346 delete [] c; 347 delete [] c2; 348 delete [] xc; 349 delete [] xc2; 350 delete [] d; 351 delete [] d2; 352 } 353 354 } -
XIOS/dev/branch_openmp/extern/remap/src/intersection_ym.cpp
r1220 r1328 13 13 14 14 namespace sphereRemap { 15 16 extern CRemapGrid srcGrid;17 #pragma omp threadprivate(srcGrid)18 19 extern CRemapGrid tgtGrid;20 #pragma omp threadprivate(tgtGrid)21 22 15 23 16 using namespace std; -
XIOS/dev/branch_openmp/extern/remap/src/libmapper.cpp
r1205 r1328 14 14 #include "mapper.hpp" 15 15 #include "cputime.hpp" // cputime 16 #include <stdio.h> 16 17 using namespace ep_lib; 17 18 18 19 using namespace sphereRemap ; … … 21 22 and deallocated during the second step (computing the weights) */ 22 23 Mapper *mapper; 23 #pragma omp threadprivate(mapper) 24 24 25 25 26 /** xxx_bounds_yyy is of length n_vert_per_cell_xxx*n_cell_xxx … … 33 34 int order, int* n_weights) 34 35 { 35 36 37 38 39 40 41 42 43 36 assert(src_bounds_lon); 37 assert(src_bounds_lat); 38 assert(n_vert_per_cell_src >= 3); 39 assert(n_cell_src >= 4); 40 assert(dst_bounds_lon); 41 assert(dst_bounds_lat); 42 assert(n_vert_per_cell_dst >= 3); 43 assert(n_cell_dst >= 4); 44 assert(1 <= order && order <= 2); 44 45 45 46 mapper = new Mapper(MPI_COMM_WORLD); … … 80 81 double tic = cputime(); 81 82 mapper = new Mapper(MPI_COMM_WORLD); 82 83 mapper->setVerbosity(PROGRESS) ; 83 84 mapper->buildSSTree(src_msh, dst_msh); 84 85 double tac = cputime(); … … 149 150 char **argv = NULL; 150 151 MPI_Init(&argc, &argv);*/ 151 //MPI_Init(NULL, NULL); 152 int provided; 153 MPI_Init_thread(NULL, NULL, 3, &provided); 154 assert(provided >= 3); 152 MPI_Init(NULL, NULL); 155 153 } 156 154 -
XIOS/dev/branch_openmp/extern/remap/src/mapper.cpp
r1220 r1328 12 12 13 13 #include "mapper.hpp" 14 using namespace ep_lib; 14 15 15 16 namespace sphereRemap { 16 17 extern CRemapGrid srcGrid;18 #pragma omp threadprivate(srcGrid)19 20 extern CRemapGrid tgtGrid;21 #pragma omp threadprivate(tgtGrid)22 23 17 24 18 /* A subdivition of an array into N sub-arays … … 72 66 void Mapper::setTargetMesh(const double* boundsLon, const double* boundsLat, int nVertex, int nbCells, const double* pole, const long int* globalId) 73 67 { 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 68 tgtGrid.pole = Coord(pole[0], pole[1], pole[2]); 69 70 int mpiRank, mpiSize; 71 MPI_Comm_rank(communicator, &mpiRank); 72 MPI_Comm_size(communicator, &mpiSize); 73 74 targetElements.reserve(nbCells); 75 targetMesh.reserve(nbCells); 76 77 targetGlobalId.resize(nbCells) ; 78 if (globalId==NULL) 79 { 80 long int offset ; 81 long int nb=nbCells ; 82 MPI_Scan(&nb,&offset,1,MPI_LONG,MPI_SUM,communicator) ; 83 offset=offset-nb ; 84 for(int i=0;i<nbCells;i++) targetGlobalId[i]=offset+i ; 85 } 86 else targetGlobalId.assign(globalId,globalId+nbCells); 87 88 for (int i = 0; i < nbCells; i++) 89 { 90 int offs = i*nVertex; 91 Elt elt(boundsLon + offs, boundsLat + offs, nVertex); 92 targetElements.push_back(elt); 93 targetMesh.push_back(Node(elt.x, cptRadius(elt), &sourceElements.back())); 94 cptEltGeom(targetElements[i], Coord(pole[0], pole[1], pole[2])); 95 } 102 96 103 97 … … 106 100 void Mapper::setSourceValue(const double* val) 107 101 { 108 109 102 int size=sourceElements.size() ; 103 for(int i=0;i<size;++i) sourceElements[i].val=val[i] ; 110 104 } 111 105 112 106 void Mapper::getTargetValue(double* val) 113 107 { 114 115 108 int size=targetElements.size() ; 109 for(int i=0;i<size;++i) val[i]=targetElements[i].val ; 116 110 } 117 111 118 112 vector<double> Mapper::computeWeights(int interpOrder, bool renormalize, bool quantity) 119 113 { 120 vector<double> timings; 121 int mpiSize, mpiRank; 122 MPI_Comm_size(communicator, &mpiSize); 123 MPI_Comm_rank(communicator, &mpiRank); 124 125 this->buildSSTree(sourceMesh, targetMesh); 126 127 if (mpiRank == 0 && verbose) cout << "Computing intersections ..." << endl; 128 double tic = cputime(); 129 computeIntersection(&targetElements[0], targetElements.size()); 130 timings.push_back(cputime() - tic); 131 132 tic = cputime(); 133 if (interpOrder == 2) { 134 if (mpiRank == 0 && verbose) cout << "Computing grads ..." << endl; 135 buildMeshTopology(); 136 computeGrads(); 137 } 138 timings.push_back(cputime() - tic); 139 140 /* Prepare computation of weights */ 141 /* compute number of intersections which for the first order case 142 corresponds to the number of edges in the remap matrix */ 143 int nIntersections = 0; 144 for (int j = 0; j < targetElements.size(); j++) 145 { 146 Elt &elt = targetElements[j]; 147 for (list<Polyg*>::iterator it = elt.is.begin(); it != elt.is.end(); it++) 148 nIntersections++; 149 } 150 /* overallocate for NMAX neighbours for each elements */ 151 remapMatrix = new double[nIntersections*NMAX]; 152 srcAddress = new int[nIntersections*NMAX]; 153 srcRank = new int[nIntersections*NMAX]; 154 dstAddress = new int[nIntersections*NMAX]; 155 sourceWeightId =new long[nIntersections*NMAX]; 156 targetWeightId =new long[nIntersections*NMAX]; 157 158 159 if (mpiRank == 0 && verbose) cout << "Remapping..." << endl; 160 tic = cputime(); 161 nWeights = remap(&targetElements[0], targetElements.size(), interpOrder, renormalize, quantity); 162 timings.push_back(cputime() - tic); 163 164 for (int i = 0; i < targetElements.size(); i++) targetElements[i].delete_intersections(); 165 166 return timings; 114 vector<double> timings; 115 int mpiSize, mpiRank; 116 MPI_Comm_size(communicator, &mpiSize); 117 MPI_Comm_rank(communicator, &mpiRank); 118 119 120 this->buildSSTree(sourceMesh, targetMesh); 121 122 if (mpiRank == 0 && verbose) cout << "Computing intersections ..." << endl; 123 double tic = cputime(); 124 computeIntersection(&targetElements[0], targetElements.size()); 125 timings.push_back(cputime() - tic); 126 127 tic = cputime(); 128 if (interpOrder == 2) { 129 if (mpiRank == 0 && verbose) cout << "Computing grads ..." << endl; 130 buildMeshTopology(); 131 computeGrads(); 132 } 133 timings.push_back(cputime() - tic); 134 135 /* Prepare computation of weights */ 136 /* compute number of intersections which for the first order case 137 corresponds to the number of edges in the remap matrix */ 138 int nIntersections = 0; 139 for (int j = 0; j < targetElements.size(); j++) 140 { 141 Elt &elt = targetElements[j]; 142 for (list<Polyg*>::iterator it = elt.is.begin(); it != elt.is.end(); it++) 143 nIntersections++; 144 } 145 /* overallocate for NMAX neighbours for each elements */ 146 remapMatrix = new double[nIntersections*NMAX]; 147 srcAddress = new int[nIntersections*NMAX]; 148 srcRank = new int[nIntersections*NMAX]; 149 dstAddress = new int[nIntersections*NMAX]; 150 sourceWeightId =new long[nIntersections*NMAX]; 151 targetWeightId =new long[nIntersections*NMAX]; 152 153 154 if (mpiRank == 0 && verbose) cout << "Remapping..." << endl; 155 tic = cputime(); 156 nWeights = remap(&targetElements[0], targetElements.size(), interpOrder, renormalize, quantity); 157 timings.push_back(cputime() - tic); 158 159 for (int i = 0; i < targetElements.size(); i++) targetElements[i].delete_intersections(); 160 161 return timings; 167 162 } 168 163 169 164 /** 170 @param elements are cells of the target grid that are distributed over CPUs171 indepentently of the distribution of the SS-tree.172 @param nbElements is the size of the elements array.173 @param order is the order of interpolaton (must be 1 or 2).174 165 @param elements are cells of the target grid that are distributed over CPUs 166 indepentently of the distribution of the SS-tree. 167 @param nbElements is the size of the elements array. 168 @param order is the order of interpolaton (must be 1 or 2). 169 */ 175 170 int Mapper::remap(Elt *elements, int nbElements, int order, bool renormalize, bool quantity) 176 171 { 177 int mpiSize, mpiRank; 178 MPI_Comm_size(communicator, &mpiSize); 179 MPI_Comm_rank(communicator, &mpiRank); 180 181 /* create list of intersections (super mesh elements) for each rank */ 182 multimap<int, Polyg *> *elementList = new multimap<int, Polyg *>[mpiSize]; 183 for (int j = 0; j < nbElements; j++) 184 { 185 Elt& e = elements[j]; 186 for (list<Polyg *>::iterator it = e.is.begin(); it != e.is.end(); it++) 187 elementList[(*it)->id.rank].insert(pair<int, Polyg *>((*it)->id.ind, *it)); 188 } 189 190 int *nbSendElement = new int[mpiSize]; 191 int **sendElement = new int*[mpiSize]; /* indices of elements required from other rank */ 192 double **recvValue = new double*[mpiSize]; 193 double **recvArea = new double*[mpiSize]; 194 Coord **recvGrad = new Coord*[mpiSize]; 195 GloId **recvNeighIds = new GloId*[mpiSize]; /* ids of the of the source neighbours which also contribute through gradient */ 196 for (int rank = 0; rank < mpiSize; rank++) 197 { 198 /* get size for allocation */ 199 int last = -1; /* compares unequal to any index */ 200 int index = -1; /* increased to starting index 0 in first iteration */ 201 for (multimap<int, Polyg *>::iterator it = elementList[rank].begin(); it != elementList[rank].end(); ++it) 202 { 203 if (last != it->first) 204 index++; 205 (it->second)->id.ind = index; 206 last = it->first; 207 } 208 nbSendElement[rank] = index + 1; 209 210 /* if size is non-zero allocate and collect indices of elements on other ranks that we intersect */ 211 if (nbSendElement[rank] > 0) 212 { 213 sendElement[rank] = new int[nbSendElement[rank]]; 214 recvValue[rank] = new double[nbSendElement[rank]]; 215 recvArea[rank] = new double[nbSendElement[rank]]; 216 if (order == 2) 217 { 218 recvNeighIds[rank] = new GloId[nbSendElement[rank]*(NMAX+1)]; 219 recvGrad[rank] = new Coord[nbSendElement[rank]*(NMAX+1)]; 220 } 221 else 222 recvNeighIds[rank] = new GloId[nbSendElement[rank]]; 223 224 last = -1; 225 index = -1; 226 for (multimap<int, Polyg *>::iterator it = elementList[rank].begin(); it != elementList[rank].end(); ++it) 227 { 228 if (last != it->first) 229 index++; 230 sendElement[rank][index] = it->first; 231 last = it->first; 232 } 233 } 234 } 235 236 /* communicate sizes of source elements to be sent (index lists and later values and gradients) */ 237 int *nbRecvElement = new int[mpiSize]; 238 MPI_Alltoall(nbSendElement, 1, MPI_INT, nbRecvElement, 1, MPI_INT, communicator); 239 240 /* communicate indices of source elements on other ranks whoes value and gradient we need (since intersection) */ 241 int nbSendRequest = 0; 242 int nbRecvRequest = 0; 243 int **recvElement = new int*[mpiSize]; 244 double **sendValue = new double*[mpiSize]; 245 double **sendArea = new double*[mpiSize]; 246 Coord **sendGrad = new Coord*[mpiSize]; 247 GloId **sendNeighIds = new GloId*[mpiSize]; 248 MPI_Request *sendRequest = new MPI_Request[4*mpiSize]; 249 MPI_Request *recvRequest = new MPI_Request[4*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 if (order == 2) 264 { 265 sendNeighIds[rank] = new GloId[nbRecvElement[rank]*(NMAX+1)]; 266 sendGrad[rank] = new Coord[nbRecvElement[rank]*(NMAX+1)]; 267 } 268 else 269 { 270 sendNeighIds[rank] = new GloId[nbRecvElement[rank]]; 271 } 272 MPI_Irecv(recvElement[rank], nbRecvElement[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 273 nbRecvRequest++; 274 } 275 } 276 277 MPI_Status *status = new MPI_Status[4*mpiSize]; 278 279 MPI_Waitall(nbSendRequest, sendRequest, status); 280 MPI_Waitall(nbRecvRequest, recvRequest, status); 281 282 /* for all indices that have been received from requesting ranks: pack values and gradients, then send */ 283 nbSendRequest = 0; 284 nbRecvRequest = 0; 285 for (int rank = 0; rank < mpiSize; rank++) 286 { 287 if (nbRecvElement[rank] > 0) 288 { 289 int jj = 0; // jj == j if no weight writing 290 for (int j = 0; j < nbRecvElement[rank]; j++) 291 { 292 sendValue[rank][j] = sstree.localElements[recvElement[rank][j]].val; 293 sendArea[rank][j] = sstree.localElements[recvElement[rank][j]].area; 294 if (order == 2) 295 { 296 sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].grad; 297 sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].src_id; 298 jj++; 299 for (int i = 0; i < NMAX; i++) 300 { 301 sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].gradNeigh[i]; 302 sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].neighId[i]; 303 jj++; 304 } 305 } 306 else 307 sendNeighIds[rank][j] = sstree.localElements[recvElement[rank][j]].src_id; 308 } 309 MPI_Issend(sendValue[rank], nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 310 nbSendRequest++; 311 MPI_Issend(sendArea[rank], nbRecvElement[rank], MPI_DOUBLE, rank, 1, communicator, &sendRequest[nbSendRequest]); 312 nbSendRequest++; 313 if (order == 2) 314 { 315 MPI_Issend(sendGrad[rank], 3*nbRecvElement[rank]*(NMAX+1), 316 MPI_DOUBLE, rank, 2, communicator, &sendRequest[nbSendRequest]); 317 nbSendRequest++; 318 MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank]*(NMAX+1), MPI_INT, rank, 3, communicator, &sendRequest[nbSendRequest]); 319 //ym --> attention taille GloId 320 nbSendRequest++; 321 } 322 else 323 { 324 MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank], MPI_INT, rank, 4, communicator, &sendRequest[nbSendRequest]); 325 //ym --> attention taille GloId 326 nbSendRequest++; 327 } 328 } 329 if (nbSendElement[rank] > 0) 330 { 331 MPI_Irecv(recvValue[rank], nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 332 nbRecvRequest++; 333 MPI_Irecv(recvArea[rank], nbSendElement[rank], MPI_DOUBLE, rank, 1, communicator, &recvRequest[nbRecvRequest]); 334 nbRecvRequest++; 335 if (order == 2) 336 { 337 MPI_Irecv(recvGrad[rank], 3*nbSendElement[rank]*(NMAX+1), 338 MPI_DOUBLE, rank, 2, communicator, &recvRequest[nbRecvRequest]); 339 nbRecvRequest++; 340 MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank]*(NMAX+1), MPI_INT, rank, 3, communicator, &recvRequest[nbRecvRequest]); 341 //ym --> attention taille GloId 342 nbRecvRequest++; 343 } 344 else 345 { 346 MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank], MPI_INT, rank, 4, communicator, &recvRequest[nbRecvRequest]); 347 //ym --> attention taille GloId 348 nbRecvRequest++; 349 } 350 } 351 } 352 353 MPI_Waitall(nbRecvRequest, recvRequest, status); 354 MPI_Waitall(nbSendRequest, sendRequest, status); 355 356 357 358 /* now that all values and gradients are available use them to computed interpolated values on target 359 and also to compute weights */ 360 int i = 0; 361 for (int j = 0; j < nbElements; j++) 362 { 363 Elt& e = elements[j]; 364 365 /* since for the 2nd order case source grid elements can contribute to a destination grid element over several "paths" 366 (step1: gradient is computed using neighbours on same grid, step2: intersection uses several elements on other grid) 367 accumulate them so that there is only one final weight between two elements */ 368 map<GloId,double> wgt_map; 369 370 /* for destination element `e` loop over all intersetions/the corresponding source elements */ 371 for (list<Polyg *>::iterator it = e.is.begin(); it != e.is.end(); it++) 372 { 373 /* it is the intersection element, so it->x and it->area are barycentre and area of intersection element (super mesh) 374 but it->id is id of the source element that it intersects */ 375 int n1 = (*it)->id.ind; 376 int rank = (*it)->id.rank; 377 double fk = recvValue[rank][n1]; 378 double srcArea = recvArea[rank][n1]; 379 double w = (*it)->area; 380 if (quantity) w/=srcArea ; 381 382 /* first order: src value times weight (weight = supermesh area), later divide by target area */ 383 int kk = (order == 2) ? n1 * (NMAX + 1) : n1; 384 GloId neighID = recvNeighIds[rank][kk]; 385 wgt_map[neighID] += w; 386 387 if (order == 2) 388 { 389 for (int k = 0; k < NMAX+1; k++) 390 { 391 int kk = n1 * (NMAX + 1) + k; 392 GloId neighID = recvNeighIds[rank][kk]; 393 if (neighID.ind != -1) wgt_map[neighID] += w * scalarprod(recvGrad[rank][kk], (*it)->x); 394 } 395 396 } 397 } 398 399 double renorm=0; 400 if (renormalize) 401 for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) renorm+=it->second / e.area; 402 else renorm=1. ; 403 404 for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) 405 { 406 if (quantity) this->remapMatrix[i] = (it->second ) / renorm; 407 else this->remapMatrix[i] = (it->second / e.area) / renorm; 408 this->srcAddress[i] = it->first.ind; 409 this->srcRank[i] = it->first.rank; 410 this->dstAddress[i] = j; 411 this->sourceWeightId[i]= it->first.globalId ; 412 this->targetWeightId[i]= targetGlobalId[j] ; 413 i++; 414 } 415 } 416 417 /* free all memory allocated in this function */ 418 for (int rank = 0; rank < mpiSize; rank++) 419 { 420 if (nbSendElement[rank] > 0) 421 { 422 delete[] sendElement[rank]; 423 delete[] recvValue[rank]; 424 delete[] recvArea[rank]; 425 if (order == 2) 426 { 427 delete[] recvGrad[rank]; 428 } 429 delete[] recvNeighIds[rank]; 430 } 431 if (nbRecvElement[rank] > 0) 432 { 433 delete[] recvElement[rank]; 434 delete[] sendValue[rank]; 435 delete[] sendArea[rank]; 436 if (order == 2) 437 delete[] sendGrad[rank]; 438 delete[] sendNeighIds[rank]; 439 } 440 } 441 delete[] status; 442 delete[] sendRequest; 443 delete[] recvRequest; 444 delete[] elementList; 445 delete[] nbSendElement; 446 delete[] nbRecvElement; 447 delete[] sendElement; 448 delete[] recvElement; 449 delete[] sendValue; 450 delete[] recvValue; 451 delete[] sendGrad; 452 delete[] recvGrad; 453 delete[] sendNeighIds; 454 delete[] recvNeighIds; 455 456 return i; 172 int mpiSize, mpiRank; 173 MPI_Comm_size(communicator, &mpiSize); 174 MPI_Comm_rank(communicator, &mpiRank); 175 176 /* create list of intersections (super mesh elements) for each rank */ 177 multimap<int, Polyg *> *elementList = new multimap<int, Polyg *>[mpiSize]; 178 for (int j = 0; j < nbElements; j++) 179 { 180 Elt& e = elements[j]; 181 for (list<Polyg *>::iterator it = e.is.begin(); it != e.is.end(); it++) 182 elementList[(*it)->id.rank].insert(pair<int, Polyg *>((*it)->id.ind, *it)); 183 } 184 185 int *nbSendElement = new int[mpiSize]; 186 int **sendElement = new int*[mpiSize]; /* indices of elements required from other rank */ 187 double **recvValue = new double*[mpiSize]; 188 double **recvArea = new double*[mpiSize]; 189 Coord **recvGrad = new Coord*[mpiSize]; 190 GloId **recvNeighIds = new GloId*[mpiSize]; /* ids of the of the source neighbours which also contribute through gradient */ 191 for (int rank = 0; rank < mpiSize; rank++) 192 { 193 /* get size for allocation */ 194 int last = -1; /* compares unequal to any index */ 195 int index = -1; /* increased to starting index 0 in first iteration */ 196 for (multimap<int, Polyg *>::iterator it = elementList[rank].begin(); it != elementList[rank].end(); ++it) 197 { 198 if (last != it->first) 199 index++; 200 (it->second)->id.ind = index; 201 last = it->first; 202 } 203 nbSendElement[rank] = index + 1; 204 205 /* if size is non-zero allocate and collect indices of elements on other ranks that we intersect */ 206 if (nbSendElement[rank] > 0) 207 { 208 sendElement[rank] = new int[nbSendElement[rank]]; 209 recvValue[rank] = new double[nbSendElement[rank]]; 210 recvArea[rank] = new double[nbSendElement[rank]]; 211 if (order == 2) 212 { 213 recvNeighIds[rank] = new GloId[nbSendElement[rank]*(NMAX+1)]; 214 recvGrad[rank] = new Coord[nbSendElement[rank]*(NMAX+1)]; 215 } 216 else 217 recvNeighIds[rank] = new GloId[nbSendElement[rank]]; 218 219 last = -1; 220 index = -1; 221 for (multimap<int, Polyg *>::iterator it = elementList[rank].begin(); it != elementList[rank].end(); ++it) 222 { 223 if (last != it->first) 224 index++; 225 sendElement[rank][index] = it->first; 226 last = it->first; 227 } 228 } 229 } 230 231 /* communicate sizes of source elements to be sent (index lists and later values and gradients) */ 232 int *nbRecvElement = new int[mpiSize]; 233 MPI_Alltoall(nbSendElement, 1, MPI_INT, nbRecvElement, 1, MPI_INT, communicator); 234 235 /* communicate indices of source elements on other ranks whoes value and gradient we need (since intersection) */ 236 int nbSendRequest = 0; 237 int nbRecvRequest = 0; 238 int **recvElement = new int*[mpiSize]; 239 double **sendValue = new double*[mpiSize]; 240 double **sendArea = new double*[mpiSize]; 241 Coord **sendGrad = new Coord*[mpiSize]; 242 GloId **sendNeighIds = new GloId*[mpiSize]; 243 MPI_Request *sendRequest = new MPI_Request[4*mpiSize]; 244 MPI_Request *recvRequest = new MPI_Request[4*mpiSize]; 245 for (int rank = 0; rank < mpiSize; rank++) 246 { 247 if (nbSendElement[rank] > 0) 248 { 249 MPI_Issend(sendElement[rank], nbSendElement[rank], MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 250 nbSendRequest++; 251 } 252 253 if (nbRecvElement[rank] > 0) 254 { 255 recvElement[rank] = new int[nbRecvElement[rank]]; 256 sendValue[rank] = new double[nbRecvElement[rank]]; 257 sendArea[rank] = new double[nbRecvElement[rank]]; 258 if (order == 2) 259 { 260 sendNeighIds[rank] = new GloId[nbRecvElement[rank]*(NMAX+1)]; 261 sendGrad[rank] = new Coord[nbRecvElement[rank]*(NMAX+1)]; 262 } 263 else 264 { 265 sendNeighIds[rank] = new GloId[nbRecvElement[rank]]; 266 } 267 MPI_Irecv(recvElement[rank], nbRecvElement[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 268 nbRecvRequest++; 269 } 270 } 271 272 MPI_Status *status = new MPI_Status[4*mpiSize]; 273 274 MPI_Waitall(nbSendRequest, sendRequest, status); 275 MPI_Waitall(nbRecvRequest, recvRequest, status); 276 277 /* for all indices that have been received from requesting ranks: pack values and gradients, then send */ 278 nbSendRequest = 0; 279 nbRecvRequest = 0; 280 for (int rank = 0; rank < mpiSize; rank++) 281 { 282 if (nbRecvElement[rank] > 0) 283 { 284 int jj = 0; // jj == j if no weight writing 285 for (int j = 0; j < nbRecvElement[rank]; j++) 286 { 287 sendValue[rank][j] = sstree.localElements[recvElement[rank][j]].val; 288 sendArea[rank][j] = sstree.localElements[recvElement[rank][j]].area; 289 if (order == 2) 290 { 291 sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].grad; 292 sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].src_id; 293 jj++; 294 for (int i = 0; i < NMAX; i++) 295 { 296 sendGrad[rank][jj] = sstree.localElements[recvElement[rank][j]].gradNeigh[i]; 297 sendNeighIds[rank][jj] = sstree.localElements[recvElement[rank][j]].neighId[i]; 298 jj++; 299 } 300 } 301 else 302 sendNeighIds[rank][j] = sstree.localElements[recvElement[rank][j]].src_id; 303 } 304 MPI_Issend(sendValue[rank], nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 305 nbSendRequest++; 306 MPI_Issend(sendArea[rank], nbRecvElement[rank], MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 307 nbSendRequest++; 308 if (order == 2) 309 { 310 MPI_Issend(sendGrad[rank], 3*nbRecvElement[rank]*(NMAX+1), 311 MPI_DOUBLE, rank, 0, communicator, &sendRequest[nbSendRequest]); 312 nbSendRequest++; 313 MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank]*(NMAX+1), MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 314 //ym --> attention taille GloId 315 nbSendRequest++; 316 } 317 else 318 { 319 MPI_Issend(sendNeighIds[rank], 4*nbRecvElement[rank], MPI_INT, rank, 0, communicator, &sendRequest[nbSendRequest]); 320 //ym --> attention taille GloId 321 nbSendRequest++; 322 } 323 } 324 if (nbSendElement[rank] > 0) 325 { 326 MPI_Irecv(recvValue[rank], nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 327 nbRecvRequest++; 328 MPI_Irecv(recvArea[rank], nbSendElement[rank], MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 329 nbRecvRequest++; 330 if (order == 2) 331 { 332 MPI_Irecv(recvGrad[rank], 3*nbSendElement[rank]*(NMAX+1), 333 MPI_DOUBLE, rank, 0, communicator, &recvRequest[nbRecvRequest]); 334 nbRecvRequest++; 335 MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank]*(NMAX+1), MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 336 //ym --> attention taille GloId 337 nbRecvRequest++; 338 } 339 else 340 { 341 MPI_Irecv(recvNeighIds[rank], 4*nbSendElement[rank], MPI_INT, rank, 0, communicator, &recvRequest[nbRecvRequest]); 342 //ym --> attention taille GloId 343 nbRecvRequest++; 344 } 345 } 346 } 347 348 MPI_Waitall(nbSendRequest, sendRequest, status); 349 MPI_Waitall(nbRecvRequest, recvRequest, status); 350 351 352 /* now that all values and gradients are available use them to computed interpolated values on target 353 and also to compute weights */ 354 int i = 0; 355 for (int j = 0; j < nbElements; j++) 356 { 357 Elt& e = elements[j]; 358 359 /* since for the 2nd order case source grid elements can contribute to a destination grid element over several "paths" 360 (step1: gradient is computed using neighbours on same grid, step2: intersection uses several elements on other grid) 361 accumulate them so that there is only one final weight between two elements */ 362 map<GloId,double> wgt_map; 363 364 /* for destination element `e` loop over all intersetions/the corresponding source elements */ 365 for (list<Polyg *>::iterator it = e.is.begin(); it != e.is.end(); it++) 366 { 367 /* it is the intersection element, so it->x and it->area are barycentre and area of intersection element (super mesh) 368 but it->id is id of the source element that it intersects */ 369 int n1 = (*it)->id.ind; 370 int rank = (*it)->id.rank; 371 double fk = recvValue[rank][n1]; 372 double srcArea = recvArea[rank][n1]; 373 double w = (*it)->area; 374 if (quantity) w/=srcArea ; 375 376 /* first order: src value times weight (weight = supermesh area), later divide by target area */ 377 int kk = (order == 2) ? n1 * (NMAX + 1) : n1; 378 GloId neighID = recvNeighIds[rank][kk]; 379 wgt_map[neighID] += w; 380 381 if (order == 2) 382 { 383 for (int k = 0; k < NMAX+1; k++) 384 { 385 int kk = n1 * (NMAX + 1) + k; 386 GloId neighID = recvNeighIds[rank][kk]; 387 if (neighID.ind != -1) wgt_map[neighID] += w * scalarprod(recvGrad[rank][kk], (*it)->x); 388 } 389 390 } 391 } 392 393 double renorm=0; 394 if (renormalize) 395 for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) renorm+=it->second / e.area; 396 else renorm=1. ; 397 398 for (map<GloId,double>::iterator it = wgt_map.begin(); it != wgt_map.end(); it++) 399 { 400 if (quantity) this->remapMatrix[i] = (it->second ) / renorm; 401 else this->remapMatrix[i] = (it->second / e.area) / renorm; 402 this->srcAddress[i] = it->first.ind; 403 this->srcRank[i] = it->first.rank; 404 this->dstAddress[i] = j; 405 this->sourceWeightId[i]= it->first.globalId ; 406 this->targetWeightId[i]= targetGlobalId[j] ; 407 i++; 408 } 409 } 410 411 /* free all memory allocated in this function */ 412 for (int rank = 0; rank < mpiSize; rank++) 413 { 414 if (nbSendElement[rank] > 0) 415 { 416 delete[] sendElement[rank]; 417 delete[] recvValue[rank]; 418 delete[] recvArea[rank]; 419 if (order == 2) 420 { 421 delete[] recvGrad[rank]; 422 } 423 delete[] recvNeighIds[rank]; 424 } 425 if (nbRecvElement[rank] > 0) 426 { 427 delete[] recvElement[rank]; 428 delete[] sendValue[rank]; 429 delete[] sendArea[rank]; 430 if (order == 2) 431 delete[] sendGrad[rank]; 432 delete[] sendNeighIds[rank]; 433 } 434 } 435 delete[] status; 436 delete[] sendRequest; 437 delete[] recvRequest; 438 delete[] elementList; 439 delete[] nbSendElement; 440 delete[] nbRecvElement; 441 delete[] sendElement; 442 delete[] recvElement; 443 delete[] sendValue; 444 delete[] recvValue; 445 delete[] sendGrad; 446 delete[] recvGrad; 447 delete[] sendNeighIds; 448 delete[] recvNeighIds; 449 return i; 457 450 } 458 451 459 452 void Mapper::computeGrads() 460 453 { 461 462 463 464 465 466 467 468 469 470 454 /* array of pointers to collect local elements and elements received from other cpu */ 455 vector<Elt*> globalElements(sstree.nbLocalElements + nbNeighbourElements); 456 int index = 0; 457 for (int i = 0; i < sstree.nbLocalElements; i++, index++) 458 globalElements[index] = &(sstree.localElements[i]); 459 for (int i = 0; i < nbNeighbourElements; i++, index++) 460 globalElements[index] = &neighbourElements[i]; 461 462 update_baryc(sstree.localElements, sstree.nbLocalElements); 463 computeGradients(&globalElements[0], sstree.nbLocalElements); 471 464 } 472 465 473 466 /** for each element of the source grid, finds all the neighbouring elements that share an edge 474 (filling array neighbourElements). This is used later to compute gradients */467 (filling array neighbourElements). This is used later to compute gradients */ 475 468 void Mapper::buildMeshTopology() 476 469 { 477 478 479 480 481 482 483 484 485 486 487 488 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 665 666 667 668 669 670 671 672 intersector of a local node */673 674 675 676 use the SS-tree to find all elements (including neighbourElements)677 who are potential neighbours because their circles intersect,678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 whoes circles intersect with the circle around `elt` (the SS intersectors)693 and check if they are neighbours in the sense that the two elements share an edge.694 If they do, save this information for elt */695 696 697 698 699 700 701 470 int mpiSize, mpiRank; 471 MPI_Comm_size(communicator, &mpiSize); 472 MPI_Comm_rank(communicator, &mpiRank); 473 474 vector<Node> *routingList = new vector<Node>[mpiSize]; 475 vector<vector<int> > routes(sstree.localTree.leafs.size()); 476 477 sstree.routeIntersections(routes, sstree.localTree.leafs); 478 479 for (int i = 0; i < routes.size(); ++i) 480 for (int k = 0; k < routes[i].size(); ++k) 481 routingList[routes[i][k]].push_back(sstree.localTree.leafs[i]); 482 routingList[mpiRank].clear(); 483 484 485 CMPIRouting mpiRoute(communicator); 486 mpiRoute.init(routes); 487 int nRecv = mpiRoute.getTotalSourceElement(); 488 489 int *nbSendNode = new int[mpiSize]; 490 int *nbRecvNode = new int[mpiSize]; 491 int *sendMessageSize = new int[mpiSize]; 492 int *recvMessageSize = new int[mpiSize]; 493 494 for (int rank = 0; rank < mpiSize; rank++) 495 { 496 nbSendNode[rank] = routingList[rank].size(); 497 sendMessageSize[rank] = 0; 498 for (size_t j = 0; j < routingList[rank].size(); j++) 499 { 500 Elt *elt = (Elt *) (routingList[rank][j].data); 501 sendMessageSize[rank] += packedPolygonSize(*elt); 502 } 503 } 504 505 MPI_Alltoall(nbSendNode, 1, MPI_INT, nbRecvNode, 1, MPI_INT, communicator); 506 MPI_Alltoall(sendMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator); 507 508 char **sendBuffer = new char*[mpiSize]; 509 char **recvBuffer = new char*[mpiSize]; 510 int *pos = new int[mpiSize]; 511 512 for (int rank = 0; rank < mpiSize; rank++) 513 { 514 if (nbSendNode[rank] > 0) sendBuffer[rank] = new char[sendMessageSize[rank]]; 515 if (nbRecvNode[rank] > 0) recvBuffer[rank] = new char[recvMessageSize[rank]]; 516 } 517 518 for (int rank = 0; rank < mpiSize; rank++) 519 { 520 pos[rank] = 0; 521 for (size_t j = 0; j < routingList[rank].size(); j++) 522 { 523 Elt *elt = (Elt *) (routingList[rank][j].data); 524 packPolygon(*elt, sendBuffer[rank], pos[rank]); 525 } 526 } 527 delete [] routingList; 528 529 530 int nbSendRequest = 0; 531 int nbRecvRequest = 0; 532 MPI_Request *sendRequest = new MPI_Request[mpiSize]; 533 MPI_Request *recvRequest = new MPI_Request[mpiSize]; 534 MPI_Status *status = new MPI_Status[mpiSize]; 535 536 for (int rank = 0; rank < mpiSize; rank++) 537 { 538 if (nbSendNode[rank] > 0) 539 { 540 MPI_Issend(sendBuffer[rank], sendMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]); 541 nbSendRequest++; 542 } 543 if (nbRecvNode[rank] > 0) 544 { 545 MPI_Irecv(recvBuffer[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]); 546 nbRecvRequest++; 547 } 548 } 549 550 MPI_Waitall(nbRecvRequest, recvRequest, status); 551 MPI_Waitall(nbSendRequest, sendRequest, status); 552 553 for (int rank = 0; rank < mpiSize; rank++) 554 if (nbSendNode[rank] > 0) delete [] sendBuffer[rank]; 555 delete [] sendBuffer; 556 557 char **sendBuffer2 = new char*[mpiSize]; 558 char **recvBuffer2 = new char*[mpiSize]; 559 560 for (int rank = 0; rank < mpiSize; rank++) 561 { 562 nbSendNode[rank] = 0; 563 sendMessageSize[rank] = 0; 564 565 if (nbRecvNode[rank] > 0) 566 { 567 set<NodePtr> neighbourList; 568 pos[rank] = 0; 569 for (int j = 0; j < nbRecvNode[rank]; j++) 570 { 571 Elt elt; 572 unpackPolygon(elt, recvBuffer[rank], pos[rank]); 573 Node node(elt.x, cptRadius(elt), &elt); 574 findNeighbour(sstree.localTree.root, &node, neighbourList); 575 } 576 nbSendNode[rank] = neighbourList.size(); 577 for (set<NodePtr>::iterator it = neighbourList.begin(); it != neighbourList.end(); it++) 578 { 579 Elt *elt = (Elt *) ((*it)->data); 580 sendMessageSize[rank] += packedPolygonSize(*elt); 581 } 582 583 sendBuffer2[rank] = new char[sendMessageSize[rank]]; 584 pos[rank] = 0; 585 586 for (set<NodePtr>::iterator it = neighbourList.begin(); it != neighbourList.end(); it++) 587 { 588 Elt *elt = (Elt *) ((*it)->data); 589 packPolygon(*elt, sendBuffer2[rank], pos[rank]); 590 } 591 } 592 } 593 for (int rank = 0; rank < mpiSize; rank++) 594 if (nbRecvNode[rank] > 0) delete [] recvBuffer[rank]; 595 delete [] recvBuffer; 596 597 598 MPI_Barrier(communicator); 599 MPI_Alltoall(nbSendNode, 1, MPI_INT, nbRecvNode, 1, MPI_INT, communicator); 600 MPI_Alltoall(sendMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator); 601 602 for (int rank = 0; rank < mpiSize; rank++) 603 if (nbRecvNode[rank] > 0) recvBuffer2[rank] = new char[recvMessageSize[rank]]; 604 605 nbSendRequest = 0; 606 nbRecvRequest = 0; 607 608 for (int rank = 0; rank < mpiSize; rank++) 609 { 610 if (nbSendNode[rank] > 0) 611 { 612 MPI_Issend(sendBuffer2[rank], sendMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]); 613 nbSendRequest++; 614 } 615 if (nbRecvNode[rank] > 0) 616 { 617 MPI_Irecv(recvBuffer2[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]); 618 nbRecvRequest++; 619 } 620 } 621 622 MPI_Waitall(nbRecvRequest, recvRequest, status); 623 MPI_Waitall(nbSendRequest, sendRequest, status); 624 625 int nbNeighbourNodes = 0; 626 for (int rank = 0; rank < mpiSize; rank++) 627 nbNeighbourNodes += nbRecvNode[rank]; 628 629 neighbourElements = new Elt[nbNeighbourNodes]; 630 nbNeighbourElements = nbNeighbourNodes; 631 632 int index = 0; 633 for (int rank = 0; rank < mpiSize; rank++) 634 { 635 pos[rank] = 0; 636 for (int j = 0; j < nbRecvNode[rank]; j++) 637 { 638 unpackPolygon(neighbourElements[index], recvBuffer2[rank], pos[rank]); 639 neighbourElements[index].id.ind = sstree.localTree.leafs.size() + index; 640 index++; 641 } 642 } 643 for (int rank = 0; rank < mpiSize; rank++) 644 { 645 if (nbRecvNode[rank] > 0) delete [] recvBuffer2[rank]; 646 if (nbSendNode[rank] > 0) delete [] sendBuffer2[rank]; 647 } 648 delete [] recvBuffer2; 649 delete [] sendBuffer2; 650 delete [] sendMessageSize; 651 delete [] recvMessageSize; 652 delete [] nbSendNode; 653 delete [] nbRecvNode; 654 delete [] sendRequest; 655 delete [] recvRequest; 656 delete [] status; 657 delete [] pos; 658 659 /* re-compute on received elements to avoid having to send this information */ 660 neighbourNodes.resize(nbNeighbourNodes); 661 setCirclesAndLinks(neighbourElements, neighbourNodes); 662 cptAllEltsGeom(neighbourElements, nbNeighbourNodes, srcGrid.pole); 663 664 /* the local SS tree must include nodes from other cpus if they are potential 665 intersector of a local node */ 666 sstree.localTree.insertNodes(neighbourNodes); 667 668 /* for every local element, 669 use the SS-tree to find all elements (including neighbourElements) 670 who are potential neighbours because their circles intersect, 671 then check all canditates for common edges to build up connectivity information 672 */ 673 for (int j = 0; j < sstree.localTree.leafs.size(); j++) 674 { 675 Node& node = sstree.localTree.leafs[j]; 676 677 /* find all leafs whoes circles that intersect node's circle and save into node->intersectors */ 678 node.search(sstree.localTree.root); 679 680 Elt *elt = (Elt *)(node.data); 681 682 for (int i = 0; i < elt->n; i++) elt->neighbour[i] = NOT_FOUND; 683 684 /* for element `elt` loop through all nodes in the SS-tree 685 whoes circles intersect with the circle around `elt` (the SS intersectors) 686 and check if they are neighbours in the sense that the two elements share an edge. 687 If they do, save this information for elt */ 688 for (list<NodePtr>::iterator it = (node.intersectors).begin(); it != (node.intersectors).end(); ++it) 689 { 690 Elt *elt2 = (Elt *)((*it)->data); 691 set_neighbour(*elt, *elt2); 692 } 693 694 } 702 695 } 703 696 … … 705 698 void Mapper::computeIntersection(Elt *elements, int nbElements) 706 699 { 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 //intersect_ym(&recvElt[j], elt2);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 700 int mpiSize, mpiRank; 701 MPI_Comm_size(communicator, &mpiSize); 702 MPI_Comm_rank(communicator, &mpiRank); 703 704 MPI_Barrier(communicator); 705 706 vector<Node> *routingList = new vector<Node>[mpiSize]; 707 708 vector<Node> routeNodes; routeNodes.reserve(nbElements); 709 for (int j = 0; j < nbElements; j++) 710 { 711 elements[j].id.ind = j; 712 elements[j].id.rank = mpiRank; 713 routeNodes.push_back(Node(elements[j].x, cptRadius(elements[j]), &elements[j])); 714 } 715 716 vector<vector<int> > routes(routeNodes.size()); 717 sstree.routeIntersections(routes, routeNodes); 718 for (int i = 0; i < routes.size(); ++i) 719 for (int k = 0; k < routes[i].size(); ++k) 720 routingList[routes[i][k]].push_back(routeNodes[i]); 721 722 if (verbose >= 2) 723 { 724 cout << " --> rank " << mpiRank << " nbElements " << nbElements << " : "; 725 for (int rank = 0; rank < mpiSize; rank++) 726 cout << routingList[rank].size() << " "; 727 cout << endl; 728 } 729 MPI_Barrier(communicator); 730 731 int *nbSendNode = new int[mpiSize]; 732 int *nbRecvNode = new int[mpiSize]; 733 int *sentMessageSize = new int[mpiSize]; 734 int *recvMessageSize = new int[mpiSize]; 735 736 for (int rank = 0; rank < mpiSize; rank++) 737 { 738 nbSendNode[rank] = routingList[rank].size(); 739 sentMessageSize[rank] = 0; 740 for (size_t j = 0; j < routingList[rank].size(); j++) 741 { 742 Elt *elt = (Elt *) (routingList[rank][j].data); 743 sentMessageSize[rank] += packedPolygonSize(*elt); 744 } 745 } 746 747 MPI_Alltoall(nbSendNode, 1, MPI_INT, nbRecvNode, 1, MPI_INT, communicator); 748 MPI_Alltoall(sentMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator); 749 750 int total = 0; 751 752 for (int rank = 0; rank < mpiSize; rank++) 753 { 754 total = total + nbRecvNode[rank]; 755 } 756 757 if (verbose >= 2) cout << "---> rank " << mpiRank << " : compute intersection : total received nodes " << total << endl; 758 char **sendBuffer = new char*[mpiSize]; 759 char **recvBuffer = new char*[mpiSize]; 760 int *pos = new int[mpiSize]; 761 762 for (int rank = 0; rank < mpiSize; rank++) 763 { 764 if (nbSendNode[rank] > 0) sendBuffer[rank] = new char[sentMessageSize[rank]]; 765 if (nbRecvNode[rank] > 0) recvBuffer[rank] = new char[recvMessageSize[rank]]; 766 } 767 768 for (int rank = 0; rank < mpiSize; rank++) 769 { 770 pos[rank] = 0; 771 for (size_t j = 0; j < routingList[rank].size(); j++) 772 { 773 Elt* elt = (Elt *) (routingList[rank][j].data); 774 packPolygon(*elt, sendBuffer[rank], pos[rank]); 775 } 776 } 777 delete [] routingList; 778 779 int nbSendRequest = 0; 780 int nbRecvRequest = 0; 781 MPI_Request *sendRequest = new MPI_Request[mpiSize]; 782 MPI_Request *recvRequest = new MPI_Request[mpiSize]; 783 MPI_Status *status = new MPI_Status[mpiSize]; 784 785 for (int rank = 0; rank < mpiSize; rank++) 786 { 787 if (nbSendNode[rank] > 0) 788 { 789 MPI_Issend(sendBuffer[rank], sentMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]); 790 nbSendRequest++; 791 } 792 if (nbRecvNode[rank] > 0) 793 { 794 MPI_Irecv(recvBuffer[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]); 795 nbRecvRequest++; 796 } 797 } 798 799 MPI_Waitall(nbRecvRequest, recvRequest, status); 800 MPI_Waitall(nbSendRequest, sendRequest, status); 801 802 char **sendBuffer2 = new char*[mpiSize]; 803 char **recvBuffer2 = new char*[mpiSize]; 804 805 double tic = cputime(); 806 for (int rank = 0; rank < mpiSize; rank++) 807 { 808 sentMessageSize[rank] = 0; 809 810 if (nbRecvNode[rank] > 0) 811 { 812 Elt *recvElt = new Elt[nbRecvNode[rank]]; 813 pos[rank] = 0; 814 for (int j = 0; j < nbRecvNode[rank]; j++) 815 { 816 unpackPolygon(recvElt[j], recvBuffer[rank], pos[rank]); 817 cptEltGeom(recvElt[j], tgtGrid.pole); 818 Node recvNode(recvElt[j].x, cptRadius(recvElt[j]), &recvElt[j]); 819 recvNode.search(sstree.localTree.root); 820 /* for a node holding an element of the target, loop throught candidates for intersecting source */ 821 for (list<NodePtr>::iterator it = (recvNode.intersectors).begin(); it != (recvNode.intersectors).end(); ++it) 822 { 823 Elt *elt2 = (Elt *) ((*it)->data); 824 /* recvElt is target, elt2 is source */ 825 // intersect(&recvElt[j], elt2); 826 intersect_ym(&recvElt[j], elt2); 827 } 828 if (recvElt[j].is.size() > 0) sentMessageSize[rank] += packIntersectionSize(recvElt[j]); 829 830 // here recvNode goes out of scope 831 } 832 833 if (sentMessageSize[rank] > 0) 834 { 835 sentMessageSize[rank] += sizeof(int); 836 sendBuffer2[rank] = new char[sentMessageSize[rank]]; 837 *((int *) sendBuffer2[rank]) = 0; 838 pos[rank] = sizeof(int); 839 for (int j = 0; j < nbRecvNode[rank]; j++) 840 { 841 packIntersection(recvElt[j], sendBuffer2[rank], pos[rank]); 842 //FIXME should be deleted: recvElt[j].delete_intersections(); // intersection areas have been packed to buffer and won't be used any more 843 } 844 } 845 delete [] recvElt; 846 } 847 } 848 delete [] pos; 849 850 for (int rank = 0; rank < mpiSize; rank++) 851 { 852 if (nbSendNode[rank] > 0) delete [] sendBuffer[rank]; 853 if (nbRecvNode[rank] > 0) delete [] recvBuffer[rank]; 854 nbSendNode[rank] = 0; 855 } 856 857 if (verbose >= 2) cout << "Rank " << mpiRank << " Compute (internal) intersection " << cputime() - tic << " s" << endl; 858 MPI_Alltoall(sentMessageSize, 1, MPI_INT, recvMessageSize, 1, MPI_INT, communicator); 859 860 for (int rank = 0; rank < mpiSize; rank++) 861 if (recvMessageSize[rank] > 0) 862 recvBuffer2[rank] = new char[recvMessageSize[rank]]; 863 864 nbSendRequest = 0; 865 nbRecvRequest = 0; 866 867 for (int rank = 0; rank < mpiSize; rank++) 868 { 869 if (sentMessageSize[rank] > 0) 870 { 871 MPI_Issend(sendBuffer2[rank], sentMessageSize[rank], MPI_CHAR, rank, 0, communicator, &sendRequest[nbSendRequest]); 872 nbSendRequest++; 873 } 874 if (recvMessageSize[rank] > 0) 875 { 876 MPI_Irecv(recvBuffer2[rank], recvMessageSize[rank], MPI_CHAR, rank, 0, communicator, &recvRequest[nbRecvRequest]); 877 nbRecvRequest++; 878 } 879 } 880 881 MPI_Waitall(nbRecvRequest, recvRequest, status); 882 MPI_Waitall(nbSendRequest, sendRequest, status); 883 884 delete [] sendRequest; 885 delete [] recvRequest; 886 delete [] status; 887 for (int rank = 0; rank < mpiSize; rank++) 888 { 889 if (nbRecvNode[rank] > 0) 890 { 891 if (sentMessageSize[rank] > 0) 892 delete [] sendBuffer2[rank]; 893 } 894 895 if (recvMessageSize[rank] > 0) 896 { 897 unpackIntersection(elements, recvBuffer2[rank]); 898 delete [] recvBuffer2[rank]; 899 } 900 } 901 delete [] sendBuffer2; 902 delete [] recvBuffer2; 903 delete [] sendBuffer; 904 delete [] recvBuffer; 905 906 delete [] nbSendNode; 907 delete [] nbRecvNode; 908 delete [] sentMessageSize; 909 delete [] recvMessageSize; 917 910 } 918 911 919 912 Mapper::~Mapper() 920 913 { 921 922 923 924 925 926 } 927 928 } 914 delete [] remapMatrix; 915 delete [] srcAddress; 916 delete [] srcRank; 917 delete [] dstAddress; 918 if (neighbourElements) delete [] neighbourElements; 919 } 920 921 } -
XIOS/dev/branch_openmp/extern/remap/src/mapper.hpp
r1134 r1328 3 3 #include "parallel_tree.hpp" 4 4 #include "mpi.hpp" 5 6 #ifdef _usingEP7 #include "ep_declaration.hpp"8 #endif9 5 10 6 namespace sphereRemap { … … 22 18 { 23 19 public: 24 Mapper(ep_lib::MPI_Comm comm =MPI_COMM_WORLD) : communicator(comm), verbose(SILENT), neighbourElements(NULL), sstree(comm) {}20 Mapper(ep_lib::MPI_Comm comm) : communicator(comm), verbose(SILENT), neighbourElements(NULL), sstree(comm) {} 25 21 ~Mapper(); 26 22 void setVerbosity(verbosity v) {verbose=v ;} -
XIOS/dev/branch_openmp/extern/remap/src/mpi_cascade.cpp
r688 r1328 1 1 #include "mpi_cascade.hpp" 2 2 #include <iostream> 3 using namespace ep_lib; 3 4 4 5 namespace sphereRemap { -
XIOS/dev/branch_openmp/extern/remap/src/mpi_cascade.hpp
r694 r1328 12 12 { 13 13 public: 14 CCascadeLevel(MPI_Comm comm) : comm(comm)15 16 17 18 19 20 14 CCascadeLevel(ep_lib::MPI_Comm comm) : comm(comm) 15 { 16 ep_lib::MPI_Comm_size(comm, &size); 17 ep_lib::MPI_Comm_rank(comm, &rank); 18 } 19 int colour() const { return rank % group_size; }; 20 int key() const { return p_colour() + rank/(p_grp_size*group_size)*p_grp_size; } 21 21 22 23 24 22 // perpendicular group 23 int p_colour() const { return (rank%group_size + rank/group_size) % p_grp_size; } 24 int p_key() const { return colour() + rank/(p_grp_size*group_size)*group_size; } 25 25 26 27 28 29 30 26 ep_lib::MPI_Comm comm, pg_comm; 27 int rank; 28 int size; 29 int group_size; // group_size and p_grp_size are interchanged?? FIXME 30 int p_grp_size; 31 31 }; 32 32 … … 34 34 { 35 35 public: 36 // 37 CMPICascade(int nodes_per_level, MPI_Comm comm); 36 CMPICascade(int nodes_per_level, ep_lib::MPI_Comm comm); 38 37 39 40 38 int num_levels; 39 std::vector<CCascadeLevel> level; 41 40 }; 42 41 -
XIOS/dev/branch_openmp/extern/remap/src/mpi_routing.cpp
r1289 r1328 5 5 #include "timerRemap.hpp" 6 6 #include <iostream> 7 #ifdef _usingEP 8 #include "ep_declaration.hpp" 9 #endif 7 using namespace ep_lib; 10 8 11 9 namespace sphereRemap { … … 125 123 CTimer::get("CMPIRouting::init(reduce_scatter)").print(); 126 124 127 MPI_Info info_null;128 129 // MPI_Alloc_mem(nbTarget *sizeof(int), info_null, &targetRank);130 // MPI_Alloc_mem(nbSource *sizeof(int), info_null, &sourceRank);131 125 MPI_Alloc_mem(nbTarget *sizeof(int), MPI_INFO_NULL, &targetRank); 132 126 MPI_Alloc_mem(nbSource *sizeof(int), MPI_INFO_NULL, &sourceRank); … … 158 152 { 159 153 #ifdef _usingEP 160 MPI_Irecv(&sourceRank[i], 1, MPI_INT, - 1, 0, communicator, &request[indexRequest]);154 MPI_Irecv(&sourceRank[i], 1, MPI_INT, -2, 0, communicator, &request[indexRequest]); 161 155 #else 162 156 MPI_Irecv(&sourceRank[i], 1, MPI_INT, MPI_ANY_SOURCE, 0, communicator, &request[indexRequest]); … … 182 176 { 183 177 #ifdef _usingEP 184 MPI_Irecv(&sourceRank[i], 1, MPI_INT, - 1, 0, communicator, &request[indexRequest]);178 MPI_Irecv(&sourceRank[i], 1, MPI_INT, -2, 0, communicator, &request[indexRequest]); 185 179 #else 186 180 MPI_Irecv(&sourceRank[i], 1, MPI_INT, MPI_ANY_SOURCE, 0, communicator, &request[indexRequest]); -
XIOS/dev/branch_openmp/extern/remap/src/mpi_routing.hpp
r694 r1328 11 11 { 12 12 13 MPI_Comm communicator;13 ep_lib::MPI_Comm communicator; 14 14 int mpiRank; 15 15 int mpiSize; … … 29 29 30 30 public: 31 CMPIRouting( MPI_Comm comm);31 CMPIRouting(ep_lib::MPI_Comm comm); 32 32 ~CMPIRouting(); 33 33 template<typename T> void init(const std::vector<T>& route, CMPICascade *cascade = NULL); … … 44 44 template <typename T> 45 45 void alltoalls_known(const std::vector<std::vector<T> >& send, std::vector<std::vector<T> >& recv, 46 const std::vector<int>& ranks, MPI_Comm communicator);46 const std::vector<int>& ranks, ep_lib::MPI_Comm communicator); 47 47 48 48 template <typename T> 49 49 void alltoalls_unknown(const std::vector<std::vector<T> >& send, std::vector<std::vector<T> >& recv, 50 const std::vector<int>& ranks, MPI_Comm communicator);50 const std::vector<int>& ranks, ep_lib::MPI_Comm communicator); 51 51 } 52 52 #endif -
XIOS/dev/branch_openmp/extern/remap/src/node.cpp
r1205 r1328 254 254 NodePtr insert(NodePtr thIs, NodePtr node) 255 255 { 256 257 258 259 260 261 262 263 264 256 int la = thIs->level; // node to be inserted 257 int lb = node->level; // node where insertation 258 assert(la < lb); // node to be inserted must have lower level then parent 259 //if (thIs->parent) assert(find_in_tree1(thIs) == true); 260 NodePtr q = NULL; 261 NodePtr chd = NULL; 262 node->move(thIs); 263 if (la == lb - 1) 264 { 265 265 node->child.push_back(thIs); 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 266 thIs->parent = node; 267 if (node->child.size() > MAX_NODE_SZ && node->tree->canSplit() ) // with us as additional child `node` is now too large :( 268 return (node->reinserted || node->parent == NULL) ? split(node) : reinsert(node); 269 } 270 else // la < lb - 1 271 { 272 chd = thIs->closest(node->child); 273 q = insert(thIs, chd); 274 } 275 if ((node->updateCount + 1) % UPDATE_EVERY == 0) 276 node->update(); 277 else 278 { 279 if (q) node->remove(q); 280 node->inflate(chd); 281 } 282 282 283 283 return q; -
XIOS/dev/branch_openmp/extern/remap/src/node.hpp
r1153 r1328 15 15 struct Circle 16 16 { 17 18 17 Coord centre; 18 double radius; 19 19 }; 20 20 … … 116 116 struct Node 117 117 { 118 119 120 121 122 123 124 125 126 127 128 129 118 int level; /* FIXME leafs are 0 and root is max level? */ 119 int leafCount; /* number of leafs that are descendants of this node (the elements in it's cycle) */ 120 Coord centre; 121 double radius; 122 NodePtr parent, ref; 123 std::vector<NodePtr> child; 124 std::list<NodePtr> intersectors; 125 bool reinserted; 126 int updateCount; // double var; 127 CBasicTree* tree; 128 void *data; 129 int route; 130 130 bool toDelete ; 131 131 132 133 134 132 Node() : level(0), leafCount(1), centre(ORIGIN), radius(0), reinserted(false), updateCount(0), toDelete(false) {} 133 Node(const Coord& centre, double radius, void *data) 134 : level(0), leafCount(1), centre(centre), radius(radius), reinserted(false), updateCount(0), data(data), toDelete(false) {} 135 135 136 136 //#ifdef DEBUG … … 178 178 //#endif 179 179 180 181 182 183 180 void move(const NodePtr node); 181 void remove(const NodePtr node); 182 void inflate(const NodePtr node); 183 void update(); 184 184 void output(std::ostream& flux, int level, int color) ; 185 186 187 188 189 190 191 192 193 185 NodePtr closest(std::vector<NodePtr>& list, int n = CLOSEST); 186 NodePtr farthest(std::vector<NodePtr>& list); 187 void findClosest(int level, NodePtr src, double& minDist, NodePtr &closest); 188 189 void search(NodePtr node); 190 bool centreInside(Node &node); 191 bool intersects(NodePtr node); 192 bool isInside(Node &node); 193 int incluCheck(); 194 194 void checkParent(void) ; 195 196 197 198 199 200 201 195 void printChildren(); 196 void assignRoute(std::vector<int>::iterator& rank, int level); 197 void assignCircleAndPropagateUp(Coord *centres, double *radia, int level); 198 void printLevel(int level); 199 void routeNode(NodePtr node, int level); 200 void routingIntersecting(std::vector<Node>* routingList, NodePtr node); 201 void routeIntersection(std::vector<int>& routes, NodePtr node); 202 202 void getNodeLevel(int level,std::list<NodePtr>& NodeList) ; 203 203 bool removeDeletedNodes(int assignLevel) ; -
XIOS/dev/branch_openmp/extern/remap/src/parallel_tree.cpp
r1295 r1328 12 12 13 13 #include "parallel_tree.hpp" 14 using namespace ep_lib; 14 15 15 16 namespace sphereRemap { 16 17 extern CRemapGrid srcGrid;18 #pragma omp threadprivate(srcGrid)19 20 extern CRemapGrid tgtGrid;21 #pragma omp threadprivate(tgtGrid)22 17 23 18 static const int assignLevel = 2; … … 292 287 { 293 288 294 295 289 int assignLevel = 2; 290 int nbSampleNodes = 2*ipow(MAX_NODE_SZ + 1, assignLevel); 296 291 297 292 … … 305 300 MPI_Comm_size(communicator,&commSize) ; 306 301 307 308 309 310 //assert( nbTot > nbSampleNodes*commSize) ;302 // make multiple of two 303 nbSampleNodes /= 2; 304 nbSampleNodes *= 2; 305 // assert( nbTot > nbSampleNodes*commSize) ; 311 306 312 307 int nbSampleNodes1 = nbSampleNodes * (nb1*commSize)/(1.*nbTot) ; … … 314 309 315 310 316 //assert(node.size() > nbSampleNodes);317 //assert(node2.size() > nbSampleNodes);318 //assert(node.size() + node2.size() > nbSampleNodes);319 320 321 322 323 324 325 326 for (int i = 0; i <nbSampleNodes1; i++) sampleNodes.push_back(Node(node[randomArray1[i%nb1]].centre, node[randomArray1[i%nb1]].radius, NULL));327 for (int i = 0; i <nbSampleNodes2; i++) sampleNodes.push_back(Node(node2[randomArray2[i%nb2]].centre, node2[randomArray2[i%nb2]].radius, NULL));311 // assert(node.size() > nbSampleNodes); 312 // assert(node2.size() > nbSampleNodes); 313 // assert(node.size() + node2.size() > nbSampleNodes); 314 vector<Node> sampleNodes; sampleNodes.reserve(nbSampleNodes1+nbSampleNodes2); 315 316 vector<int> randomArray1(node.size()); 317 randomizeArray(randomArray1); 318 vector<int> randomArray2(node2.size()); 319 randomizeArray(randomArray2); 320 321 for (int i = 0; i <nbSampleNodes1; i++) sampleNodes.push_back(Node(node[randomArray1[i%nb1]].centre, node[randomArray1[i%nb1]].radius, NULL)); 322 for (int i = 0; i <nbSampleNodes2; i++) sampleNodes.push_back(Node(node2[randomArray2[i%nb2]].centre, node2[randomArray2[i%nb2]].radius, NULL)); 328 323 329 324 CTimer::get("buildParallelSampleTree").resume(); … … 336 331 CTimer::get("parallelRouteNode").resume(); 337 332 vector<int> route(node.size()); 333 cout<<"node.size = "<<node.size()<<endl; 338 334 routeNodes(route /*out*/, node); 339 335 CTimer::get("parallelRouteNode").suspend(); -
XIOS/dev/branch_openmp/extern/remap/src/parallel_tree.hpp
r1134 r1328 6 6 #include "mpi_cascade.hpp" 7 7 #include "mpi.hpp" 8 #ifdef _usingEP9 #include "ep_declaration.hpp"10 #endif11 12 8 namespace sphereRemap { 13 9 … … 15 11 { 16 12 public: 17 18 13 CParallelTree(ep_lib::MPI_Comm comm); 14 ~CParallelTree(); 19 15 20 16 void build(vector<Node>& node, vector<Node>& node2); 21 17 22 23 18 void routeNodes(vector<int>& route, vector<Node>& nodes, int level = 0); 19 void routeIntersections(vector<vector<int> >& route, vector<Node>& nodes, int level = 0); 24 20 25 26 21 int nbLocalElements; 22 Elt* localElements; 27 23 28 24 CTree localTree; 29 25 30 26 private: 31 32 33 34 27 void updateCirclesForRouting(Coord rootCentre, double rootRadius, int level = 0); 28 void buildSampleTreeCascade(vector<Node>& sampleNodes, int level = 0); 29 void buildLocalTree(const vector<Node>& node, const vector<int>& route); 30 void buildRouteTree(); 35 31 36 37 38 39 32 //CSampleTree sampleTree; 33 vector<CSampleTree> treeCascade; // first for sample tree, then for routing tree 34 CMPICascade cascade; 35 ep_lib::MPI_Comm communicator ; 40 36 41 37 }; -
XIOS/dev/branch_openmp/extern/remap/src/polyg.cpp
r1289 r1328 7 7 8 8 #include "polyg.hpp" 9 #include <stdio.h>10 9 11 10 namespace sphereRemap { … … 46 45 Coord barycentre(const Coord *x, int n) 47 46 { 48 49 47 if (n == 0) return ORIGIN; 50 48 Coord bc = ORIGIN; 51 49 for (int i = 0; i < n; i++) 52 {53 50 bc = bc + x[i]; 54 }55 51 /* both distances can be equal down to roundoff when norm(bc) < mashineepsilon 56 52 which can occur when weighted with tiny area */ 57 58 59 //assert(squaredist(bc, proj(bc)) <= squaredist(bc, proj(bc * (-1.0)))); 53 54 assert(squaredist(bc, proj(bc)) <= squaredist(bc, proj(bc * (-1.0)))); 55 //if (squaredist(bc, proj(bc)) > squaredist(bc, proj(bc * (-1.0)))) return proj(bc * (-1.0)); 60 56 61 57 return proj(bc); … … 165 161 { 166 162 if (N < 3) 167 163 return 0; /* polygons with less than three vertices have zero area */ 168 164 Coord t[3]; 169 165 t[0] = barycentre(x, N); … … 177 173 t[1] = x[i]; 178 174 t[2] = x[ii]; 179 double sc=scalarprod(crossprod(t[1] - t[0], t[2] - t[0]), t[0]) ;175 double sc=scalarprod(crossprod(t[1] - t[0], t[2] - t[0]), t[0]) ; 180 176 assert(sc >= -1e-10); // Error: tri a l'env (wrong orientation) 181 177 double area_gc = triarea(t[0], t[1], t[2]); 182 //if(area_gc<=0) printf("area_gc = %e\n", area_gc);183 178 double area_sc_gc_moon = 0; 184 179 if (d[i]) /* handle small circle case */ … … 188 183 char sgl = (mext > 0) ? -1 : 1; 189 184 area_sc_gc_moon = sgl * alun(arcdist(t[1], t[2]), fabs(scalarprod(t[1], pole))); 190 //if(area_sc_gc_moon<=0) printf("area_sc_gc_moon = %e\n", area_sc_gc_moon);191 185 gg_exact = gg_exact + sc_gc_moon_normalintegral(t[1], t[2], pole); 192 186 } 193 187 area += area_gc + area_sc_gc_moon; /* for "spherical circle segment" sum triangular part (at) and "small moon" and => account for small circle */ 194 188 g[i] = barycentre(t, 3) * (area_gc + area_sc_gc_moon); 195 //printf("g[%d] = (%e,%e,%e) * (%e+%e) = (%e,%e,%e) norm = %e\n", i, barycentre(t, 3).x, barycentre(t, 3).y, barycentre(t, 3).z, area_gc, area_sc_gc_moon, g[i].x, g[i].y, g[i].z, norm(g[i]));196 189 } 197 190 gg = barycentre(g, N); -
XIOS/dev/branch_openmp/extern/remap/src/timerRemap.cpp
r1146 r1328 4 4 #include <map> 5 5 #include <iostream> 6 using namespace ep_lib; 6 7 7 8 namespace sphereRemap { … … 9 10 using namespace std; 10 11 11 map<string,CTimer*> *CTimer::allTimer = 0; 12 //map<string,CTimer*> CTimer::allTimer; 13 map<string,CTimer*> *CTimer::allTimer_ptr = 0; 12 14 13 15 CTimer::CTimer(const string& name_) : name(name_) … … 55 57 CTimer& CTimer::get(const string name) 56 58 { 57 if(allTimer == 0) allTimer = new map<string,CTimer*>;58 59 map<string,CTimer*>::iterator it; 59 it=(*allTimer).find(name); 60 if (it==(*allTimer).end()) it=(*allTimer).insert(pair<string,CTimer*>(name,new CTimer(name))).first; 60 if(allTimer_ptr == 0) allTimer_ptr = new map<string,CTimer*>; 61 //it=allTimer.find(name); 62 it=allTimer_ptr->find(name); 63 //if (it==allTimer.end()) it=allTimer.insert(pair<string,CTimer*>(name,new CTimer(name))).first; 64 if (it==allTimer_ptr->end()) it=allTimer_ptr->insert(pair<string,CTimer*>(name,new CTimer(name))).first; 61 65 return *(it->second); 62 66 } -
XIOS/dev/branch_openmp/extern/remap/src/timerRemap.hpp
r1146 r1328 27 27 void print(void); 28 28 //static map<string,CTimer*> allTimer; 29 static map<string,CTimer*> *allTimer; 30 #pragma omp threadprivate(allTimer) 31 29 static map<string,CTimer*> *allTimer_ptr; 32 30 static double getTime(void); 33 31 static CTimer& get(string name); -
XIOS/dev/branch_openmp/extern/remap/src/tree.cpp
r1172 r1328 29 29 void CBasicTree::routeNodes(vector<int>& route, vector<Node>& nodes, int assignLevel) 30 30 { 31 32 33 34 35 31 for (int i = 0; i < nodes.size(); i++) 32 { 33 root->routeNode(&nodes[i], assignLevel); 34 route[i] = nodes[i].route; 35 } 36 36 } 37 37 38 38 void CBasicTree::routeIntersections(vector<vector<int> >& routes, vector<Node>& nodes) 39 39 { 40 41 40 for (int i = 0; i < nodes.size(); i++) 41 root->routeIntersection(routes[i], &nodes[i]); 42 42 } 43 43 44 44 void CBasicTree::build(vector<Node>& nodes) 45 45 { 46 47 46 newRoot(1); 47 insertNodes(nodes); 48 48 } 49 49 50 50 void CBasicTree::output(ostream& flux, int level) 51 51 { 52 52 root->output(flux,level,0) ; 53 53 } 54 54 void CBasicTree::slim(int nbIts) 55 55 { 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 56 for (int i = 0; i < nbIts; i++) 57 { 58 for (int level = root->level - 1; level > 0; level--) 59 { 60 slim2(root, level); 61 ri = 0; 62 emptyPool(); 63 } 64 65 for (int level = 2; level < root->level; level++) 66 { 67 slim2(root, level); 68 ri = 0; 69 emptyPool(); 70 } 71 } 72 72 } 73 73 … … 76 76 void CBasicTree::insertNode(NodePtr node) 77 77 { 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 78 node->tree = this; 79 increaseLevelSize(0); 80 push_back(node); 81 82 NodePtr q; 83 while (pool.size()) 84 { 85 q = pool.front(); 86 pool.pop_front(); 87 q = insert(q, root); 88 if (ri) 89 { 90 delete q; 91 ri = 0; 92 } 93 } 94 94 } 95 95 96 96 void CBasicTree::emptyPool(void) 97 97 { 98 99 100 101 102 103 104 105 106 107 108 98 while (pool.size()) 99 { 100 NodePtr q = pool.front(); 101 pool.pop_front(); 102 q = insert(q, root); 103 if (ri) 104 { 105 delete q; 106 ri = 0; 107 } 108 } 109 109 } 110 110 … … 142 142 root->parent = 0; 143 143 root->leafCount = 0; 144 // initialize root node on the sphere 145 root->centre.x=1 ; 146 root->centre.y=0 ; 147 root->centre.z=0 ; 144 // initialize root node on the sphere 145 root->centre.x=1 ; root->centre.y=0 ; root->centre.z=0 ; 148 146 root->radius = 0.; 149 147 root->reinserted = false; -
XIOS/dev/branch_openmp/extern/remap/src/tree.hpp
r1172 r1328 14 14 class CBasicTree 15 15 { 16 16 public: 17 17 18 19 20 21 22 18 NodePtr root; /* The main tree is stored as Nodes which can be reached through traversal starting here */ 19 NodePtr ref; // FIXME this reference, set by a node is odd, try to remove 20 int ri; /** this is set to one by a node in case of reinsertion */ 21 vector<int> levelSize; /** e.g. levelSize[0] == leafs.size() */ 22 vector<Node> leafs; /** leafs are stored in vector for easy access and rest of the tree nodes as separate allocations, only reachable through tree traversal */ 23 23 24 25 26 27 28 24 CBasicTree() : ri(0), levelSize(MAX_LEVEL_SIZE), root(NULL), isAssignedLevel(false), okSplit(true), isActiveOkSplit(false) {} 25 ~CBasicTree(); 26 void build(vector<Node>& nodes); 27 void slim(int nbIts = 1); 28 virtual void insertNodes(vector<Node>& node) = 0; 29 29 30 31 30 void routeNodes(vector<int>& route, vector<Node>& nodes, int assignLevel); 31 void routeIntersections(vector<vector<int> >& route, vector<Node>& nodes); 32 32 33 34 35 36 37 38 33 void push_back(NodePtr node); 34 void push_front(NodePtr node); 35 void increaseLevelSize(int level); 36 void decreaseLevelSize(int level); 37 void newRoot(int level); 38 void insertNode(NodePtr node); 39 39 void output(ostream& flux, int level) ; 40 40 41 41 int keepNodes; 42 42 bool isAssignedLevel ; 43 43 int assignLevel; … … 50 50 51 51 52 53 52 private: 53 deque<NodePtr > pool; 54 54 55 55 bool okSplit ; 56 56 57 57 protected: 58 58 void emptyPool(); 59 CBasicTree(int keepNodes_, int assignLevel_) : ri(0), levelSize(MAX_LEVEL_SIZE), root(NULL), keepNodes(keepNodes_), assignLevel(assignLevel_), isAssignedLevel(true), 60 okSplit(true), isActiveOkSplit(false) {} 59 CBasicTree(int keepNodes_, int assignLevel_) : ri(0), levelSize(MAX_LEVEL_SIZE), root(NULL), keepNodes(keepNodes_), assignLevel(assignLevel_), isAssignedLevel(true), okSplit(true), isActiveOkSplit(false) {} 61 60 }; 62 61 63 62 class CTree : public CBasicTree 64 63 { 65 66 64 public: 65 void insertNodes(vector<Node>& nodes); 67 66 }; 68 67 … … 70 69 { 71 70 72 73 71 public: 72 CSampleTree(int keepNodes_, int assignLevel_) : CBasicTree(keepNodes_,assignLevel_) {} 74 73 void slimAssignedLevel() ; 75 74 void removeExtraNode(void) ; 76 75 void insertNodes(vector<Node>& nodes); 77 76 }; 78 77 -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_fortran.cpp
r1287 r1328 31 31 { 32 32 fc_comm_map.insert(std::make_pair( std::make_pair( fint, omp_get_thread_num()) , comm)); 33 //printf("EP_Comm_c2f : MAP %p insert: %d, %d, %p\n", &fc_comm_map, fint, omp_get_thread_num(), comm.ep_comm_ptr);33 printf("EP_Comm_c2f : MAP %p insert: %d, %d, %p\n", &fc_comm_map, fint, omp_get_thread_num(), comm.ep_comm_ptr); 34 34 } 35 35 } … … 54 54 MPI_Comm comm_ptr; 55 55 comm_ptr = it->second; 56 //printf("EP_Comm_f2c : MAP %p find: %d, %d, %p\n", &fc_comm_map, it->first.first, it->first.second, comm_ptr.ep_comm_ptr);56 printf("EP_Comm_f2c : MAP %p find: %d, %d, %p\n", &fc_comm_map, it->first.first, it->first.second, comm_ptr.ep_comm_ptr); 57 57 return comm_ptr; 58 58 } -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_intercomm.cpp
r1287 r1328 47 47 48 48 MPI_Waitall(2, request, status); 49 50 //MPI_Send(&leader_ranks[0], 3, static_cast< ::MPI_Datatype>(MPI_INT), remote_leader, tag, peer_comm); 51 //MPI_Recv(&leader_ranks[3], 3, static_cast< ::MPI_Datatype>(MPI_INT), remote_leader, tag, peer_comm, &status[1]); 49 52 } 50 53 -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_lib.cpp
r1295 r1328 10 10 #pragma omp threadprivate(EP_PendingRequests) 11 11 12 13 12 14 namespace ep_lib 13 15 { 16 bool MPI_Comm::is_null() 17 { 18 if(!this->is_intercomm) 19 return this->mpi_comm == MPI_COMM_NULL.mpi_comm; 20 else 21 return this->ep_comm_ptr->intercomm->mpi_inter_comm == MPI_COMM_NULL.mpi_comm; 22 } 14 23 15 24 int tag_combine(int real_tag, int src, int dest) -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_send.cpp
r1295 r1328 19 19 if(!comm.is_ep) 20 20 return ::MPI_Send(buf, count, to_mpi_type(datatype), dest, tag, to_mpi_comm(comm.mpi_comm)); 21 22 MPI_Request request; 23 MPI_Status status; 24 MPI_Isend(buf, count, datatype, dest, tag, comm, &request); 25 MPI_Wait(&request, &status); 26 21 if(comm.is_intercomm) 22 { 23 MPI_Request request; 24 MPI_Status status; 25 MPI_Isend(buf, count, datatype, dest, tag, comm, &request); 26 MPI_Wait(&request, &status); 27 } 28 else 29 { 30 int ep_src_loc = comm.ep_comm_ptr->size_rank_info[1].first; 31 int ep_dest_loc = comm.ep_comm_ptr->comm_list->rank_map->at(dest).first; 32 int mpi_tag = tag_combine(tag, ep_src_loc, ep_dest_loc); 33 int mpi_dest = comm.ep_comm_ptr->comm_list->rank_map->at(dest).second; 34 35 ::MPI_Send(buf, count, to_mpi_type(datatype), mpi_dest, mpi_tag, to_mpi_comm(comm.mpi_comm)); 36 //printf("call mpi_send for intracomm, dest = %d, tag = %d\n", dest, tag); 37 } 27 38 //check_sum_send(buf, count, datatype, dest, tag, comm); 28 39 -
XIOS/dev/branch_openmp/extern/src_ep_dev/ep_type.hpp
r1295 r1328 344 344 } 345 345 346 bool is_null(); 347 346 348 }; 347 349
Note: See TracChangeset
for help on using the changeset viewer.