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 deallocate( sidmapg ) 319 320c 321c.... complete the averaging, adjust iBC, adjust BCtmp 322c 323 wnrm(:,1:3)=zero 324 do nn = 1, numnp ! loop over all nodes 325c leave nodes without wall-modeled surfaces unchanged 326 if(nsurf(nn).eq.0) cycle 327c 328c mark this as a node with comp3 329c 330 ikp=0 331 if(ibits(iBC(nn),3,3).eq.7 ) ikp=1 332c if(ibits(ibc(nn),3,3).eq.7 .and. BCtmp(nn,7).eq.zero) cycle 333c find out which velocity BC's were set by the user, and cancel them 334 ixset=ibits(iBC(nn),3,1) 335 iyset=ibits(iBC(nn),4,1) 336 izset=ibits(iBC(nn),5,1) 337 iBC(nn)=iBC(nn)-ixset*8-iyset*16-izset*32 338c 339c save this stripped iBC for later un-fixing 340c 341 iBCSAV=iBC(nn) 342c 343 if(abs(itwmod).eq.1) then ! slip velocity wall-model 344c If we're using the slip-velocity wall-model, then the velocity 345c boundary condition for this node will depend on how many modeled 346c surfaces share it... 347 if(nsurf(nn).eq.1) then ! 1 modeled wall 348c If this node is part of only one modeled wall, then the component 349c of the velocity normal to the surface is set to zero. In this case 350c only the first vector of BCtmp received normal contributions 351 iBC(nn)=iBC(nn)+8 ! assume normal is x-dominated first 352 wnrm(nn,:)=BCtmp(nn,4:6) 353 if(abs(wnrm(nn,3)).gt.abs(wnrm(nn,2)))then ! z beats y 354 if(abs(wnrm(nn,3)).gt.abs(wnrm(nn,1)))then ! z beats x 355 iBC(nn)=iBC(nn)+24 356 endif ! z beats x 357 endif ! z beats y 358 if(abs(wnrm(nn,2)).ge.abs(wnrm(nn,3)))then ! y beats z 359 if(abs(wnrm(nn,2)).gt.abs(wnrm(nn,1)))then ! y beats x 360 iBC(nn)=iBC(nn)+8 361 endif ! y beats x 362 endif ! y beats z !(adjusted iBC) 363 BCtmp(nn,7)=zero 364 else if(nsurf(nn).eq.2) then 365c If this node is shared by two modeled walls, then each wall 366c provides a normal vector along which the velocity must be zero. 367c This leaves only one "free" direction, along which the flow is nonzero. 368c The two normal vectors define a plane. Any pair of non-parallel 369c vectors in this plane can also be specified, leaving the same free 370c direction. Area-weighted average normal vectors for the two surfaces 371c sharing this node have been stored in BCtmp(4:6) and BCtmp(8:10) 372c We normalize the first of these, and then orthogonalize the second 373c against the first. After then normalizing the second, we choose which 374c cartesian direction dominates each of the vectors, and adjust iBC 375c and BCtmp accordingly 376 e1=BCtmp(nn,4:6) 377 wmag=e1(1)*e1(1)+e1(2)*e1(2)+e1(3)*e1(3) 378 wmag=1./sqrt(wmag) 379 e1=e1*wmag ! first vector is normalized 380c 381 e2=BCtmp(nn,8:10) 382 wmag=e2(1)*e1(1)+e2(2)*e1(2)+e2(3)*e1(3) ! wmag=e2.n1 383 e2=e2-wmag*e1 ! second vector is orthogonalized 384c 385 wmag=e2(1)*e2(1)+e2(2)*e2(2)+e2(3)*e2(3) 386 wmag=1./sqrt(wmag) 387 e2=e2*wmag ! second vector is normalized 388c 389 ! idom tells us which component is currently dominant 390 ! ivec(n) tells us which vector is dominated by comp n 391 ivec(:)=0 ! initialize 392 idom=1 ! assume x-comp dominates e1 393 bigcomp=abs(e1(1)) 394 ivec(idom)=1 395 do i = 2, nsd 396 if(abs(e1(i)).gt.bigcomp) then 397 ivec(idom)=0 398 bigcomp=abs(e1(i)) 399 idom=i 400 ivec(i)=1 401 endif 402 enddo 403 if(idom.ne.1) then 404 idom=1 ! assume x-comp dominates e2 405 else 406 idom=3 ! assume z-comp dominates e2 407 endif 408 bigcomp=abs(e2(idom)) 409 410 ivec(idom)=2 411 do i = 1, nsd 412 if(abs(e2(i)).gt.bigcomp) then 413 if(ivec(i).eq.1) cycle 414 ivec(idom)=0 415 bigcomp=abs(e2(i)) 416 idom=i 417 ivec(i)=2 418 endif 419 enddo 420 ! now we know which components dominate each vector 421 ixset=0 ! 422 iyset=0 ! initialize 423 izset=0 ! 424 if(ivec(1).ne.0) ixset=1 425 if(ivec(2).ne.0) iyset=1 426 if(ivec(3).ne.0) izset=1 427 ncomp=ixset+iyset+izset 428 if(ncomp.ne.2) write(*,*) 'WARNING: TROUBLE IN GENBC' 429 BCvecs(1,1:3)=e1(:) 430 BCvecs(2,1:3)=e2(:) 431 if((ixset.eq.1).and.(iyset.eq.1)) then 432 BCtmp(nn,4:6)=BCvecs(ivec(1),1:3) 433 BCtmp(nn,8:10)=BCvecs(ivec(2),1:3) 434 else if((ixset.eq.1).and.(izset.eq.1)) then 435 BCtmp(nn,4:6)=BCvecs(ivec(1),1:3) 436 BCtmp(nn,8:10)=BCvecs(ivec(3),1:3) 437 else if((iyset.eq.1).and.(izset.eq.1)) then 438 BCtmp(nn,4:6)=BCvecs(ivec(2),1:3) 439 BCtmp(nn,8:10)=BCvecs(ivec(3),1:3) 440 endif 441 iBC(nn)=iBC(nn)+ixset*8+iyset*16+izset*32 442 BCtmp(nn,7)=zero 443 BCtmp(nn,11)=zero 444 wnrm(nn,:)=BCtmp(nn,4:6)+BCtmp(nn,8:10) 445 else if(nsurf(nn).gt.2) then 446c If this node is shared by more than two modeled walls, then 447c it is a corner node and it will be treated as having no slip 448c The normals for all surfaces beyond the first two were 449c contributed to the first vector of BCtmp 450 wnrm(nn,:)=BCtmp(nn,4:6)+BCtmp(nn,8:10) 451 iBC(nn)=iBC(nn)+56 452 BCtmp(nn,7)=zero 453 endif ! curved surface 454 else if(abs(itwmod).eq.2) then ! effective viscosity wall-model 455c Otherwise, we're using the effective-viscosity wall-model and the 456c nodes on modeled surfaces are treated as no-slip 457 iBC(nn)=iBC(nn)+56 ! set 3-comp 458 if(itwmod.eq.-2) BCtmp(nn,7)=zero ! no slip 459 if(nsurf(nn).eq.1) 460 & wnrm(nn,:)=BCtmp(nn,4:6) 461 if(nsurf(nn).ge.2) 462 & wnrm(nn,:)=BCtmp(nn,4:6)+BCtmp(nn,8:10) 463 endif 464c Now normalize the wall normal for this node 465 if(itwmod.ne.0) then 466 wmag=sqrt(wnrm(nn,1)*wnrm(nn,1) 467 & +wnrm(nn,2)*wnrm(nn,2)+wnrm(nn,3)*wnrm(nn,3)) 468 wnrm(nn,:)=wnrm(nn,:)/wmag 469 endif 470c 471c.... put back the comp3 info for bctmp 472c 473 if(ikp.eq.1) then 474 iBC(nn)=iBCSAV+56 ! put it back to a comp3 475 BCtmp(nn,:)=BCtmpSAV(nn,:) ! ditto 476 endif 477 enddo ! loop over all nodes 478 endif ! end "if there are any surfID's" 479c 480c If you are using the turbulence wall with axisymmetry you 481c need to modify the axisymmetry angle to account for the discrete 482c normals at the wall being different than the exact normals 483c 484c find the my normal, my partners normal and correct the angle 485c$$$ 486c$$$ do i=1,numnp 487c$$$ wmag = wnrm(i,1) * wnrm(i,1) 488c$$$ & + wnrm(i,2) * wnrm(i,2) 489c$$$ & + wnrm(i,3) * wnrm(i,3) 490c$$$c 491c$$$c only "fix" wall nodes...other nodes still have a zero in wnrm 492c$$$c 493c$$$ if((btest(iBC(i),12)).and.(wmag.ne.zero)) then 494c$$$ BCtmp(i,1)=acos( wnrm(i,1) * wnrm(iper(i),1) 495c$$$ & + wnrm(i,2) * wnrm(iper(i),2) 496c$$$ & + wnrm(i,3) * wnrm(iper(i),3) ) 497c$$$ endif 498c$$$ enddo 499c 500c.... return 501c 502 return 503c 504 end 505 506 507 subroutine gensidcount(nsidg) 508c--------------------------------------------------------------------- 509c 510c This routine counts up the total number of surface ID's across 511c all processors and makes a list of them 512c 513c Inputs: 514c iBCB natural boundary condition switches and surfIDs 515c 516c Outputs: 517c nsidg number of surface ID's globally (including all procs) 518c sidmapg global list of surface ID's, lowest to highest 519c 520c--------------------------------------------------------------------- 521 use pointer_data ! access to miBCB 522 use turbSA ! access to sidmapg 523 include "common.h" 524 include "mpif.h" 525c 526 integer newflag, i, j 527 integer nsidl ! number of surface IDs on-proc 528 integer nsidt ! sum of all nsidl's 529 integer nsidg ! number of surface IDs globally 530 integer nsid(numpe) ! array of nsidl's 531 integer idisp(numpe) ! needed by mpi: see note below 532 type llnod ! structure for one node in a linked list 533 integer :: value 534 type (llnod), pointer :: next 535 end type llnod 536 type (llnod), pointer :: sidlist ! points to first elt of linked list 537 type (llnod), pointer :: sidelt ! points to generic elt of linked list 538 type (llnod), pointer :: nextelt ! points to generic elt of linked list 539 integer, allocatable :: sidmapl(:) ! list of surfID's on-proc 540 integer, allocatable :: temp(:) ! temp space 541c integer iBCB(numelb,ndiBCB) ! column 1: naturalBC switches 542 ! column 2: surface ID's 543c Since we don't know a priori how many surface ID's there are, 544c on-processor or globally, we will store the ID's as an open-ended 545c link list while we determine the total number of distinct ID's 546 allocate (sidlist) ! allocate the first element of the sid 547 ! linked list and point sidlist to it 548 nsidl=0 ! initialize # of sids on this processor 549 nsidg=0 550 nullify(sidlist % next) ! next does not exist yet 551 do iblk=1, nelblb ! loop over boundary element blocks 552 npro = lcblkb(1,iblk+1)-lcblkb(1,iblk) 553 do i = 1, npro ! loop over boundary elements (belts) 554 iBCB2=miBCB(iblk)%p(i,2) 555 if(iBCB2.ne.zero) then ! if a sid is given for this belt 556 if(nsidl.eq.0) then ! if this is the first sid we've seen 557 nsidl=1 ! increment our count and 558 sidlist % value = iBCB2 ! record its value 559 nullify(sidlist % next) ! next does not exist yet 560 else ! if this isn't the first sid 561 newflag = 1 ! assume this is a new sid 562 sidelt => sidlist ! starting with the first sid 563 do j = 1, nsidl ! check the assumption 564 if((sidelt % value).eq.iBCB2) then 565 newflag=0 ! ... 566 endif 567 if(j.ne.nsidl) then ! ... 568 sidelt => sidelt % next 569 endif! ... 570 enddo 571 if(newflag.eq.1) then! if it really is new to us 572 nsidl = nsidl + 1 ! increment our count 573 allocate (sidelt % next)! tack a new element to the end 574 sidelt => sidelt % next! point to the new element 575 sidelt % value = iBCB2 ! record the new sid 576 nullify(sidelt % next) ! next does not exist yet 577 endif ! (really is a new sid) 578 endif ! (first sid) 579 endif ! (belt has a sid) 580 enddo ! (loop over belts) 581 enddo ! (loop over boundary element blocks) 582c Copy the data from the linked list to a more easily-used array form 583 if(nsidl.gt.0) then 584 allocate( sidmapl(nsidl) ) 585 sidelt => sidlist ! starting with the first sid 586 do j = 1, nsidl 587 sidmapl(j)=sidelt%value 588 if(j.ne.nsidl) sidelt => sidelt%next 589 enddo 590 else 591 allocate( sidmapl(1)) ! some compilers/MPI don't like to send unallocated arrays 592 endif 593! Deallocate the link list 594! http://stackoverflow.com/questions/9184675/how-does-fortran-deallocate-linked-lists 595 sidelt => sidlist 596 nextelt => sidelt % next 597 do 598 deallocate(sidelt) 599 if( .not. associated(nextelt) ) exit 600 sidelt => nextelt 601 nextelt => sidelt % next 602 enddo 603 604c Gather the number of surface ID's on each processor 605 if (numpe.gt.1) then ! multiple processors 606c write the number of surfID's on the jth processor to slot j of nsid 607 call MPI_ALLGATHER(nsidl,1,MPI_INTEGER,nsid,1, 608 & MPI_INTEGER,MPI_COMM_WORLD,ierr) 609c count up the total number of surfID's among all processes 610 nsidt=0 611 do j=1,numpe 612 nsidt=nsidt+nsid(j) 613 enddo 614 else ! single processor 615c the local information is the global information for single-processor 616 nsid=nsidl 617 nsidt=nsidl 618 endif ! if-else for multiple processors 619 if(nsidt.gt.0) then 620c 621c Make all-processor surfID collage 622c 623c there will be some duplicate surface ID's when we gather, so we 624c will use a temporary array 625 allocate (temp(nsidt)) 626 if (numpe.gt.1) then ! multiple processors 627c we will gather surfID's from local on-proc sets to a global set 628c we will stack each processor's surfID list atop that of the previous 629c processor. If the list for processor i is called sidmapi, then our 630c global coordinate list sidmap will look like this: 631c --------------------------------------------------------------- 632c | sidmap1 | sidmap2 | sidmap3 | ... | 633c --------------------------------------------------------------- 634c <---nsid(1)---> <-----nsid(2)-----> <-nsid(3)-> 635c <------------------------nsidt-----------------------...----> 636c To accomplish this with MPI, we use MPI_ALLGATHERV, summarized as: 637cMPI_ALLGATHERV(sendbuf,sendcount,sendtype,recvbuf,recvcount,disp,recvtype,comm) 638c[ IN sendbuf] starting address of send buffer (choice) 639c[ IN sendcount] number of elements in send buffer (integer) 640c[ IN sendtype] data type of send buffer elements (handle) 641c[ OUT recvbuf] address of receive buffer (choice) 642c[ IN recvcount] number of elements received from each process (int array) 643c[ IN disp] displacement array 644c[ IN recvtype] data type of receive buffer elements (handle) 645c[ IN comm] communicator (handle) 646c The displacement array disp is an array of integers in which the jth 647c entry indicates which slot of sidmap marks beginning of sidmapj 648c So, first we will build this displacement array 649 idisp(:)=0 ! starting with zero, since MPI likes C-numbering 650 do j=2,numpe 651 idisp(j)=idisp(j-1)+nsid(j-1) ! see diagram above 652 enddo 653c Now, we gather the data 654 call MPI_ALLGATHERV(sidmapl(:),nsidl, 655 & MPI_INTEGER,temp(:),nsid,idisp, 656 & MPI_INTEGER,MPI_COMM_WORLD,ierr) 657c sort surfID's, lowest to highest 658 isorted = 0 659 do while (isorted.eq.0) ! while the list isn't sorted 660 isorted = 1 ! assume the list is sorted this time 661 do j = 2, nsidt ! loop over ID's 662 if(temp(j).lt.temp(j-1)) then ! ID exceeds predecessor 663 itmp=temp(j-1) 664 temp(j-1)=temp(j) 665 temp(j)=itmp 666 isorted=0 ! not sorted yet 667 endif 668 enddo !loop over ID's 669 enddo ! while not sorted 670c determine the total number of distinct surfID's, globally 671 nsidg=nsidt ! assume there are no duplicate SurfID's 672 do j = 2, nsidt 673 if(temp(j).eq.temp(j-1)) nsidg = nsidg - 1 ! correction 674 enddo 675c create duplicate-free surfID list 676 allocate( sidmapg(nsidg) ) 677 sidmapg(1)=temp(1) 678 nsidg = 1 679 do j = 2, nsidt 680 if(temp(j).ne.temp(j-1)) then 681 nsidg = nsidg + 1 682 sidmapg(nsidg)=temp(j) 683 endif 684 enddo 685 deallocate( temp ) 686 else ! single-processor 687c global data is local data in single processor case 688 nsidg=nsidl 689 allocate( sidmapg(nsidg) ) 690 sidmapg=sidmapl 691 deallocate(sidmapl) 692 endif ! if-else multiple processor 693c 694 endif ! if at least one surfid 695c 696 return 697c 698 end 699 700 701 702 subroutine genotwn(x, BCtmp, iBC, nsurf) 703c 704c---------------------------------------------------------------------- 705c This routine determines which node to use as the first node off the 706c wall in near-wall modeling traction calculations for each wall node. 707c Each wall node has a normal, calculated from the wall elements to 708c which that node belongs. We seek the node within the boundary 709c element that lies closest to the line defined by the normal vector. 710c We create normalized vectors pointing from the wall node in 711c question to each of the nodes in the boundary element. The vector 712c that has the largest projection onto the wall node's normal vector 713c points to the node we want. Nodes that are not on a wall point to 714c themselves as their own off-the-wall node. 715c 716c input: 717c x (nshg,3) : : nodal position vectors 718c wnrm (nshg,3) : (mod) : normal vectors for each node 719c iBCB5 (nshg) : (file) : wall flag for belt 720c ienb (numelb,nen) : (file) : connectivity for belts 721c 722c output: 723c otwn (nshg) : (mod) : off-the-wall nodes for each node 724c 725c 726c Kenneth Jansen, Summer 2000 (algorithm) 727c Michael Yaworski, Summer 2000 (code) 728c---------------------------------------------------------------------- 729c 730 use pointer_data ! used for mienb, miBCB 731 use turbSA 732 include "common.h" 733c 734 integer iel, nod, can 735 real*8 vec(3), leng, dp, bigdp, lil 736 real*8 x(numnp,nsd),BCtmp(nshg,ndof+7) 737 integer iBC(nshg), nsurf(nshg) 738 integer gbits 739 integer, allocatable :: ienb(:) 740 741 allocate( otwn(nshg) ) 742c 743c initialize otwn to oneself 744c 745 do nod = 1, nshg 746 otwn(nod)=nod 747 enddo 748c 749c determine otwn for each wall node 750c 751 do iblk=1, nelblb ! loop over boundary element blocks 752 npro = lcblkb(1,iblk+1)-lcblkb(1,iblk) 753 nenl = lcblkb(5,iblk) 754 nenbl = lcblkb(6,iblk) 755 nshl = lcblkb(9,iblk) 756 allocate( ienb(nshl) ) 757 do i = 1, npro ! loop over boundary elements 758 iBCB1=miBCB(iblk)%p(i,1) 759 iBCB2=miBCB(iblk)%p(i,2) 760 ienb(1:nshl)=mienb(iblk)%p(i,1:nshl) 761 if (btest(iBCB1,4)) then ! (on the wall) 762 do nod = 1, nenbl ! for each wall node in this belt 763 bigdp = zero ! initialize largest dot product 764 do can=nenbl+1,nenl ! loop over off-the-wall candidates 765 nn=ienb(can) 766 ! don't bother with wall nodes 767 if(nsurf(nn).ne.0) cycle 768 ! don't bother with no-slip nodes 769 if(ibits(iBC(nn),3,3).eq.7 .and. 770 & BCtmp(nn,7).eq.zero) cycle 771c ! calculate candidate vector 772 vec(:)=x(ienb(can),:)-x(ienb(nod),:) 773 leng=sqrt(vec(1)**2+vec(2)**2+vec(3)**2) 774 vec(:)=vec(:)/leng 775c ! calculate dot product with wnrm 776 vec(:)=vec(:)*wnrm(ienb(nod),:) 777 dp=vec(1)+vec(2)+vec(3) 778c ! set candidate as otwn if necessary 779c ! (wnrm points into fluid, so 780c ! we want the most positive dp) 781 if (dp.gt.bigdp) then 782 otwn(ienb(nod)) = ienb(can) 783 bigdp=dp 784 endif 785 enddo !(loop over off-the-wall candidates) 786 enddo !(loop over wall nodes in current belt) 787 endif 788 enddo !(loop over elts in block) 789 deallocate(ienb) 790 enddo !(loop over boundary element blocks) 791 do nn = 1, nshg 792 if((otwn(nn).eq.nn).and.(nsurf(nn).ne.0)) then ! if a node on a 793 ! modeled surface 794 ! didn't find a node 795 ! off the wall 796 lil = 1.0e32 ! distance to current closest prospect 797 do can = 1, nshg ! loop over nodes 798 if(nsurf(can).eq.0) then ! if this candidate is off the 799 ! wall 800 if(ibits(iBC(can),3,3).eq.7 .and. 801 & BCtmp(can,7).eq.zero) then ! no slip nodes not allowed 802 else ! not a forbidden node 803 vec(:)=x(nn,:)-x(can,:) 804 leng=vec(1)**2+vec(2)**2+vec(3)**2 805 if(leng.lt.lil) then ! this candidate is closest so far 806 lil=leng 807 otwn(nn)=can 808 endif ! closest so far 809 endif ! end of no slip nodes 810 endif ! off-the-wall check 811 enddo ! loop over nodes 812 endif ! lonely wall-node check 813 enddo ! loop over nodes 814c 815 return 816c 817 end 818 819