1 /* ---------------------------------------------------------------------- 2 This is the 3 4 ██╗ ██╗ ██████╗ ██████╗ ██████╗ ██╗ ██╗████████╗███████╗ 5 ██║ ██║██╔════╝ ██╔════╝ ██╔════╝ ██║ ██║╚══██╔══╝██╔════╝ 6 ██║ ██║██║ ███╗██║ ███╗██║ ███╗███████║ ██║ ███████╗ 7 ██║ ██║██║ ██║██║ ██║██║ ██║██╔══██║ ██║ ╚════██║ 8 ███████╗██║╚██████╔╝╚██████╔╝╚██████╔╝██║ ██║ ██║ ███████║ 9 ╚══════╝╚═╝ ╚═════╝ ╚═════╝ ╚═════╝ ╚═╝ ╚═╝ ╚═╝ ╚══════╝® 10 11 DEM simulation engine, released by 12 DCS Computing Gmbh, Linz, Austria 13 http://www.dcs-computing.com, office@dcs-computing.com 14 15 LIGGGHTS® is part of CFDEM®project: 16 http://www.liggghts.com | http://www.cfdem.com 17 18 Core developer and main author: 19 Christoph Kloss, christoph.kloss@dcs-computing.com 20 21 LIGGGHTS® is open-source, distributed under the terms of the GNU Public 22 License, version 2 or later. It is distributed in the hope that it will 23 be useful, but WITHOUT ANY WARRANTY; without even the implied warranty 24 of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. You should have 25 received a copy of the GNU General Public License along with LIGGGHTS®. 26 If not, see http://www.gnu.org/licenses . See also top-level README 27 and LICENSE files. 28 29 LIGGGHTS® and CFDEM® are registered trade marks of DCS Computing GmbH, 30 the producer of the LIGGGHTS® software and the CFDEM®coupling software 31 See http://www.cfdem.com/terms-trademark-policy for details. 32 33 ------------------------------------------------------------------------- 34 Contributing author and copyright for this file: 35 (if not contributing author is listed, this file has been contributed 36 by the core developer) 37 38 Copyright 2012- DCS Computing GmbH, Linz 39 Copyright 2009-2012 JKU Linz 40 ------------------------------------------------------------------------- */ 41 42 #ifndef LMP_MULTI_NODE_MESH_PARALLEL_I_H 43 #define LMP_MULTI_NODE_MESH_PARALLEL_I_H 44 45 #define BIG_MNMP 1.0e20 46 #define BUFFACTOR_MNMP 1.5 47 #define BUFMIN_MNMP 2000 48 #define BUFEXTRA_MNMP 2000 49 50 /* ---------------------------------------------------------------------- 51 consturctors 52 ------------------------------------------------------------------------- */ 53 54 template<int NUM_NODES> MultiNodeMeshParallel(LAMMPS * lmp)55 MultiNodeMeshParallel<NUM_NODES>::MultiNodeMeshParallel(LAMMPS *lmp) 56 : MultiNodeMesh<NUM_NODES>(lmp), 57 doParallellization_(true), 58 nLocal_(0), nGhost_(0), nGlobal_(0), nGlobalOrig_(0), 59 isParallel_(false), 60 isInsertionMesh_(false), 61 maxsend_(0), maxrecv_(0), 62 buf_send_(0), buf_recv_(0), 63 half_atom_cut_(0.), 64 size_exchange_(0), 65 size_forward_(0), 66 size_border_(0), 67 maxforward_(0),maxreverse_(0), 68 nswap_(0), 69 maxswap_(0), 70 sendnum_(0),recvnum_(0), 71 firstrecv_(0), 72 sendproc_(0),recvproc_(0), 73 size_forward_recv_(0), 74 size_reverse_recv_(0), 75 slablo_(0),slabhi_(0), 76 sendlist_(0), 77 sendwraplist_(0), 78 maxsendlist_(0), 79 pbc_flag_(0), 80 pbc_(0) 81 { 82 // initialize comm buffers & exchange memory 83 84 maxsend_ = BUFMIN_MNMP; 85 this->memory->create(buf_send_,maxsend_+BUFEXTRA_MNMP,"MultiNodeMeshParallel:buf_send"); 86 maxrecv_ = BUFMIN_MNMP; 87 this->memory->create(buf_recv_,maxrecv_,"MultiNodeMeshParallel:buf_recv"); 88 89 maxswap_ = 6; 90 allocate_swap(maxswap_); 91 92 sendlist_ = (int **) this->memory->smalloc(maxswap_*sizeof(int *),"MultiNodeMeshParallel:sendlist"); 93 sendwraplist_ = (int **) this->memory->smalloc(maxswap_*sizeof(int *),"MultiNodeMeshParallel:sendlist"); 94 this->memory->create(maxsendlist_,maxswap_,"MultiNodeMeshParallel:maxsendlist"); 95 for (int i = 0; i < maxswap_; i++) { 96 maxsendlist_[i] = BUFMIN_MNMP; 97 this->memory->create(sendlist_[i],BUFMIN_MNMP,"MultiNodeMeshParallel:sendlist[i]"); 98 this->memory->create(sendwraplist_[i],BUFMIN_MNMP,"MultiNodeMeshParallel:sendlist[i]"); 99 } 100 } 101 102 /* ---------------------------------------------------------------------- 103 destructor 104 ------------------------------------------------------------------------- */ 105 106 template<int NUM_NODES> ~MultiNodeMeshParallel()107 MultiNodeMeshParallel<NUM_NODES>::~MultiNodeMeshParallel() 108 { 109 free_swap(); 110 111 if (sendlist_) 112 for (int i = 0; i < maxswap_; i++) 113 this->memory->destroy(sendlist_[i]); 114 if (sendwraplist_) 115 for (int i = 0; i < maxswap_; i++) 116 this->memory->destroy(sendwraplist_[i]); 117 118 this->memory->sfree(sendlist_); 119 this->memory->sfree(sendwraplist_); 120 this->memory->destroy(maxsendlist_); 121 122 this->memory->destroy(buf_send_); 123 this->memory->destroy(buf_recv_); 124 } 125 126 /* ---------------------------------------------------------------------- 127 add and delete elements 128 ------------------------------------------------------------------------- */ 129 130 template<int NUM_NODES> addElement(double ** nodeToAdd)131 bool MultiNodeMeshParallel<NUM_NODES>::addElement(double **nodeToAdd) 132 { 133 134 if(MultiNodeMesh<NUM_NODES>::addElement(nodeToAdd)) 135 { 136 nLocal_++; 137 return true; 138 } 139 return false; 140 } 141 142 template<int NUM_NODES> deleteElement(int n)143 void MultiNodeMeshParallel<NUM_NODES>::deleteElement(int n) 144 { 145 146 if(n < nLocal_ && nGhost_ != 0) 147 this->error->one(FLERR,"Illegal call to MultiNodeMeshParallel<NUM_NODES>::deleteElement"); 148 149 MultiNodeMesh<NUM_NODES>::deleteElement(n); 150 151 if(n >= nLocal_) 152 nGhost_--; 153 else 154 nLocal_--; 155 } 156 157 /* ---------------------------------------------------------------------- 158 recalculate properties on setup (on start and during simulation) 159 ------------------------------------------------------------------------- */ 160 161 template<int NUM_NODES> refreshOwned(int setupFlag)162 void MultiNodeMeshParallel<NUM_NODES>::refreshOwned(int setupFlag) 163 { 164 MultiNodeMesh<NUM_NODES>::refreshOwned(setupFlag); 165 } 166 167 template<int NUM_NODES> refreshGhosts(int setupFlag)168 void MultiNodeMeshParallel<NUM_NODES>::refreshGhosts(int setupFlag) 169 { 170 MultiNodeMesh<NUM_NODES>::refreshGhosts(setupFlag); 171 } 172 173 /* ---------------------------------------------------------------------- 174 completely clear ghosts - called in borders() 175 ------------------------------------------------------------------------- */ 176 177 template<int NUM_NODES> clearGhosts()178 void MultiNodeMeshParallel<NUM_NODES>::clearGhosts() 179 { 180 // delete ghost data from container classes 181 182 while(nGhost_ > 0) 183 { 184 185 deleteElement(nLocal_); 186 } 187 } 188 189 /* ---------------------------------------------------------------------- 190 clear ghost data that is communicated via forward comm - called in forw comm 191 ------------------------------------------------------------------------- */ 192 193 template<int NUM_NODES> clearGhostForward(bool scale,bool translate,bool rotate)194 void MultiNodeMeshParallel<NUM_NODES>::clearGhostForward(bool scale,bool translate,bool rotate) 195 { 196 // delete ghost data from container classes 197 // delete only data that is communicated afterwards 198 199 for(int i = this->sizeLocal()+this->sizeGhost()-1; i >= this->sizeLocal(); i--) 200 { 201 // clear ghost data that belongs to this class 202 // must match push/pop implementation for forward comm in this class 203 if(translate || rotate || scale) 204 { 205 this->node_.del(i); 206 this->center_.del(i); 207 } 208 if(scale) 209 this->rBound_.del(i); 210 } 211 } 212 213 /* ---------------------------------------------------------------------- 214 check if all elements are in domain 215 ------------------------------------------------------------------------- */ 216 217 template<int NUM_NODES> allNodesInsideSimulationBox()218 bool MultiNodeMeshParallel<NUM_NODES>::allNodesInsideSimulationBox() 219 { 220 int flag = 0; 221 for(int i=0;i<sizeLocal();i++) 222 for(int j=0;j<NUM_NODES;j++) 223 { 224 225 if(!this->domain->is_in_domain(this->node_(i)[j])) 226 { 227 flag = 1; 228 break; 229 } 230 } 231 232 MPI_Max_Scalar(flag,this->world); 233 if(flag) return false; 234 else return true; 235 } 236 237 /* ---------------------------------------------------------------------- 238 set flag if used as insertion mesh 239 ------------------------------------------------------------------------- */ 240 241 template<int NUM_NODES> useAsInsertionMesh(bool parallelflag)242 void MultiNodeMeshParallel<NUM_NODES>::useAsInsertionMesh(bool parallelflag) 243 { 244 245 isInsertionMesh_ = true; 246 247 if(!parallelflag) 248 { 249 if(isParallel()) 250 this->error->all(FLERR,"If a run command is between the fix mesh/surface and the " 251 "fix insert command, you have to use fix mesh/surface/planar for " 252 "the insertion mesh"); 253 doParallellization_ = false; 254 } 255 } 256 257 /* ---------------------------------------------------------------------- 258 setup of communication 259 ------------------------------------------------------------------------- */ 260 261 template<int NUM_NODES> setup()262 void MultiNodeMeshParallel<NUM_NODES>::setup() 263 { 264 if(!doParallellization_) return; 265 266 double sublo[3],subhi[3], extent_acc; 267 double rBound_max, cut_ghost; 268 double **sublo_all, **subhi_all; 269 270 int nprocs = this->comm->nprocs; 271 int myloc[3], loc_dim, nextproc, need_this; 272 273 // get required size of communication per element 274 bool scale = this->isScaling(); 275 bool translate = this->isTranslating(); 276 bool rotate = this->isRotating(); 277 278 size_exchange_ = elemBufSize(OPERATION_COMM_EXCHANGE, NULL, scale,translate,rotate) + 1; 279 size_border_ = elemBufSize(OPERATION_COMM_BORDERS, NULL, scale,translate,rotate); 280 size_forward_ = elemBufSize(OPERATION_COMM_FORWARD, NULL, scale,translate,rotate); 281 size_reverse_ = elemBufSize(OPERATION_COMM_REVERSE, NULL, scale,translate,rotate); 282 283 // maxforward = # of datums in largest forward communication 284 // maxreverse = # of datums in largest reverse communication 285 286 maxforward_ = MathExtraLiggghts::max(size_exchange_,size_border_,size_forward_); 287 maxreverse_ = size_reverse_; 288 289 // copy comm and domain data 290 vectorCopy3D(this->comm->myloc,myloc); 291 vectorCopy3D(this->domain->sublo,sublo); 292 vectorCopy3D(this->domain->subhi,subhi); 293 294 this->memory->create(sublo_all,nprocs,3,"MultiNodeMeshParallel::setup() sublo_all"); 295 this->memory->create(subhi_all,nprocs,3,"MultiNodeMeshParallel::setup() subhi_all"); 296 297 // ghost elements are for computing interaction with owned particles 298 // so need to aquire ghost elements that overlap my subbox extened by 299 // half neigh cutoff 300 301 half_atom_cut_ = this->neighbor->cutneighmax / 2.; 302 303 if(this->isMoving()) 304 half_atom_cut_+= this->neighbor->skin / 2.; 305 306 // calculate maximum bounding radius of elements across all procs 307 rBound_max = 0.; 308 for(int i = 0; i < sizeLocal(); i++) 309 rBound_max = std::max(this->rBound_(i),rBound_max); 310 MPI_Max_Scalar(rBound_max,this->world); 311 312 // mesh element ghost cutoff is element bounding radius plus half atom neigh cut 313 cut_ghost = rBound_max + half_atom_cut_; 314 315 // set up maxneed_, sendneed_ 316 // account for non-uniform boundaries due to load-balancing 317 // so aquire sub-box bounds from all processors 318 319 MPI_Allgather(sublo,3,MPI_DOUBLE,&(sublo_all[0][0]),3,MPI_DOUBLE,this->world); 320 MPI_Allgather(subhi,3,MPI_DOUBLE,&(subhi_all[0][0]),3,MPI_DOUBLE,this->world); 321 322 // set up maxneed_ and sendneed_ 323 // assume element with max bound radius is in my subbox 324 325 for(int dim = 0; dim < 3; dim++) 326 { 327 bool is_x = dim == 0 ? true : false; 328 bool is_y = dim == 1 ? true : false; 329 bool is_z = dim == 2 ? true : false; 330 331 // go each direction (N-S-E-W-UP-DN) 332 333 maxneed_[dim] = 0; 334 for(int way = -1; way <= 1; way += 2) 335 { 336 // start from location of myself 337 // reset accumulated extent 338 loc_dim = myloc[dim]; 339 extent_acc = 0.; 340 need_this = 0; 341 sendneed_[dim][way == -1 ? 0 : 1] = 0; 342 343 while(extent_acc < cut_ghost) 344 { 345 // increase or decrease location 346 loc_dim += way; 347 348 // break if at dead end and non-pbc 349 if( (loc_dim < 0 && !this->domain->periodicity[dim]) || 350 (loc_dim > this->comm->procgrid[dim]-1 && !this->domain->periodicity[dim]) ) 351 break; 352 353 // wrap around PBCs 354 if(loc_dim < 0 && this->domain->periodicity[dim]) 355 loc_dim = this->comm->procgrid[dim]-1; 356 357 if(loc_dim > this->comm->procgrid[dim]-1) 358 loc_dim = 0; 359 360 // increase counters 361 need_this++; 362 sendneed_[dim][way == -1 ? 0 : 1]++; 363 364 // go to next proc in proc grid and add its extent 365 nextproc = this->comm->grid2proc[is_x ? loc_dim : myloc[0]] 366 [is_y ? loc_dim : myloc[1]] 367 [is_z ? loc_dim : myloc[2]]; 368 extent_acc += subhi_all[nextproc][dim] - sublo_all[nextproc][dim]; 369 } 370 371 maxneed_[dim] = std::max(maxneed_[dim],need_this); 372 } 373 374 // limit maxneed for non-pbc 375 376 if(maxneed_[dim] > this->comm->procgrid[dim]-1 && !this->domain->periodicity[dim]) 377 maxneed_[dim] = this->comm->procgrid[dim]-1; 378 } 379 380 // maxneed_ summed accross all processors 381 MPI_Max_Vector(maxneed_,3,this->world); 382 383 destroy(sublo_all); 384 destroy(subhi_all); 385 386 // allocate comm memory 387 388 nswap_ = 2 * (maxneed_[0]+maxneed_[1]+maxneed_[2]); 389 if (nswap_ > maxswap_) grow_swap(nswap_); 390 391 // setup parameters for each exchange: 392 // slablo_/slabhi_ = boundaries for slab of elements to send at each swap 393 // use -BIG/midpt/BIG to insure all elements included even if round-off occurs 394 // if round-off, atoms elements across PBC can be < or > than subbox boundary 395 // note that borders() only loops over subset of elements during each swap 396 397 // treat all as PBC here, non-PBC is handled in borders() via r/s need[][] 398 // pbc_flag_: 0 = nothing across a boundary, 1 = something across a boundary 399 // pbc_ = -1/0/1 for PBC factor in each of 3/6 orthogonal/triclinic dirs 400 // 1st part of if statement is sending to the west/south/down 401 // 2nd part of if statement is sending to the east/north/up 402 403 int dim,ineed; 404 405 int iswap = 0; 406 for (dim = 0; dim < 3; dim++) 407 { 408 for (ineed = 0; ineed < 2*maxneed_[dim]; ineed++) 409 { 410 pbc_flag_[iswap] = 0; 411 vectorZeroizeN(pbc_[iswap],6); 412 413 // send left, receive right 414 if (ineed % 2 == 0) 415 { 416 sendproc_[iswap] = this->comm->procneigh[dim][0]; 417 recvproc_[iswap] = this->comm->procneigh[dim][1]; 418 419 if (ineed < 2) slablo_[iswap] = -BIG_MNMP; 420 else slablo_[iswap] = 0.5 * (this->domain->sublo[dim] + this->domain->subhi[dim]); 421 422 // use half cut here, since rBound is used (added) in checkBorderElement() 423 slabhi_[iswap] = this->domain->sublo[dim] + half_atom_cut_; 424 425 if (myloc[dim] == 0) 426 { 427 pbc_flag_[iswap] = 1; 428 pbc_[iswap][dim] = 1; 429 } 430 } 431 // send right, receive left 432 else 433 { 434 sendproc_[iswap] = this->comm->procneigh[dim][1]; 435 recvproc_[iswap] = this->comm->procneigh[dim][0]; 436 437 // use half cut here, since rBound is used (added) in checkBorderElement() 438 slablo_[iswap] = this->domain->subhi[dim] - half_atom_cut_; 439 if (ineed < 2) slabhi_[iswap] = BIG_MNMP; 440 else slabhi_[iswap] = 0.5 * (this->domain->sublo[dim] + this->domain->subhi[dim]); 441 442 if (myloc[dim] == this->comm->procgrid[dim]-1) 443 { 444 pbc_flag_[iswap] = 1; 445 pbc_[iswap][dim] = -1; 446 } 447 } 448 iswap++; 449 } 450 } 451 } 452 453 /* ---------------------------------------------------------------------- 454 realloc the buffers needed for communication and swaps 455 ------------------------------------------------------------------------- */ 456 457 template<int NUM_NODES> grow_swap(int n)458 void MultiNodeMeshParallel<NUM_NODES>::grow_swap(int n) 459 { 460 free_swap(); 461 allocate_swap(n); 462 463 sendlist_ = (int **) 464 this->memory->srealloc(sendlist_,n*sizeof(int *),"MultiNodeMeshParallel:sendlist_"); 465 sendwraplist_ = (int **) 466 this->memory->srealloc(sendwraplist_,n*sizeof(int *),"MultiNodeMeshParallel:sendwraplist_"); 467 this->memory->grow(maxsendlist_,n,"MultiNodeMeshParallel:maxsendlist_"); 468 for (int i = maxswap_; i < n; i++) 469 { 470 maxsendlist_[i] = BUFMIN_MNMP; 471 this->memory->create(sendlist_[i],BUFMIN_MNMP,"MultiNodeMeshParallel:sendlist_[i]"); 472 this->memory->create(sendwraplist_[i],BUFMIN_MNMP,"MultiNodeMeshParallel:sendwraplist_[i]"); 473 } 474 maxswap_ = n; 475 } 476 477 template<int NUM_NODES> allocate_swap(int n)478 void MultiNodeMeshParallel<NUM_NODES>::allocate_swap(int n) 479 { 480 this->memory->create(sendnum_,n,"MultiNodeMeshParallel:sendnum_"); 481 this->memory->create(recvnum_,n,"MultiNodeMeshParallel:recvnum_"); 482 this->memory->create(sendproc_,n,"MultiNodeMeshParallel:sendproc_"); 483 this->memory->create(recvproc_,n,"MultiNodeMeshParallel:recvproc_"); 484 this->memory->create(size_forward_recv_,n,"MultiNodeMeshParallel:size"); 485 this->memory->create(size_reverse_recv_,n,"MultiNodeMeshParallel:size"); 486 this->memory->create(slablo_,n,"MultiNodeMeshParallel:slablo_"); 487 this->memory->create(slabhi_,n,"MultiNodeMeshParallel:slabhi_"); 488 this->memory->create(firstrecv_,n,"MultiNodeMeshParallel:firstrecv"); 489 this->memory->create(pbc_flag_,n,"MultiNodeMeshParallel:pbc_flag_"); 490 this->memory->create(pbc_,n,6,"MultiNodeMeshParallel:pbc_"); 491 } 492 493 template<int NUM_NODES> free_swap()494 void MultiNodeMeshParallel<NUM_NODES>::free_swap() 495 { 496 this->memory->destroy(sendnum_); 497 this->memory->destroy(recvnum_); 498 this->memory->destroy(sendproc_); 499 this->memory->destroy(recvproc_); 500 this->memory->destroy(size_forward_recv_); 501 this->memory->destroy(size_reverse_recv_); 502 this->memory->destroy(slablo_); 503 this->memory->destroy(slabhi_); 504 this->memory->destroy(firstrecv_); 505 this->memory->destroy(pbc_flag_); 506 this->memory->destroy(pbc_); 507 } 508 509 /* ---------------------------------------------------------------------- 510 realloc the size of the send buffer as needed with BUFFACTOR & BUFEXTRA 511 if flag = 1, realloc 512 if flag = 0, don't need to realloc with copy, just free/malloc 513 ------------------------------------------------------------------------- */ 514 515 template<int NUM_NODES> grow_send(int n,int flag)516 void MultiNodeMeshParallel<NUM_NODES>::grow_send(int n, int flag) 517 { 518 maxsend_ = static_cast<int> (BUFFACTOR_MNMP * n); 519 520 if (flag) 521 this->memory->grow(buf_send_,(maxsend_+BUFEXTRA_MNMP),"MultiNodeMeshParallel:buf_send"); 522 else { 523 this->memory->destroy(buf_send_); 524 this->memory->create(buf_send_,maxsend_+BUFEXTRA_MNMP,"MultiNodeMeshParallel:buf_send"); 525 } 526 } 527 528 /* ---------------------------------------------------------------------- 529 free/malloc the size of the recv buffer as needed with BUFFACTOR 530 ------------------------------------------------------------------------- */ 531 532 template<int NUM_NODES> grow_recv(int n)533 void MultiNodeMeshParallel<NUM_NODES>::grow_recv(int n) 534 { 535 maxrecv_ = static_cast<int> (BUFFACTOR_MNMP * n); 536 this->memory->destroy(buf_recv_); 537 this->memory->create(buf_recv_,maxrecv_,"MultiNodeMeshParallel:buf_recv"); 538 } 539 540 /* ---------------------------------------------------------------------- 541 realloc the size of the iswap sendlist as needed with BUFFACTOR 542 ------------------------------------------------------------------------- */ 543 544 template<int NUM_NODES> grow_list(int iswap,int n)545 void MultiNodeMeshParallel<NUM_NODES>::grow_list(int iswap, int n) 546 { 547 maxsendlist_[iswap] = static_cast<int> (BUFFACTOR_MNMP * n)+1; 548 this->memory->grow(sendlist_[iswap],maxsendlist_[iswap],"MultiNodeMeshParallel:sendlist[iswap]"); 549 this->memory->grow(sendwraplist_[iswap],maxsendlist_[iswap],"MultiNodeMeshParallel:sendlist[iswap]"); 550 } 551 552 /* ---------------------------------------------------------------------- 553 parallelization - 554 initially, all processes have read the whole data 555 ------------------------------------------------------------------------- */ 556 557 template<int NUM_NODES> initialSetup()558 void MultiNodeMeshParallel<NUM_NODES>::initialSetup() 559 { 560 nGlobalOrig_ = sizeLocal(); 561 562 // check for possible round-off isues 563 564 double span = this->node_.max_scalar()-this->node_.min_scalar(); 565 if(span < 1e-4) 566 this->error->all(FLERR,"Mesh error - root causes: (a) mesh empty or (b) dimensions too small - use different unit system"); 567 568 double comBefore[3]; 569 this->center_of_mass(comBefore); 570 571 // delete all elements that do not belong to this processor 572 573 deleteUnowned(); 574 575 if(sizeGlobal() != sizeGlobalOrig()) 576 { 577 578 char errstr[1024]; 579 580 if(0 == sizeGlobal()) 581 { 582 sprintf(errstr,"Mesh (id %s): All %d mesh elements have been lost / left the domain. \n" 583 "Please use 'boundary m m m' or scale/translate/rotate the mesh or change its dynamics\n" 584 "FYI: center of mass of mesh including scale/tranlate/rotate is %f / %f / %f\n" 585 " simulation box x from %f to %f y from %f to %f z from %f to %f\n" 586 " (gives indication about changes in scale/tranlate/rotate necessary to make simulation run)\n", 587 this->mesh_id_,sizeGlobalOrig()-sizeGlobal(),comBefore[0],comBefore[1],comBefore[2], 588 this->domain->boxlo[0],this->domain->boxhi[0],this->domain->boxlo[1],this->domain->boxhi[1],this->domain->boxlo[2],this->domain->boxhi[2]); 589 } 590 else 591 { 592 double comAfter[3]; 593 this->center_of_mass(comAfter); 594 595 sprintf(errstr,"Mesh (id %s): %d mesh elements have been lost / left the domain. \n" 596 "Please use 'boundary m m m' or scale/translate/rotate the mesh or change its dynamics\n" 597 "FYI: center of mass of mesh including scale/tranlate/rotate before cutting out elements is %f / %f / %f\n" 598 " simulation box x from %f to %f y from %f to %f z from %f to %f\n" 599 " center of mass of mesh after cutting out elements outside simulation box is is %f / %f / %f\n" 600 " (gives indication about changes in scale/tranlate/rotate necessary to make simulation run)\n", 601 this->mesh_id_,sizeGlobalOrig()-sizeGlobal(),comBefore[0],comBefore[1],comBefore[2], 602 this->domain->boxlo[0],this->domain->boxhi[0],this->domain->boxlo[1],this->domain->boxhi[1],this->domain->boxlo[2],this->domain->boxhi[2], 603 comAfter[0],comAfter[1],comAfter[2]); 604 } 605 this->error->all(FLERR,errstr); 606 } 607 608 // perform operations that should be done before initial setup 609 610 preInitialSetup(); 611 612 // set-up mesh parallelism 613 614 setup(); 615 616 // re-calculate properties for owned particles 617 618 refreshOwned(1); 619 620 // identify elements that are near borders 621 // forward communicate them 622 623 borders(); 624 625 // re-calculate properties for ghost particles 626 627 refreshGhosts(1); 628 629 // build mesh topology and neigh list 630 631 buildNeighbours(); 632 633 // perform quality check on the mesh 634 635 qualityCheck(); 636 637 if(doParallellization_) isParallel_ = true; 638 639 postInitialSetup(); 640 641 // stuff that should be done before resuming simulation 642 643 postBorders(); 644 645 } 646 647 /* ---------------------------------------------------------------------- 648 parallelization - aggregates pbc, exchange and borders 649 ------------------------------------------------------------------------- */ 650 651 template<int NUM_NODES> pbcExchangeBorders(int setupFlag)652 void MultiNodeMeshParallel<NUM_NODES>::pbcExchangeBorders(int setupFlag) 653 { 654 // need not do this during simulation for non-moving mesh and non-changing simulation box 655 656 if(setupFlag) this->reset_stepLastReset(); 657 658 // perform operations that should be done before setting up parallellism and exchanging elements 659 preSetup(); 660 661 if(!setupFlag && !this->isMoving() && !this->isDeforming() && !this->domain->box_change) return; 662 663 // set-up mesh parallelism 664 setup(); 665 666 // enforce pbc 667 pbc(); 668 669 // communicate particles 670 exchange(); 671 672 if(sizeGlobal() != sizeGlobalOrig()) 673 { 674 675 //this->error->all(FLERR,"Mesh elements have been lost"); 676 char errstr[500]; 677 sprintf(errstr,"Mesh (id %s): Mesh elements have been lost / left the domain. Please use " 678 "'boundary m m m' or scale/translate/rotate the mesh or change its dynamics", 679 this->mesh_id_); 680 this->error->all(FLERR,errstr); 681 } 682 683 // re-calculate properties for owned particles 684 685 refreshOwned(setupFlag); 686 687 // identify elements that are near borders 688 // forward communicate them 689 690 borders(); 691 692 // re-calculate properties for ghosts 693 refreshGhosts(setupFlag); 694 695 // stuff that should be done before resuming simulation 696 postBorders(); 697 698 } 699 700 /* ---------------------------------------------------------------------- 701 parallelization - clear data of reverse comm properties 702 ------------------------------------------------------------------------- */ 703 704 template<int NUM_NODES> clearReverse()705 void MultiNodeMeshParallel<NUM_NODES>::clearReverse() 706 { 707 // nothing to do here 708 } 709 710 /* ---------------------------------------------------------------------- 711 delete all particles which are not owned on this proc 712 ------------------------------------------------------------------------- */ 713 714 template<int NUM_NODES> deleteUnowned()715 void MultiNodeMeshParallel<NUM_NODES>::deleteUnowned() 716 { 717 718 int i = 0; 719 720 if(doParallellization_) 721 { 722 723 while(i < nLocal_) 724 { 725 if(!this->domain->is_in_subdomain(this->center_(i))) 726 this->deleteElement(i); 727 else i++; 728 } 729 730 // calculate nGlobal for the first time 731 MPI_Sum_Scalar(nLocal_,nGlobal_,this->world); 732 } 733 else 734 nGlobal_ = nLocal_; 735 736 } 737 738 /* ---------------------------------------------------------------------- 739 enforce periodic boundary conditions 740 ------------------------------------------------------------------------- */ 741 742 template<int NUM_NODES> pbc()743 void MultiNodeMeshParallel<NUM_NODES>::pbc() 744 { 745 if(!doParallellization_) return; 746 747 double centerNew[3], delta[3]; 748 749 for(int i = 0; i < this->sizeLocal(); i++) 750 { 751 vectorCopy3D(this->center_(i),centerNew); 752 this->domain->remap(centerNew); 753 vectorSubtract3D(centerNew,this->center_(i),delta); 754 755 // move element i incremental 756 if(vectorMag3DSquared(delta) > 1e-9) 757 this->moveElement(i,delta); 758 } 759 } 760 761 /* ---------------------------------------------------------------------- 762 exchange elements with nearby processors 763 ------------------------------------------------------------------------- */ 764 765 template<int NUM_NODES> exchange()766 void MultiNodeMeshParallel<NUM_NODES>::exchange() 767 { 768 if(!doParallellization_) return; 769 770 int nrecv, nsend = 0; 771 int nrecv1,nrecv2; 772 double *buf; 773 MPI_Request request; 774 MPI_Status status; 775 MPI_Comm world = this->world; 776 777 //int nprocs = this->comm->nprocs; 778 int *procgrid = this->comm->procgrid; 779 int procneigh[3][2]; 780 781 // clear global->local map for owned and ghost atoms 782 783 clearMap(); 784 785 // clear old ghosts 786 787 clearGhosts(); 788 789 // copy procneigh 790 for (int i = 0; i < 3; i++) 791 for( int j = 0; j < 2; j++) 792 procneigh[i][j] = this->comm->procneigh[i][j]; 793 794 for (int dim = 0; dim < 3; dim++) 795 { 796 // push data to buffer 797 798 nsend = pushExchange(dim); 799 800 // send/recv in both directions 801 // if 1 proc in dimension, no send/recv, set recv buf to send buf 802 // if 2 procs in dimension, single send/recv 803 // if more than 2 procs in dimension, send/recv to both neighbors 804 805 if (procgrid[dim] == 1) 806 { 807 nrecv = nsend; 808 buf = buf_send_; 809 } 810 else 811 { 812 MPI_Sendrecv(&nsend,1,MPI_INT,procneigh[dim][0],0,&nrecv1,1,MPI_INT,procneigh[dim][1],0,world,&status); 813 nrecv = nrecv1; 814 815 if (this->comm->procgrid[dim] > 2) 816 { 817 MPI_Sendrecv(&nsend,1,MPI_INT,procneigh[dim][1],0,&nrecv2,1,MPI_INT,procneigh[dim][0],0,world,&status); 818 nrecv += nrecv2; 819 } 820 821 if (nrecv > maxrecv_) grow_recv(nrecv); 822 823 MPI_Irecv(buf_recv_,nrecv1,MPI_DOUBLE,procneigh[dim][1],0,world,&request); 824 MPI_Send(buf_send_,nsend,MPI_DOUBLE,procneigh[dim][0],0,world); 825 MPI_Wait(&request,&status); 826 827 if (procgrid[dim] > 2) 828 { 829 MPI_Irecv(&buf_recv_[nrecv1],nrecv2,MPI_DOUBLE,procneigh[dim][0],0,world,&request); 830 MPI_Send(buf_send_,nsend,MPI_DOUBLE,procneigh[dim][1],0,world); 831 MPI_Wait(&request,&status); 832 } 833 834 buf = buf_recv_; 835 } 836 837 // check incoming elements to see if they are in my box 838 // if so, add on this proc 839 840 popExchange(nrecv,dim, buf); 841 842 } 843 844 // re-calculate nGlobal as some element might have been lost 845 MPI_Sum_Scalar(nLocal_,nGlobal_,world); 846 } 847 848 /* ---------------------------------------------------------------------- 849 generate ghost elements, refresh global map 850 ------------------------------------------------------------------------- */ 851 852 template<int NUM_NODES> borders()853 void MultiNodeMeshParallel<NUM_NODES>::borders() 854 { 855 if(doParallellization_) 856 { 857 int iswap, twoneed, nfirst, nlast, n, nsend, nrecv, smax, rmax; 858 bool sendflag, dummy = false; 859 double lo,hi; 860 MPI_Request request; 861 MPI_Status status; 862 863 nfirst = 0; 864 iswap = 0; 865 smax = rmax = 0; 866 867 for (int dim = 0; dim < 3; dim++) 868 { 869 nlast = 0; 870 871 // need to go left and right in each dim 872 twoneed = 2*maxneed_[dim]; 873 for (int ineed = 0; ineed < twoneed; ineed++) 874 { 875 lo = slablo_[iswap]; 876 hi = slabhi_[iswap]; 877 878 // find elements within slab boundaries lo/hi using <= and >= 879 880 if (ineed % 2 == 0) 881 { 882 nfirst = nlast; 883 nlast = sizeLocal() + sizeGhost(); 884 } 885 886 nsend = 0; 887 888 // sendflag = 0 if I do not send on this swap 889 890 sendflag = true; 891 int wrap = 0; 892 893 if(ineed % 2 == 0 && this->comm->myloc[dim] == 0) 894 { 895 if(this->domain->periodicity[dim] && !this->domain->triclinic && !dynamic_cast<DomainWedge*>(this->domain)) 896 wrap = 1; 897 else 898 sendflag = false; 899 } 900 901 if(ineed % 2 == 1 && this->comm->myloc[dim] == this->comm->procgrid[dim]-1) 902 { 903 if(this->domain->periodicity[dim] && !this->domain->triclinic && !dynamic_cast<DomainWedge*>(this->domain)) 904 wrap = -1; 905 else 906 sendflag = false; 907 } 908 909 // find send elements 910 if(sendflag) 911 { 912 913 for (int i = nfirst; i < nlast; i++) 914 { 915 int type = checkBorderElement(ineed, i, dim, lo, hi); 916 if(type != NOT_GHOST) 917 { 918 if (nsend >= maxsendlist_[iswap]) 919 grow_list(iswap,nsend); 920 sendlist_[iswap][nsend] = i; 921 if (wrap == 1) 922 { 923 switch (dim) 924 { 925 case 0: 926 type = IS_GHOST_WRAP_DIM_0_POS; 927 break; 928 case 1: 929 type = IS_GHOST_WRAP_DIM_1_POS; 930 break; 931 case 2: 932 type = IS_GHOST_WRAP_DIM_2_POS; 933 break; 934 } 935 } 936 else if (wrap == -1) 937 { 938 switch (dim) 939 { 940 case 0: 941 type = IS_GHOST_WRAP_DIM_0_NEG; 942 break; 943 case 1: 944 type = IS_GHOST_WRAP_DIM_1_NEG; 945 break; 946 case 2: 947 type = IS_GHOST_WRAP_DIM_2_NEG; 948 break; 949 } 950 } 951 sendwraplist_[iswap][nsend] = type; 952 nsend++; 953 954 } 955 } 956 } 957 958 // pack up list of border elements 959 960 if(nsend*size_border_ > maxsend_) 961 grow_send(nsend*size_border_,0); 962 963 n = pushElemListToBuffer(nsend, sendlist_[iswap], sendwraplist_[iswap], buf_send_, OPERATION_COMM_BORDERS, NULL, this->domain->boxlo, this->domain->boxhi,dummy,dummy,dummy); 964 965 // swap atoms with other proc 966 // no MPI calls except SendRecv if nsend/nrecv = 0 967 // put incoming ghosts at end of my atom arrays 968 // if swapping with self, simply copy, no messages 969 970 double *buf = NULL; 971 972 if (sendproc_[iswap] != this->comm->me) 973 { 974 MPI_Sendrecv(&nsend,1,MPI_INT,sendproc_[iswap],0,&nrecv,1,MPI_INT,recvproc_[iswap],0,this->world,&status); 975 if (nrecv*size_border_ > maxrecv_) 976 grow_recv(nrecv*size_border_); 977 if (nrecv) 978 MPI_Irecv(buf_recv_,nrecv*size_border_,MPI_DOUBLE,recvproc_[iswap],0,this->world,&request); 979 980 if (n) 981 MPI_Send(buf_send_,n,MPI_DOUBLE,sendproc_[iswap],0,this->world); 982 983 if (nrecv) 984 MPI_Wait(&request,&status); 985 986 buf = buf_recv_; 987 } 988 else 989 { 990 nrecv = nsend; 991 buf = buf_send_; 992 } 993 994 // unpack buffer 995 996 n = popElemListFromBuffer(nLocal_+nGhost_, nrecv, buf, OPERATION_COMM_BORDERS, NULL, dummy,dummy,dummy); 997 998 // set pointers & counters 999 1000 smax = MAX(smax,nsend); 1001 rmax = MAX(rmax,nrecv); 1002 sendnum_[iswap] = nsend; 1003 recvnum_[iswap] = nrecv; 1004 size_forward_recv_[iswap] = nrecv*size_forward_; 1005 size_reverse_recv_[iswap] = nsend*size_reverse_; 1006 firstrecv_[iswap] = nLocal_+nGhost_; 1007 nGhost_ += nrecv; 1008 iswap++; 1009 } 1010 } 1011 1012 // insure send/recv buffers are long enough for all forward & reverse comm 1013 int max = MAX(maxforward_*smax,maxreverse_*rmax); 1014 if (max > maxsend_) grow_send(max,0); 1015 max = MAX(maxforward_*rmax,maxreverse_*smax); 1016 if (max > maxrecv_) grow_recv(max); 1017 } 1018 1019 // build global-local map 1020 this->generateMap(); 1021 } 1022 1023 /* ---------------------------------------------------------------------- 1024 check if element qualifies as ghost 1025 ------------------------------------------------------------------------- */ 1026 1027 template<int NUM_NODES> checkBorderElement(const int ineed,const int i,const int dim,const double lo,const double hi)1028 inline int MultiNodeMeshParallel<NUM_NODES>::checkBorderElement(const int ineed, const int i, const int dim, const double lo, const double hi) const 1029 { 1030 if (ineed % 2 == 0) 1031 return checkBorderElementLeft(i,dim,lo,hi); 1032 else 1033 return checkBorderElementRight(i,dim,lo,hi); 1034 } 1035 1036 template<int NUM_NODES> checkBorderElementLeft(const int i,const int dim,const double lo,const double hi)1037 inline int MultiNodeMeshParallel<NUM_NODES>::checkBorderElementLeft(const int i, const int dim, const double lo, const double hi) const 1038 { 1039 1040 // center of triangle 1041 const double pos = this->center_(i)[dim]; 1042 // hi is extended by the bounding radius 1043 const double hi_extended = hi + this->rBound_(i); 1044 // check whether center is inside interval 1045 if (pos >= lo && pos <= hi_extended) 1046 return IS_GHOST; 1047 1048 return NOT_GHOST; 1049 } 1050 1051 template<int NUM_NODES> checkBorderElementRight(const int i,const int dim,const double lo,const double hi)1052 inline int MultiNodeMeshParallel<NUM_NODES>::checkBorderElementRight(const int i, const int dim, const double lo, const double hi) const 1053 { 1054 // center of triangle 1055 const double pos = this->center_(i)[dim]; 1056 // lo is extended by the bounding radius 1057 const double lo_extended = lo - this->rBound_(i); 1058 // check whether center is inside interval 1059 if (pos >= lo_extended && pos <= hi) 1060 return IS_GHOST; 1061 1062 return NOT_GHOST; 1063 } 1064 1065 /* ---------------------------------------------------------------------- 1066 communicate properties to ghost elements 1067 ------------------------------------------------------------------------- */ 1068 1069 template<int NUM_NODES> forwardComm(std::string property)1070 void MultiNodeMeshParallel<NUM_NODES>::forwardComm(std::string property) 1071 { 1072 std::list<std::string> properties (1, property); 1073 forwardComm(&properties); 1074 } 1075 1076 template<int NUM_NODES> forwardComm(std::list<std::string> * properties)1077 void MultiNodeMeshParallel<NUM_NODES>::forwardComm(std::list<std::string> * properties) 1078 { 1079 int n; 1080 MPI_Request request; 1081 MPI_Status status; 1082 int me = this->comm->me; 1083 1084 bool scale = this->isScaling(); 1085 bool translate = this->isTranslating(); 1086 bool rotate = this->isRotating(); 1087 1088 // exit here if no forward communication at all 1089 1090 if(size_forward_ == 0) 1091 return; 1092 1093 const int size_this = properties ? elemBufSize(OPERATION_COMM_REVERSE, properties, scale, translate, rotate) : 1; 1094 1095 // exchange data with another proc 1096 // if other proc is self, just copy 1097 1098 for (int iswap = 0; iswap < nswap_; iswap++) 1099 { 1100 if (sendproc_[iswap] != me) 1101 { 1102 if (size_forward_recv_[iswap] && size_this) 1103 { 1104 int nrecv = size_forward_recv_[iswap]; 1105 if (properties) 1106 { 1107 // size forward is the size of all forward buffers 1108 nrecv /= size_forward_; 1109 // size_this is the size of the buffers listed in properties 1110 nrecv *= size_this; 1111 } 1112 MPI_Irecv(buf_recv_, nrecv, MPI_DOUBLE,recvproc_[iswap],0,this->world,&request); 1113 } 1114 1115 n = pushElemListToBuffer(sendnum_[iswap],sendlist_[iswap], sendwraplist_[iswap],buf_send_,OPERATION_COMM_FORWARD, properties, this->domain->boxlo, this->domain->boxhi,scale,translate,rotate); 1116 1117 if (n) 1118 MPI_Send(buf_send_,n,MPI_DOUBLE,sendproc_[iswap],0,this->world); 1119 1120 if (size_forward_recv_[iswap] && size_this) 1121 MPI_Wait(&request,&status); 1122 1123 n = popElemListFromBuffer(firstrecv_[iswap],recvnum_[iswap],buf_recv_,OPERATION_COMM_FORWARD, properties, scale,translate,rotate); 1124 1125 } 1126 else 1127 { 1128 n = pushElemListToBuffer(sendnum_[iswap], sendlist_[iswap], sendwraplist_[iswap], buf_send_, OPERATION_COMM_FORWARD, properties, this->domain->boxlo, this->domain->boxhi, scale, translate, rotate); 1129 1130 // note buf_recv_ not used in this case (just use buf_send_ as receive buffer 1131 n = popElemListFromBuffer(firstrecv_[iswap], recvnum_[iswap], buf_send_, OPERATION_COMM_FORWARD, properties, scale, translate, rotate); 1132 } 1133 } 1134 } 1135 1136 /* ---------------------------------------------------------------------- 1137 reverse communication of properties on atoms every timestep 1138 ------------------------------------------------------------------------- */ 1139 1140 template<int NUM_NODES> reverseComm(std::string property)1141 void MultiNodeMeshParallel<NUM_NODES>::reverseComm(std::string property) 1142 { 1143 std::list<std::string> properties (1, property); 1144 reverseComm(&properties); 1145 } 1146 1147 template<int NUM_NODES> reverseComm(std::list<std::string> * properties)1148 void MultiNodeMeshParallel<NUM_NODES>::reverseComm(std::list<std::string> * properties) 1149 { 1150 int n; 1151 MPI_Request request; 1152 MPI_Status status; 1153 int me = this->comm->me; 1154 1155 bool scale = this->isScaling(); 1156 bool translate = this->isTranslating(); 1157 bool rotate = this->isRotating(); 1158 1159 const int size_this = properties ? elemBufSize(OPERATION_COMM_REVERSE, properties, scale, translate, rotate) : 1; 1160 1161 // exchange data with another proc 1162 // if other proc is self, just copy 1163 1164 for (int iswap = nswap_-1; iswap >= 0; iswap--) 1165 { 1166 if (sendproc_[iswap] != me) 1167 { 1168 if (size_reverse_recv_[iswap] && size_this) 1169 { 1170 int nrecv = size_reverse_recv_[iswap]; 1171 if (properties) 1172 { 1173 // size reverse is the size of all reverse buffers 1174 nrecv /= size_reverse_; 1175 // size_this is the size of the buffers listed in properties 1176 nrecv *= size_this; 1177 } 1178 MPI_Irecv(buf_recv_, nrecv, MPI_DOUBLE, sendproc_[iswap], 0, this->world, &request); 1179 } 1180 1181 n = pushElemListToBufferReverse(firstrecv_[iswap], recvnum_[iswap], buf_send_, OPERATION_COMM_REVERSE, properties, scale, translate, rotate); 1182 1183 if (n) 1184 MPI_Send(buf_send_,n,MPI_DOUBLE,recvproc_[iswap],0,this->world); 1185 1186 if (size_reverse_recv_[iswap] && size_this) 1187 MPI_Wait(&request,&status); 1188 1189 n = popElemListFromBufferReverse(sendnum_[iswap], sendlist_[iswap], buf_recv_, OPERATION_COMM_REVERSE, properties, scale, translate, rotate); 1190 } 1191 else 1192 { 1193 n = pushElemListToBufferReverse(firstrecv_[iswap], recvnum_[iswap], buf_send_, OPERATION_COMM_REVERSE, properties, scale, translate, rotate); 1194 // note buf_recv_ not used in this case (just use buf_send_ as receive buffer 1195 n = popElemListFromBufferReverse(sendnum_[iswap], sendlist_[iswap], buf_send_, OPERATION_COMM_REVERSE, properties, scale, translate, rotate); 1196 } 1197 } 1198 } 1199 1200 #endif 1201