xref: /phasta/phSolver/common/slpwBC.f (revision 712d3df0b59ebebaaeaea358162c8d2c043c6e08) 
1       subroutine slpwBC(nslpw,ihis,iBCg,BCg,BCtmpg,iBC,BC,wlnorm)
2       include "common.h"
3       include "mpif.h"
4
5       integer nslpw,ihis,iBCg,iBC
6       real*8 BCg(ndofBC)
7       real*8 BCtmpg(ndof+7)
8       real*8 BC(ndofBC)
9       real*8 wlnorm(3,9)
10       real*8 c4,c5,c6,c7
11       real*8 c8,c9,c10,c11
12       real*8 c12,c13,c14,c15
13       real*8 ck1,ck2,ck3
14       real*8 det,det3x3,tmp,u,v,w
15       integer nmax
16       real*8 wn1(3),wn2(3),wn3(3)
17
18       icd      = ibits(iBCg,3,3)
19       ixset    = ibits(iBCg,3,1)
20       iyset    = ibits(iBCg,4,1)
21       izset    = ibits(iBCg,5,1)
22C
23C...   release the previous setting for velocities
24C
25       iBC      = iBCg-(ixset*8+iyset*16+izset*32)
26
27
28ccc case1: one slipwall and no user-specified velocity BC
29       if(nslpw.eq.1.and.icd.eq.0)then
30        do islpw=1,9
31           if(btest(ihis,islpw))exit
32        enddo
33        wn1(:)=wlnorm(:,islpw)
34        ii=nmax(wn1(1),wn1(2),wn1(3))
35        if(ii.eq.1)then  ! u has the largest constrain, u should be constrained
36         iBC  = iBC+8
37         BC(3)= 0
38         BC(4)= wn1(2)/wn1(1)
39         BC(5)= wn1(3)/wn1(1)
40        elseif(ii.eq.2)then
41         iBC=iBC+16
42         BC(3)=0
43         BC(4)=wn1(1)/wn1(2)
44         BC(5)=wn1(3)/wn1(2)
45        elseif(ii.eq.3)then
46         iBC=iBC+32
47         BC(3)=0
48         BC(4)=wn1(1)/wn1(3)
49         BC(5)=wn1(2)/wn1(3)
50        endif
51
52ccc case2: two slipwall and no user-specified velocity BC
53       elseif(nslpw.eq.2.and.icd.eq.0)then
54        do islpw=1,9
55           if(btest(ihis,islpw))exit
56        enddo
57        do jslpw=islpw+1,9
58           if(btest(ihis,jslpw))exit
59        enddo
60        wn1(:)=wlnorm(:,islpw)
61        wn2(:)=wlnorm(:,jslpw)
62        c4=wn1(1)
63        c5=wn1(2)
64        c6=wn1(3)
65        c7=0
66        c8=wn2(1)
67        c9=wn2(2)
68        c10=wn2(3)
69        c11=0
70C
71C... ck is the cross product of wn1 and wn2
72C
73        ck1=c5*c10 - c9*c6
74        ck2=c6*c8  - c4*c10
75        ck3=c4*c9  - c8*c5
76
77        kk=nmax(ck1,ck2,ck3)
78
79        if(kk.eq.1)then  ! u has the largest freedom, v and w should be constrained
80C
81          iBC=iBC + 48
82          det=c5*c10-c9*c6
83          BC(3)=(c10*c7-c6*c11)/det
84          BC(4)=(c10*c4-c6*c8)/det
85          BC(5)=(c11*c5-c9*c7)/det
86          BC(6)=(c5*c8-c9*c4)/det
87        elseif(kk.eq.2)then
88          iBC=iBC+40
89          det=c4*c10-c8*c6
90          BC(3)=(c10*c7-c6*c11)/det
91          BC(4)=(c10*c5-c6*c9)/det
92          BC(5)=(c11*c4-c8*c7)/det
93          BC(6)=(c4*c9-c5*c8)/det
94        elseif(kk.eq.3)then
95          iBC=iBC+24
96          det=c4*c9-c8*c5
97          BC(3)=(c9*c7-c5*c11)/det
98          BC(4)=(c9*c6-c5*c10)/det
99          BC(5)=(c4*c11-c8*c7)/det
100          BC(6)=(c4*c10-c8*c6)/det
101        endif
102
103ccc case3: one slipwall and one user-specified velocity BC
104       elseif(nslpw.eq.1.and.(icd.eq.1.or.icd.eq.2.or.icd.eq.4))then
105        do islpw=1,9
106           if(btest(ihis,islpw))exit
107        enddo
108        wn1(:)=wlnorm(:,islpw)
109        c4=BCtmpg(4)
110        c5=BCtmpg(5)
111        c6=BCtmpg(6)
112        c7=BCtmpg(7)
113        c8=wn1(1)
114        c9=wn1(2)
115        c10=wn1(3)
116        c11=0
117
118        ck1=c5*c10 - c9*c6
119        ck2=c6*c8  - c4*c10
120        ck3=c4*c9  - c8*c5
121
122        kk=nmax(ck1,ck2,ck3)
123        if(kk.eq.1)then  ! u has the largest freedom
124          iBC=iBC+48
125          det=c5*c10-c9*c6
126          BC(3)=(c10*c7-c6*c11)/det
127          BC(4)=(c10*c4-c6*c8)/det
128          BC(5)=(c11*c5-c9*c7)/det
129          BC(6)=(c5*c8-c9*c4)/det
130        elseif(kk.eq.2)then
131          iBC=iBC+40
132          det=c4*c10-c8*c6
133          BC(3)=(c10*c7-c6*c11)/det
134          BC(4)=(c10*c5-c6*c9)/det
135          BC(5)=(c11*c4-c8*c7)/det
136          BC(6)=(c4*c9-c5*c8)/det
137        elseif(kk.eq.3)then
138          iBC=iBC+24
139          det=c4*c9-c8*c5
140          BC(3)=(c9*c7-c5*c11)/det
141          BC(4)=(c9*c6-c5*c10)/det
142          BC(5)=(c4*c11-c8*c7)/det
143          BC(6)=(c4*c10-c8*c6)/det
144        endif
145
146ccc case4: two slipwalls and one user-specified velocity BC
147       elseif(nslpw.eq.2.and.(icd.eq.1.or.icd.eq.2.or.icd.eq.4))then
148        iBC=iBC+56
149        c4=BCtmpg(4)
150        c5=BCtmpg(5)
151        c6=BCtmpg(6)
152        c7=BCtmpg(7)
153        do islpw=1,9
154           if(btest(ihis,islpw))exit
155        enddo
156        do jslpw=islpw+1,9
157           if(btest(ihis,jslpw))exit
158        enddo
159        wn1(:)=wlnorm(:,islpw)
160        wn2(:)=wlnorm(:,jslpw)
161        c8=wn1(1)
162        c9=wn1(2)
163        c10=wn1(3)
164        c11=0
165        c12=wn2(1)
166        c13=wn2(2)
167        c14=wn2(3)
168        c15=0
169        tmp=det3x3(c4,c5,c6,c8,c9,c10,c12,c13,c14)
170        BC(4)=det3x3(c4,c7,c6,c8,c11,c10,c12,c15,c14)/tmp
171        BC(3)=det3x3(c7,c5,c6,c11,c9,c10,c15,c13,c14)/tmp
172        BC(5)=det3x3(c4,c5,c7,c8,c9,c11,c12,c13,c15)/tmp
173
174ccc case5: one slipwall and two user-specified velocity BC
175       elseif(nslpw.eq.1.and.(icd.eq.3.or.icd.eq.5.or.icd.eq.6))then
176        iBC=iBC+56
177        c4=BCtmpg(4)
178        c5=BCtmpg(5)
179        c6=BCtmpg(6)
180        c7=BCtmpg(7)
181        c8=BCtmpg(8)
182        c9=BCtmpg(9)
183        c10=BCtmpg(10)
184        c11=BCtmpg(11)
185        do islpw=1,9
186           if(btest(ihis,islpw))exit
187        enddo
188        wn1(:)=wlnorm(:,islpw)
189        c12=wn1(1)
190        c13=wn1(2)
191        c14=wn1(3)
192        c15=0
193        tmp=det3x3(c4,c5,c6,c8,c9,c10,c12,c13,c14)
194        BC(4)=det3x3(c4,c7,c6,c8,c11,c10,c12,c15,c14)/tmp
195        BC(3)=det3x3(c7,c5,c6,c11,c9,c10,c15,c13,c14)/tmp
196        BC(5)=det3x3(c4,c5,c7,c8,c9,c11,c12,c13,c15)/tmp
197
198ccc case6: two slipwalls and two user-specified velocity BC
199
200ccc case7: one slipwall and three user-specified velocity BC
201       elseif(nslpw.eq.1.and.icd.eq.7)then
202        iBC=iBC+56
203        u=BC(3)
204        v=BC(4)
205        w=BC(5)
206        do islpw=1,9
207           if(btest(ihis,islpw))exit
208        enddo
209        wn1(:)=wlnorm(:,islpw)
210        ii=nmax(wn1(1),wn1(2),wn1(3))
211        if(ii.eq.1)then  ! u has the largest constrain
212         BC(3)=-wn1(2)/wn1(1)*v-wn1(3)/wn1(1)*w  ! u is determined by v and w
213        elseif(ii.eq.2)then ! v has the largest constrain
214         BC(4)=-wn1(1)/wn1(2)*u-wn1(3)/wn1(2)*w  ! v is determined by u and w
215        elseif(ii.eq.3)then ! w has the largest constrain
216         BC(5)=-wn1(1)/wn1(3)*u-wn1(2)/wn1(3)*v  ! w is determined by u and v
217        endif
218
219ccc case8: two slipwall and three user-specified velocity BC
220       elseif(nslpw.eq.2.and.icd.eq.7)then
221        iBC=iBC+56
222        u=BC(3)
223        v=BC(4)
224        w=BC(5)
225        do islpw=1,9
226           if(btest(ihis,islpw))exit
227        enddo
228        do jslpw=islpw+1,9
229           if(btest(ihis,jslpw))exit
230        enddo
231        wn1(:)=wlnorm(:,islpw)
232        wn2(:)=wlnorm(:,jslpw)
233        c4=wn1(1)
234        c5=wn1(2)
235        c6=wn1(3)
236        c7=0
237        c8=wn2(1)
238        c9=wn2(2)
239        c10=wn2(3)
240        c11=0
241
242        ck1=c5*c10 - c9*c6     ! ck is the cross product of wn1 and wn2
243        ck2=c6*c8  - c4*c10
244        ck3=c4*c9  - c8*c5
245
246        kk=nmax(ck1,ck2,ck3)
247        if(kk.eq.1)then  ! u has the largest freedom
248          det=c5*c10-c9*c6
249          BC(4)=(c10*c7-c6*c11)/det-(c10*c4-c6*c8)/det*u ! v is determined by u
250          BC(5)=(c11*c5-c9*c7)/det-(c5*c8-c9*c4)/det*u   ! w is determined by u
251        elseif(kk.eq.2)then ! v has the largest freedom
252          det=c4*c10-c8*c6
253          BC(3)=(c10*c7-c6*c11)/det-(c10*c5-c6*c9)/det*v ! u is determined by v
254          BC(5)=(c11*c4-c8*c7)/det-(c4*c9-c5*c8)/det*v   ! w is determined by v
255        elseif(kk.eq.3)then  ! w has the largest freedom
256          det=c4*c9-c8*c5
257          BC(3)=(c9*c7-c5*c11)/det-(c9*c6-c5*c10)/det*w  ! u is determined by w
258          BC(4)=(c4*c11-c8*c7)/det-(c4*c10-c8*c6)/det*w  ! v is determined by w
259        endif
260
261ccc case9: more than three slipwalls is acutally non-slip wall, regardless of user-specified velocity BC
262       elseif(nslpw.ge.3)then
263          iBC=iBC+56
264          BC(3:5)=0
265ccc
266
267       endif
268
269       return
270       end
271
272C-------------------------------------------------
273C function nmax
274C     returns the index at which the input number
275C     has the largest absolute value.
276C     For example nmax( 0.1, -3.4, 1.2  ) = 2
277C-------------------------------------------------
278
279       integer function nmax(nx,ny,nz)
280       real*8 nx,ny,nz
281       real*8 lnx,lny,lnz
282          lnx=abs(nx)
283          lny=abs(ny)
284          lnz=abs(nz)
285          if(lnx.gt.lny.and.lnx.gt.lnz)nmax=1
286          if(lny.gt.lnx.and.lny.gt.lnz)nmax=2
287          if(lnz.gt.lnx.and.lnz.gt.lny)nmax=3
288          return
289       end
290
291C-----------------------------------------------------------
292C function det3x3 computes the determinant of a 3x3 matrix
293C-----------------------------------------------------------
294
295       real*8 function det3x3(a1,a2,a3,a4,a5,a6,a7,a8,a9)
296       real*8 a1,a2,a3,a4,a5,a6,a7,a8,a9
297         det3x3=a1*a5*a9+a4*a8*a3+a7*a6*a2-a7*a5*a3-a8*a6*a1-a9*a4*a2
298       return
299       end
300
301
302
303
304
305
306
307
308