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, 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 con (npro) : conductivity 218c g1yti (npro) : grad-Sclr in direction 1 219c g2yti (npro) : grad-Sclr in direction 2 220c g3yti (npro) : grad-Sclr in direction 3 221c 222c output: 223c rti (npro,4) : components of residual at intpt 224c rLyti (npro) : GLS stabilization 225c 226c--------------------------------------------------------------------- 227c 228 use turbSA 229 include "common.h" 230c 231 dimension Sclr (npro), ycl(npro,nshl,ndof), 232 & dist2w (npro), shape(npro,nshl), 233 & vort (npro), rho (npro), 234 & rmu (npro), con (npro), 235 & g1yti (npro), g2yti (npro), 236 & g3yti (npro), u1 (npro), 237 & u2 (npro), u3 (npro) 238c 239 dimension rti (npro,4), rLyti (npro) 240c 241 dimension ft1 (npro), 242c unfix later -- pieces used in acusim: 243 & srcrat (npro), vdgn (npro), 244c & term1 (npro), term2 (npro), 245c & term3 (npro), 246c 247 & chi (npro), fv1 (npro), 248 & fv2 (npro), Stilde (npro), 249 & r (npro), g (npro), 250 & fw (npro), ft2 (npro), 251 & fv1p (npro), fv2p (npro), 252 & stp (npro), rp (npro), 253 & gp (npro), fwp (npro), 254 & bf (npro), srcp (npro), 255 & gp6 (npro), tmp (npro), 256 & tmp1 (npro), fwog (npro) 257 real*8 elDwl(npro) ! local quadrature point DES dvar 258 real*8 sclrm(npro) ! modified for non-negativity 259c... for levelset 260 real*8 sign_levelset(npro), sclr_ls(npro), mytmp(npro), 261 & xl(npro,nenl,nsd) 262 263c 264 if(iRANS.lt.0) then ! spalart almaras model 265 sclrm=max(rmu/100.0,Sclr) 266 if(iles.lt.0) then 267 do i=1,npro 268 dx=maxval(xl(i,:,1))-minval(xl(i,:,1)) 269 dy=maxval(xl(i,:,2))-minval(xl(i,:,2)) 270 dz=maxval(xl(i,:,3))-minval(xl(i,:,3)) 271 dmax=max(dx,max(dy,dz)) 272c.... limit edge length for DES based on SA model 273c.... (only useful when DES_SA_hmin is greater than 0.0 as element lengths are positive) 274 dmax=max(DES_SA_hmin,dmax) 275 dmax=0.65d0*dmax 276 dist2w(i)=min(dmax,dist2w(i)) 277 enddo 278 endif 279 280 elDwl(:)=elDwl(:)+dist2w(:) 281c 282c determine chi 283 chi = rho*sclrm/rmu 284c determine f_v1 285 fv1 = chi**3/(chi**3+saCv1**3) 286c determine f_v2 287 fv2 = one - chi/(one+chi*fv1) 288c determine Stilde 289 Stilde = vort + sclrm*fv2/(saKappa*dist2w)**2!unfix 290c determine r 291 where(Stilde(:).ne.zero) 292 r(:) = sclrm(:)/Stilde(:)/(saKappa*dist2w(:))**2 293 elsewhere 294 r(:) = 1.0d32 295 endwhere 296c determine g 297 saCw3l=saCw3 298 g = r + saCw2*(r**6-r) 299 sixth = 1.0/6.0 300c gp = rp * (tmp + 5 * saCw2 * rP5) 301c 302c gP6 = (g * g * g) ** 2 303c tmp = 1 / (gP6 + (saCw3*saCw3*saCw3)**2) 304c tmp1 = ( (1 + (saCw3*saCw3*saCw3)**2) * tmp ) ** sixth 305c fw = g * tmp1 306c fwp = gp * tmp1 * (saCw3*saCw3*saCw3)**2 * tmp 307c determine f_w and f_w/g 308 fwog = ((one+saCw3**6)/(g**6+saCw3**6))**sixth 309 fw = g*fwog 310c determine f_t2 311c ft2 = ct3*exp(-ct4*chi**2) 312 ft2 = zero 313 314c srcrat=saCb1*(one-ft2)*Stilde*sclrm 315c & -(saCw1*fw - saCb1*ft2/saKappa**2)*(sclrm/dist2w)**2 316c srcrat=srcrat/sclrm 317 318 srcrat=saCb1*(one-ft2)*Stilde 319 & -(saCw1*fw - saCb1*ft2/saKappa**2)*(sclrm/dist2w/dist2w) 320 321c 322c term1=saCb1*(one-ft2)*Stilde*sclrm 323c term2=saCb2*saSigmaInv*(g1yti**2+g2yti**2+g3yti**2) 324c term3=-(saCw1*fw - saCb1*ft2/saKappa**2)*(sclrm/dist2w)**2 325c determine d()/d(sclrm) 326 fv1p = 3*(saCv1**3)*(chi**2)*rho 327 fv1p = fv1p/(rmu*(chi**3+saCv1**3)**2) 328 fv2p = (chi**2)*fv1p-(one/rmu) 329 fv2p = fv2p/(one+chi*fv1)**2 330 stp = fv2 + sclrm*fv2p 331 stp = stp/(saKappa*dist2w)**2 332 where(Stilde(:).ne.zero) 333 rp(:) = Stilde(:) - sclrm(:)*stp(:) 334 rp(:) = rp(:)/(saKappa*dist2w(:)*Stilde(:))**2 335 elsewhere 336 rp(:) = 1.0d32 337 endwhere 338 gp = one+saCw2*(6*r**5 - one) 339 gp = gp*rp 340 fwp = (saCw3**6)*fwog 341 fwp = fwp*gp/(g**6+saCw3**6) 342c determine source term 343 bf = saCb2*saSigmaInv*(g1yti**2+g2yti**2+g3yti**2) 344 & +saCb1*(one-ft2)*Stilde*sclrm 345 & -(saCw1*fw - saCb1*ft2/saKappa**2)*(sclrm/dist2w)**2 346 bf = bf * rho 347c determine d(source)/d(sclrm) 348 srcp = rho*saCb1*(sclrm*stp+Stilde) 349 & -rho*saCw1*(fwp*sclrm**2 + 2*sclrm*fw)/dist2w**2 350 do i=1, npro 351 if(srcp(i).le.zero .and. srcp(i).le.srcrat(i)) then 352 srcp(i)=srcp(i) 353 else if(srcrat(i).lt.zero) then 354 srcp(i)=srcrat(i) 355 else 356 srcp(i)=zero 357 endif 358 enddo 359c 360c==========================>> IRES = 1 or 3 <<======================= 361c 362c if ((ires .eq. 1) .or. (ires .eq. 3)) then 363 rti (:,4) = rti (:,4) - bf(:) 364c endif !ires 365 366c rmti (:,4) = rmti (:,4) - bf(:) 367 rLyti(:) = rLyti(:) - bf(:) 368c 369 elseif (iLSet.ne.0) then 370 if (isclr.eq.1) then 371 srcp = zero 372 373 elseif (isclr.eq.2) then !we are redistancing level-sets 374 375 sclr_ls = zero !zero out temp variable 376 377 do ii=1,npro 378 379 do jj = 1, nshl ! first find the value of levelset at point ii 380 381 sclr_ls(ii) = sclr_ls(ii) 382 & + shape(ii,jj) * ycl(ii,jj,6) 383 384 enddo 385 386 if (sclr_ls(ii) .lt. - epsilon_ls)then 387 388 sign_levelset(ii) = - one 389 390 elseif (abs(sclr_ls(ii)) .le. epsilon_ls)then 391c sign_levelset(ii) = zero 392c 393 sign_levelset(ii) =sclr_ls(ii)/epsilon_ls 394 & + sin(pi*sclr_ls(ii)/epsilon_ls)/pi 395 396 397 elseif (sclr_ls(ii) .gt. epsilon_ls) then 398 399 sign_levelset(ii) = one 400 401 endif 402 srcp(ii) = sign_levelset(ii) 403 404 enddo 405c 406c ad The redistancing equation can be written in the following form 407c ad 408c ad d_{,t} + sign(phi)*( d_{,i}/|d_{,i}| )* d_{,i} = sign(phi) 409c ad 410c ad This is rewritten in the form 411c ad 412c ad d_{,t} + u * d_{,i} = sign(phi) 413c ad 414 415c$$$ CAD For the redistancing equation the "pseudo velocity" term is 416c$$$ CAD calculated as follows 417 418 419 420 mytmp = srcp(:) / sqrt( g1yti(:) * g1yti(:) 421 & + g2yti(:) * g2yti(:) 422 & + g3yti(:) * g3yti(:) ) 423 424 u1 = mytmp(:) * g1yti(:) 425 u2 = mytmp(:) * g2yti(:) 426 u3 = mytmp(:) * g3yti(:) 427c 428c==========================>> IRES = 1 or 3 <<======================= 429c 430c if ((ires .eq. 1) .or. (ires .eq. 3)) then 431 rti (:,4) = rti (:,4) - srcp(:) 432c endif !ires 433 434c rmti (:,4) = rmti (:,4) - srcp(:) 435 rLyti(:) = rLyti(:) - srcp(:) 436c 437 endif ! close of scalar 2 of level set 438 439 else ! NOT turbulence and NOT level set so this is a simple 440 ! scalar. If your scalar equation has a source term 441 ! then add your own like the above but leave an unforced case 442 ! as an option like you see here 443 444 srcp = zero 445 endif 446 447 448c 449c.... Return and end 450c 451 return 452 end 453 454