1 subroutine bc3Res (y, iBC, BC, res, iper, ilwork) 2c 3c---------------------------------------------------------------------- 4c 5c This routine satisfies the BC of the residual vector for 3D elements. 6c 7c input: 8c y (nshg,ndof) : Y variables 9c iBC (nshg) : Boundary Condition Code 10c BC (nshg,ndofBC) : the boundary condition constraint parameters 11c res (nshg,nflow) : residual before BC is applied 12c 13c output: 14c res (nshg,nflow) : residual after satisfaction of BC 15c 16c 17c Thuc Bui, Winter 1989. 18c Zdenek Johan, Winter 1991. (Fortran 90) 19c---------------------------------------------------------------------- 20c 21 include "common.h" 22c 23 dimension y(nshg,ndof), iBC(nshg), 24 & BC(nshg,ndofBC), 25 & res(nshg,nflow), ilwork(nlwork), 26 & iper(nshg) 27c 28c.... density 29c 30 where (btest(iBC,0)) 31 res(:,5) = res(:,5) + BC(:,1)*Rgas* res(:,1) !IDEAL GAS ASSUMED 32 res(:,1) = zero 33 endwhere 34c 35c.... pressure 36c 37 if(EntropyPressure.eq.1) then 38 where (btest(iBC,2)) 39 40c Thought this would be correct if W was in the tangent space of V 41c as described in Shakib's thesis it does not work for primitive 42c variables 43c 44c Instead we use the entropy variable W 45c 46 res(:,2) = res(:,2) -y(:,1)*res(:,1) 47 res(:,3) = res(:,3) -y(:,2)*res(:,1) 48 res(:,4) = res(:,4) -y(:,3)*res(:,1) 49 res(:,5) = res(:,5) - 50 & (gamma*Rgas/gamma1*y(:,5) 51 & + pt5 * ( y(:,1)**2 + y(:,2)**2 + y(:,3)**2))*res(:,1) 52 res(:,1) = zero 53 endwhere 54 else 55 where (btest(iBC,2)) 56 res(:,1) = zero 57 endwhere 58 endif 59c 60c.... velocities 61c 62c ibits(n1,n2,n3) extracts bits n2+1 through n2+n3 (extending to the left 63c as is traditional in binary) of the integer n1 64c and returns the base 10 integer. In examples below x y z a b can 65c be 1 or zero without any effect. 66c 67c.... x1-velocity 68c 69c if iBC=4 bits of ibc =00000100 => ibits(4,3,3)=0 70c if iBC=40 bits of ibc =00101000 => ibits(4,3,3)=5 71c if iBC=40 bits of ibc =00101000 => ibits(4,3,2)=1 72c 73 where (ibits(iBC,3,3) .eq. 1) ! bits of iBC= xy001zab 74c 75c notice that the extracted 3 bits form the number 1. below 76c you will see the combinations which make up 2-7, all of the 77c possible velocity combinations 78c 79 res(:,3) = res(:,3) - BC(:,4) * res(:,2) 80 res(:,4) = res(:,4) - BC(:,5) * res(:,2) 81 res(:,2) = zero 82 endwhere 83c 84c.... x2-velocity 85c 86 where (ibits(iBC,3,3) .eq. 2) ! bits of iBC= xy010zab 87 res(:,2) = res(:,2) - BC(:,4) * res(:,3) 88 res(:,4) = res(:,4) - BC(:,5) * res(:,3) 89 res(:,3) = zero 90 endwhere 91c 92c.... x1-velocity and x2-velocity 93c 94 where (ibits(iBC,3,3) .eq. 3) ! bits of iBC= xy011zab 95 res(:,4) = res(:,4) - BC(:,4) * res(:,2) - BC(:,6) * res(:,3) 96 res(:,2) = zero 97 res(:,3) = zero 98 endwhere 99c 100c.... x3-velocity 101c 102 where (ibits(iBC,3,3) .eq. 4) ! bits of iBC= xy100zab 103 res(:,2) = res(:,2) - BC(:,4) * res(:,4) 104 res(:,3) = res(:,3) - BC(:,5) * res(:,4) 105 res(:,4) = zero 106 endwhere 107c 108c.... x1-velocity and x3-velocity 109c 110 where (ibits(iBC,3,3) .eq. 5) ! bits of iBC= xy101zab 111 res(:,3) = res(:,3) - BC(:,4) * res(:,2) - BC(:,6) * res(:,4) 112 res(:,2) = zero 113 res(:,4) = zero 114 endwhere 115c 116c.... x2-velocity and x3-velocity 117c 118 where (ibits(iBC,3,3) .eq. 6) ! bits of iBC= xy110zab 119 res(:,2) = res(:,2) - BC(:,4) * res(:,3) - BC(:,6) * res(:,4) 120 res(:,3) = zero 121 res(:,4) = zero 122 endwhere 123c 124c.... x1-velocity, x2-velocity and x3-velocity 125c 126 where (ibits(iBC,3,3) .eq. 7) ! bits of iBC= xy111zab 127 res(:,2) = zero 128 res(:,3) = zero 129 res(:,4) = zero 130 endwhere 131c 132c.... scaled plane extraction boundary condition 133c 134 if(intpres.eq.1) then ! interpolating pressure so zero continuity res 135 where (btest(iBC,11)) 136 res(:,1) = zero 137 res(:,2) = zero 138 res(:,3) = zero 139 res(:,4) = zero 140 res(:,5) = zero ! added to correspond to genscale (Elaine) 141 endwhere 142 else ! leave residual in continuity equation 143 where (btest(iBC,11)) 144 res(:,2) = zero 145 res(:,3) = zero 146 res(:,4) = zero 147 res(:,5) = zero ! added to correspond to genscale (Elaine) 148 endwhere 149 endif 150c 151c.... temperature 152c 153 where (btest(iBC,1)) res(:,5) = zero 154c 155c.... local periodic boundary conditions (no communications) 156c 157 do j = 1,nshg 158 if (btest(iBC(j),10)) then 159 i = iper(j) 160 res(i,:) = res(i,:) + res(j,:) 161 res(j,:) = zero 162 endif 163 enddo 164c 165c.... periodic slaves get the residual values of the masters 166c 167c do i = 1,nshg 168c if (btest(iBC(i),10)) then 169c res(i,:) = res(iper(i),:) 170c endif 171c enddo 172 if(numpe.gt.1) then 173c 174c.... nodes treated on another processor are eliminated 175c 176 numtask = ilwork(1) 177 itkbeg = 1 178 179 do itask = 1, numtask 180 181 iacc = ilwork (itkbeg + 2) 182 numseg = ilwork (itkbeg + 4) 183 184 if (iacc .eq. 0) then 185 do is = 1,numseg 186 isgbeg = ilwork (itkbeg + 3 + 2*is) 187 lenseg = ilwork (itkbeg + 4 + 2*is) 188 isgend = isgbeg + lenseg - 1 189 res(isgbeg:isgend,:) = zero 190 enddo 191 endif 192 193 itkbeg = itkbeg + 4 + 2*numseg 194 195 enddo 196 endif 197c 198c.... return 199c 200 return 201 end 202c 203c 204c 205 subroutine bc3ResSclr (y, iBC, BC, rest, 206 & iper, ilwork) 207c 208c---------------------------------------------------------------------- 209c 210c This routine satisfies the BC of the residual vector for 3D elements. 211c 212c input: 213c Y (nshg,ndof) : Y Variables 214c iBC (nshg) : Boundary Condition Code 215c BC (nshg,ndofBC) : the boundary condition constraint parameters 216c rest (nshg) : residual before BC is applied 217c 218c output: 219c rest (nshg) : residual after satisfaction of BC 220c 221c 222c Thuc Bui, Winter 1989. 223c Zdenek Johan, Winter 1991. (Fortran 90) 224c---------------------------------------------------------------------- 225c 226 include "common.h" 227c 228 dimension y(nshg,ndof), iBC(nshg), 229 & BC(nshg,ndofBC), 230 & rest(nshg), ilwork(nlwork), 231 & iper(nshg) 232c 233c 234 id = isclr+5 235c.... Scalar Variable 236c 237 where (btest(iBC,id)) 238 rest(:) = zero 239 endwhere 240c 241c.... local periodic boundary conditions (no communications) 242c 243 do j = 1,nshg 244 if (btest(iBC(j),10)) then 245 i = iper(j) 246 rest(i) = rest(i) + rest(j) 247 rest(j) = zero !changed 248 endif 249 enddo 250c 251c.... periodic slaves get the residual values of the masters 252c 253c$$$ do i = 1,nshg 254c$$$ if (btest(iBC(i),10)) then 255c$$$ rest(i) = rest(iper(i)) 256c$$$ endif 257c$$$ enddo 258c 259c removed for impl=4 as we have set the loops over ntopsh 260c 261 if(numpe.gt.1 ) then 262c 263c.... nodes treated on another processor are eliminated 264c 265 numtask = ilwork(1) 266 itkbeg = 1 267 268 do itask = 1, numtask 269 270 iacc = ilwork (itkbeg + 2) 271 numseg = ilwork (itkbeg + 4) 272 273 if (iacc .eq. 0) then 274 do is = 1,numseg 275 isgbeg = ilwork (itkbeg + 3 + 2*is) 276 lenseg = ilwork (itkbeg + 4 + 2*is) 277 isgend = isgbeg + lenseg - 1 278 rest(isgbeg:isgend) = zero 279 enddo 280 endif 281 282 itkbeg = itkbeg + 4 + 2*numseg 283 284 enddo 285 endif 286c 287c 288c.... return 289c 290 return 291 end 292 293