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