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) 270c 271 dimension rti (npro,4), rLyti (npro) 272c 273 dimension ft1 (npro), 274c unfix later -- pieces used in acusim: 275 & srcrat (npro), vdgn (npro), 276c & term1 (npro), term2 (npro), 277c & term3 (npro), 278c 279 & chi (npro), fv1 (npro), 280 & fv2 (npro), Stilde (npro), 281 & r (npro), g (npro), 282 & fw (npro), ft2 (npro), 283 & fv1p (npro), fv2p (npro), 284 & stp (npro), rp (npro), 285 & gp (npro), fwp (npro), 286 & bf (npro), srcp (npro), 287 & gp6 (npro), tmp (npro), 288 & tmp1 (npro), fwog (npro) 289 real*8 elDwl(npro) ! local quadrature point DES dvar 290 real*8 sclrm(npro) ! modified for non-negativity 291 real*8 saCb1Scale(npro) !Hack to change the production term and BL thickness 292 real*8 xl_xbar(npro) !Hack to store mean x location of element. 293c... for levelset 294 real*8 sign_levelset(npro), sclr_ls(npro), mytmp(npro), 295 & xl(npro,nenl,nsd) 296 297c 298 if(iRANS.lt.0) then ! spalart almaras model 299 sclrm=max(rmu/100.0,Sclr) 300 if(iles.lt.0) then 301 do i=1,npro 302 dx=maxval(xl(i,:,1))-minval(xl(i,:,1)) 303 dy=maxval(xl(i,:,2))-minval(xl(i,:,2)) 304 dz=maxval(xl(i,:,3))-minval(xl(i,:,3)) 305 dmax=max(dx,max(dy,dz)) 306 dmax=0.65d0*dmax 307 if( iles.eq.-1) then !original DES97 308 dist2w(i)=min(dmax,dist2w(i)) 309 elseif(iles.eq.-2) then ! DDES 310 rd=sclrm(i)*saKappaP2Inv/(dist2w(i)**2*gVnrm(i)+1.0d-12) 311 fd=one-tanh((8.0000000000000000d0*rd)**3) 312 dist2w(i)=dist2w(i)-fd*max(zero,dist2w(i)-dmax) 313 endif 314 enddo 315 endif 316 317 elDwl(:)=elDwl(:)+dist2w(:) 318c 319c determine chi 320 chi = rho*sclrm/rmu 321c determine f_v1 322 fv1 = chi**3/(chi**3+saCv1**3) 323c determine f_v2 324 fv2 = one - chi/(one+chi*fv1) 325c determine Stilde 326 Stilde = vort + sclrm*fv2/(saKappa*dist2w)**2!unfix 327c determine r 328 where(Stilde(:).ne.zero) 329 r(:) = sclrm(:)/Stilde(:)/(saKappa*dist2w(:))**2 330 elsewhere 331 r(:) = 1.0d32 332 endwhere 333c determine g 334 saCw3l=saCw3 335 g = r + saCw2*(r**6-r) 336 sixth = 1.0/6.0 337c gp = rp * (tmp + 5 * saCw2 * rP5) 338c 339c gP6 = (g * g * g) ** 2 340c tmp = 1 / (gP6 + (saCw3*saCw3*saCw3)**2) 341c tmp1 = ( (1 + (saCw3*saCw3*saCw3)**2) * tmp ) ** sixth 342c fw = g * tmp1 343c fwp = gp * tmp1 * (saCw3*saCw3*saCw3)**2 * tmp 344c determine f_w and f_w/g 345 fwog = ((one+saCw3**6)/(g**6+saCw3**6))**sixth 346 fw = g*fwog 347c determine f_t2 348c ft2 = ct3*exp(-ct4*chi**2) 349 ft2 = zero 350 351c srcrat=saCb1*(one-ft2)*Stilde*sclrm 352c & -(saCw1*fw - saCb1*ft2/saKappa**2)*(sclrm/dist2w)**2 353c srcrat=srcrat/sclrm 354 355!---------------------------------------------------------------------------- 356!HACK: lower the EV production rate within a region to decrease BL thickness. 357! Appear NM was not finished yet if(scrScaleEnable) then 358 if(one.eq.zero) then 359 do i = 1,nenl !average the x-locations 360 xl_xbar(:) = xl_xbar(:) + xl(:,i,1) 361 enddo 362 xl_xbar = xl_xbar/nenl 363 364 saCb1Scale = one 365 where(xl_xbar < saCb1alterXmin .and. xl_xbar > saCb1alterXmax) 366 saCb1Scale(:) = seCb1alter 367 endwhere 368 369 srcrat = saCb1Scale*saCb1*(one-ft2)*Stilde 370 & -(saCw1*fw - saCb1*ft2/saKappa**2)*(sclrm/dist2w/dist2w) 371 else 372 srcrat=saCb1*(one-ft2)*Stilde 373 & -(saCw1*fw - saCb1*ft2/saKappa**2)*(sclrm/dist2w/dist2w) 374 endif 375 376!Original: 377! srcrat=saCb1*(one-ft2)*Stilde 378! & -(saCw1*fw - saCb1*ft2/saKappa**2)*(sclrm/dist2w/dist2w) 379!End Hack 380!---------------------------------------------------------------------------- 381 382c 383c term1=saCb1*(one-ft2)*Stilde*sclrm 384c term2=saCb2*saSigmaInv*(g1yti**2+g2yti**2+g3yti**2) 385c term3=-(saCw1*fw - saCb1*ft2/saKappa**2)*(sclrm/dist2w)**2 386c determine d()/d(sclrm) 387 fv1p = 3*(saCv1**3)*(chi**2)*rho 388 fv1p = fv1p/(rmu*(chi**3+saCv1**3)**2) 389 fv2p = (chi**2)*fv1p-(one/rmu) 390 fv2p = fv2p/(one+chi*fv1)**2 391 stp = fv2 + sclrm*fv2p 392 stp = stp/(saKappa*dist2w)**2 393 where(Stilde(:).ne.zero) 394 rp(:) = Stilde(:) - sclrm(:)*stp(:) 395 rp(:) = rp(:)/(saKappa*dist2w(:)*Stilde(:))**2 396 elsewhere 397 rp(:) = 1.0d32 398 endwhere 399 gp = one+saCw2*(6*r**5 - one) 400 gp = gp*rp 401 fwp = (saCw3**6)*fwog 402 fwp = fwp*gp/(g**6+saCw3**6) 403c determine source term 404 bf = saCb2*saSigmaInv*(g1yti**2+g2yti**2+g3yti**2) 405 & +saCb1*(one-ft2)*Stilde*sclrm 406 & -(saCw1*fw - saCb1*ft2/saKappa**2)*(sclrm/dist2w)**2 407 bf = bf * rho 408c determine d(source)/d(sclrm) 409 srcp = rho*saCb1*(sclrm*stp+Stilde) 410 & -rho*saCw1*(fwp*sclrm**2 + 2*sclrm*fw)/dist2w**2 411 do i=1, npro 412 if(srcp(i).le.zero .and. srcp(i).le.srcrat(i)) then 413 srcp(i)=srcp(i) 414 else if(srcrat(i).lt.zero) then 415 srcp(i)=srcrat(i) 416 else 417 srcp(i)=zero 418 endif 419 enddo 420c 421c==========================>> IRES = 1 or 3 <<======================= 422c 423c if ((ires .eq. 1) .or. (ires .eq. 3)) then 424 rti (:,4) = rti (:,4) - bf(:) 425c endif !ires 426 427c rmti (:,4) = rmti (:,4) - bf(:) 428 rLyti(:) = rLyti(:) - bf(:) 429c 430 elseif (iLSet.ne.0) then 431 if (isclr.eq.1) then 432 srcp = zero 433 434 elseif (isclr.eq.2) then !we are redistancing level-sets 435 436 sclr_ls = zero !zero out temp variable 437 438 do ii=1,npro 439 440 do jj = 1, nshl ! first find the value of levelset at point ii 441 442 sclr_ls(ii) = sclr_ls(ii) 443 & + shape(ii,jj) * ycl(ii,jj,6) 444 445 enddo 446 447 if (sclr_ls(ii) .lt. - epsilon_ls)then 448 449 sign_levelset(ii) = - one 450 451 elseif (abs(sclr_ls(ii)) .le. epsilon_ls)then 452c sign_levelset(ii) = zero 453c 454 sign_levelset(ii) =sclr_ls(ii)/epsilon_ls 455 & + sin(pi*sclr_ls(ii)/epsilon_ls)/pi 456 457 458 elseif (sclr_ls(ii) .gt. epsilon_ls) then 459 460 sign_levelset(ii) = one 461 462 endif 463 srcp(ii) = sign_levelset(ii) 464 465 enddo 466c 467c ad The redistancing equation can be written in the following form 468c ad 469c ad d_{,t} + sign(phi)*( d_{,i}/|d_{,i}| )* d_{,i} = sign(phi) 470c ad 471c ad This is rewritten in the form 472c ad 473c ad d_{,t} + u * d_{,i} = sign(phi) 474c ad 475 476c$$$ CAD For the redistancing equation the "pseudo velocity" term is 477c$$$ CAD calculated as follows 478 479 480 481 mytmp = srcp(:) / sqrt( g1yti(:) * g1yti(:) 482 & + g2yti(:) * g2yti(:) 483 & + g3yti(:) * g3yti(:) ) 484 485 u1 = mytmp(:) * g1yti(:) 486 u2 = mytmp(:) * g2yti(:) 487 u3 = mytmp(:) * g3yti(:) 488c 489c==========================>> IRES = 1 or 3 <<======================= 490c 491c if ((ires .eq. 1) .or. (ires .eq. 3)) then 492 rti (:,4) = rti (:,4) - srcp(:) 493c endif !ires 494 495c rmti (:,4) = rmti (:,4) - srcp(:) 496 rLyti(:) = rLyti(:) - srcp(:) 497c 498 endif ! close of scalar 2 of level set 499 500 else ! NOT turbulence and NOT level set so this is a simple 501 ! scalar. If your scalar equation has a source term 502 ! then add your own like the above but leave an unforced case 503 ! as an option like you see here 504 505 srcp = zero 506 endif 507 508 509c 510c.... Return and end 511c 512 return 513 end 514 515