1 subroutine e3source (ri, rmi, rlyi, 2 & rho, u1, u2, 3 & u3, pres, sforce, 4 & dui, dxidx, ytargetl, 5 & xl, shpfun, bcool) 6c 7c---------------------------------------------------------------------- 8c 9c This routine calculates the contribution of the bodyforce and surface 10c tension force operator to the RHS vector and LHS tangent matrix. The 11c temporary results are put in ri. 12c 13c u1 (npro) : x1-velocity component 14c u2 (npro) : x2-velocity component 15c u3 (npro) : x3-velocity component 16c ri (npro,nflow*(nsd+1)) : partial residual 17c rmi (npro,nflow*(nsd+1)) : partial modified residual 18c rLyi (npro,nflow) : least-squares residual vector 19c shape (npro,nshl) : element shape functions 20c g1yti (npro) : grad-Sclr in direction 1 at intpt 21c g2yti (npro) : grad-Sclr in direction 2 at intpt 22c g3yti (npro) : grad-Sclr in direction 3 at intpt 23c 24 use turbSA 25 use specialBC 26 include "common.h" 27c 28 dimension ri(npro,nflow*(nsd+1)), rmi(npro,nflow*(nsd+1)), 29 & u1(npro), u2(npro), 30 & u3(npro), rho(npro), 31 & pres(npro), 32 & rLyi(npro,nflow), sforce(npro,3), 33 & shpfun(npro,nshl), 34 & xl(npro,nenl,3), xx(npro,3) 35 36 real*8 ytargeti(npro,nflow), ytargetl(npro,nshl,nflow) 37 38 real*8 src(npro,nflow), bcool(npro), 39 & dui(npro,nflow), duitarg(npro,nflow), 40 & dxidx( npro, nsd, nsd), xfind( npro ), delta(npro), rat 41c 42c......contribution of body force 43c 44 bcool=zero 45 src=zero 46c 47 48 49 if(matflg(5,1).eq.1) then ! usual case 50 src(:,1) = zero 51 src(:,2) = rho(:) * datmat(1,5,1) 52 src(:,3) = rho(:) * datmat(2,5,1) 53 src(:,4) = rho(:) * datmat(3,5,1) 54 src(:,5) = u1*src(:,2) + u2*src(:,3) + u3*src(:,4) 55 else if(matflg(5,1).eq.3) then ! user supplied white noise 56 57 xsor = 18 58c ampl = spamp(lstep+1) 59c rat = Delt(1)/0.1 60 ampl = 0.002*exp(-(0.1248222*(lstep)-2.9957323)**2) 61c if((myrank.eq.zero).and.(intp.eq.ngauss)) write(*,*) ampl 62 delta(:) = 0.5*sqrt(dxidx(:,1,1)*dxidx(:,1,1) ! 1/dx 63 . +dxidx(:,2,1)*dxidx(:,2,1) 64 . +dxidx(:,3,1)*dxidx(:,3,1)) 65 do i=1,npro 66 xfind(i) = (xsor-minval(xl(i,:,1))) 67 & *(maxval(xl(i,:,1))-xsor) 68 enddo 69 70 where ( xfind .ge. 0. ) 71 src(:,2) = rho(:) * ampl * delta 72c scaling by element size is removed not to mess up 73c refinement studies 74c src(:,2) = rho(:) * ampl 75 src(:,5) = u1*src(:,2) 76 endwhere 77 78 else if(matflg(5,1).ge.4) then ! cool case (sponge outside of a 79 ! revolved box defined from zinSponge to 80 ! zoutSponge in axial extent and 0 81 ! to radSponge in radial extent for 82 ! all theta) 83 84c determine coordinates of quadrature pt 85 xx=zero 86 do n = 1,nenl 87 xx(:,1) = xx(:,1) + shpfun(:,n) * xl(:,n,1) 88 xx(:,2) = xx(:,2) + shpfun(:,n) * xl(:,n,2) 89 xx(:,3) = xx(:,3) + shpfun(:,n) * xl(:,n,3) 90 enddo 91 ytargeti=zero 92 do j=1,nflow 93 do n=1,nshl 94 ytargeti(:,j) = ytargeti(:,j) 95 & + shpfun(:,n)*ytargetl(:,n,j) 96 enddo 97 enddo 98 if (1.eq.1) then ! bringing in x BL sponge 99 do id=1,npro 100 if((xx(id,1).gt.zoutSponge)) then 101 bcool(id)=grthOSponge*(xx(id,1)-zoutSponge)**2 102 bcool(id)=min(bcool(id),betamax) 103c Determine the resulting density and energies 104 den = ytargeti(id,1) / (Rgas * ytargeti(id,5)) 105 ei = ytargeti(id,5) * ( Rgas / gamma1 ) 106 rk = pt5 * ( ytargeti(id,2)**2+ytargeti(id,3)**2 107 & +ytargeti(id,4)**2 ) 108c Determine the resulting conservation variables 109 duitarg(id,1) = den 110 duitarg(id,2) = den * ytargeti(id,2) 111 duitarg(id,3) = den * ytargeti(id,3) 112 duitarg(id,4) = den * ytargeti(id,4) 113 duitarg(id,5) = den * (ei + rk) 114c Apply the sponge 115 if(spongeContinuity.eq.1) 116 & src(id,1) = -bcool(id)*(dui(id,1) - duitarg(id,1)) 117 if(spongeMomentum1.eq.1) 118 & src(id,2) = -bcool(id)*(dui(id,2) - duitarg(id,2)) 119 if(spongeMomentum2.eq.1) 120 & src(id,3) = -bcool(id)*(dui(id,3) - duitarg(id,3)) 121 if(spongeMomentum3.eq.1) 122 & src(id,4) = -bcool(id)*(dui(id,4) - duitarg(id,4)) 123 if(spongeEnergy.eq.1) 124 & src(id,5) = -bcool(id)*(dui(id,5) - duitarg(id,5)) 125 endif 126 enddo 127 else !keep the original sponge below but 128 129 130c we=3.0*29./682. 131 rsteep=3.0 132 src=zero 133 radsts=radSponge*radSponge 134 CoefRatioI2O = grthISponge/grthOSponge 135 do id=1,npro 136 radsqr=xx(id,2)**2+xx(id,1)**2 137 if(xx(id,3).lt. zinSponge) then ! map this into big outflow 138 ! sponge to keep logic 139 ! below simple 140 141 xx(id,3)=(zinSponge-xx(id,3))*CoefRatioI2O 142 & + zoutSponge 143 ! 144 ! CoefRatioI2O is the ratio of the inlet quadratic coefficient to the 145 ! outlet quadratic coeficient (basically how much faster sponge 146 ! coefficient grows in inlet region relative to outlet region) 147 ! 148 endif 149 if((xx(id,3).gt.zoutSponge).or.(radsqr.gt.radsts)) then 150 rad=sqrt(radsqr) 151 radc=max(rad,radSponge) 152 zval=max(xx(id,3),zoutSponge) 153 bcool(id)=grthOSponge*((zval-zoutSponge)**2 154 & +(radc-radSponge)**2) 155 bcool(id)=min(bcool(id),betamax) 156c Determine the resulting density and energies 157 den = ytargeti(id,1) / (Rgas * ytargeti(id,5)) 158 ei = ytargeti(id,5) * ( Rgas / gamma1 ) 159 rk = pt5 * ( ytargeti(id,2)**2+ytargeti(id,3)**2 160 & +ytargeti(id,4)**2 ) 161c Determine the resulting conservation variables 162 duitarg(id,1) = den 163 duitarg(id,2) = den * ytargeti(id,2) 164 duitarg(id,3) = den * ytargeti(id,3) 165 duitarg(id,4) = den * ytargeti(id,4) 166 duitarg(id,5) = den * (ei + rk) 167c Apply the sponge 168 if(spongeContinuity.eq.1) 169 & src(id,1) = -bcool(id)*(dui(id,1) - duitarg(id,1)) 170 if(spongeMomentum1.eq.1) 171 & src(id,2) = -bcool(id)*(dui(id,2) - duitarg(id,2)) 172 if(spongeMomentum2.eq.1) 173 & src(id,3) = -bcool(id)*(dui(id,3) - duitarg(id,3)) 174 if(spongeMomentum3.eq.1) 175 & src(id,4) = -bcool(id)*(dui(id,4) - duitarg(id,4)) 176 if(spongeEnergy.eq.1) 177 & src(id,5) = -bcool(id)*(dui(id,5) - duitarg(id,5)) 178 endif 179 enddo 180 endif ! end of initial sponge circumvented for BL 181 else 182 if(isurf .ne. 1) then 183 write(*,*) 'only vector (1) and cooling (4) implemented' 184 stop 185 endif 186 endif 187 188 if (isurf .eq. 1) then ! add the surface tension force 189 src(:,2) = src(:,2) + rho(:)*sforce(:,1) 190 src(:,3) = src(:,3) + rho(:)*sforce(:,2) 191 src(:,4) = src(:,4) + rho(:)*sforce(:,3) 192 src(:,5) = src(:,5) + (u1*sforce(:,1)+u2*sforce(:,2) 193 & + u3*sforce(:,3))*rho(:) 194 endif 195 196c 197c==========================>> IRES = 1 or 3 <<======================= 198c 199 if (ivart.gt.1) then 200 rLyi(:,1) = rLyi(:,1) - src(:,1) 201 rLyi(:,2) = rLyi(:,2) - src(:,2) 202 rLyi(:,3) = rLyi(:,3) - src(:,3) 203 rLyi(:,4) = rLyi(:,4) - src(:,4) 204 rLyi(:,5) = rLyi(:,5) - src(:,5) 205 endif 206 207 if ((ires .eq. 1) .or. (ires .eq. 3)) then ! we need ri built 208 ri (:,16) = ri (:,16) - src(:,1) 209 ri (:,17) = ri (:,17) - src(:,2) 210 ri (:,18) = ri (:,18) - src(:,3) 211 ri (:,19) = ri (:,19) - src(:,4) 212 ri (:,20) = ri (:,20) - src(:,5) 213 214 endif 215 216 if ((ires.eq.2) .or. (ires.eq.3)) then ! we need rmi built 217 rmi (:,16) = rmi (:,16) - src(:,1) 218 rmi (:,17) = rmi (:,17) - src(:,2) 219 rmi (:,18) = rmi (:,18) - src(:,3) 220 rmi (:,19) = rmi (:,19) - src(:,4) 221 rmi (:,20) = rmi (:,20) - src(:,5) 222 endif 223c 224 return 225 end 226c 227c 228c 229 subroutine e3sourceSclr(Sclr, rho, rmu, 230 & dist2w, vort, gVnrm, con, 231 & g1yti, g2yti, g3yti, 232 & rti, rLyti, srcp, 233 & ycl, shape, u1, 234 & u2, u3, xl, elDwl) 235c 236c--------------------------------------------------------------------- 237c 238c This routine calculates the source term indicated in the Spalart- 239c Allmaras eddy viscosity model. After term is stored in rti(:,4), 240c for later use by e3wmltSclr, and in rLyti(:) for later use by e3lsSclr. 241c 242c input: 243c Sclr (npro) : working turbulence variable 244c rho (npro) : density at intpt 245c rmu (npro) : molecular viscosity 246c dist2w (npro) : distance from intpt to the nearest wall 247c vort (npro) : magnitude of the vorticity 248c gVnrm (npro) : magnitude of the velocity gradient 249c con (npro) : conductivity 250c g1yti (npro) : grad-Sclr in direction 1 251c g2yti (npro) : grad-Sclr in direction 2 252c g3yti (npro) : grad-Sclr in direction 3 253c 254c output: 255c rti (npro,4) : components of residual at intpt 256c rLyti (npro) : GLS stabilization 257c 258c--------------------------------------------------------------------- 259c 260 use turbSA 261 include "common.h" 262c 263 dimension Sclr (npro), ycl(npro,nshl,ndof), 264 & dist2w (npro), shape(npro,nshl), 265 & vort (npro), gVnrm(npro), rho (npro), 266 & rmu (npro), con (npro), 267 & g1yti (npro), g2yti (npro), 268 & g3yti (npro), u1 (npro), 269 & u2 (npro), u3 (npro), 270 & rnu (npro) 271c 272 dimension rti (npro,4), rLyti (npro) 273c 274 dimension ft1 (npro), 275c unfix later -- pieces used in acusim: 276 & srcrat (npro), vdgn (npro), 277c & term1 (npro), term2 (npro), 278c & term3 (npro), 279c 280 & chi (npro), fv1 (npro), 281 & fv2 (npro), Stilde (npro), 282 & r (npro), g (npro), 283 & fw (npro), ft2 (npro), 284 & fv1p (npro), fv2p (npro), 285 & stp (npro), rp (npro), 286 & gp (npro), fwp (npro), 287 & bf (npro), srcp (npro), 288 & gp6 (npro), tmp (npro), 289 & tmp1 (npro), fwog (npro) 290 real*8 elDwl(npro) ! local quadrature point DES dvar 291 real*8 sclrm(npro) ! modified for non-negativity 292 real*8 saCb1Scale(npro) !Hack to change the production term and BL thickness 293 real*8 xl_xbar(npro) !Hack to store mean x location of element. 294c... for levelset 295 real*8 sign_levelset(npro), sclr_ls(npro), mytmp(npro), 296 & xl(npro,nenl,nsd) 297 298c 299 rnu=rmu/rho ! SA variable is nu_til not mu so nu needed in lots of places where xi=nu_til/nu 300 if(iRANS.lt.0) then ! spalart almaras model 301 sclrm=max(rnu/100.0,Sclr) ! sets a floor on SCLR 302 if(iles.lt.0) then 303 do i=1,npro 304 dx=maxval(xl(i,:,1))-minval(xl(i,:,1)) 305 dy=maxval(xl(i,:,2))-minval(xl(i,:,2)) 306 dz=maxval(xl(i,:,3))-minval(xl(i,:,3)) 307 dmax=max(dx,max(dy,dz)) 308 dmax=0.65d0*dmax 309 if( iles.eq.-1) then !original DES97 310 dist2w(i)=min(dmax,dist2w(i)) 311 elseif(iles.eq.-2) then ! DDES 312 rd=sclrm(i)*saKappaP2Inv/(dist2w(i)**2*gVnrm(i)+1.0d-12) 313 fd=one-tanh((8.0000000000000000d0*rd)**3) 314 dist2w(i)=dist2w(i)-fd*max(zero,dist2w(i)-dmax) 315 endif 316 enddo 317 endif 318 319 elDwl(:)=elDwl(:)+dist2w(:) 320c 321c determine chi 322 chi = sclrm/rnu 323c determine f_v1 324 fv1 = chi**3/(chi**3+saCv1**3) 325c determine f_v2 326 fv2 = one - chi/(one+chi*fv1) 327c determine Stilde 328 Stilde = vort + sclrm*fv2/(saKappa*dist2w)**2!unfix 329c determine r 330 where(Stilde(:).ne.zero) 331 r(:) = sclrm(:)/Stilde(:)/(saKappa*dist2w(:))**2 332 elsewhere 333 r(:) = 1.0d32 334 endwhere 335c determine g 336 saCw3l=saCw3 337 g = r + saCw2*(r**6-r) 338 sixth = 1.0/6.0 339c gp = rp * (tmp + 5 * saCw2 * rP5) 340c 341c gP6 = (g * g * g) ** 2 342c tmp = 1 / (gP6 + (saCw3*saCw3*saCw3)**2) 343c tmp1 = ( (1 + (saCw3*saCw3*saCw3)**2) * tmp ) ** sixth 344c fw = g * tmp1 345c fwp = gp * tmp1 * (saCw3*saCw3*saCw3)**2 * tmp 346c determine f_w and f_w/g 347 fwog = ((one+saCw3**6)/(g**6+saCw3**6))**sixth 348 fw = g*fwog 349c determine f_t2 350c ft2 = ct3*exp(-ct4*chi**2) 351 ft2 = zero 352 353c srcrat=saCb1*(one-ft2)*Stilde*sclrm 354c & -(saCw1*fw - saCb1*ft2/saKappa**2)*(sclrm/dist2w)**2 355c srcrat=srcrat/sclrm 356 357!---------------------------------------------------------------------------- 358!HACK: lower the EV production rate within a region to decrease BL thickness. 359! Appear NM was not finished yet if(scrScaleEnable) then 360 if(one.eq.zero) then 361 do i = 1,nenl !average the x-locations 362 xl_xbar(:) = xl_xbar(:) + xl(:,i,1) 363 enddo 364 xl_xbar = xl_xbar/nenl 365 366 saCb1Scale = one 367 where(xl_xbar < saCb1alterXmin .and. xl_xbar > saCb1alterXmax) 368 saCb1Scale(:) = seCb1alter 369 endwhere 370 371 srcrat = saCb1Scale*saCb1*(one-ft2)*Stilde 372 & -(saCw1*fw - saCb1*ft2/saKappa**2)*(sclrm/dist2w/dist2w) 373 else 374 srcrat=saCb1*(one-ft2)*Stilde 375 & -(saCw1*fw - saCb1*ft2/saKappa**2)*(sclrm/dist2w/dist2w) 376 endif 377 378!Original: 379! srcrat=saCb1*(one-ft2)*Stilde 380! & -(saCw1*fw - saCb1*ft2/saKappa**2)*(sclrm/dist2w/dist2w) 381!End Hack 382!---------------------------------------------------------------------------- 383 384c 385c term1=saCb1*(one-ft2)*Stilde*sclrm 386c term2=saCb2*saSigmaInv*(g1yti**2+g2yti**2+g3yti**2) 387c term3=-(saCw1*fw - saCb1*ft2/saKappa**2)*(sclrm/dist2w)**2 388c determine d()/d(sclrm) 389 fv1p = 3*(saCv1**3)*(chi**2) 390! fv1p = 3*(saCv1**3)*(chi**2) 391! rho stays as chi=nutil/nu = rho nutil/mu -> dxi/dnutil=rho/rmu 392 fv1p = fv1p/(rnu*(chi**3+saCv1**3)**2) 393 fv2p = (chi**2)*fv1p-(one/rnu) 394 fv2p = fv2p/(one+chi*fv1)**2 395 stp = fv2 + sclrm*fv2p 396 stp = stp/(saKappa*dist2w)**2 397 where(Stilde(:).ne.zero) 398 rp(:) = Stilde(:) - sclrm(:)*stp(:) 399 rp(:) = rp(:)/(saKappa*dist2w(:)*Stilde(:))**2 400 elsewhere 401 rp(:) = 1.0d32 402 endwhere 403 gp = one+saCw2*(6*r**5 - one) 404 gp = gp*rp 405 fwp = (saCw3**6)*fwog 406 fwp = fwp*gp/(g**6+saCw3**6) 407c determine source term 408 bf = saCb2*saSigmaInv*(g1yti**2+g2yti**2+g3yti**2) 409 & +saCb1*(one-ft2)*Stilde*sclrm 410 & -(saCw1*fw - saCb1*ft2/saKappa**2)*(sclrm/dist2w)**2 411c determine d(source)/d(sclrm) 412 srcp = saCb1*(sclrm*stp+Stilde) 413 & -saCw1*(fwp*sclrm**2 + 2*sclrm*fw)/dist2w**2 414 do i=1, npro 415 if(srcp(i).le.zero .and. srcp(i).le.srcrat(i)) then 416 srcp(i)=srcp(i) 417 else if(srcrat(i).lt.zero) then 418 srcp(i)=srcrat(i) 419 else 420 srcp(i)=zero 421 endif 422 enddo 423c 424c==========================>> IRES = 1 or 3 <<======================= 425c 426c if ((ires .eq. 1) .or. (ires .eq. 3)) then 427 rti (:,4) = rti (:,4) - bf(:) 428c endif !ires 429 430c rmti (:,4) = rmti (:,4) - bf(:) 431 rLyti(:) = rLyti(:) - bf(:) 432c 433 elseif (iLSet.ne.0) then 434 if (isclr.eq.1) then 435 srcp = zero 436 437 elseif (isclr.eq.2) then !we are redistancing level-sets 438 439 sclr_ls = zero !zero out temp variable 440 441 do ii=1,npro 442 443 do jj = 1, nshl ! first find the value of levelset at point ii 444 445 sclr_ls(ii) = sclr_ls(ii) 446 & + shape(ii,jj) * ycl(ii,jj,6) 447 448 enddo 449 450 if (sclr_ls(ii) .lt. - epsilon_ls)then 451 452 sign_levelset(ii) = - one 453 454 elseif (abs(sclr_ls(ii)) .le. epsilon_ls)then 455c sign_levelset(ii) = zero 456c 457 sign_levelset(ii) =sclr_ls(ii)/epsilon_ls 458 & + sin(pi*sclr_ls(ii)/epsilon_ls)/pi 459 460 461 elseif (sclr_ls(ii) .gt. epsilon_ls) then 462 463 sign_levelset(ii) = one 464 465 endif 466 srcp(ii) = sign_levelset(ii) 467 468 enddo 469c 470c ad The redistancing equation can be written in the following form 471c ad 472c ad d_{,t} + sign(phi)*( d_{,i}/|d_{,i}| )* d_{,i} = sign(phi) 473c ad 474c ad This is rewritten in the form 475c ad 476c ad d_{,t} + u * d_{,i} = sign(phi) 477c ad 478 479c$$$ CAD For the redistancing equation the "pseudo velocity" term is 480c$$$ CAD calculated as follows 481 482 483 484 mytmp = srcp(:) / sqrt( g1yti(:) * g1yti(:) 485 & + g2yti(:) * g2yti(:) 486 & + g3yti(:) * g3yti(:) ) 487 488 u1 = mytmp(:) * g1yti(:) 489 u2 = mytmp(:) * g2yti(:) 490 u3 = mytmp(:) * g3yti(:) 491c 492c==========================>> IRES = 1 or 3 <<======================= 493c 494c if ((ires .eq. 1) .or. (ires .eq. 3)) then 495 rti (:,4) = rti (:,4) - srcp(:) 496c endif !ires 497 498c rmti (:,4) = rmti (:,4) - srcp(:) 499 rLyti(:) = rLyti(:) - srcp(:) 500c 501 endif ! close of scalar 2 of level set 502 503 else ! NOT turbulence and NOT level set so this is a simple 504 ! scalar. If your scalar equation has a source term 505 ! then add your own like the above but leave an unforced case 506 ! as an option like you see here 507 508 srcp = zero 509 endif 510 511 512c 513c.... Return and end 514c 515 return 516 end 517 518