1 module wallData !was suctionDuct 2 logical wnormCalced 3 !integer nsuctionface 4 !integer, allocatable :: suctionf(:) 5 6 !Threshold for determining whether a node is on a corner 7 real*8 :: cornerCosThresh 8 9 !Wall normal array of size (nshg, 3) 10 real*8, allocatable :: wnorm(:,:) 11 contains 12 13 subroutine findWallNorm(x,iBC,ilwork,iper) 14 15! use wallData 16 use pointer_data 17 include "common.h" 18 include "mpif.h" 19 20 integer iBC(nshg) 21 real*8 x(numnp,nsd) 22 integer iper(nshg) 23 integer ilwork(nlwork) 24 integer, allocatable :: ienb(:) 25 real*8, allocatable :: wNormTmp(:,:) 26 real*8 eNorm(nsd),en(nsd),wn(nsd),aedg(nsd),bedg(nsd) 27 real*8 mag 28! integer isfID 29 30 real*8 wallBCout(nshg, 6) 31 real*8 BC(nshg,ndofBC) 32 33 cornerCosThresh = 0.75d00 34 !The number 0.75 is somewhat arbitrary. The thinking is that this 35 !code is also used for the normal on curved surfaces which will have 36 !a non unit dot product from elements surrounding a wall vertex but 37 !acos(0.75) is a pretty big angle. On coarse meshes this number may 38 !have to be dialed down to get normal on curved walls correct. 39 40 allocate(wNorm(nshg,nsd)) 41 wNorm = 0 42 43 do iblk=1,nelblb !loop over all the bounday element blocks at current process 44 npro=lcblkb(1,iblk+1)-lcblkb(1,iblk) 45 nenbl=lcblkb(6,iblk) 46 nshl=lcblkb(9,iblk) 47 allocate(ienb(nshl)) 48 do i=1,npro ! loop over boundary elements 49! isfID=miBCB(iblk)%p(i,2) ! miBCB(2) is the surfID 50 ienb(1:nshl)=mienb(iblk)%p(i,1:nshl) 51 do j=1,nenbl ! loop over boundary nodes 52 nn=ienb(j) !nn is the global index of suction node 53 if(j.ne.1.and.j.ne.nenbl)then 54 aedg(:)=x(ienb(j+1),:)-x(nn,:) 55 bedg(:)=x(ienb(j-1),:)-x(nn,:) 56 elseif(j.eq.1)then 57 aedg(:)=x(ienb(j+1),:)-x(nn,:) 58 bedg(:)=x(ienb(nenbl),:)-x(nn,:) 59 elseif(j.eq.nenbl)then 60 aedg(:)=x(ienb(1),:)-x(nn,:) 61 bedg(:)=x(ienb(j-1),:)-x(nn,:) 62 endif 63 eNorm(1)=aedg(2)*bedg(3)-aedg(3)*bedg(2) 64 eNorm(2)=aedg(3)*bedg(1)-aedg(1)*bedg(3) 65 eNorm(3)=aedg(1)*bedg(2)-aedg(2)*bedg(1) 66 67 call addElemNormal(wNorm, eNorm, nn) 68 69 enddo !end do j = 1, nenbl 70 enddo !end do i = 1, npro 71 deallocate(ienb) 72 enddo !end do iblk = 1, nelblk 73 74 !Now communicate between parts. Note: there still may be a 75 !conflict if the boundary edge is on a corner, so a separate array 76 !needs to be used to allow us to repeat the z-norm comparison. 77 if(numpe.gt.1) then 78 allocate(wNormTmp(nshg,nsd)) 79 wNormTmp = wNorm 80 call commu(wNormTmp(:,:), ilwork, 3, 'in ') 81 wNormTmp = wNormTmp - wNorm !undo the addition in commu to just 82 !get contributions from other processors 83 84 do i = 1,nshg !loop over each node and add the off 85 wn(:) = wNormTmp(i,:) !processor contribution. If this is 86 call addElemNormal(wNorm, wn, i) !an interior node, zero gets 87 enddo !added to zero. 88 endif 89 90 !It is unclear how to treat periodic boundary conditions, i.e. 91 !should periodic walls have zero norm, a correct normal for each 92 !side, or a correct normal for the master. With bc3per commented, 93 !the wall normal should be calculated as though there is no 94 !periodic BC. 95! call bc3per(iBC,wNorm(:,:),iper,ilwork,3) 96 97 if(numpe.gt.1) call commu(wNorm(:,:),ilwork,3,'out') 98 99 do nn = 1,nshg 100 wn(:) = wNorm(nn,:) 101 mag = (wn(1)**2+wn(2)**2+wn(3)**2)**0.5 102 103 if(mag > 0) wNorm(nn,:) = wn(:)/mag 104 enddo 105 106 zwall=0.0762d00 107 do nn=1,nshg ! normalize wNorm 108 wn(:)=wNorm(nn,:) 109 mag=(wn(1)**2+wn(2)**2+wn(3)**2)**0.5 110 if(mag.ne.0) then 111 if(((dabs(dabs(x(nn,3))-zwall)).lt.0.000000001d00) 112 & .and. (dabs(wn(3)/mag) .lt.0.99)) then 113! logic to trap wall normals in z above failed, possibley in commu? so 114! report it and "fix it" 115 write(*,*) 116 & "Warning: still screwing up wall Norm, proc ", 117 & myrank, ", node ", nn, 118 & ", x = (", x(nn,1), x(nn,2), x(nn,3), ")", 119 & ", n = (", wn(1), wn(2), wn(3), ")" 120! write(myrank+1000,1234) 121! & nn,x(nn,1),x(nn,2),x(nn,3),wn(1),wn(2),wn(3) 122 wn(1)=0 123 wn(2)=0 124 wn(3)=x(nn,3) 125 mag=abs(wn(3)) 126 endif 127 wNorm(nn,:)=wn(:)/mag 128 endif 129 enddo 130 1311234 format(i5,6(2x,e12.5)) 132!!DEBUG 133!! dimension u(nshg,n) 134! do j = 1,3 135! do i = 1,nshg 136! u(i,j)=9.876543e21 137! wallBCout(i, j ) = BC( i, j+2) 138! wallBCout(i, j+3) = wNorm(i, j ) 139! enddo 140! enddo 141! call write_restart(myrank,100002,nshg,6,wallBCout,wallBCout) 142!!DEBUG 143 144 wNormCalced = .true. 145 146 return 147 end 148 149 subroutine addElemNormal(wNorm, eNorm, nn) 150 !Adds contributions from the wall element normal en into the wall 151 !normal array wNorm at index nn. 152 153 include "common.h" 154 155 real*8 :: wNorm(nshg, nsd) 156 real*8 :: eNorm(nsd), en(nsd), wn(nsd) 157 real*8 :: mag, oDotn 158 integer :: nn 159 160 wn(:)=wNorm(nn,:) 161 mag=(wn(1)**2+wn(2)**2+wn(3)**2)**0.5 162 163 if(mag < 1.0e-14) then ! no direction yet so give it this 164 wNorm(nn,:) = eNorm(:) 165 166 else ! wNorm already has a direction 167 wn(:)=wn(:)/mag 168 169 mag=(eNorm(1)**2+eNorm(2)**2+eNorm(3)**2)**0.5 170 en(:)=eNorm(:)/mag 171 172 ! when we are in corners we need to establish a precidence since 173 ! averaging will produce irregular results. In this first pass, 174 ! let's let the vector with the largest z-normal take precedence. 175 oDotn=en(1)*wn(1)+en(2)*wn(2)+en(3)*wn(3) 176 if(abs(oDotn) < cornerCosThresh) then ! this is a corner-Don't average 177 if(abs(en(3)) > abs(wn(3))) then 178 wNorm(nn,:) = eNorm(:) ! overwrite wNorm 179 endif ! since no else, wNorm unaffected by weaker z 180 ! e.g., for the else, eNorm is ignored 181 else ! not a corner so average 182 wNorm(nn,:)=wNorm(nn,:) + eNorm(:) 183 endif !end if(abs(oDotn) < cornerCosThresh) 184 endif !end mag < 1e-14 185 186 end subroutine 187 188 end module 189 190