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