xref: /phasta/phSolver/common/setSuction_Duct3.f90 (revision 712d3df0b59ebebaaeaea358162c8d2c043c6e08)
1c========================================================================
2c Dynamically controll the suction level based on main contraction mdot,
3c called by itrdrv.f
4c========================================================================
5        subroutine setSuction_Duct3(x,BC,y, ilwork)
6
7        use wallData ! wnorm
8        include "common.h"
9        include "mpif.h"
10        include "auxmpi.h"
11
12        integer i, j, k, nn, nL
13        integer ilwork(nlwork)
14        integer nSuctionNode
15        integer, allocatable :: suctionNodeMap(:)
16!        logical inPatch
17        real*8 xcoor,ycoor,zcoor
18        real*8 BC(nshg,ndofBC)
19        real*8 BC3(numnp, 5)                      !used for hack to get suction right on part boundaries
20        real*8 x(nshg,nsd), y(nshg,ndof)
21!       real*8 xmin_t, xmax_t, xmin_b, xmax_b     !ramp geometry controls
22!       real*8 xmin_su, xmax_su, xmin_sl, xmax_sl !side wall patch limits
23!       real*8 t_upper                            !upper thickness of suction area
24!       real*8 ymin_su, ymax_sl
25        real*8 nvel_s, nvel_t, nvel_b             !normal velocity on the side walls, top, and bottom
26
27!        real*8 delX, delY, X1, Y1     !Variables used to map upper surface or trumpet
28!        real*8 Xline(30,2)
29!        integer suctionIDs(4)
30
31!        real*8 wallBCout(nshg, 6)
32
33        if(myrank.eq.0)write(*,*)'Setting side suctions'
34
35        nvel_b     = suctionVbottom      !normal velocity at the bottom patch
36        nvel_sl    = suctionVside_lower  !normal velocity of lower side wall patches
37        nvel_su    = suctionVside_upper  !normal velocity of upper side wall patches
38        nvel_t     = suctionVtop         !normal velocity at the top patch
39
40        allocate(suctionNodeMap(nshg))
41        call sfID2np(isetSuctionID_Duct, nSuctionNode, suctionNodeMap)  !to figure out whether it's the top, bottom, or side surface.
42
43        smallNum=1.0e-5
44        xsideBotMax=0.0428625+smallNum
45        xsideTopMax=0.130175+smallNum
46
47        do i = 1, nSuctionNode
48          nn = suctionNodeMap(i) ! map face node index to global node
49          xcoor = x(nn,1)
50          ycoor = x(nn,2)
51          !zcoor = x(nn,3)
52
53          !wnormx = wnorm(nn, 1)
54          wnormy = wnorm(nn, 2)
55          wnormz = wnorm(nn, 3)
56
57          !Test if the point lies in one of the suction patches
58
59          if((abs(nvel_sl) >  0  .or. abs(nvel_su) > 0) .and.  !if either upper or lower side walls are on...
60     &       (abs(wnormz) >=  0.999)) then                     ! and the node normal is in the z direction
61            if((ycoor > 0) .and. (xcoor .le. xsideTopMax)) then !hack to prevent edges from turning on.
62              BC(nn, 3) = nvel_su*wnorm(nn, 1) ! set xVel
63              BC(nn, 4) = nvel_su*wnorm(nn, 2) ! set yVel
64              BC(nn, 5) = nvel_su*wnorm(nn, 3) ! set zVel
65            elseif (xcoor .le. xsideBotMax) then
66              BC(nn, 3) = nvel_sl*wnorm(nn, 1) ! set xVel
67              BC(nn, 4) = nvel_sl*wnorm(nn, 2) ! set yVel
68              BC(nn, 5) = nvel_sl*wnorm(nn, 3) ! set zVel
69            endif
70          else if(abs(wnormz) > 0.1) then      ! side walls should have a larger znormal than this - it's ambiguous what to do so do nothing
71            cycle
72          else if(abs(nvel_t) > 0      .and. !top suction patch
73     &                 wnormy > 0.001) then  !face is pointing up (in the +y direction)
74            BC(nn, 3) = nvel_t*wnorm(nn, 1) ! set xVel using wall normal
75            BC(nn, 4) = nvel_t*wnorm(nn, 2) ! set yVel
76            BC(nn, 5) = nvel_t*wnorm(nn, 3) ! set zVel
77
78          else if(abs(nvel_b) > 0       .and. !bottom suction patches
79     &                wnormy  < -0.001) then  !face is pointing down (in the -y direction)
80            BC(nn, 3) = nvel_b*wnorm(nn, 1) ! set xVel using wall normal
81            BC(nn, 4) = nvel_b*wnorm(nn, 2) ! set yVel
82            BC(nn, 5) = nvel_b*wnorm(nn, 3) ! set zVel
83          endif
84
85        enddo
86
87        BC3=zero
88        if(numpe.gt.1) then
89            BC3(:,1:3)=BC(:,3:5)
90            do i=1,nshg
91                bc3mag=BC3(i,1)**2+BC3(i,2)**2+BC3(i,3)**2
92                if(bc3mag.gt.0) BC3(i,4)=one
93            enddo
94            call commu(BC3,ilwork,5,'in ')
95! This accumulated and failed                 call commu(BC(:,3:5),ilwork,3,'in ')
96            do i=1,nshg
97                bcmag=BC(i,3)**2+BC(i,4)**2+BC(i,5)**2
98                if((bcmag.eq.0).and.(BC3(i,4).ne.0)) then !node is slave and has not been correctly set.
99                     BC(i,3:5)=BC3(i,1:3)/BC3(i,4)  !division by BC3 is necessary to account for more than two blocks sharing the same node
100                endif
101            enddo
102       endif
103
104        !Debugging loop
105!        if(nSuctionNode.gt.0) then
106!          icount=0
107!          do i=1,nSuctionNode
108!             nn=suctionNodeMap(i)
109!             dz=dabs(dabs(x(nn,3))-0.0762)
110!             if(dz.lt.0.0001) then
111!               icount=icount+1
112!               if(icount.eq.1)
113!     & write(myrank+1000,*) 'writing wnorm and BC setSuctionNG3'
114!               write(myrank+1000,1234) nn,i,
115!     &           x(nn,1), x(nn,2), x(nn,3),
116!     &           wnorm(nn,1),wnorm(nn,2),wnorm(nn,3),
117!     &           BC(nn,3),BC(nn,4),BC(nn,5)
118!             endif
119!          enddo
120!         if(icount.gt.0) close(myrank+1000)
121!        endif
122!1234           format(i5,2x,i5,9(2x,e14.7))
123
124      end
125
126
127
128
129
130
131