1 subroutine e3ivar (yl, acl, shpfun, 2 & shgl, xl, 3 & aci, g1yi, g2yi, 4 & g3yi, shg, dxidx, 5 & WdetJ, rho, pres, 6 & u1, u2, u3, 7 & ql, rLui, src, 8 & rerrl, rlsl, rlsli, 9 & dwl) 10c 11c---------------------------------------------------------------------- 12c 13c This routine computes the variables at integration point. 14c 15c input: 16c yl (npro,nshl,ndof) : primitive variables 17c acl (npro,nshl,ndof) : prim.var. accel. 18c shp (nen) : element shape-functions 19c shgl (nsd,nen) : element local-grad-shape-functions 20c xl (npro,nenl,nsd) : nodal coordinates at current step 21c ql (npro,nshl,nsd*nsd) : diffusive flux vector 22c rlsl (npro,nshl,6) : resolved Leonard stresses 23c 24c output: 25c aci (npro,3) : primvar accel. variables 26c g1yi (npro,ndof) : grad-y in direction 1 27c g2yi (npro,ndof) : grad-y in direction 2 28c g3yi (npro,ndof) : grad-y in direction 3 29c shg (npro,nshl,nsd) : element global grad-shape-functions 30c dxidx (npro,nsd,nsd) : inverse of deformation gradient 31c WdetJ (npro) : weighted Jacobian 32c rho (npro) : density 33c pres (npro) : pressure 34c u1 (npro) : x1-velocity component 35c u2 (npro) : x2-velocity component 36c u3 (npro) : x3-velocity component 37c rLui (npro,nsd) : xi-momentum residual 38c src (npro,nsd) : body force term (not density weighted) 39c rlsli (npro,6) : resolved Leonard stresses at quad pt 40c 41c locally calculated and used 42c divqi (npro,nsd+isurf) : divergence of reconstructed quantity 43c 44c Zdenek Johan, Summer 1990. (Modified from e2ivar.f) 45c Zdenek Johan, Winter 1991. (Fortran 90) 46c Kenneth Jansen, Winter 1997. Primitive Variables 47c Christian Whiting, Winter 1999. (uBar formulation) 48c 49c---------------------------------------------------------------------- 50c 51 include "common.h" 52c 53c passed arrays 54c 55 dimension yl(npro,nshl,ndof), dwl(npro,nenl), 56 & acl(npro,nshl,ndof), shpfun(npro,nshl), 57 & shgl(npro,nsd,nshl), xl(npro,nenl,nsd), 58 & aci(npro,nsd), g1yi(npro,ndof), 59 & g2yi(npro,ndof), g3yi(npro,ndof), 60 & shg(npro,nshl,nsd), dxidx(npro,nsd,nsd), 61 & WdetJ(npro), 62 & rho(npro), pres(npro), 63 & u1(npro), u2(npro), 64 & u3(npro), divqi(npro,nflow-1+isurf), 65 & ql(npro,nshl,idflx), rLui(npro,nsd), 66 & src(npro,nsd), Temp(npro),xx(npro,nsd) 67c 68 dimension tmp(npro), tmp1(npro), dkei(npro), dist2w(npro) 69c 70 dimension rlsl(npro,nshl,6), rlsli(npro,6) 71c 72 real*8 rerrl(npro,nshl,6), omega(3), divu(npro) 73 dimension gyti(npro,nsd), gradh(npro,nsd), 74 & sforce(npro,3), weber(npro), 75 & Sclr(npro) 76c 77c.... -------------> Primitive variables at int. point <-------------- 78c 79c.... compute primitive variables 80c 81 pres = zero 82 u1 = zero 83 u2 = zero 84 u3 = zero 85c 86 do n = 1, nshl 87 pres = pres + shpfun(:,n) * yl(:,n,1) 88 u1 = u1 + shpfun(:,n) * yl(:,n,2) 89 u2 = u2 + shpfun(:,n) * yl(:,n,3) 90 u3 = u3 + shpfun(:,n) * yl(:,n,4) 91 enddo 92 if(matflg(5,1).eq.2) then ! boussinesq body force 93 Temp = zero 94 do n = 1, nshl 95 Temp = Temp + shpfun(:,n) * yl(:,n,5) 96 enddo 97 endif 98 if(matflg(5,1).eq.3.or.matflg(6,1).eq.1) then 99c user-specified body force or coriolis force specified 100 xx = zero 101 do n = 1,nenl 102 xx(:,1) = xx(:,1) + shpfun(:,n) * xl(:,n,1) 103 xx(:,2) = xx(:,2) + shpfun(:,n) * xl(:,n,2) 104 xx(:,3) = xx(:,3) + shpfun(:,n) * xl(:,n,3) 105 enddo 106 endif 107c 108 if(iRANS.eq.-2) then ! kay-epsilon 109 dist2w = zero 110 do n = 1, nenl 111 dist2w = dist2w + shpfun(:,n) * dwl(:,n) 112 enddo 113 endif 114c 115 116 if( (iLES.gt.10).and.(iLES.lt.20)) then ! bardina 117 rlsli = zero 118 do n = 1, nshl 119 120 rlsli(:,1) = rlsli(:,1) + shpfun(:,n) * rlsl(:,n,1) 121 rlsli(:,2) = rlsli(:,2) + shpfun(:,n) * rlsl(:,n,2) 122 rlsli(:,3) = rlsli(:,3) + shpfun(:,n) * rlsl(:,n,3) 123 rlsli(:,4) = rlsli(:,4) + shpfun(:,n) * rlsl(:,n,4) 124 rlsli(:,5) = rlsli(:,5) + shpfun(:,n) * rlsl(:,n,5) 125 rlsli(:,6) = rlsli(:,6) + shpfun(:,n) * rlsl(:,n,6) 126 127 enddo 128 else 129 rlsli = zero 130 endif 131c 132c.... -----------------------> accel. at int. point <---------------------- 133c 134 aci = zero 135 do n = 1, nshl 136 aci(:,1) = aci(:,1) + shpfun(:,n) * acl(:,n,2) 137 aci(:,2) = aci(:,2) + shpfun(:,n) * acl(:,n,3) 138 aci(:,3) = aci(:,3) + shpfun(:,n) * acl(:,n,4) 139 enddo 140c 141c.... ---------------------> Element Metrics <----------------------- 142c 143 call e3metric( xl, shgl, dxidx, 144 & shg, WdetJ) 145c 146c.... compute the global gradient of u and P 147c 148c 149 g1yi = zero 150 g2yi = zero 151 g3yi = zero 152 do n = 1, nshl 153 g1yi(:,1) = g1yi(:,1) + shg(:,n,1) * yl(:,n,1) 154 g1yi(:,2) = g1yi(:,2) + shg(:,n,1) * yl(:,n,2) 155 g1yi(:,3) = g1yi(:,3) + shg(:,n,1) * yl(:,n,3) 156 g1yi(:,4) = g1yi(:,4) + shg(:,n,1) * yl(:,n,4) 157c 158 g2yi(:,1) = g2yi(:,1) + shg(:,n,2) * yl(:,n,1) 159 g2yi(:,2) = g2yi(:,2) + shg(:,n,2) * yl(:,n,2) 160 g2yi(:,3) = g2yi(:,3) + shg(:,n,2) * yl(:,n,3) 161 g2yi(:,4) = g2yi(:,4) + shg(:,n,2) * yl(:,n,4) 162c 163 g3yi(:,1) = g3yi(:,1) + shg(:,n,3) * yl(:,n,1) 164 g3yi(:,2) = g3yi(:,2) + shg(:,n,3) * yl(:,n,2) 165 g3yi(:,3) = g3yi(:,3) + shg(:,n,3) * yl(:,n,3) 166 g3yi(:,4) = g3yi(:,4) + shg(:,n,3) * yl(:,n,4) 167 enddo 168 169 divqi = zero 170 idflow = 3 171 if ( idiff >= 1 .or. isurf==1 ) then 172c 173c.... compute divergence of diffusive flux vector, qi,i 174c 175 if(idiff >= 1) then 176 do n=1, nshl 177 divqi(:,1) = divqi(:,1) + shg(:,n,1)*ql(:,n,1 ) 178 & + shg(:,n,2)*ql(:,n,4 ) 179 & + shg(:,n,3)*ql(:,n,7 ) 180 181 divqi(:,2) = divqi(:,2) + shg(:,n,1)*ql(:,n,2 ) 182 & + shg(:,n,2)*ql(:,n,5 ) 183 & + shg(:,n,3)*ql(:,n,8) 184 185 divqi(:,3) = divqi(:,3) + shg(:,n,1)*ql(:,n,3 ) 186 & + shg(:,n,2)*ql(:,n,6 ) 187 & + shg(:,n,3)*ql(:,n,9 ) 188 189 enddo 190 191 endif !end of idiff 192c 193 if (isurf .eq. 1) then 194c .... divergence of normal calculation (curvature) 195 do n=1, nshl 196 divqi(:,idflow+1) = divqi(:,idflow+1) 197 & + shg(:,n,1)*ql(:,n,idflx-2) 198 & + shg(:,n,2)*ql(:,n,idflx-1) 199 & + shg(:,n,3)*ql(:,n,idflx) 200 enddo 201c .... initialization of some variables 202 Sclr = zero 203 gradh= zero 204 gyti = zero 205 sforce=zero 206 do i = 1, npro 207 do n = 1, nshl 208 Sclr(i) = Sclr(i) + shpfun(i,n) * yl(i,n,6) !scalar 209c 210c .... compute the global gradient of Scalar variable 211c 212 gyti(i,1) = gyti(i,1) + shg(i,n,1) * yl(i,n,6) 213 gyti(i,2) = gyti(i,2) + shg(i,n,2) * yl(i,n,6) 214 gyti(i,3) = gyti(i,3) + shg(i,n,3) * yl(i,n,6) 215c 216 enddo 217 218 if (abs (sclr(i)) .le. epsilon_ls) then 219 gradh(i,1) = 0.5/epsilon_ls * (1.0 220 & + cos(pi*Sclr(i)/epsilon_ls)) * gyti(i,1) 221 gradh(i,2) = 0.5/epsilon_ls * (1.0 222 & + cos(pi*Sclr(i)/epsilon_ls)) * gyti(i,2) 223 gradh(i,3) = 0.5/epsilon_ls * (1.0 224 & + cos(pi*Sclr(i)/epsilon_ls)) * gyti(i,3) 225 endif 226 enddo !end of the loop over npro 227c 228c .. surface tension force calculation 229c .. divide by density now as it gets multiplied in e3res.f, as surface 230c tension force is already in the form of force per unti volume 231c 232 weber(:) = Bo 233 sforce(:,1) = -(1.0/weber(:)) * divqi(:,idflow+1) !x-direction 234 & *gradh(:,1) /rho(:) 235 sforce(:,2) = -(1.0/weber(:)) * divqi(:,idflow+1) !y-direction 236 & *gradh(:,2) /rho(:) 237 sforce(:,3) = -(1.0/weber(:)) * divqi(:,idflow+1) !z-direction 238 & *gradh(:,3) /rho(:) 239c 240 endif ! end of the surface tension force calculation 241 endif ! diffusive flux computation 242c 243c Calculate strong form of pde as well as the source term 244c 245 call e3resStrongPDE( 246 & aci, u1, u2, u3, Temp, rho, xx, 247 & g1yi, g2yi, g3yi, 248 & rLui, src, divqi) 249c 250c.... take care of the surface tension force term here 251c 252 if (isurf .eq. 1) then ! note multiplied by density in e3res.f 253 src(:,1) = src(:,1) + sforce(:,1) 254 src(:,2) = src(:,2) + sforce(:,2) 255 src(:,3) = src(:,3) + sforce(:,3) 256 endif 257c 258c.... -------------------> error calculation <----------------- 259c 260 if((ierrcalc.eq.1).and.(nitr.eq.iter)) then 261 do ia=1,nshl 262 tmp=shpfun(:,ia)*WdetJ(:) 263 tmp1=shpfun(:,ia)*Qwt(lcsyst,intp) 264 rerrl(:,ia,1) = rerrl(:,ia,1) + 265 & tmp1(:)*rLui(:,1)*rLui(:,1) 266 rerrl(:,ia,2) = rerrl(:,ia,2) + 267 & tmp1(:)*rLui(:,2)*rLui(:,2) 268 rerrl(:,ia,3) = rerrl(:,ia,3) + 269 & tmp1(:)*rLui(:,3)*rLui(:,3) 270 271 rerrl(:,ia,4) = rerrl(:,ia,4) + 272 & tmp(:)*divqi(:,1)*divqi(:,1) 273 rerrl(:,ia,5) = rerrl(:,ia,5) + 274 & tmp(:)*divqi(:,2)*divqi(:,2) 275 rerrl(:,ia,6) = rerrl(:,ia,6) + 276 & tmp(:)*divqi(:,3)*divqi(:,3) 277 enddo 278 endif 279 distcalc=0 ! return to 1 if you want to compute T-S instability 280 if(distcalc.eq.1) then 281c 282c.... -----------------------> dist. kin energy at int. point <-------------- 283c 284 285 if (ires .ne. 2 .and. iter.eq.1) then !only do at beginning of step 286c 287c calc exact velocity for a channel at quadrature points. 288c 289 dkei=0.0 290c 291 do n = 1, nenl 292 dkei = dkei + shpfun(:,n) * (1.0-xl(:,n,2)**2) !u_ex^~ (in FEM space) 293 enddo 294 dkei = (u1(:)-dkei)**2 +u2(:)**2 ! u'^2+v'^2 295 dkei = dkei*WdetJ ! mult function*W*det of jacobian to 296c get this quadrature point contribution 297 dke = dke+sum(dkei) ! we move the sum over elements inside of the 298c sum over quadrature to save memory (we want 299c a scalar only) 300 endif 301 endif 302c 303c.... return 304c 305 return 306 end 307 308c----------------------------------------------------------------------- 309c 310c Calculate the variables for the scalar advection-diffusion 311c equation. 312c 313c----------------------------------------------------------------------- 314 subroutine e3ivarSclr (yl, acl, shpfun, 315 & shgl, xl, xmudmi, 316 & Sclr, Sdot, gradS, 317 & shg, dxidx, WdetJ, 318 & u1, u2, u3, 319 & ql, rLS , SrcR, 320 & SrcL, uMod, dwl, 321 & diffus, srcRat) 322c 323 include "common.h" 324c 325c passed arrays 326c 327 dimension yl(npro,nshl,ndof), acl(npro,nshl,ndof), 328 & Sclr(npro), Sdot(npro), 329 & gradS(npro,nsd), shpfun(npro,nshl), 330 & shgl(npro,nsd,nshl), xl(npro,nenl,nsd), 331 & shg(npro,nshl,nsd), dxidx(npro,nsd,nsd), 332 & WdetJ(npro), 333 & u1(npro), u2(npro), 334 & u3(npro), divS(npro), 335 & ql(npro,nshl,nsd), rLS(npro), 336 & SrcR(npro), SrcL(npro), 337 & dwl(npro,nshl), diffus(npro), 338 & umod(npro,nsd), Temp(npro),xx(npro,nsd), 339 & divqi(npro) 340c 341 dimension tmp(npro), srcRat(npro) 342 real*8 rLui(npro,nsd), aci(npro,nsd), 343 & g1yi(npro,nflow), g2yi(npro,nflow), 344 & g3yi(npro,nflow), 345 & src(npro,nsd), rho(npro), 346 & rmu(npro) 347 real*8 uBar(npro,nsd), xmudmi(npro,ngauss) 348 349c 350c.... -------------> Primitive variables at int. point <-------------- 351c 352c.... compute primitive variables 353c 354 u1 = zero 355 u2 = zero 356 u3 = zero 357 Sclr = zero 358c 359 id=isclr+5 360 do n = 1, nshl 361 u1 = u1 + shpfun(:,n) * yl(:,n,2) 362 u2 = u2 + shpfun(:,n) * yl(:,n,3) 363 u3 = u3 + shpfun(:,n) * yl(:,n,4) 364 Sclr = Sclr + shpfun(:,n) * yl(:,n,id) 365 enddo 366c 367c 368c.... -----------------------> dS/dt at int. point <---------------------- 369c 370 Sdot = zero 371 do n = 1, nshl 372 Sdot = Sdot + shpfun(:,n) * acl(:,n,id) 373 enddo 374c 375c.... ---------------------> Element Metrics <----------------------- 376c 377 378 call e3metric( xl, shgl, dxidx, 379 & shg, WdetJ) 380 381c 382c.... compute the global gradient of u and P 383c 384c 385 gradS = zero 386 do n = 1, nshl 387 gradS(:,1) = gradS(:,1) + shg(:,n,1) * yl(:,n,id) 388 gradS(:,2) = gradS(:,2) + shg(:,n,2) * yl(:,n,id) 389 gradS(:,3) = gradS(:,3) + shg(:,n,3) * yl(:,n,id) 390 enddo 391 392 divS = zero 393 if ( idiff >= 1 ) then 394c 395c.... compute divergence of diffusive flux vector, qi,i 396c 397 do n=1, nshl 398 divS(:) = divS(:) + shg(:,n,1)*ql(:,n,1 ) 399 & + shg(:,n,2)*ql(:,n,2 ) 400 & + shg(:,n,3)*ql(:,n,3 ) 401 enddo 402 endif ! diffusive flux computation 403 404 if(consrv_sclr_conv_vel.eq.1) then 405c Calculate uBar = u - TauM*L, where TauM is the momentum 406c stabilization factor and L is the momentum residual 407 408 if(matflg(5,1).eq.2) then ! boussinesq body force 409 Temp = zero 410 do n = 1, nshl 411 Temp = Temp + shpfun(:,n) * yl(:,n,5) 412 enddo 413 endif 414 if(matflg(5,1).eq.3.or.matflg(6,1).eq.1) then 415c user-specified body force or coriolis force specified 416 xx = zero 417 do n = 1,nenl 418 xx(:,1) = xx(:,1) + shpfun(:,n) * xl(:,n,1) 419 xx(:,2) = xx(:,2) + shpfun(:,n) * xl(:,n,2) 420 xx(:,3) = xx(:,3) + shpfun(:,n) * xl(:,n,3) 421 enddo 422 endif 423 aci = zero 424 do n = 1, nshl 425 aci(:,1) = aci(:,1) + shpfun(:,n) * acl(:,n,2) 426 aci(:,2) = aci(:,2) + shpfun(:,n) * acl(:,n,3) 427 aci(:,3) = aci(:,3) + shpfun(:,n) * acl(:,n,4) 428 enddo 429 g1yi = zero 430 g2yi = zero 431 g3yi = zero 432 do n = 1, nshl 433 g1yi(:,1) = g1yi(:,1) + shg(:,n,1) * yl(:,n,1) 434 g1yi(:,2) = g1yi(:,2) + shg(:,n,1) * yl(:,n,2) 435 g1yi(:,3) = g1yi(:,3) + shg(:,n,1) * yl(:,n,3) 436 g1yi(:,4) = g1yi(:,4) + shg(:,n,1) * yl(:,n,4) 437c 438 g2yi(:,1) = g2yi(:,1) + shg(:,n,2) * yl(:,n,1) 439 g2yi(:,2) = g2yi(:,2) + shg(:,n,2) * yl(:,n,2) 440 g2yi(:,3) = g2yi(:,3) + shg(:,n,2) * yl(:,n,3) 441 g2yi(:,4) = g2yi(:,4) + shg(:,n,2) * yl(:,n,4) 442c 443 g3yi(:,1) = g3yi(:,1) + shg(:,n,3) * yl(:,n,1) 444 g3yi(:,2) = g3yi(:,2) + shg(:,n,3) * yl(:,n,2) 445 g3yi(:,3) = g3yi(:,3) + shg(:,n,3) * yl(:,n,3) 446 g3yi(:,4) = g3yi(:,4) + shg(:,n,3) * yl(:,n,4) 447 enddo 448c 449 if (iLSet .eq. 0)then 450 rho = datmat(1,1,1) 451 rmu = datmat(1,2,1) 452 else 453 write(*,*) 'Not sure if we can handle level set with K-E' 454 write(*,*) '(different uMods? correct value of rho?)' 455 endif 456 divqi=zero ! until we reconstruct q_flow for scalar solve 457 call e3resStrongPDE( 458 & aci, u1, u2, u3, Temp, rho, x, 459 & g1yi, g2yi, g3yi, 460 & rLui, src, divqi) 461 src(:,1)=u1 ! 462 src(:,2)=u2 ! store u in src memory 463 src(:,3)=u3 ! 464c e3uBar calculates Tau_M and assembles uBar 465 call getdiff(dwl, yl, shpfun, xmudmi, xl, rmu, rho) 466 call e3uBar(rho, src, dxidx, rLui, rmu, uBar) 467 u1=ubar(:,1) ! the entire scalar residual 468 u2=ubar(:,2) ! is based on the modified 469 u3=ubar(:,3) ! velocity for conservation 470 endif 471c 472c.... Initialize uMod, the modified velocity uMod 473c We initialize it to u_i and then calculate 474c the correction in e3sourcesclr 475c 476 477 umod(:,1) = u1 478 umod(:,2) = u2 479 umod(:,3) = u3 480c 481c.... compute source terms 482c 483cad 484cad if we are solving the redistancing equation, the umod(:,:) are 485CAD modified in e3sourceSclr. 486CAD 487CAD if we are redistancing levelset variable we want to use a use the 488CAD convective term from the equation. 489 490 491 if(nosource.ne.1) then 492 call e3sourceSclr ( Sclr, Sdot, gradS, dwl, 493 & shpfun, shg, yl, dxidx, 494 & diffus, u1, u2, u3, 495 & xl, srcR, srcL, uMod, 496 & srcRat) 497 else 498 srcRat = zero 499 srcR = zero 500 srcL = zero 501 endif 502c 503c.... -------------------> Scalar residual <----------------- 504c 505 506 rLS(:) = ( Sdot(:) + (u1*gradS(:,1) + 507 & u2*gradS(:,2) + 508 & u3*gradS(:,3)) ) 509 & - divS(:) 510 511c 512c.... return 513c 514 return 515 end 516 517