1 subroutine e3b (yl, ycl, iBCB, BCB, shpb, shglb, 2 & xlb, rl, rml, sgn) 3c 4c---------------------------------------------------------------------- 5c 6c This routine calculates the 3D RHS residual of the fluid boundary 7c elements. 8c 9c input: 10c yl (npro,nshl,nflow) : Y variables (perturbed, no scalars) 11c ycl (npro,nshl,ndof) : Y variables 12c iBCB (npro,ndiBCB) : boundary condition code (iBCB(:,1) is 13c a bit tested boundary integral flag i.e. 14c if set to value of BCB if set to floating value 15c iBCB(:,1) : convective flux * 1 0 (ditto to all below) 16c pressure flux * 2 17c viscous flux * 4 18c heat flux * 8 19c turbulence wall * 16 20c scalarI flux * 16*2^I 21c (where I is the scalar number) 22c 23c iBCB(:,2) is the srfID given by the user in MGI that we will 24c collect integrated fluxes for. 25c 26c BCB (npro,nshlb,ndBCB) : Boundary Condition values 27c BCB (1) : mass flux 28c BCB (2) : pressure 29c BCB (3) : viscous flux in x1-direc. 30c BCB (4) : viscous flux in x2-direc. 31c BCB (5) : viscous flux in x3-direc. 32c BCB (6) : heat flux 33c shpb (nshl,ngaussb) : boundary element shape-functions 34c shglb (nsd,nshl,ngaussb) : boundary element grad-shape-functions 35c xlb (npro,nenl,nsd) : nodal coordinates at current step 36c 37c output: 38c rl (npro,nshl,nflow) : element residual 39c rml (npro,nshl,nflow) : element modified residual 40c 41c 42c Note: Always the first side of the element is on the boundary. 43c However, note that for higher-order elements the nodes on 44c the boundary side are not the first nshlb nodes, see the 45c array lnode. 46c 47c 48c Zdenek Johan, Summer 1990. (Modified from e2b.f) 49c Zdenek Johan, Winter 1991. (Fortran 90) 50c Anilkumar Karanam Spring 2000 (Modified for Hierarchic Hexes) 51c---------------------------------------------------------------------- 52c 53 include "common.h" 54c 55 dimension yl(npro,nshl,nflow), iBCB(npro,ndiBCB), 56 & ycl(npro,nshl,ndof), 57 & BCB(npro,nshlb,ndBCB), shpb(nshl,ngaussb), 58 & shglb(nsd,nshl,ngaussb), 59 & xlb(npro,nenl,nsd), 60 & rl(npro,nshl,nflow), rml(npro,nshl,nflow) 61c 62 dimension g1yi(npro,nflow), g2yi(npro,nflow), 63 & g3yi(npro,nflow), WdetJb(npro), 64 & bnorm(npro,nsd) 65c 66 dimension un(npro), rk(npro), 67 & u1(npro), u2(npro), 68 & u3(npro), 69 & rho(npro), pres(npro), 70 & T(npro), ei(npro), 71 & cp(npro) 72c 73 dimension rou(npro), p(npro), 74 & F1(npro), F2(npro), 75 & F3(npro), F4(npro), 76 & F5(npro), Fv2(npro), 77 & Fv3(npro), Fv4(npro), 78 & Fv5(npro), Fh5(npro) 79c 80 dimension rmu(npro), rlm(npro), 81 & rlm2mu(npro), con(npro), 82 & tau1n(npro), 83 & tau2n(npro), tau3n(npro), 84 & heat(npro) 85c 86 dimension lnode(27), sgn(npro,nshl), 87 & shape(npro,nshl), shdrv(npro,nsd,nshl) 88c 89 dimension xmudum(npro,ngauss) 90 91 ttim(40) = ttim(40) - secs(0.0) 92 93c 94c.... compute the nodes which lie on the boundary 95c 96 call getbnodes(lnode) 97 98c.... loop through the integration points 99 100 if(lcsyst.eq.3.or.lcsyst.eq.4) then 101 ngaussb = nintb(lcsyst) 102 else 103 ngaussb = nintb(lcsyst) 104 endif 105 106 do intp = 1, ngaussb 107c 108c.... if Det. .eq. 0, do not include this point 109c 110 if (Qwtb(lcsyst,intp) .eq. zero) cycle ! precaution 111c 112c.... create a matrix of shape functions (and derivatives) for each 113c element at this quadrature point. These arrays will contain 114c the correct signs for the hierarchic basis 115c 116c 117 call getshpb(shpb, shglb, sgn, 118 & shape, shdrv) 119c 120c.... calculate the integration variables 121c 122 call e3bvar (yl, ycl, BCB, shape, 123 & shdrv, xlb, 124 & lnode, g1yi, 125 & g2yi, g3yi, WdetJb, 126 & bnorm, pres, T, 127 & u1, u2, u3, 128 & rho, ei, cp, 129 & rk, rou, p, 130 & Fv2, Fv3, Fv4, 131 & Fh5) 132c 133c.... ires = 1 or 3 134c 135 if ((ires .eq. 1) .or. (ires .eq. 3)) then 136c 137c.... clear some variables 138c 139 tau1n = zero 140 tau2n = zero 141 tau3n = zero 142 heat = zero 143c 144c.... -------------------------> convective <------------------------ 145c 146c 147 where (.not.btest(iBCB(:,1),0) ) 148 un = bnorm(:,1) * u1 + bnorm(:,2) * u2 + bnorm(:,3) * u3 149 rou = rho * ( un ) 150 elsewhere 151 un = (rou / rho) 152 endwhere 153c 154c.... -------------------------> pressure <-------------------------- 155c 156c.... use one-point quadrature in time 157c 158 where (.not.btest(iBCB(:,1),1)) p = pres 159c 160c.... add the Euler contribution 161c 162 F1 = rou 163 F2 = rou * u1 + bnorm(:,1) * p 164 F3 = rou * u2 + bnorm(:,2) * p 165 F4 = rou * u3 + bnorm(:,3) * p 166 F5 = rou * (ei + rk) + un * p 167c 168c.... flop count 169c 170 flops = flops + 23*npro 171c 172c.... end of ires = 1 or 3 173c 174 endif 175c 176c.... -----------------------> Navier-Stokes <----------------------- 177c 178 if (Navier .eq. 1) then 179 180 xmudum = zero 181 182c 183c.... get the material properties 184c 185 call getDiff (T, cp, rho, ycl, 186 & rmu, rlm, rlm2mu, con, shape, 187 & xmudum, xl) 188c 189c.... ------------------------> viscous flux <------------------------ 190c 191c.... floating viscous flux 192c 193 tau1n = bnorm(:,1) * (rlm2mu* g1yi(:,2) + rlm *g2yi(:,3) 194 & + rlm *g3yi(:,4)) 195 & + bnorm(:,2) * (rmu *(g2yi(:,2) + g1yi(:,3))) 196 & + bnorm(:,3) * (rmu *(g3yi(:,2) + g1yi(:,4))) 197 tau2n = bnorm(:,1) * (rmu *(g2yi(:,2) + g1yi(:,3))) 198 & + bnorm(:,2) * (rlm * g1yi(:,2) + rlm2mu*g2yi(:,3) 199 & + rlm *g3yi(:,4)) 200 & + bnorm(:,3) * (rmu *(g3yi(:,3) + g2yi(:,4))) 201 tau3n = bnorm(:,1) * (rmu *(g3yi(:,2) + g1yi(:,4))) 202 & + bnorm(:,2) * (rmu *(g3yi(:,3) + g2yi(:,4))) 203 & + bnorm(:,3) * (rlm * g1yi(:,2) + rlm *g2yi(:,3) 204 & + rlm2mu*g3yi(:,4)) 205c 206 where (.not.btest(iBCB(:,1),2) ) 207 Fv2 = tau1n ! wherever traction is not set, use the 208 Fv3 = tau2n ! viscous flux calculated from the field 209 Fv4 = tau3n ! 210 endwhere 211c 212 Fv5 = u1 * Fv2 + u2 * Fv3 + u3 * Fv4 213c 214c.... --------------------------> heat flux <------------------------- 215c 216c.... floating heat flux 217c 218 heat = con * ( bnorm(:,1) * g1yi(:,5) + 219 & bnorm(:,2) * g2yi(:,5) + 220 & bnorm(:,3) * g3yi(:,5) ) 221c 222 where (.not.btest(iBCB(:,1),3) ) Fh5 = heat 223c 224c.... add the Navier-Stokes contribution 225c 226 F2 = F2 - Fv2 227 F3 = F3 - Fv3 228 F4 = F4 - Fv4 229 F5 = F5 - Fv5 + Fh5 230c 231c.... flop count 232c 233 flops = flops + 27*npro 234c 235c.... end of Navier Stokes part 236c 237 endif 238c 239c.... -------------------------> Residual <-------------------------- 240c 241c.... add the flux to the residual 242c 243 if ((ires .eq. 1) .or. (ires .eq. 3)) then 244c 245c 246 do n = 1, nshlb 247 nodlcl = lnode(n) 248c 249 rl(:,nodlcl,1) = rl(:,nodlcl,1) 250 & + WdetJb * shape(:,nodlcl) * F1 251 rl(:,nodlcl,2) = rl(:,nodlcl,2) 252 & + WdetJb * shape(:,nodlcl) * F2 253 rl(:,nodlcl,3) = rl(:,nodlcl,3) 254 & + WdetJb * shape(:,nodlcl) * F3 255 rl(:,nodlcl,4) = rl(:,nodlcl,4) 256 & + WdetJb * shape(:,nodlcl) * F4 257 rl(:,nodlcl,5) = rl(:,nodlcl,5) 258 & + WdetJb * shape(:,nodlcl) * F5 259 enddo 260c 261 flops = flops + 12*nshlb*npro 262c 263 endif 264c 265c.... add the flux to the modified residual 266c 267 if (((ires .eq. 2) .or. (ires .eq. 3)) 268 & .and. (Navier .eq. 1) .and. (Jactyp .eq. 1)) then 269c 270 do n = 1, nshlb 271 nodlcl = lnode(n) 272c 273 rml(:,nodlcl,2) = rml(:,nodlcl,2) - WdetJb * 274 & shape(:,nodlcl) * Fv2 275 rml(:,nodlcl,3) = rml(:,nodlcl,3) - WdetJb * 276 & shape(:,nodlcl) * Fv3 277 rml(:,nodlcl,4) = rml(:,nodlcl,4) - WdetJb * 278 & shape(:,nodlcl) * Fv4 279 rml(:,nodlcl,5) = rml(:,nodlcl,5) - WdetJb * 280 & shape(:,nodlcl) * (Fv5 - Fh5) 281 enddo 282c 283 flops = flops + 11*nenbl*npro 284c 285 endif 286c 287c uncomment and run 1 step to get estimate of wall shear vs z 288c 289c do i=1,npro 290c tnorm= sqrt(tau1n(i)*tau1n(i) 291c & +tau2n(i)*tau2n(i)+tau3n(i)*tau3n(i)) 292c 293c write(700+myrank,*) xlb(i,1,3),tnorm 294c enddo 295 296 297 do iel = 1, npro 298c 299c if we have a nonzero value then 300c calculate the fluxes through this surface 301c 302 iface = abs(iBCB(iel,2)) 303 if (iface .ne. 0 .and. ires.ne.2) then 304 flxID(1,iface) = flxID(1,iface) + WdetJb(iel)! measure area too 305c flxID(2,iface) = flxID(2,iface) - WdetJb(iel) * un(iel) 306 flxID(2,iface) = flxID(2,iface) - WdetJb(iel) * rou(iel) 307 flxID(3,iface) = flxID(3,iface) 308 & - ( tau1n(iel) - bnorm(iel,1)*pres(iel)) 309 & * WdetJb(iel) 310 flxID(4,iface) = flxID(4,iface) 311 & - ( tau2n(iel) - bnorm(iel,2)*pres(iel)) 312 & * WdetJb(iel) 313 flxID(5,iface) = flxID(5,iface) 314 & - ( tau3n(iel) - bnorm(iel,3)*pres(iel)) 315 & * WdetJb(iel) 316 317 endif 318 enddo 319 320c 321c.... --------------------> Aerodynamic Forces <--------------------- 322c 323 if ((ires .ne. 2) .and. (iter .eq. nitr)) then 324c 325c.... compute the forces on the body 326c 327 where (.not.btest(iBCB(:,1),0) ) 328 tau1n = ( pres * bnorm(:,1) - tau1n ) * WdetJb 329 tau2n = ( pres * bnorm(:,2) - tau2n ) * WdetJb 330 tau3n = ( pres * bnorm(:,3) - tau3n ) * WdetJb 331 heat = - heat * WdetJb 332 elsewhere 333 tau1n = zero 334 tau2n = zero 335 tau3n = zero 336 heat = zero 337 endwhere 338c 339 Force(1) = Force(1) + sum(tau1n) 340 Force(2) = Force(2) + sum(tau2n) 341 Force(3) = Force(3) + sum(tau3n) 342 HFlux = HFlux + sum(heat) 343c 344 endif 345c 346c.... end of integration loop 347c 348 enddo 349 350 ttim(40) = ttim(40) + secs(0.0) 351c 352c.... return 353c 354 return 355 end 356c 357c 358c 359 subroutine e3bSclr (ycl, iBCB, BCB, 360 & shpb, shglb, sgn, 361 & xlb, rtl, rmtl) 362c 363c---------------------------------------------------------------------- 364c 365c This routine calculates the 3D RHS residual of the fluid boundary 366c elements. 367c 368c input: 369c ycl (npro,nshl,ndof) : Y variables 370c iBCB (npro,ndiBCB) : boundary condition code & surfID 371c BCB (npro,nenbl,ndBCB) : Boundary Condition values 372c BCB (1) : mass flux 373c BCB (2) : pressure 374c BCB (3) : viscous flux in x1-direc. 375c BCB (4) : viscous flux in x2-direc. 376c BCB (5) : viscous flux in x3-direc. 377c BCB (6) : heat flux 378c BCB (7) : eddy visc flux 379c shpb (nen,nintg) : boundary element shape-functions 380c shglb (nsd,nen,nintg) : boundary element grad-shape-functions 381c xlb (npro,nenl,nsd) : nodal coordinates at current step 382c 383c output: 384c rtl (npro,nenl,nflow) : element residual 385c rmtl (npro,nenl,nflow) : element modified residual 386c 387c 388c Note: Always the first side of the element is on the boundary. 389c However, note that for higher-order elements the nodes on 390c the boundary side are not the first nenbl nodes, see the 391c array mnodeb. 392c 393c 394c Zdenek Johan, Summer 1990. (Modified from e2b.f) 395c Zdenek Johan, Winter 1991. (Fortran 90) 396c---------------------------------------------------------------------- 397c 398 use turbSA ! for saSigma 399 include "common.h" 400c 401 dimension ycl(npro,nshl,ndof), iBCB(npro,ndiBCB), 402 & BCB(npro,nshlb,ndBCB), shpb(nshl,ngaussb), 403 & shglb(nsd,nshl,ngaussb), sgn(npro,nshl), 404 & xlb(npro,nenl,nsd), shape(npro,nshl), 405 & rtl(npro,nshl), rmtl(npro,nshl), 406 & shdrv(npro,nsd,nshl) 407c 408 dimension u1(npro), u2(npro), 409 & u3(npro), 410 & g1yti(npro), g2yti(npro), 411 & g3yti(npro), WdetJb(npro), 412 & bnorm(npro,nsd) 413c 414 dimension rho(npro), Sclr(npro), 415 & T(npro), cp(npro) 416c 417 dimension rou(npro), F(npro), 418 & un(npro), Sclrn(npro) 419c 420 dimension rmu(npro), rlm(npro), 421 & rlm2mu(npro), con(npro), 422 & heat(npro), srcp(npro) 423c 424 dimension lnode(27) 425 real*8 sign_levelset(npro), sclr_ls(npro), mytmp(npro) 426 ttim(40) = ttim(40) - tmr() 427c 428c.... get the boundary nodes 429c 430 id = isclr + 5 431 ib = isclr + 4 432 ibb = isclr + 6 433 call getbnodes(lnode) 434c 435c.... loop through the integration points 436c 437 ngaussb = nintb(lcsyst) 438c 439 do intp = 1, ngaussb 440c 441c.... if Det. .eq. 0, do not include this point 442c 443 if (Qwtb(lcsyst,intp) .eq. zero) cycle ! precaution 444c 445c.... create a matrix of shape functions (and derivatives) for each 446c element at this quadrature point. These arrays will contain 447c the correct signs for the hierarchic basis 448c 449 call getshpb(shpb, shglb, sgn, 450 & shape, shdrv) 451c 452c.... calculate the integraton variables 453c 454 call e3bvarSclr (ycl, BCB, 455 & shape, shdrv, 456 & xlb, lnode, u1, 457 & u2, u3, g1yti, 458 & g2yti, g3yti, WdetJb, 459 & bnorm, T, rho, 460 & cp, rou, Sclr, 461 & F) 462c.......********************modification for Ilset=2********************** 463 if (ilset.eq.2 .and. isclr.eq.2) then !we are redistancing level-sets 464 465CAD If Sclr(:,1).gt.zero, result of sign_term function 1 466CAD If Sclr(:,1).eq.zero, result of sign_term function 0 467CAD If Sclr(:,1).lt.zero, result of sign_term function -1 468 469 sclr_ls = zero !zero out temp variable 470 471 do ii=1,npro 472 473 do jj = 1, nshl ! first find the value of levelset at point ii 474 475 sclr_ls(ii) = sclr_ls(ii) 476 & + shape(ii,jj) * ycl(ii,jj,6) 477 478 enddo 479 if (sclr_ls(ii) .lt. - epsilon_ls)then 480 481 sign_levelset(ii) = - one 482 483 elseif (abs(sclr_ls(ii)) .le. epsilon_ls)then 484c 485 sign_levelset(ii) =sclr_ls(ii)/epsilon_ls 486 & + sin(pi*sclr_ls(ii)/epsilon_ls)/pi 487 488 elseif (sclr_ls(ii) .gt. epsilon_ls) then 489 490 sign_levelset(ii) = one 491 492 endif 493c 494 srcp(ii) = sign_levelset(ii) 495c 496 enddo 497c 498cad The redistancing equation can be written in the following form 499cad 500cad d_{,t} + sign(phi)*( d_{,i}/|d_{,i}| )* d_{,i} = sign(phi) 501cad 502cad This is rewritten in the form 503cad 504cad d_{,t} + u * d_{,i} = sign(phi) 505cad 506 507c$$$CAD For the redistancing equation the "pseudo velocity" term is 508c$$$CAD calculated as follows 509 510 511 512 mytmp = srcp / sqrt ( g1yti * g1yti 513 & + g2yti * g2yti 514 & + g3yti * g3yti) 515 516 u1 = mytmp * g1yti 517 u2 = mytmp * g2yti 518 u3 = mytmp * g3yti 519 endif 520 521c 522c.... ires = 1 or 3 523c 524 if ((ires .eq. 1) .or. (ires .eq. 3)) then 525 526c.... -------------------------> convective <------------------------ 527c 528 where (.not.btest(iBCB(:,1),0) ) 529 un = bnorm(:,1)*u1 + bnorm(:,2)*u2 + bnorm(:,3)*u3 530 rou = rho * ( un ) 531 elsewhere 532 un = (rou / rho) 533 endwhere 534c 535c.... calculate flux where unconstrained 536c 537 where (.not.btest(iBCB(:,1),ib) ) 538 F = Sclr *rou 539 endwhere 540c 541c.... get the material properties 542c 543 544 call getDiffSclr (T, cp, rmu, 545 & rlm, rlm2mu, con, rho, Sclr) 546 547c 548c.... ----------> DiffFlux for Scalar Variable <-------- 549c 550 if (ilset.ne.2) then 551 552 where (.not.btest(iBCB(:,1),ib) ) 553 Sclrn = bnorm(:,1) * g1yti(:) 554 & + bnorm(:,2) * g2yti(:) 555 & + bnorm(:,3) * g3yti(:) 556C 557 558c F = F + rmu*Sclrn !!!! CHECK THIS 559 560 F = F + saSigmaInv*rho*((rmu/rho)+Sclr)*Sclrn!confirm the modificationc in getdiffsclr 561 562c.....this modification of viscosity goes in getdiffsclr 563 564 endwhere 565 endif 566c 567c.... end of ires = 1 or 3 568c 569 endif 570c 571c.... -------------------------> Residual <-------------------------- 572c 573c.... add the flux to the residual 574c 575 if ((ires .eq. 1) .or. (ires .eq. 3)) then 576 if (iconvsclr.eq.1) then !conservative boundary integral 577 do n = 1, nshlb 578 nodlcl = lnode(n) 579 rtl(:,nodlcl) = rtl(:,nodlcl) 580 & + WdetJb * shape(:,nodlcl) * F 581 enddo 582 flops = flops + 12*nshlb*npro 583 endif 584 endif 585c 586c.... add the flux to the modified residual 587c 588c if (((ires .eq. 2) .or. (ires .eq. 3)) 589c & .and. (Navier .eq. 1) .and. (Jactyp .eq. 1)) then 590c 591c do n = 1, nenbl 592c nodlcl = lnode(n) 593c 594c rml(:,nodlcl,2) = rml(:,nodlcl,2) - WdetJb * 595c & shpb(nodlcl,intp) * Fv2 596c rml(:,nodlcl,3) = rml(:,nodlcl,3) - WdetJb * 597c & shpb(nodlcl,intp) * Fv3 598c rml(:,nodlcl,4) = rml(:,nodlcl,4) - WdetJb * 599c & shpb(nodlcl,intp) * Fv4 600c rml(:,nodlcl,5) = rml(:,nodlcl,5) - WdetJb * 601c & shpb(nodlcl,intp) * (Fv5 - Fh5) 602c enddo 603c 604c flops = flops + 11*nenbl*npro 605c 606c endif 607c 608 609c 610c.... end of integration loop 611c 612 enddo 613 614 ttim(40) = ttim(40) + tmr() 615c 616c.... return 617c 618 return 619 end 620 621