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