1 subroutine bc3BDg (y, iBC, BC, BDiag, iper, ilwork) 2c 3c---------------------------------------------------------------------- 4c 5c This routine satisfies the BC of the block-diagonal preconditioning 6c matrix for 3D elements. 7c 8c input: 9c y (nshg,ndof) : Y variables 10c iBC (nshg) : boundary condition code 11c BC (nshg,ndofBC) : Dirichlet BC constraint parameters 12c BDiag (nshg,nflow,nflow) : preconditionning matrix before BC 13c (only upper part) 14c 15c output: 16c BDiag (nshg,nflow,nflow) : preconditionning matrix after BC 17c is satisfied 18c 19c 20c Zdenek Johan, Summer 1990. (Modified from g3bce.f) 21c Zdenek Johan, Winter 1991. (Fortran 90) 22c---------------------------------------------------------------------- 23c 24 include "common.h" 25c 26 dimension y(nshg,ndof), iBC(nshg), 27 & BC(nshg,ndofBC), 28 & BDiag(nshg,nflow,nflow), ilwork(nlwork), 29 & iper(nshg) 30c 31 real*8 a5(nshg) 32c 33c.... density 34c 35 do i = 1, nshg 36 a5(i) = - y(i,5) * (Rgas * gamma / gamma1) !IDEAL GAS ASSUMED 37 end do 38 39 where (btest(iBC,0)) 40c 41c.... engbc was replaced for a5 by following 42 43 BDiag(:,5,5) = BDiag(:,5,5) + a5 * a5 * BDiag(:,1,1) + 44 & a5 * BDiag(:,1,5) + 45 & a5 * BDiag(:,5,1) 46 BDiag(:,4,5) = BDiag(:,4,5) + a5 * BDiag(:,4,1) 47 BDiag(:,3,5) = BDiag(:,3,5) + a5 * BDiag(:,3,1) 48 BDiag(:,2,5) = BDiag(:,2,5) + a5 * BDiag(:,2,1) 49 BDiag(:,5,4) = BDiag(:,5,4) + a5 * BDiag(:,1,4) 50 BDiag(:,5,3) = BDiag(:,5,3) + a5 * BDiag(:,1,3) 51 BDiag(:,5,2) = BDiag(:,5,2) + a5 * BDiag(:,1,2) 52 BDiag(:,1,2) = zero 53 BDiag(:,1,3) = zero 54 BDiag(:,1,4) = zero 55 BDiag(:,1,5) = zero 56 BDiag(:,2,1) = zero 57 BDiag(:,3,1) = zero 58 BDiag(:,4,1) = zero 59 BDiag(:,5,1) = zero 60 BDiag(:,1,1) = one 61 endwhere 62 63c where (btest(iBC,11)) ! pressure to deactivate 64 where (btest(iBC,2)) ! pressure 65 66 BDiag(:,1,2) = zero 67 BDiag(:,1,3) = zero 68 BDiag(:,1,4) = zero 69 BDiag(:,1,5) = zero 70 BDiag(:,2,1) = zero 71 BDiag(:,3,1) = zero 72 BDiag(:,4,1) = zero 73 BDiag(:,5,1) = zero 74 BDiag(:,1,1) = one 75 endwhere 76 77c 78c.... velocities 79c 80c.... x1-velocity 81c 82 where (ibits(iBC,3,3) .eq. 1) 83 BDiag(:,5,4) = BDiag(:,5,4) - BC(:,5) * BDiag(:,5,2) 84 BDiag(:,5,3) = BDiag(:,5,3) - BC(:,4) * BDiag(:,5,2) 85 86 BDiag(:,4,5) = BDiag(:,4,5) - BC(:,5) * BDiag(:,2,5) 87 BDiag(:,3,5) = BDiag(:,3,5) - BC(:,4) * BDiag(:,2,5) 88 89 BDiag(:,4,1) = BDiag(:,4,1) - BC(:,5) * BDiag(:,2,1) 90 BDiag(:,3,1) = BDiag(:,3,1) - BC(:,4) * BDiag(:,2,1) 91 92 BDiag(:,1,4) = BDiag(:,1,4) - BC(:,5) * BDiag(:,1,2) 93 BDiag(:,1,3) = BDiag(:,1,3) - BC(:,4) * BDiag(:,1,2) 94 95 BDiag(:,4,4) = BDiag(:,4,4) + BC(:,5) * BC(:,5) * BDiag(:,2,2) 96 & - BC(:,5) * BDiag(:,2,4) 97 & - BC(:,5) * BDiag(:,4,2) 98 BDiag(:,3,4) = BDiag(:,3,4) + BC(:,4) * BC(:,5) * BDiag(:,2,2) 99 & - BC(:,5) * BDiag(:,3,2) 100 & - BC(:,4) * BDiag(:,2,4) 101 BDiag(:,4,3) = BDiag(:,4,3) + BC(:,4) * BC(:,5) * BDiag(:,2,2) 102 & - BC(:,5) * BDiag(:,2,3) 103 & - BC(:,4) * BDiag(:,4,2) 104 BDiag(:,3,3) = BDiag(:,3,3) + BC(:,4) * BC(:,4) * BDiag(:,2,2) 105 & - BC(:,4) * BDiag(:,2,3) 106 & - BC(:,4) * BDiag(:,3,2) 107 BDiag(:,2,1) = zero 108 BDiag(:,1,2) = zero 109 BDiag(:,2,3) = zero 110 BDiag(:,2,4) = zero 111 BDiag(:,2,5) = zero 112 BDiag(:,3,2) = zero 113 BDiag(:,4,2) = zero 114 BDiag(:,5,2) = zero 115 BDiag(:,2,2) = one 116 endwhere 117c 118c.... x2-velocity 119c 120 where (ibits(iBC,3,3) .eq. 2) 121 BDiag(:,5,4) = BDiag(:,5,4) - BC(:,5) * BDiag(:,5,3) 122 BDiag(:,5,2) = BDiag(:,5,2) - BC(:,4) * BDiag(:,5,3) 123 124 BDiag(:,4,5) = BDiag(:,4,5) - BC(:,5) * BDiag(:,3,5) 125 BDiag(:,2,5) = BDiag(:,2,5) - BC(:,4) * BDiag(:,3,5) 126 127 BDiag(:,4,1) = BDiag(:,4,1) - BC(:,5) * BDiag(:,3,1) 128 BDiag(:,2,1) = BDiag(:,2,1) - BC(:,4) * BDiag(:,3,1) 129 130 BDiag(:,1,4) = BDiag(:,1,4) - BC(:,5) * BDiag(:,1,3) 131 BDiag(:,1,2) = BDiag(:,1,2) - BC(:,4) * BDiag(:,1,3) 132 133 BDiag(:,4,4) = BDiag(:,4,4) + BC(:,5) * BC(:,5) * BDiag(:,3,3) 134 & - BC(:,5) * BDiag(:,3,4) 135 & - BC(:,5) * BDiag(:,4,3) 136 BDiag(:,2,4) = BDiag(:,2,4) + BC(:,4) * BC(:,5) * BDiag(:,3,3) 137 & - BC(:,5) * BDiag(:,2,3) 138 & - BC(:,4) * BDiag(:,3,4) 139 BDiag(:,4,2) = BDiag(:,4,2) + BC(:,4) * BC(:,5) * BDiag(:,3,3) 140 & - BC(:,5) * BDiag(:,3,2) 141 & - BC(:,4) * BDiag(:,4,3) 142 BDiag(:,2,2) = BDiag(:,2,2) + BC(:,4) * BC(:,4) * BDiag(:,3,3) 143 & - BC(:,4) * BDiag(:,2,3) 144 & - BC(:,4) * BDiag(:,3,2) 145 BDiag(:,3,1) = zero 146 BDiag(:,3,2) = zero 147 BDiag(:,3,4) = zero 148 BDiag(:,3,5) = zero 149 BDiag(:,1,3) = zero 150 BDiag(:,2,3) = zero 151 BDiag(:,4,3) = zero 152 BDiag(:,5,3) = zero 153 BDiag(:,3,3) = one 154 endwhere 155c 156c.... x1-velocity and x2-velocity 157c 158 where (ibits(iBC,3,3) .eq. 3) 159 BDiag(:,4,4) = BDiag(:,4,4) + BC(:,4) * BC(:,4) * BDiag(:,2,2) 160 & + BC(:,6) * BC(:,6) * BDiag(:,3,3) 161 & + BC(:,4) * BC(:,6) * ( BDiag(:,2,3) * BDiag(:,3,2)) 162 & - BC(:,6) * ( BDiag(:,4,3) * BDiag(:,3,4)) 163 & - BC(:,4) * ( BDiag(:,4,2) * BDiag(:,2,4)) 164 BDiag(:,1,4) = BDiag(:,1,4) - BC(:,4) * BDiag(:,1,2) 165 & - BC(:,6) * BDiag(:,1,3) 166 BDiag(:,4,1) = BDiag(:,4,1) - BC(:,4) * BDiag(:,2,1) 167 & - BC(:,6) * BDiag(:,3,1) 168 BDiag(:,5,4) = BDiag(:,5,4) - BC(:,4) * BDiag(:,5,2) 169 & - BC(:,6) * BDiag(:,5,3) 170 BDiag(:,4,5) = BDiag(:,4,5) - BC(:,4) * BDiag(:,2,5) 171 & - BC(:,6) * BDiag(:,3,5) 172 BDiag(:,2,1) = zero 173 BDiag(:,2,3) = zero 174 BDiag(:,2,4) = zero 175 BDiag(:,2,5) = zero 176 BDiag(:,3,1) = zero 177 BDiag(:,3,2) = zero 178 BDiag(:,3,4) = zero 179 BDiag(:,3,5) = zero 180 BDiag(:,1,2) = zero 181 BDiag(:,4,2) = zero 182 BDiag(:,5,2) = zero 183 BDiag(:,1,3) = zero 184 BDiag(:,4,3) = zero 185 BDiag(:,5,3) = zero 186 BDiag(:,3,3) = one 187 BDiag(:,2,2) = one 188 endwhere 189c 190c.... x3-velocity 191c 192 where (ibits(iBC,3,3) .eq. 4) 193 BDiag(:,5,3) = BDiag(:,5,3) - BC(:,5) * BDiag(:,5,4) 194 BDiag(:,5,2) = BDiag(:,5,2) - BC(:,4) * BDiag(:,5,4) 195 196 BDiag(:,3,5) = BDiag(:,3,5) - BC(:,5) * BDiag(:,4,5) 197 BDiag(:,2,5) = BDiag(:,2,5) - BC(:,4) * BDiag(:,4,5) 198 199 BDiag(:,3,1) = BDiag(:,3,1) - BC(:,5) * BDiag(:,4,1) 200 BDiag(:,2,1) = BDiag(:,2,1) - BC(:,4) * BDiag(:,4,1) 201 202 BDiag(:,1,3) = BDiag(:,1,3) - BC(:,5) * BDiag(:,1,4) 203 BDiag(:,1,2) = BDiag(:,1,2) - BC(:,4) * BDiag(:,1,4) 204 205 BDiag(:,3,3) = BDiag(:,3,3) + BC(:,5) * BC(:,5) * BDiag(:,4,4) 206 & - BC(:,5) * BDiag(:,3,4) 207 & - BC(:,5) * BDiag(:,4,3) 208 BDiag(:,2,3) = BDiag(:,2,3) + BC(:,4) * BC(:,5) * BDiag(:,4,4) 209 & - BC(:,5) * BDiag(:,2,4) 210 & - BC(:,4) * BDiag(:,4,3) 211 BDiag(:,3,2) = BDiag(:,3,2) + BC(:,4) * BC(:,5) * BDiag(:,4,4) 212 & - BC(:,5) * BDiag(:,4,2) 213 & - BC(:,4) * BDiag(:,3,4) 214 BDiag(:,2,2) = BDiag(:,2,2) + BC(:,4) * BC(:,4) * BDiag(:,4,4) 215 & - BC(:,4) * BDiag(:,2,4) 216 & - BC(:,4) * BDiag(:,4,2) 217 BDiag(:,4,1) = zero 218 BDiag(:,4,2) = zero 219 BDiag(:,4,3) = zero 220 BDiag(:,4,5) = zero 221 BDiag(:,1,4) = zero 222 BDiag(:,2,4) = zero 223 BDiag(:,3,4) = zero 224 BDiag(:,5,4) = zero 225 BDiag(:,4,4) = one 226 endwhere 227c 228c.... x1-velocity and x3-velocity 229c 230 where (ibits(iBC,3,3) .eq. 5) 231 BDiag(:,3,3) = BDiag(:,3,3) + BC(:,4) * BC(:,4) * BDiag(:,2,2) 232 & + BC(:,6) * BC(:,6) * BDiag(:,4,4) 233 & + BC(:,4) * BC(:,6) *(BDiag(:,2,4) + BDiag(:,4,2)) 234 & - BC(:,4) *(BDiag(:,2,3) + BDiag(:,3,2)) 235 & - BC(:,6) *(BDiag(:,4,3) + BDiag(:,3,4)) 236 BDiag(:,1,3) = BDiag(:,1,3) - BC(:,4) * BDiag(:,1,2) 237 & - BC(:,6) * BDiag(:,1,4) 238 BDiag(:,3,1) = BDiag(:,3,1) - BC(:,4) * BDiag(:,2,1) 239 & - BC(:,6) * BDiag(:,4,1) 240 BDiag(:,5,3) = BDiag(:,5,3) - BC(:,4) * BDiag(:,5,2) 241 & - BC(:,6) * BDiag(:,5,4) 242 BDiag(:,3,5) = BDiag(:,3,5) - BC(:,4) * BDiag(:,2,5) 243 & - BC(:,6) * BDiag(:,4,5) 244 BDiag(:,2,1) = zero 245 BDiag(:,2,3) = zero 246 BDiag(:,2,4) = zero 247 BDiag(:,2,5) = zero 248 BDiag(:,4,1) = zero 249 BDiag(:,4,2) = zero 250 BDiag(:,4,3) = zero 251 BDiag(:,4,5) = zero 252 BDiag(:,1,2) = zero 253 BDiag(:,4,2) = zero 254 BDiag(:,5,2) = zero 255 BDiag(:,1,4) = zero 256 BDiag(:,3,4) = zero 257 BDiag(:,5,4) = zero 258 BDiag(:,4,4) = one 259 BDiag(:,2,2) = one 260 endwhere 261c 262c.... x2-velocity and x3-velocity 263c 264 where (ibits(iBC,3,3) .eq. 6) 265 BDiag(:,2,2) = BDiag(:,2,2) + BC(:,4) * BC(:,4) * BDiag(:,3,3) 266 & + BC(:,6) * BC(:,6) * BDiag(:,4,4) 267 & + BC(:,4) * BC(:,6) * (BDiag(:,3,4) + BDiag(:,4,3)) 268 & - BC(:,4) *(BDiag(:,2,3) + BDiag(:,3,2)) 269 & - BC(:,6) *(BDiag(:,4,2) + BDiag(:,2,4)) 270 BDiag(:,1,2) = BDiag(:,1,2) - BC(:,4) * BDiag(:,1,3) 271 & - BC(:,6) * BDiag(:,1,4) 272 BDiag(:,2,1) = BDiag(:,2,1) - BC(:,4) * BDiag(:,3,1) 273 & - BC(:,6) * BDiag(:,4,1) 274 BDiag(:,5,2) = BDiag(:,5,2) - BC(:,4) * BDiag(:,5,3) 275 & - BC(:,6) * BDiag(:,5,4) 276 BDiag(:,2,5) = BDiag(:,2,5) - BC(:,4) * BDiag(:,3,5) 277 & - BC(:,6) * BDiag(:,4,5) 278 BDiag(:,3,1) = zero 279 BDiag(:,3,2) = zero 280 BDiag(:,3,4) = zero 281 BDiag(:,3,5) = zero 282 BDiag(:,4,1) = zero 283 BDiag(:,4,2) = zero 284 BDiag(:,4,3) = zero 285 BDiag(:,4,5) = zero 286 BDiag(:,1,3) = zero 287 BDiag(:,2,3) = zero 288 BDiag(:,5,3) = zero 289 BDiag(:,1,4) = zero 290 BDiag(:,2,4) = zero 291 BDiag(:,5,4) = zero 292 BDiag(:,4,4) = one 293 BDiag(:,3,3) = one 294 endwhere 295c 296c.... x1-velocity and x2-velocity and x3-velocity 297c 298 where (ibits(iBC,3,3) .eq. 7) 299 BDiag(:,2,1) = zero 300 BDiag(:,2,3) = zero 301 BDiag(:,2,4) = zero 302 BDiag(:,2,5) = zero 303 BDiag(:,3,1) = zero 304 BDiag(:,3,2) = zero 305 BDiag(:,3,4) = zero 306 BDiag(:,3,5) = zero 307 BDiag(:,4,1) = zero 308 BDiag(:,4,2) = zero 309 BDiag(:,4,3) = zero 310 BDiag(:,4,5) = zero 311 BDiag(:,1,2) = zero 312 BDiag(:,5,2) = zero 313 BDiag(:,1,3) = zero 314 BDiag(:,5,3) = zero 315 BDiag(:,1,4) = zero 316 BDiag(:,5,4) = zero 317 BDiag(:,2,2) = one 318 BDiag(:,3,3) = one 319 BDiag(:,4,4) = one 320 endwhere 321c 322c.... temperature 323c 324 where (btest(iBC,1)) 325 BDiag(:,5,5) = one 326 BDiag(:,1,5) = zero 327 BDiag(:,2,5) = zero 328 BDiag(:,3,5) = zero 329 BDiag(:,4,5) = zero 330 BDiag(:,5,1) = zero 331 BDiag(:,5,2) = zero 332 BDiag(:,5,3) = zero 333 BDiag(:,5,4) = zero 334 endwhere 335c 336c.... local periodic boundary conditions (no communications) 337c 338 339 do j = 1,nshg 340 if (btest(iBC(j),10)) then 341 i = iper(j) 342 BDiag(i, :,:) = BDiag(i,:,:) + BDiag(j,:,:) 343 endif 344 enddo 345c 346c.... periodic slaves get the residual values of the masters 347c 348 do j = 1,nshg 349 if (btest(iBC(j),10)) then 350 i=iper(j) 351 BDiag(j,:,:) = BDiag(i,:,:) 352 endif 353 enddo 354c$$$ endif 355 356 if(numpe.gt.1) then 357c 358c.... nodes treated on another processor are eliminated 359c 360 numtask = ilwork(1) 361 itkbeg = 1 362 363 do itask = 1, numtask 364 365 iacc = ilwork (itkbeg + 2) 366 numseg = ilwork (itkbeg + 4) 367 368 if (iacc .eq. 0) then 369 do is = 1,numseg 370 isgbeg = ilwork (itkbeg + 3 + 2*is) 371 lenseg = ilwork (itkbeg + 4 + 2*is) 372 isgend = isgbeg + lenseg - 1 373 BDiag(isgbeg:isgend,:,:) = zero 374 BDiag(isgbeg:isgend,1,1) = one 375 BDiag(isgbeg:isgend,2,2) = one 376 BDiag(isgbeg:isgend,3,3) = one 377 BDiag(isgbeg:isgend,4,4) = one 378 BDiag(isgbeg:isgend,5,5) = one 379 enddo 380 endif 381 382 itkbeg = itkbeg + 4 + 2*numseg 383 384 enddo 385 endif 386c 387c.... return 388c 389 return 390 end 391c 392c 393c 394 subroutine bc3BDgSclr (iBC, Diag, iper, ilwork) 395c 396c---------------------------------------------------------------------- 397c 398c This routine satisfies the BC of the block-diagonal preconditioning 399c matrix for 3D elements. 400c 401c input: 402c iBC (numnp) : boundary condition code 403c BC (numnp,ndofBC) : Dirichlet BC constraint parameters 404c Diag (numnp) : preconditionning diagonal matrix before BC 405c 406c output: 407c Diag (numnp) : preconditionning matrix after BC 408c is satisfied 409c 410c 411c Zdenek Johan, Summer 1990. (Modified from g3bce.f) 412c Zdenek Johan, Winter 1991. (Fortran 90) 413c---------------------------------------------------------------------- 414c 415 include "common.h" 416c 417 dimension iBC(nshg), 418 & Diag(nshg), ilwork(nlwork), 419 & iper(nshg) 420c 421c 422 id = isclr+5 423c 424c.... scalar variable 425c 426 where (btest(iBC,id)) 427 Diag(:) = one 428 endwhere 429c 430c.... local periodic boundary conditions (no communications) 431c 432 do j = 1,nshg 433 if (btest(iBC(j),10)) then 434 i = iper(j) 435 Diag(i) = Diag(i) + Diag(j) 436 endif 437 enddo 438c 439c.... periodic slaves get the residual values of the masters 440c 441 do i = 1,nshg 442 if (btest(iBC(i),10)) then 443 Diag(i) = Diag(iper(i)) 444 endif 445 enddo 446c 447c.... nodes treated on another processor are eliminated 448c 449 if(numpe.gt.1) then 450 numtask = ilwork(1) 451 itkbeg = 1 452 453 do itask = 1, numtask 454 455 iacc = ilwork (itkbeg + 2) 456 numseg = ilwork (itkbeg + 4) 457 458 if (iacc .eq. 0) then 459 do is = 1,numseg 460 isgbeg = ilwork (itkbeg + 3 + 2*is) 461 lenseg = ilwork (itkbeg + 4 + 2*is) 462 isgend = isgbeg + lenseg - 1 463 Diag(isgbeg:isgend) = one 464 enddo 465 endif 466 467 itkbeg = itkbeg + 4 + 2*numseg 468 469 enddo 470 endif 471c 472c.... return 473c 474 return 475 end 476 477 478 479