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