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 173 if(usingPETSc.eq.0) then !kill this code for petsc 174c 175c.... nodes treated on another processor are eliminated 176c 177 numtask = ilwork(1) 178 itkbeg = 1 179 180 do itask = 1, numtask 181 182 iacc = ilwork (itkbeg + 2) 183 numseg = ilwork (itkbeg + 4) 184 185 if (iacc .eq. 0) then 186 do is = 1,numseg 187 isgbeg = ilwork (itkbeg + 3 + 2*is) 188 lenseg = ilwork (itkbeg + 4 + 2*is) 189 isgend = isgbeg + lenseg - 1 190 res(isgbeg:isgend,:) = zero 191 enddo 192 endif 193 194 itkbeg = itkbeg + 4 + 2*numseg 195 196 enddo 197 endif 198 endif 199c 200c.... return 201c 202 return 203 end 204c 205c 206c 207 subroutine bc3ResSclr (y, iBC, BC, rest, 208 & iper, ilwork) 209c 210c---------------------------------------------------------------------- 211c 212c This routine satisfies the BC of the residual vector for 3D elements. 213c 214c input: 215c Y (nshg,ndof) : Y Variables 216c iBC (nshg) : Boundary Condition Code 217c BC (nshg,ndofBC) : the boundary condition constraint parameters 218c rest (nshg) : residual before BC is applied 219c 220c output: 221c rest (nshg) : residual after satisfaction of BC 222c 223c 224c Thuc Bui, Winter 1989. 225c Zdenek Johan, Winter 1991. (Fortran 90) 226c---------------------------------------------------------------------- 227c 228 include "common.h" 229c 230 dimension y(nshg,ndof), iBC(nshg), 231 & BC(nshg,ndofBC), 232 & rest(nshg), ilwork(nlwork), 233 & iper(nshg) 234c 235c 236 id = isclr+5 237c.... Scalar Variable 238c 239 where (btest(iBC,id)) 240 rest(:) = zero 241 endwhere 242c 243c.... local periodic boundary conditions (no communications) 244c 245 do j = 1,nshg 246 if (btest(iBC(j),10)) then 247 i = iper(j) 248 rest(i) = rest(i) + rest(j) 249 rest(j) = zero !changed 250 endif 251 enddo 252c 253c.... periodic slaves get the residual values of the masters 254c 255c$$$ do i = 1,nshg 256c$$$ if (btest(iBC(i),10)) then 257c$$$ rest(i) = rest(iper(i)) 258c$$$ endif 259c$$$ enddo 260c 261c removed for impl=4 as we have set the loops over ntopsh 262c 263 if(numpe.gt.1 ) then 264 if(usingPETSc.eq.0) then !kill this code for petsc 265c 266c.... nodes treated on another processor are eliminated 267c 268 numtask = ilwork(1) 269 itkbeg = 1 270 271 do itask = 1, numtask 272 273 iacc = ilwork (itkbeg + 2) 274 numseg = ilwork (itkbeg + 4) 275 276 if (iacc .eq. 0) then 277 do is = 1,numseg 278 isgbeg = ilwork (itkbeg + 3 + 2*is) 279 lenseg = ilwork (itkbeg + 4 + 2*is) 280 isgend = isgbeg + lenseg - 1 281 rest(isgbeg:isgend) = zero 282 enddo 283 endif 284 285 itkbeg = itkbeg + 4 + 2*numseg 286 287 enddo 288 endif 289 endif 290c 291c 292c.... return 293c 294 return 295 end 296 297