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 99 100c we=3.0*29./682. 101 rsteep=3.0 102 src=zero 103 radsts=radSponge*radSponge 104 CoefRatioI2O = grthISponge/grthOSponge 105 do id=1,npro 106 radsqr=xx(id,2)**2+xx(id,1)**2 107 if(xx(id,3).lt. zinSponge) then ! map this into big outflow 108 ! sponge to keep logic 109 ! below simple 110 111 xx(id,3)=(zinSponge-xx(id,3))*CoefRatioI2O 112 & + zoutSponge 113 ! 114 ! CoefRatioI2O is the ratio of the inlet quadratic coefficient to the 115 ! outlet quadratic coeficient (basically how much faster sponge 116 ! coefficient grows in inlet region relative to outlet region) 117 ! 118 endif 119 if((xx(id,3).gt.zoutSponge).or.(radsqr.gt.radsts)) then 120 rad=sqrt(radsqr) 121 radc=max(rad,radSponge) 122 zval=max(xx(id,3),zoutSponge) 123 bcool(id)=grthOSponge*((zval-zoutSponge)**2 124 & +(radc-radSponge)**2) 125 bcool(id)=min(bcool(id),betamax) 126c Determine the resulting density and energies 127 den = ytargeti(id,1) / (Rgas * ytargeti(id,5)) 128 ei = ytargeti(id,5) * ( Rgas / gamma1 ) 129 rk = pt5 * ( ytargeti(id,2)**2+ytargeti(id,3)**2 130 & +ytargeti(id,4)**2 ) 131c Determine the resulting conservation variables 132 duitarg(id,1) = den 133 duitarg(id,2) = den * ytargeti(id,2) 134 duitarg(id,3) = den * ytargeti(id,3) 135 duitarg(id,4) = den * ytargeti(id,4) 136 duitarg(id,5) = den * (ei + rk) 137c Apply the sponge 138 if(spongeContinuity.eq.1) 139 & src(id,1) = -bcool(id)*(dui(id,1) - duitarg(id,1)) 140 if(spongeMomentum1.eq.1) 141 & src(id,2) = -bcool(id)*(dui(id,2) - duitarg(id,2)) 142 if(spongeMomentum2.eq.1) 143 & src(id,3) = -bcool(id)*(dui(id,3) - duitarg(id,3)) 144 if(spongeMomentum3.eq.1) 145 & src(id,4) = -bcool(id)*(dui(id,4) - duitarg(id,4)) 146 if(spongeEnergy.eq.1) 147 & src(id,5) = -bcool(id)*(dui(id,5) - duitarg(id,5)) 148 endif 149 enddo 150 else 151 if(isurf .ne. 1) then 152 write(*,*) 'only vector (1) and cooling (4) implemented' 153 stop 154 endif 155 endif 156 157 if (isurf .eq. 1) then ! add the surface tension force 158 src(:,2) = src(:,2) + rho(:)*sforce(:,1) 159 src(:,3) = src(:,3) + rho(:)*sforce(:,2) 160 src(:,4) = src(:,4) + rho(:)*sforce(:,3) 161 src(:,5) = src(:,5) + (u1*sforce(:,1)+u2*sforce(:,2) 162 & + u3*sforce(:,3))*rho(:) 163 endif 164 165c 166c==========================>> IRES = 1 or 3 <<======================= 167c 168 if (ivart.gt.1) then 169 rLyi(:,1) = rLyi(:,1) - src(:,1) 170 rLyi(:,2) = rLyi(:,2) - src(:,2) 171 rLyi(:,3) = rLyi(:,3) - src(:,3) 172 rLyi(:,4) = rLyi(:,4) - src(:,4) 173 rLyi(:,5) = rLyi(:,5) - src(:,5) 174 endif 175 176 if ((ires .eq. 1) .or. (ires .eq. 3)) then ! we need ri built 177 ri (:,16) = ri (:,16) - src(:,1) 178 ri (:,17) = ri (:,17) - src(:,2) 179 ri (:,18) = ri (:,18) - src(:,3) 180 ri (:,19) = ri (:,19) - src(:,4) 181 ri (:,20) = ri (:,20) - src(:,5) 182 183 endif 184 185 if ((ires.eq.2) .or. (ires.eq.3)) then ! we need rmi built 186 rmi (:,16) = rmi (:,16) - src(:,1) 187 rmi (:,17) = rmi (:,17) - src(:,2) 188 rmi (:,18) = rmi (:,18) - src(:,3) 189 rmi (:,19) = rmi (:,19) - src(:,4) 190 rmi (:,20) = rmi (:,20) - src(:,5) 191 endif 192c 193 return 194 end 195c 196c 197c 198 subroutine e3sourceSclr(Sclr, rho, rmu, 199 & dist2w, vort, gVnrm, con, 200 & g1yti, g2yti, g3yti, 201 & rti, rLyti, srcp, 202 & ycl, shape, u1, 203 & u2, u3, xl, elDwl) 204c 205c--------------------------------------------------------------------- 206c 207c This routine calculates the source term indicated in the Spalart- 208c Allmaras eddy viscosity model. After term is stored in rti(:,4), 209c for later use by e3wmltSclr, and in rLyti(:) for later use by e3lsSclr. 210c 211c input: 212c Sclr (npro) : working turbulence variable 213c rho (npro) : density at intpt 214c rmu (npro) : molecular viscosity 215c dist2w (npro) : distance from intpt to the nearest wall 216c vort (npro) : magnitude of the vorticity 217c gVnrm (npro) : magnitude of the velocity gradient 218c con (npro) : conductivity 219c g1yti (npro) : grad-Sclr in direction 1 220c g2yti (npro) : grad-Sclr in direction 2 221c g3yti (npro) : grad-Sclr in direction 3 222c 223c output: 224c rti (npro,4) : components of residual at intpt 225c rLyti (npro) : GLS stabilization 226c 227c--------------------------------------------------------------------- 228c 229 use turbSA 230 include "common.h" 231c 232 dimension Sclr (npro), ycl(npro,nshl,ndof), 233 & dist2w (npro), shape(npro,nshl), 234 & vort (npro), gVnrm(npro), rho (npro), 235 & rmu (npro), con (npro), 236 & g1yti (npro), g2yti (npro), 237 & g3yti (npro), u1 (npro), 238 & u2 (npro), u3 (npro) 239c 240 dimension rti (npro,4), rLyti (npro) 241c 242 dimension ft1 (npro), 243c unfix later -- pieces used in acusim: 244 & srcrat (npro), vdgn (npro), 245c & term1 (npro), term2 (npro), 246c & term3 (npro), 247c 248 & chi (npro), fv1 (npro), 249 & fv2 (npro), Stilde (npro), 250 & r (npro), g (npro), 251 & fw (npro), ft2 (npro), 252 & fv1p (npro), fv2p (npro), 253 & stp (npro), rp (npro), 254 & gp (npro), fwp (npro), 255 & bf (npro), srcp (npro), 256 & gp6 (npro), tmp (npro), 257 & tmp1 (npro), fwog (npro) 258 real*8 elDwl(npro) ! local quadrature point DES dvar 259 real*8 sclrm(npro) ! modified for non-negativity 260 real*8 saCb1Scale(npro) !Hack to change the production term and BL thickness 261 real*8 xl_xbar(npro) !Hack to store mean x location of element. 262c... for levelset 263 real*8 sign_levelset(npro), sclr_ls(npro), mytmp(npro), 264 & xl(npro,nenl,nsd) 265 266c 267 if(iRANS.lt.0) then ! spalart almaras model 268 sclrm=max(rmu/100.0,Sclr) 269 if(iles.lt.0) then 270 do i=1,npro 271 dx=maxval(xl(i,:,1))-minval(xl(i,:,1)) 272 dy=maxval(xl(i,:,2))-minval(xl(i,:,2)) 273 dz=maxval(xl(i,:,3))-minval(xl(i,:,3)) 274 dmax=max(dx,max(dy,dz)) 275 dmax=0.65d0*dmax 276 if( iles.eq.-1) then !original DES97 277 dist2w(i)=min(dmax,dist2w(i)) 278 elseif(iles.eq.-2) then ! DDES 279 rd=sclrm(i)*saKappaP2Inv/(dist2w(i)**2*gVnrm(i)+1.0d-12) 280 fd=one-tanh((8.0000000000000000d0*rd)**3) 281 dist2w(i)=dist2w(i)-fd*max(zero,dist2w(i)-dmax) 282 endif 283 enddo 284 endif 285 286 elDwl(:)=elDwl(:)+dist2w(:) 287c 288c determine chi 289 chi = rho*sclrm/rmu 290c determine f_v1 291 fv1 = chi**3/(chi**3+saCv1**3) 292c determine f_v2 293 fv2 = one - chi/(one+chi*fv1) 294c determine Stilde 295 Stilde = vort + sclrm*fv2/(saKappa*dist2w)**2!unfix 296c determine r 297 where(Stilde(:).ne.zero) 298 r(:) = sclrm(:)/Stilde(:)/(saKappa*dist2w(:))**2 299 elsewhere 300 r(:) = 1.0d32 301 endwhere 302c determine g 303 saCw3l=saCw3 304 g = r + saCw2*(r**6-r) 305 sixth = 1.0/6.0 306c gp = rp * (tmp + 5 * saCw2 * rP5) 307c 308c gP6 = (g * g * g) ** 2 309c tmp = 1 / (gP6 + (saCw3*saCw3*saCw3)**2) 310c tmp1 = ( (1 + (saCw3*saCw3*saCw3)**2) * tmp ) ** sixth 311c fw = g * tmp1 312c fwp = gp * tmp1 * (saCw3*saCw3*saCw3)**2 * tmp 313c determine f_w and f_w/g 314 fwog = ((one+saCw3**6)/(g**6+saCw3**6))**sixth 315 fw = g*fwog 316c determine f_t2 317c ft2 = ct3*exp(-ct4*chi**2) 318 ft2 = zero 319 320c srcrat=saCb1*(one-ft2)*Stilde*sclrm 321c & -(saCw1*fw - saCb1*ft2/saKappa**2)*(sclrm/dist2w)**2 322c srcrat=srcrat/sclrm 323 324!---------------------------------------------------------------------------- 325!HACK: lower the EV production rate within a region to decrease BL thickness. 326! Appear NM was not finished yet if(scrScaleEnable) then 327 if(one.eq.zero) then 328 do i = 1,nenl !average the x-locations 329 xl_xbar(:) = xl_xbar(:) + xl(:,i,1) 330 enddo 331 xl_xbar = xl_xbar/nenl 332 333 saCb1Scale = one 334 where(xl_xbar < saCb1alterXmin .and. xl_xbar > saCb1alterXmax) 335 saCb1Scale(:) = seCb1alter 336 endwhere 337 338 srcrat = saCb1Scale*saCb1*(one-ft2)*Stilde 339 & -(saCw1*fw - saCb1*ft2/saKappa**2)*(sclrm/dist2w/dist2w) 340 else 341 srcrat=saCb1*(one-ft2)*Stilde 342 & -(saCw1*fw - saCb1*ft2/saKappa**2)*(sclrm/dist2w/dist2w) 343 endif 344 345!Original: 346! srcrat=saCb1*(one-ft2)*Stilde 347! & -(saCw1*fw - saCb1*ft2/saKappa**2)*(sclrm/dist2w/dist2w) 348!End Hack 349!---------------------------------------------------------------------------- 350 351c 352c term1=saCb1*(one-ft2)*Stilde*sclrm 353c term2=saCb2*saSigmaInv*(g1yti**2+g2yti**2+g3yti**2) 354c term3=-(saCw1*fw - saCb1*ft2/saKappa**2)*(sclrm/dist2w)**2 355c determine d()/d(sclrm) 356 fv1p = 3*(saCv1**3)*(chi**2)*rho 357 fv1p = fv1p/(rmu*(chi**3+saCv1**3)**2) 358 fv2p = (chi**2)*fv1p-(one/rmu) 359 fv2p = fv2p/(one+chi*fv1)**2 360 stp = fv2 + sclrm*fv2p 361 stp = stp/(saKappa*dist2w)**2 362 where(Stilde(:).ne.zero) 363 rp(:) = Stilde(:) - sclrm(:)*stp(:) 364 rp(:) = rp(:)/(saKappa*dist2w(:)*Stilde(:))**2 365 elsewhere 366 rp(:) = 1.0d32 367 endwhere 368 gp = one+saCw2*(6*r**5 - one) 369 gp = gp*rp 370 fwp = (saCw3**6)*fwog 371 fwp = fwp*gp/(g**6+saCw3**6) 372c determine source term 373 bf = saCb2*saSigmaInv*(g1yti**2+g2yti**2+g3yti**2) 374 & +saCb1*(one-ft2)*Stilde*sclrm 375 & -(saCw1*fw - saCb1*ft2/saKappa**2)*(sclrm/dist2w)**2 376 bf = bf * rho 377c determine d(source)/d(sclrm) 378 srcp = rho*saCb1*(sclrm*stp+Stilde) 379 & -rho*saCw1*(fwp*sclrm**2 + 2*sclrm*fw)/dist2w**2 380 do i=1, npro 381 if(srcp(i).le.zero .and. srcp(i).le.srcrat(i)) then 382 srcp(i)=srcp(i) 383 else if(srcrat(i).lt.zero) then 384 srcp(i)=srcrat(i) 385 else 386 srcp(i)=zero 387 endif 388 enddo 389c 390c==========================>> IRES = 1 or 3 <<======================= 391c 392c if ((ires .eq. 1) .or. (ires .eq. 3)) then 393 rti (:,4) = rti (:,4) - bf(:) 394c endif !ires 395 396c rmti (:,4) = rmti (:,4) - bf(:) 397 rLyti(:) = rLyti(:) - bf(:) 398c 399 elseif (iLSet.ne.0) then 400 if (isclr.eq.1) then 401 srcp = zero 402 403 elseif (isclr.eq.2) then !we are redistancing level-sets 404 405 sclr_ls = zero !zero out temp variable 406 407 do ii=1,npro 408 409 do jj = 1, nshl ! first find the value of levelset at point ii 410 411 sclr_ls(ii) = sclr_ls(ii) 412 & + shape(ii,jj) * ycl(ii,jj,6) 413 414 enddo 415 416 if (sclr_ls(ii) .lt. - epsilon_ls)then 417 418 sign_levelset(ii) = - one 419 420 elseif (abs(sclr_ls(ii)) .le. epsilon_ls)then 421c sign_levelset(ii) = zero 422c 423 sign_levelset(ii) =sclr_ls(ii)/epsilon_ls 424 & + sin(pi*sclr_ls(ii)/epsilon_ls)/pi 425 426 427 elseif (sclr_ls(ii) .gt. epsilon_ls) then 428 429 sign_levelset(ii) = one 430 431 endif 432 srcp(ii) = sign_levelset(ii) 433 434 enddo 435c 436c ad The redistancing equation can be written in the following form 437c ad 438c ad d_{,t} + sign(phi)*( d_{,i}/|d_{,i}| )* d_{,i} = sign(phi) 439c ad 440c ad This is rewritten in the form 441c ad 442c ad d_{,t} + u * d_{,i} = sign(phi) 443c ad 444 445c$$$ CAD For the redistancing equation the "pseudo velocity" term is 446c$$$ CAD calculated as follows 447 448 449 450 mytmp = srcp(:) / sqrt( g1yti(:) * g1yti(:) 451 & + g2yti(:) * g2yti(:) 452 & + g3yti(:) * g3yti(:) ) 453 454 u1 = mytmp(:) * g1yti(:) 455 u2 = mytmp(:) * g2yti(:) 456 u3 = mytmp(:) * g3yti(:) 457c 458c==========================>> IRES = 1 or 3 <<======================= 459c 460c if ((ires .eq. 1) .or. (ires .eq. 3)) then 461 rti (:,4) = rti (:,4) - srcp(:) 462c endif !ires 463 464c rmti (:,4) = rmti (:,4) - srcp(:) 465 rLyti(:) = rLyti(:) - srcp(:) 466c 467 endif ! close of scalar 2 of level set 468 469 else ! NOT turbulence and NOT level set so this is a simple 470 ! scalar. If your scalar equation has a source term 471 ! then add your own like the above but leave an unforced case 472 ! as an option like you see here 473 474 srcp = zero 475 endif 476 477 478c 479c.... Return and end 480c 481 return 482 end 483 484