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