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 integer, allocatable :: sidmapl(:) ! list of surfID's on-proc 538 integer, allocatable :: temp(:) ! temp space 539c integer iBCB(numelb,ndiBCB) ! column 1: naturalBC switches 540 ! column 2: surface ID's 541c Since we don't know a priori how many surface ID's there are, 542c on-processor or globally, we will store the ID's as an open-ended 543c link list while we determine the total number of distinct ID's 544 allocate (sidlist) ! allocate the first element of the sid 545 ! linked list and point sidlist to it 546 nsidl=0 ! initialize # of sids on this processor 547 nsidg=0 548 do iblk=1, nelblb ! loop over boundary element blocks 549 npro = lcblkb(1,iblk+1)-lcblkb(1,iblk) 550 do i = 1, npro ! loop over boundary elements (belts) 551 iBCB2=miBCB(iblk)%p(i,2) 552 if(iBCB2.ne.zero) then ! if a sid is given for this belt 553 if(nsidl.eq.0) then ! if this is the first sid we've seen 554 nsidl=1 ! increment our count and 555 sidlist % value = iBCB2 ! record its value 556 else ! if this isn't the first sid 557 newflag = 1 ! assume this is a new sid 558 sidelt => sidlist ! starting with the first sid 559 do j = 1, nsidl ! check the assumption 560 if((sidelt % value).eq.iBCB2) then 561 newflag=0 ! ... 562 endif 563 if(j.ne.nsidl) then ! ... 564 sidelt => sidelt % next 565 endif! ... 566 enddo 567 if(newflag.eq.1) then! if it really is new to us 568 nsidl = nsidl + 1 ! increment our count 569 allocate (sidelt % next)! tack a new element to the end 570 sidelt => sidelt % next! point to the new element 571 sidelt % value = iBCB2 ! record the new sid 572 endif ! (really is a new sid) 573 endif ! (first sid) 574 endif ! (belt has a sid) 575 enddo ! (loop over belts) 576 enddo ! (loop over boundary element blocks) 577c Copy the data from the linked list to a more easily-used array form 578 if(nsidl.gt.0) then 579 allocate( sidmapl(nsidl) ) 580 sidelt => sidlist ! starting with the first sid 581 do j = 1, nsidl 582 sidmapl(j)=sidelt%value 583 if(j.ne.nsidl) sidelt => sidelt%next 584 enddo 585 else 586 allocate( sidmapl(1)) ! some compilers/MPI don't like to send unallocated arrays 587 endif 588c Gather the number of surface ID's on each processor 589 if (numpe.gt.1) then ! multiple processors 590c write the number of surfID's on the jth processor to slot j of nsid 591 call MPI_ALLGATHER(nsidl,1,MPI_INTEGER,nsid,1, 592 & MPI_INTEGER,MPI_COMM_WORLD,ierr) 593c count up the total number of surfID's among all processes 594 nsidt=0 595 do j=1,numpe 596 nsidt=nsidt+nsid(j) 597 enddo 598 else ! single processor 599c the local information is the global information for single-processor 600 nsid=nsidl 601 nsidt=nsidl 602 endif ! if-else for multiple processors 603 if(nsidt.gt.0) then 604c 605c Make all-processor surfID collage 606c 607c there will be some duplicate surface ID's when we gather, so we 608c will use a temporary array 609 allocate (temp(nsidt)) 610 if (numpe.gt.1) then ! multiple processors 611c we will gather surfID's from local on-proc sets to a global set 612c we will stack each processor's surfID list atop that of the previous 613c processor. If the list for processor i is called sidmapi, then our 614c global coordinate list sidmap will look like this: 615c --------------------------------------------------------------- 616c | sidmap1 | sidmap2 | sidmap3 | ... | 617c --------------------------------------------------------------- 618c <---nsid(1)---> <-----nsid(2)-----> <-nsid(3)-> 619c <------------------------nsidt-----------------------...----> 620c To accomplish this with MPI, we use MPI_ALLGATHERV, summarized as: 621cMPI_ALLGATHERV(sendbuf,sendcount,sendtype,recvbuf,recvcount,disp,recvtype,comm) 622c[ IN sendbuf] starting address of send buffer (choice) 623c[ IN sendcount] number of elements in send buffer (integer) 624c[ IN sendtype] data type of send buffer elements (handle) 625c[ OUT recvbuf] address of receive buffer (choice) 626c[ IN recvcount] number of elements received from each process (int array) 627c[ IN disp] displacement array 628c[ IN recvtype] data type of receive buffer elements (handle) 629c[ IN comm] communicator (handle) 630c The displacement array disp is an array of integers in which the jth 631c entry indicates which slot of sidmap marks beginning of sidmapj 632c So, first we will build this displacement array 633 idisp(:)=0 ! starting with zero, since MPI likes C-numbering 634 do j=2,numpe 635 idisp(j)=idisp(j-1)+nsid(j-1) ! see diagram above 636 enddo 637c Now, we gather the data 638 call MPI_ALLGATHERV(sidmapl(:),nsidl, 639 & MPI_INTEGER,temp(:),nsid,idisp, 640 & MPI_INTEGER,MPI_COMM_WORLD,ierr) 641c sort surfID's, lowest to highest 642 isorted = 0 643 do while (isorted.eq.0) ! while the list isn't sorted 644 isorted = 1 ! assume the list is sorted this time 645 do j = 2, nsidt ! loop over ID's 646 if(temp(j).lt.temp(j-1)) then ! ID exceeds predecessor 647 itmp=temp(j-1) 648 temp(j-1)=temp(j) 649 temp(j)=itmp 650 isorted=0 ! not sorted yet 651 endif 652 enddo !loop over ID's 653 enddo ! while not sorted 654c determine the total number of distinct surfID's, globally 655 nsidg=nsidt ! assume there are no duplicate SurfID's 656 do j = 2, nsidt 657 if(temp(j).eq.temp(j-1)) nsidg = nsidg - 1 ! correction 658 enddo 659c create duplicate-free surfID list 660 allocate( sidmapg(nsidg) ) 661 sidmapg(1)=temp(1) 662 nsidg = 1 663 do j = 2, nsidt 664 if(temp(j).ne.temp(j-1)) then 665 nsidg = nsidg + 1 666 sidmapg(nsidg)=temp(j) 667 endif 668 enddo 669 deallocate( temp ) 670 else ! single-processor 671c global data is local data in single processor case 672 nsidg=nsidl 673 allocate( sidmapg(nsidg) ) 674 sidmapg=sidmapl 675 endif ! if-else multiple processor 676c 677 endif ! if at least one surfid 678c 679 return 680c 681 end 682 683 684 685 subroutine genotwn(x, BCtmp, iBC, nsurf) 686c 687c---------------------------------------------------------------------- 688c This routine determines which node to use as the first node off the 689c wall in near-wall modeling traction calculations for each wall node. 690c Each wall node has a normal, calculated from the wall elements to 691c which that node belongs. We seek the node within the boundary 692c element that lies closest to the line defined by the normal vector. 693c We create normalized vectors pointing from the wall node in 694c question to each of the nodes in the boundary element. The vector 695c that has the largest projection onto the wall node's normal vector 696c points to the node we want. Nodes that are not on a wall point to 697c themselves as their own off-the-wall node. 698c 699c input: 700c x (nshg,3) : : nodal position vectors 701c wnrm (nshg,3) : (mod) : normal vectors for each node 702c iBCB5 (nshg) : (file) : wall flag for belt 703c ienb (numelb,nen) : (file) : connectivity for belts 704c 705c output: 706c otwn (nshg) : (mod) : off-the-wall nodes for each node 707c 708c 709c Kenneth Jansen, Summer 2000 (algorithm) 710c Michael Yaworski, Summer 2000 (code) 711c---------------------------------------------------------------------- 712c 713 use pointer_data ! used for mienb, miBCB 714 use turbSA 715 include "common.h" 716c 717 integer iel, nod, can 718 real*8 vec(3), leng, dp, bigdp, lil 719 real*8 x(numnp,nsd),BCtmp(nshg,ndof+7) 720 integer iBC(nshg), nsurf(nshg) 721 integer gbits 722 integer, allocatable :: ienb(:) 723 724 allocate( otwn(nshg) ) 725c 726c initialize otwn to oneself 727c 728 do nod = 1, nshg 729 otwn(nod)=nod 730 enddo 731c 732c determine otwn for each wall node 733c 734 do iblk=1, nelblb ! loop over boundary element blocks 735 npro = lcblkb(1,iblk+1)-lcblkb(1,iblk) 736 nenl = lcblkb(5,iblk) 737 nenbl = lcblkb(6,iblk) 738 nshl = lcblkb(9,iblk) 739 allocate( ienb(nshl) ) 740 do i = 1, npro ! loop over boundary elements 741 iBCB1=miBCB(iblk)%p(i,1) 742 iBCB2=miBCB(iblk)%p(i,2) 743 ienb(1:nshl)=mienb(iblk)%p(i,1:nshl) 744 if (btest(iBCB1,4)) then ! (on the wall) 745 do nod = 1, nenbl ! for each wall node in this belt 746 bigdp = zero ! initialize largest dot product 747 do can=nenbl+1,nenl ! loop over off-the-wall candidates 748 nn=ienb(can) 749 ! don't bother with wall nodes 750 if(nsurf(nn).ne.0) cycle 751 ! don't bother with no-slip nodes 752 if(ibits(iBC(nn),3,3).eq.7 .and. 753 & BCtmp(nn,7).eq.zero) cycle 754c ! calculate candidate vector 755 vec(:)=x(ienb(can),:)-x(ienb(nod),:) 756 leng=sqrt(vec(1)**2+vec(2)**2+vec(3)**2) 757 vec(:)=vec(:)/leng 758c ! calculate dot product with wnrm 759 vec(:)=vec(:)*wnrm(ienb(nod),:) 760 dp=vec(1)+vec(2)+vec(3) 761c ! set candidate as otwn if necessary 762c ! (wnrm points into fluid, so 763c ! we want the most positive dp) 764 if (dp.gt.bigdp) then 765 otwn(ienb(nod)) = ienb(can) 766 bigdp=dp 767 endif 768 enddo !(loop over off-the-wall candidates) 769 enddo !(loop over wall nodes in current belt) 770 endif 771 enddo !(loop over elts in block) 772 deallocate(ienb) 773 enddo !(loop over boundary element blocks) 774 do nn = 1, nshg 775 if((otwn(nn).eq.nn).and.(nsurf(nn).ne.0)) then ! if a node on a 776 ! modeled surface 777 ! didn't find a node 778 ! off the wall 779 lil = 1.0e32 ! distance to current closest prospect 780 do can = 1, nshg ! loop over nodes 781 if(nsurf(can).eq.0) then ! if this candidate is off the 782 ! wall 783 if(ibits(iBC(can),3,3).eq.7 .and. 784 & BCtmp(can,7).eq.zero) then ! no slip nodes not allowed 785 else ! not a forbidden node 786 vec(:)=x(nn,:)-x(can,:) 787 leng=vec(1)**2+vec(2)**2+vec(3)**2 788 if(leng.lt.lil) then ! this candidate is closest so far 789 lil=leng 790 otwn(nn)=can 791 endif ! closest so far 792 endif ! end of no slip nodes 793 endif ! off-the-wall check 794 enddo ! loop over nodes 795 endif ! lonely wall-node check 796 enddo ! loop over nodes 797c 798 return 799c 800 end 801 802