1 subroutine genBC (iBC, BC, x, ilwork, iper) 2c 3c---------------------------------------------------------------------- 4c This routine generates the essential prescribed boundary conditions. 5c 6c input: 7c iBC (nshg) : boundary condition code 8c nBC (nshg) : boundary condition mapping array 9c 10c output: 11c BC (nshg,ndofBC) : The constraint data for prescribed BC 12c 13c 14c Note: genBC1 reduces the input data for the velocity. In the 15c case of varying velocity direction in the generation, the 16c results may not be correct. (since a linearity assumption is 17c made in the generation). 18c 19c 20c Farzin Shakib, Spring 1986. 21c Zdenek Johan, Winter 1991. (Fortran 90) 22c---------------------------------------------------------------------- 23c 24 use slpw 25 use readarrays ! used to access BCinp, nBC 26 use specialBC ! filling acs here 27 include "common.h" 28c 29 dimension iBC(nshg), nsurf(nshg), 30 & BC(nshg,ndofBC), 31 & x(numnp,nsd), ilwork(nlwork), 32 & iper(nshg) 33c 34c BCinp for each point has: 35c D T P c11 c12 c13 M1 c21 c22 c23 M2 theta S1 S2 S3... 36c 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15... 37c Remember, ndof=nsd+2+nsclr 38c 39c Arrays in the following 1 line are now dimensioned in readnblk 40c dimension BCinp(numpbc,ndof+7) 41c 42 dimension BCtmp(nshg,ndof+7) 43c 44c ndof+7= 3(thermos) + (nsd-1)*(nsd+1) + nscalars + 1 (theta) 45c #vect *(vec dir +mag) 46c 47c.... ---------------------------> Input <--------------------------- 48c 49c.... convert boundary condition data 50c 51 BCtmp = zero 52c 53 if(numpbc.ne.0) then 54 do i = 1, ndof+7 55 where (nBC(:) .ne. 0) BCtmp(:,i) = BCinp(nBC(:),i) 56 enddo 57 deallocate(BCinp) 58 endif 59 60c 61 if(any(BCtmp(:,12).ne.0)) then 62 iabc=1 63 allocate (acs(nshg,2)) 64 where (btest(iBC,10)) 65 acs(:,1) = cos(BCtmp(:,12)) 66 acs(:,2) = sin(BCtmp(:,12)) 67 endwhere 68 endif 69 70c 71c.... ----------------------> Wall Normals <-------------------------- 72c (calculate the normal and adjust BCinp to the true normal as needed) 73c 74 if(navier.eq.1)then 75 call genwnm (iBC, BCtmp, x, ilwork, iper, nsurf) 76 endif 77c determine the first point off the wall for each wall node 78 if(navier.eq.1)then 79 call genotwn (x,BCtmp, iBC, nsurf) 80 endif 81c.... ------------------------> Conversion <------------------------- 82c 83c.... convert the input boundary conditions to condensed version 84c 85 BC = zero 86c 87 if(navier.eq.1)then ! zero navier means Euler simulation 88 call genBC1 (BCtmp, iBC, BC) 89 else ! navier equals zero 90 allocate(BCtmpg(nshg,ndof+7)) 91 allocate(BCg(nshg,ndofBC)) 92 allocate(iBCg(nshg)) 93 BCtmpg=BCtmp 94 iBCg=iBC 95 call genBC1 (BCtmp, iBC, BC) 96c... genBC1 convert BCtmp to BC 97 BCg=BC 98 icdg=icd 99c... find slip wall 100 call findslpw(x,ilwork,iper,iBC) 101c... apply slip wall condition to wall nodes 102 do i=1,nslpwnd 103 nn=idxslpw(i) 104 call slpwBC(mpslpw(nn,1),mpslpw(nn,2),iBCg(nn), 105 & BCg(nn,:), BCtmpg(nn,:), 106 & iBC(nn), BC(nn,:), 107 & wlnorm(nn,:,:) ) 108 enddo 109 icd=icdg 110 endif 111 112 113c 114c.... ---------------------------> Echo <---------------------------- 115c 116c.... echo the input data 117c 118 if (necho .lt. 3) then 119 nn = 0 120 do n = 1, nshg 121 if (nBC(n) .ne. 0) then 122 nn = nn + 1 123 if(mod(nn,50).eq.1) 124 & write(iecho,1000)ititle,(j,j=1,ndofBC) 125 write (iecho,1100) n, (BC(n,i),i=1,ndofBC) 126 endif 127 enddo 128 endif 129c 130c.... return 131c 132 return 133c 134 1000 format(a80,//, 135 &' P r e s c r i b e d B o u n d a r y C o n d i t i o n s',//, 136 & ' Node ',/, 137 & ' Number ',5x,6('BC',i1,:,10x)) 138 1100 format(1p,2x,i5,3x,6(e12.5,1x)) 139c 140 end 141 142 143 144 145 146 147 148 149 150 151 subroutine genwnm (iBC, BCtmp, x, ilwork, iper, nsurf) 152c---------------------------------------------------------------------- 153c This routine generates the normal to a wall 154c 155c input: 156c iBC (nshg) : boundary condition code 157c nBC (nshg) : boundary condition mapping array 158c 159c output: 160c BCtmp (nshg,ndof+6) : The constraint data for prescribed BC 161c 162c 163c---------------------------------------------------------------------- 164c 165 use turbSA 166 use pointer_data ! used for mienb, mibcb 167 include "common.h" 168 include "mpif.h" 169c 170 character*20 fname1, fmt1 171 character*5 cname 172 dimension iBC(nshg), iper(nshg), 173 & x(numnp,nsd), ilwork(nlwork) 174c 175 dimension BCtmpSAV(nshg,ndof+7) 176 dimension BCtmp(nshg,ndof+7), fBC(nshg,ndofBC), 177 & e1(3), e2(3), 178 & elnrm(3), asum(numnp) 179c 180 integer sid, nsidg 181 integer nsurf(nshg), ivec(nsd) 182 logical :: firstvisit(nshg) 183 real*8 BCvecs(2,nsd) 184 integer, allocatable :: ienb(:) 185 dimension wnmdb(nshg,nsd) 186c 187c wnrm is dimensioned nshg but the math is only done for straight sided 188c elements at this point so wnrm will not be calculated for the hierarchic 189c modes. Note that the wall model creates a p.w. linear representation 190c only at this time. 191c 192 allocate ( wnrm(nshg,3) ) 193c 194c.... ----------------------> Wall Normals <-------------------------- 195c (calculate the normal and adjust BCinp to the true normal as needed) 196c 197c 198 asum = zero 199 wnrm = zero 200 201c 202c.... Save a copy of BCtmp so that after we calculate the normals we 203c can recover the comp3 information. 204c 205 BCtmpSAV=BCtmp 206c 207c Count out the number of surface ID's on-processor, and map them 208 call gensidcount(nsidg) 209 210c 211 if(nsidg.gt.0) then ! if there are any surfID's 212 nsurf(:) = 0 213 do k = 1, nsidg ! loop over Surface ID's 214 sid = sidmapg(k) 215 firstvisit(:)=.true. 216 wnrm(:,1:3)=zero 217 do iblk=1, nelblb ! loop over boundary element blocks 218 npro = lcblkb(1,iblk+1)-lcblkb(1,iblk) 219 nenbl = lcblkb(6,iblk) 220 nshl = lcblkb(9,iblk) 221 allocate( ienb(nshl) ) 222 do i = 1, npro ! loop over boundary elements 223 iBCB1=miBCB(iblk)%p(i,1) 224 iBCB2=miBCB(iblk)%p(i,2) 225 ienb(1:nshl)=mienb(iblk)%p(i,1:nshl) 226c don't bother with elements that aren't on modeled surfaces 227c if ( not (wall set and traction set) ) 228 if (.not.(btest(iBCB1,2).and.btest(iBCB1,4))) 229 & cycle 230c don't bother with elements that don't lie on the current surface 231 if (iBCB2.ne.sid) cycle 232c 233c.... calculate this element's area-weighted normal vector 234c 235 e1 = x(ienb(2),:)-x(ienb(1),:) 236 e2 = x(ienb(3),:)-x(ienb(1),:) 237 elnrm(1) = e1(2)*e2(3)-e1(3)*e2(2) 238 elnrm(2) = e1(3)*e2(1)-e1(1)*e2(3) 239 elnrm(3) = e1(1)*e2(2)-e1(2)*e2(1) 240c Tetrahedral elements have negative volumes in phastaI, so 241c the normals calculated from the boundary faces must be inverted 242c to point into the fluid 243 if(nenbl.eq.3) elnrm(:)=-elnrm(:) 244c 245c.... add area-weighted normals to the nodal tallies for this surface 246c 247 do j = 1, nenbl ! loop over elt boundary nodes 248 nn=ienb(j) ! global node number 249 if(firstvisit(nn)) then 250 firstvisit(nn)=.false. 251 nsurf(nn)=nsurf(nn)+1 252 if(nsurf(nn).eq.1) BCtmp(nn,4:6)=zero 253 if(nsurf(nn).eq.2) BCtmp(nn,8:10)=zero 254 endif 255 wnrm(nn,:)=wnrm(nn,:)+elnrm 256 enddo ! loop over elt boundary nodes 257 enddo ! end loop over boundary elements in block 258 deallocate(ienb) 259 enddo ! end loop over boundary element blocks 260c Now we have all of this surface's contributions to wall normals 261c for all nodes, along with an indication of how many surfaces 262c each node has encountered so far. Contributions from other processors 263c should now be accumulated for this surface 264c 265c axisymm BC's need BC (and we have not built it yet) so we need to create 266c the entries it needs. 267c 268 if(iabc==1) then 269 where (btest(iBC,10)) 270 fBC(:,1) = cos(BCtmp(:,12)) 271 fBC(:,2) = sin(BCtmp(:,12)) 272 endwhere 273 274c before the commu we need to rotate the residual vector for axisymmetric 275c boundary conditions (so that off processor periodicity is a dof add instead 276c of a dof combination). Take care of all nodes now so periodicity, like 277c commu is a simple dof add. 278 call rotabc(wnrm, iBC, 'in ') 279 endif 280 281 if (numpe.gt.1) then 282c.... add areas and norms contributed from other processors 283 call commu (wnrm, ilwork, 3, 'in ') 284 endif 285c.... account for periodicity and axisymmetry 286 call bc3per(iBC,wnrm, iper,ilwork, 3) 287c at this point the master has all the information (slaves are zeroed and 288c waiting to be told what the master has...lets do that). 289c 290c.... local periodic (and axisymmetric) boundary conditions (no communications) 291 wnmdb=wnrm 292 do i = 1,numnp ! only use the vertices in the normal calc 293 wnrm(i,:) = wnrm(iper(i),:) 294 enddo 295 wnmdb=wnrm 296 if (numpe.gt.1) then 297c.... tell other processors the resulting (and now complete) sums 298 call commu (wnrm, ilwork, 3, 'out') 299 endif 300c Axisymmetric slaves need to be rotated back 301 if(iabc==1) then !are there any axisym bc's 302 call rotabc(wnrm, iBC, 'out') 303 endif 304c Now all nodes have all the normal contributions for this surface, 305c including those from off-processor and periodic nodes. 306c We distribute these summed vectors to the proper slots in BCtmp, 307c according to how many surfaces each node has seen so far 308 do nn = 1, nshg 309 if(nsurf(nn).eq.1) 310 & BCtmp(nn,4:6)=BCtmp(nn,4:6)+wnrm(nn,:) 311 if(nsurf(nn).eq.2) 312 & BCtmp(nn,8:10)=BCtmp(nn,8:10)+wnrm(nn,:) 313 if(nsurf(nn).gt.2) 314 & BCtmp(nn,4:6)=BCtmp(nn,4:6)+wnrm(nn,:) 315 enddo 316c That's all for this surface; move on to the next 317 enddo ! end loop over surface ID's 318 319c 320c.... complete the averaging, adjust iBC, adjust BCtmp 321c 322 wnrm(:,1:3)=zero 323 do nn = 1, numnp ! loop over all nodes 324c leave nodes without wall-modeled surfaces unchanged 325 if(nsurf(nn).eq.0) cycle 326c 327c mark this as a node with comp3 328c 329 ikp=0 330 if(ibits(iBC(nn),3,3).eq.7 ) ikp=1 331c if(ibits(ibc(nn),3,3).eq.7 .and. BCtmp(nn,7).eq.zero) cycle 332c find out which velocity BC's were set by the user, and cancel them 333 ixset=ibits(iBC(nn),3,1) 334 iyset=ibits(iBC(nn),4,1) 335 izset=ibits(iBC(nn),5,1) 336 iBC(nn)=iBC(nn)-ixset*8-iyset*16-izset*32 337c 338c save this stripped iBC for later un-fixing 339c 340 iBCSAV=iBC(nn) 341c 342 if(abs(itwmod).eq.1) then ! slip velocity wall-model 343c If we're using the slip-velocity wall-model, then the velocity 344c boundary condition for this node will depend on how many modeled 345c surfaces share it... 346 if(nsurf(nn).eq.1) then ! 1 modeled wall 347c If this node is part of only one modeled wall, then the component 348c of the velocity normal to the surface is set to zero. In this case 349c only the first vector of BCtmp received normal contributions 350 iBC(nn)=iBC(nn)+8 ! assume normal is x-dominated first 351 wnrm(nn,:)=BCtmp(nn,4:6) 352 if(abs(wnrm(nn,3)).gt.abs(wnrm(nn,2)))then ! z beats y 353 if(abs(wnrm(nn,3)).gt.abs(wnrm(nn,1)))then ! z beats x 354 iBC(nn)=iBC(nn)+24 355 endif ! z beats x 356 endif ! z beats y 357 if(abs(wnrm(nn,2)).ge.abs(wnrm(nn,3)))then ! y beats z 358 if(abs(wnrm(nn,2)).gt.abs(wnrm(nn,1)))then ! y beats x 359 iBC(nn)=iBC(nn)+8 360 endif ! y beats x 361 endif ! y beats z !(adjusted iBC) 362 BCtmp(nn,7)=zero 363 else if(nsurf(nn).eq.2) then 364c If this node is shared by two modeled walls, then each wall 365c provides a normal vector along which the velocity must be zero. 366c This leaves only one "free" direction, along which the flow is nonzero. 367c The two normal vectors define a plane. Any pair of non-parallel 368c vectors in this plane can also be specified, leaving the same free 369c direction. Area-weighted average normal vectors for the two surfaces 370c sharing this node have been stored in BCtmp(4:6) and BCtmp(8:10) 371c We normalize the first of these, and then orthogonalize the second 372c against the first. After then normalizing the second, we choose which 373c cartesian direction dominates each of the vectors, and adjust iBC 374c and BCtmp accordingly 375 e1=BCtmp(nn,4:6) 376 wmag=e1(1)*e1(1)+e1(2)*e1(2)+e1(3)*e1(3) 377 wmag=1./sqrt(wmag) 378 e1=e1*wmag ! first vector is normalized 379c 380 e2=BCtmp(nn,8:10) 381 wmag=e2(1)*e1(1)+e2(2)*e1(2)+e2(3)*e1(3) ! wmag=e2.n1 382 e2=e2-wmag*e1 ! second vector is orthogonalized 383c 384 wmag=e2(1)*e2(1)+e2(2)*e2(2)+e2(3)*e2(3) 385 wmag=1./sqrt(wmag) 386 e2=e2*wmag ! second vector is normalized 387c 388 ! idom tells us which component is currently dominant 389 ! ivec(n) tells us which vector is dominated by comp n 390 ivec(:)=0 ! initialize 391 idom=1 ! assume x-comp dominates e1 392 bigcomp=abs(e1(1)) 393 ivec(idom)=1 394 do i = 2, nsd 395 if(abs(e1(i)).gt.bigcomp) then 396 ivec(idom)=0 397 bigcomp=abs(e1(i)) 398 idom=i 399 ivec(i)=1 400 endif 401 enddo 402 if(idom.ne.1) then 403 idom=1 ! assume x-comp dominates e2 404 else 405 idom=3 ! assume z-comp dominates e2 406 endif 407 bigcomp=abs(e2(idom)) 408 409 ivec(idom)=2 410 do i = 1, nsd 411 if(abs(e2(i)).gt.bigcomp) then 412 if(ivec(i).eq.1) cycle 413 ivec(idom)=0 414 bigcomp=abs(e2(i)) 415 idom=i 416 ivec(i)=2 417 endif 418 enddo 419 ! now we know which components dominate each vector 420 ixset=0 ! 421 iyset=0 ! initialize 422 izset=0 ! 423 if(ivec(1).ne.0) ixset=1 424 if(ivec(2).ne.0) iyset=1 425 if(ivec(3).ne.0) izset=1 426 ncomp=ixset+iyset+izset 427 if(ncomp.ne.2) write(*,*) 'WARNING: TROUBLE IN GENBC' 428 BCvecs(1,1:3)=e1(:) 429 BCvecs(2,1:3)=e2(:) 430 if((ixset.eq.1).and.(iyset.eq.1)) then 431 BCtmp(nn,4:6)=BCvecs(ivec(1),1:3) 432 BCtmp(nn,8:10)=BCvecs(ivec(2),1:3) 433 else if((ixset.eq.1).and.(izset.eq.1)) then 434 BCtmp(nn,4:6)=BCvecs(ivec(1),1:3) 435 BCtmp(nn,8:10)=BCvecs(ivec(3),1:3) 436 else if((iyset.eq.1).and.(izset.eq.1)) then 437 BCtmp(nn,4:6)=BCvecs(ivec(2),1:3) 438 BCtmp(nn,8:10)=BCvecs(ivec(3),1:3) 439 endif 440 iBC(nn)=iBC(nn)+ixset*8+iyset*16+izset*32 441 BCtmp(nn,7)=zero 442 BCtmp(nn,11)=zero 443 wnrm(nn,:)=BCtmp(nn,4:6)+BCtmp(nn,8:10) 444 else if(nsurf(nn).gt.2) then 445c If this node is shared by more than two modeled walls, then 446c it is a corner node and it will be treated as having no slip 447c The normals for all surfaces beyond the first two were 448c contributed to the first vector of BCtmp 449 wnrm(nn,:)=BCtmp(nn,4:6)+BCtmp(nn,8:10) 450 iBC(nn)=iBC(nn)+56 451 BCtmp(nn,7)=zero 452 endif ! curved surface 453 else if(abs(itwmod).eq.2) then ! effective viscosity wall-model 454c Otherwise, we're using the effective-viscosity wall-model and the 455c nodes on modeled surfaces are treated as no-slip 456 iBC(nn)=iBC(nn)+56 ! set 3-comp 457 if(itwmod.eq.-2) BCtmp(nn,7)=zero ! no slip 458 if(nsurf(nn).eq.1) 459 & wnrm(nn,:)=BCtmp(nn,4:6) 460 if(nsurf(nn).ge.2) 461 & wnrm(nn,:)=BCtmp(nn,4:6)+BCtmp(nn,8:10) 462 endif 463c Now normalize the wall normal for this node 464 if(itwmod.ne.0) then 465 wmag=sqrt(wnrm(nn,1)*wnrm(nn,1) 466 & +wnrm(nn,2)*wnrm(nn,2)+wnrm(nn,3)*wnrm(nn,3)) 467 wnrm(nn,:)=wnrm(nn,:)/wmag 468 endif 469c 470c.... put back the comp3 info for bctmp 471c 472 if(ikp.eq.1) then 473 iBC(nn)=iBCSAV+56 ! put it back to a comp3 474 BCtmp(nn,:)=BCtmpSAV(nn,:) ! ditto 475 endif 476 enddo ! loop over all nodes 477 endif ! end "if there are any surfID's" 478c 479c If you are using the turbulence wall with axisymmetry you 480c need to modify the axisymmetry angle to account for the discrete 481c normals at the wall being different than the exact normals 482c 483c find the my normal, my partners normal and correct the angle 484c$$$ 485c$$$ do i=1,numnp 486c$$$ wmag = wnrm(i,1) * wnrm(i,1) 487c$$$ & + wnrm(i,2) * wnrm(i,2) 488c$$$ & + wnrm(i,3) * wnrm(i,3) 489c$$$c 490c$$$c only "fix" wall nodes...other nodes still have a zero in wnrm 491c$$$c 492c$$$ if((btest(iBC(i),12)).and.(wmag.ne.zero)) then 493c$$$ BCtmp(i,1)=acos( wnrm(i,1) * wnrm(iper(i),1) 494c$$$ & + wnrm(i,2) * wnrm(iper(i),2) 495c$$$ & + wnrm(i,3) * wnrm(iper(i),3) ) 496c$$$ endif 497c$$$ enddo 498c 499c.... return 500c 501 return 502c 503 end 504 505 506 subroutine gensidcount(nsidg) 507c--------------------------------------------------------------------- 508c 509c This routine counts up the total number of surface ID's across 510c all processors and makes a list of them 511c 512c Inputs: 513c iBCB natural boundary condition switches and surfIDs 514c 515c Outputs: 516c nsidg number of surface ID's globally (including all procs) 517c sidmapg global list of surface ID's, lowest to highest 518c 519c--------------------------------------------------------------------- 520 use pointer_data ! access to miBCB 521 use turbSA ! access to sidmapg 522 include "common.h" 523 include "mpif.h" 524c 525 integer newflag, i, j 526 integer nsidl ! number of surface IDs on-proc 527 integer nsidt ! sum of all nsidl's 528 integer nsidg ! number of surface IDs globally 529 integer nsid(numpe) ! array of nsidl's 530 integer idisp(numpe) ! needed by mpi: see note below 531 type llnod ! structure for one node in a linked list 532 integer :: value 533 type (llnod), pointer :: next 534 end type llnod 535 type (llnod), pointer :: sidlist ! points to first elt of linked list 536 type (llnod), pointer :: sidelt ! points to generic elt of linked list 537 type (llnod), pointer :: nextelt ! points to generic elt of linked list 538 integer, allocatable :: sidmapl(:) ! list of surfID's on-proc 539 integer, allocatable :: temp(:) ! temp space 540c integer iBCB(numelb,ndiBCB) ! column 1: naturalBC switches 541 ! column 2: surface ID's 542c Since we don't know a priori how many surface ID's there are, 543c on-processor or globally, we will store the ID's as an open-ended 544c link list while we determine the total number of distinct ID's 545 allocate (sidlist) ! allocate the first element of the sid 546 ! linked list and point sidlist to it 547 nsidl=0 ! initialize # of sids on this processor 548 nsidg=0 549 nullify(sidlist % next) ! next does not exist yet 550 do iblk=1, nelblb ! loop over boundary element blocks 551 npro = lcblkb(1,iblk+1)-lcblkb(1,iblk) 552 do i = 1, npro ! loop over boundary elements (belts) 553 iBCB2=miBCB(iblk)%p(i,2) 554 if(iBCB2.ne.zero) then ! if a sid is given for this belt 555 if(nsidl.eq.0) then ! if this is the first sid we've seen 556 nsidl=1 ! increment our count and 557 sidlist % value = iBCB2 ! record its value 558 nullify(sidlist % next) ! next does not exist yet 559 else ! if this isn't the first sid 560 newflag = 1 ! assume this is a new sid 561 sidelt => sidlist ! starting with the first sid 562 do j = 1, nsidl ! check the assumption 563 if((sidelt % value).eq.iBCB2) then 564 newflag=0 ! ... 565 endif 566 if(j.ne.nsidl) then ! ... 567 sidelt => sidelt % next 568 endif! ... 569 enddo 570 if(newflag.eq.1) then! if it really is new to us 571 nsidl = nsidl + 1 ! increment our count 572 allocate (sidelt % next)! tack a new element to the end 573 sidelt => sidelt % next! point to the new element 574 sidelt % value = iBCB2 ! record the new sid 575 nullify(sidelt % next) ! next does not exist yet 576 endif ! (really is a new sid) 577 endif ! (first sid) 578 endif ! (belt has a sid) 579 enddo ! (loop over belts) 580 enddo ! (loop over boundary element blocks) 581c Copy the data from the linked list to a more easily-used array form 582 if(nsidl.gt.0) then 583 allocate( sidmapl(nsidl) ) 584 sidelt => sidlist ! starting with the first sid 585 do j = 1, nsidl 586 sidmapl(j)=sidelt%value 587 if(j.ne.nsidl) sidelt => sidelt%next 588 enddo 589 else 590 allocate( sidmapl(1)) ! some compilers/MPI don't like to send unallocated arrays 591 endif 592! Deallocate the link list 593! http://stackoverflow.com/questions/9184675/how-does-fortran-deallocate-linked-lists 594 sidelt => sidlist 595 nextelt => sidelt % next 596 do 597 deallocate(sidelt) 598 if( .not. associated(nextelt) ) exit 599 sidelt => nextelt 600 nextelt => sidelt % next 601 enddo 602 603c Gather the number of surface ID's on each processor 604 if (numpe.gt.1) then ! multiple processors 605c write the number of surfID's on the jth processor to slot j of nsid 606 call MPI_ALLGATHER(nsidl,1,MPI_INTEGER,nsid,1, 607 & MPI_INTEGER,MPI_COMM_WORLD,ierr) 608c count up the total number of surfID's among all processes 609 nsidt=0 610 do j=1,numpe 611 nsidt=nsidt+nsid(j) 612 enddo 613 else ! single processor 614c the local information is the global information for single-processor 615 nsid=nsidl 616 nsidt=nsidl 617 endif ! if-else for multiple processors 618 if(nsidt.gt.0) then 619c 620c Make all-processor surfID collage 621c 622c there will be some duplicate surface ID's when we gather, so we 623c will use a temporary array 624 allocate (temp(nsidt)) 625 if (numpe.gt.1) then ! multiple processors 626c we will gather surfID's from local on-proc sets to a global set 627c we will stack each processor's surfID list atop that of the previous 628c processor. If the list for processor i is called sidmapi, then our 629c global coordinate list sidmap will look like this: 630c --------------------------------------------------------------- 631c | sidmap1 | sidmap2 | sidmap3 | ... | 632c --------------------------------------------------------------- 633c <---nsid(1)---> <-----nsid(2)-----> <-nsid(3)-> 634c <------------------------nsidt-----------------------...----> 635c To accomplish this with MPI, we use MPI_ALLGATHERV, summarized as: 636cMPI_ALLGATHERV(sendbuf,sendcount,sendtype,recvbuf,recvcount,disp,recvtype,comm) 637c[ IN sendbuf] starting address of send buffer (choice) 638c[ IN sendcount] number of elements in send buffer (integer) 639c[ IN sendtype] data type of send buffer elements (handle) 640c[ OUT recvbuf] address of receive buffer (choice) 641c[ IN recvcount] number of elements received from each process (int array) 642c[ IN disp] displacement array 643c[ IN recvtype] data type of receive buffer elements (handle) 644c[ IN comm] communicator (handle) 645c The displacement array disp is an array of integers in which the jth 646c entry indicates which slot of sidmap marks beginning of sidmapj 647c So, first we will build this displacement array 648 idisp(:)=0 ! starting with zero, since MPI likes C-numbering 649 do j=2,numpe 650 idisp(j)=idisp(j-1)+nsid(j-1) ! see diagram above 651 enddo 652c Now, we gather the data 653 call MPI_ALLGATHERV(sidmapl(:),nsidl, 654 & MPI_INTEGER,temp(:),nsid,idisp, 655 & MPI_INTEGER,MPI_COMM_WORLD,ierr) 656c sort surfID's, lowest to highest 657 isorted = 0 658 do while (isorted.eq.0) ! while the list isn't sorted 659 isorted = 1 ! assume the list is sorted this time 660 do j = 2, nsidt ! loop over ID's 661 if(temp(j).lt.temp(j-1)) then ! ID exceeds predecessor 662 itmp=temp(j-1) 663 temp(j-1)=temp(j) 664 temp(j)=itmp 665 isorted=0 ! not sorted yet 666 endif 667 enddo !loop over ID's 668 enddo ! while not sorted 669c determine the total number of distinct surfID's, globally 670 nsidg=nsidt ! assume there are no duplicate SurfID's 671 do j = 2, nsidt 672 if(temp(j).eq.temp(j-1)) nsidg = nsidg - 1 ! correction 673 enddo 674c create duplicate-free surfID list 675 allocate( sidmapg(nsidg) ) 676 sidmapg(1)=temp(1) 677 nsidg = 1 678 do j = 2, nsidt 679 if(temp(j).ne.temp(j-1)) then 680 nsidg = nsidg + 1 681 sidmapg(nsidg)=temp(j) 682 endif 683 enddo 684 deallocate( temp ) 685 else ! single-processor 686c global data is local data in single processor case 687 nsidg=nsidl 688 allocate( sidmapg(nsidg) ) 689 sidmapg=sidmapl 690 deallocate(sidmapl) 691 endif ! if-else multiple processor 692c 693 endif ! if at least one surfid 694c 695 return 696c 697 end 698 699 700 701 subroutine genotwn(x, BCtmp, iBC, nsurf) 702c 703c---------------------------------------------------------------------- 704c This routine determines which node to use as the first node off the 705c wall in near-wall modeling traction calculations for each wall node. 706c Each wall node has a normal, calculated from the wall elements to 707c which that node belongs. We seek the node within the boundary 708c element that lies closest to the line defined by the normal vector. 709c We create normalized vectors pointing from the wall node in 710c question to each of the nodes in the boundary element. The vector 711c that has the largest projection onto the wall node's normal vector 712c points to the node we want. Nodes that are not on a wall point to 713c themselves as their own off-the-wall node. 714c 715c input: 716c x (nshg,3) : : nodal position vectors 717c wnrm (nshg,3) : (mod) : normal vectors for each node 718c iBCB5 (nshg) : (file) : wall flag for belt 719c ienb (numelb,nen) : (file) : connectivity for belts 720c 721c output: 722c otwn (nshg) : (mod) : off-the-wall nodes for each node 723c 724c 725c Kenneth Jansen, Summer 2000 (algorithm) 726c Michael Yaworski, Summer 2000 (code) 727c---------------------------------------------------------------------- 728c 729 use pointer_data ! used for mienb, miBCB 730 use turbSA 731 include "common.h" 732c 733 integer iel, nod, can 734 real*8 vec(3), leng, dp, bigdp, lil 735 real*8 x(numnp,nsd),BCtmp(nshg,ndof+7) 736 integer iBC(nshg), nsurf(nshg) 737 integer gbits 738 integer, allocatable :: ienb(:) 739 740 allocate( otwn(nshg) ) 741c 742c initialize otwn to oneself 743c 744 do nod = 1, nshg 745 otwn(nod)=nod 746 enddo 747c 748c determine otwn for each wall node 749c 750 do iblk=1, nelblb ! loop over boundary element blocks 751 npro = lcblkb(1,iblk+1)-lcblkb(1,iblk) 752 nenl = lcblkb(5,iblk) 753 nenbl = lcblkb(6,iblk) 754 nshl = lcblkb(9,iblk) 755 allocate( ienb(nshl) ) 756 do i = 1, npro ! loop over boundary elements 757 iBCB1=miBCB(iblk)%p(i,1) 758 iBCB2=miBCB(iblk)%p(i,2) 759 ienb(1:nshl)=mienb(iblk)%p(i,1:nshl) 760 if (btest(iBCB1,4)) then ! (on the wall) 761 do nod = 1, nenbl ! for each wall node in this belt 762 bigdp = zero ! initialize largest dot product 763 do can=nenbl+1,nenl ! loop over off-the-wall candidates 764 nn=ienb(can) 765 ! don't bother with wall nodes 766 if(nsurf(nn).ne.0) cycle 767 ! don't bother with no-slip nodes 768 if(ibits(iBC(nn),3,3).eq.7 .and. 769 & BCtmp(nn,7).eq.zero) cycle 770c ! calculate candidate vector 771 vec(:)=x(ienb(can),:)-x(ienb(nod),:) 772 leng=sqrt(vec(1)**2+vec(2)**2+vec(3)**2) 773 vec(:)=vec(:)/leng 774c ! calculate dot product with wnrm 775 vec(:)=vec(:)*wnrm(ienb(nod),:) 776 dp=vec(1)+vec(2)+vec(3) 777c ! set candidate as otwn if necessary 778c ! (wnrm points into fluid, so 779c ! we want the most positive dp) 780 if (dp.gt.bigdp) then 781 otwn(ienb(nod)) = ienb(can) 782 bigdp=dp 783 endif 784 enddo !(loop over off-the-wall candidates) 785 enddo !(loop over wall nodes in current belt) 786 endif 787 enddo !(loop over elts in block) 788 deallocate(ienb) 789 enddo !(loop over boundary element blocks) 790 do nn = 1, nshg 791 if((otwn(nn).eq.nn).and.(nsurf(nn).ne.0)) then ! if a node on a 792 ! modeled surface 793 ! didn't find a node 794 ! off the wall 795 lil = 1.0e32 ! distance to current closest prospect 796 do can = 1, nshg ! loop over nodes 797 if(nsurf(can).eq.0) then ! if this candidate is off the 798 ! wall 799 if(ibits(iBC(can),3,3).eq.7 .and. 800 & BCtmp(can,7).eq.zero) then ! no slip nodes not allowed 801 else ! not a forbidden node 802 vec(:)=x(nn,:)-x(can,:) 803 leng=vec(1)**2+vec(2)**2+vec(3)**2 804 if(leng.lt.lil) then ! this candidate is closest so far 805 lil=leng 806 otwn(nn)=can 807 endif ! closest so far 808 endif ! end of no slip nodes 809 endif ! off-the-wall check 810 enddo ! loop over nodes 811 endif ! lonely wall-node check 812 enddo ! loop over nodes 813c 814 return 815c 816 end 817 818