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