1 subroutine e3ivar (yl, ycl, acl, 2 & Sclr, shape, shgl, 3 & xl, dui, aci, 4 & g1yi, g2yi, g3yi, 5 & shg, dxidx, WdetJ, 6 & rho, pres, T, 7 & ei, h, alfap, 8 & betaT, cp, rk, 9 & u1, u2, u3, 10 & ql, divqi, sgn, tmp, 11 & rmu, rlm, rlm2mu, 12 & con, rlsl, rlsli, 13 & xmudmi, sforce, cv) 14c 15c---------------------------------------------------------------------- 16c 17c This routine computes the variables at integration point. 18c 19c input: 20c yl (npro,nshl,nflow) : primitive variables (NO SCALARS) 21c ycl (npro,nshl,ndof) : primitive variables at current step 22c acl (npro,nshl,ndof) : prim.var. accel. at the current step 23c shape (npro,nshl) : element shape-functions 24c shgl (nsd,nen) : element local-grad-shape-functions 25c xl (npro,nenl,nsd) : nodal coordinates at current step 26c ql (npro,nshl,(nflow-1)*nsd) : diffusive flux vector 27c rlsl (npro,nshl,6) : resolved Leonard stresses 28c sgn (npro,nshl) : signs of shape functions 29c 30c output: 31c dui (npro,nflow) : delta U variables at current step 32c aci (npro,nflow) : primvar accel. variables at current step 33c g1yi (npro,nflow) : grad-y in direction 1 34c g2yi (npro,nflow) : grad-y in direction 2 35c g3yi (npro,nflow) : grad-y in direction 3 36c shg (npro,nshl,nsd) : element global grad-shape-functions 37c dxidx (npro,nsd,nsd) : inverse of deformation gradient 38c WdetJ (npro) : weighted Jacobian 39c rho (npro) : density 40c pres (npro) : pressure 41c T (npro) : temperature 42c ei (npro) : internal energy 43c h (npro) : enthalpy 44c alfap (npro) : expansivity 45c betaT (npro) : isothermal compressibility 46c cp (npro) : specific heat at constant pressure 47c rk (npro) : kinetic energy 48c u1 (npro) : x1-velocity component 49c u2 (npro) : x2-velocity component 50c u3 (npro) : x3-velocity component 51c divqi (npro,nflow-1) : divergence of diffusive flux 52c rlsli (npro,6) : resolved Leonard stresses at quad pt 53c 54c Zdenek Johan, Summer 1990. (Modified from e2ivar.f) 55c Zdenek Johan, Winter 1991. (Fortran 90) 56c Kenneth Jansen, Winter 1997. Primitive Variables 57c---------------------------------------------------------------------- 58c 59 include "common.h" 60c 61c passed arrays 62c 63 dimension yl(npro,nshl,nflow), ycl(npro,nshl,ndof), 64 & acl(npro,nshl,ndof), 65 & shape(npro,nshl), 66 & shgl(npro,nsd,nshl), xl(npro,nenl,nsd), 67 & dui(npro,nflow), 68 & aci(npro,nflow), g1yi(npro,nflow), 69 & g2yi(npro,nflow), g3yi(npro,nflow), 70 & shg(npro,nshl,nsd), dxidx(npro,nsd,nsd), 71 & WdetJ(npro), 72 & rho(npro), pres(npro), 73 & T(npro), ei(npro), 74 & h(npro), alfap(npro), 75 & betaT(npro), cp(npro), 76 & rk(npro), cv(npro), 77 & u1(npro), u2(npro), 78 & u3(npro), divqi(npro,nflow), 79 & ql(npro,nshl,idflx), 80 & rmu(npro), rlm(npro), 81 & rlm2mu(npro), con(npro), 82 & Sclr(npro) 83c 84 dimension tmp(npro), dxdxi(npro,nsd,nsd), sgn(npro,nshl) 85 dimension rlsl(npro,nshl,6), rlsli(npro,6), 86 & xmudmi(npro,ngauss) 87 dimension gyti(npro,nsd), gradh(npro,nsd), 88 & sforce(npro,3), weber(npro) 89 ttim(20) = ttim(20) - secs(0.0) 90 91c 92 ttim(10) = ttim(10) - secs(0.0) 93 94 dui = zero 95c 96 do n = 1, nshl 97 dui(:,1) = dui(:,1) + shape(:,n) * yl(:,n,1) ! p 98 dui(:,2) = dui(:,2) + shape(:,n) * yl(:,n,2) ! u1 99 dui(:,3) = dui(:,3) + shape(:,n) * yl(:,n,3) ! u2 100 dui(:,4) = dui(:,4) + shape(:,n) * yl(:,n,4) ! u3 101 dui(:,5) = dui(:,5) + shape(:,n) * yl(:,n,5) ! T 102 enddo 103c 104! flops = flops + 10*nshl*npro 105c 106c 107c.... compute conservative variables 108c 109 rk = pt5 * (dui(:,2)**2 + dui(:,3)**2 + dui(:,4)**2) 110c 111 if (iLSet .ne. 0)then 112 Sclr = zero 113 isc=abs(iRANS)+6 114 do n = 1, nshl 115 Sclr = Sclr + shape(:,n) * ycl(:,n,isc) 116 enddo 117 endif 118 119 ithm = 6 120 call getthm (dui(:,1), dui(:,5), Sclr, 121 & rk, rho, ei, 122 & tmp, tmp, tmp, 123 & tmp, tmp, tmp, 124 & tmp, tmp) 125c 126c 127 dui(:,1) = rho 128 dui(:,2) = rho * dui(:,2) 129 dui(:,3) = rho * dui(:,3) 130 dui(:,4) = rho * dui(:,4) 131 dui(:,5) = rho * (ei + rk) 132 133 134 ttim(10) = ttim(10) + secs(0.0) 135c 136c.... -------------> Primitive variables at int. point <-------------- 137c 138c.... compute primitive variables 139c 140 ttim(11) = ttim(11) - secs(0.0) 141 142 pres = zero 143 u1 = zero 144 u2 = zero 145 u3 = zero 146 T = zero 147 do n = 1, nshl 148c 149c y(int)=SUM_{a=1}^nshl (N_a(int) Ya) 150c 151 pres = pres + shape(:,n) * ycl(:,n,1) 152 u1 = u1 + shape(:,n) * ycl(:,n,2) 153 u2 = u2 + shape(:,n) * ycl(:,n,3) 154 u3 = u3 + shape(:,n) * ycl(:,n,4) 155 T = T + shape(:,n) * ycl(:,n,5) 156 enddo 157 158 if( (iLES.gt.10).and.(iLES.lt.20)) then ! bardina 159 rlsli = zero 160 do n = 1, nshl 161 162 rlsli(:,1) = rlsli(:,1) + shape(:,n) * rlsl(:,n,1) 163 rlsli(:,2) = rlsli(:,2) + shape(:,n) * rlsl(:,n,2) 164 rlsli(:,3) = rlsli(:,3) + shape(:,n) * rlsl(:,n,3) 165 rlsli(:,4) = rlsli(:,4) + shape(:,n) * rlsl(:,n,4) 166 rlsli(:,5) = rlsli(:,5) + shape(:,n) * rlsl(:,n,5) 167 rlsli(:,6) = rlsli(:,6) + shape(:,n) * rlsl(:,n,6) 168 169 enddo 170 else 171 rlsli = zero 172 endif 173c 174 175 ttim(11) = ttim(11) + secs(0.0) 176 177c 178c.... -----------------------> accel. at int. point <---------------------- 179c 180 181c if (ires .ne. 2) then 182 ttim(12) = ttim(12) - secs(0.0) 183c 184c.... compute primitive variables 185c 186 aci = zero 187c 188 do n = 1, nshl 189 aci(:,1) = aci(:,1) + shape(:,n) * acl(:,n,1) 190 aci(:,2) = aci(:,2) + shape(:,n) * acl(:,n,2) 191 aci(:,3) = aci(:,3) + shape(:,n) * acl(:,n,3) 192 aci(:,4) = aci(:,4) + shape(:,n) * acl(:,n,4) 193 aci(:,5) = aci(:,5) + shape(:,n) * acl(:,n,5) 194 enddo 195c 196! flops = flops + 10*nshl*npro 197 ttim(12) = ttim(12) + secs(0.0) 198c endif !ires .ne. 2 199c 200c.... -----------------> Thermodynamic Properties <------------------- 201c 202c.... compute the kinetic energy 203c 204 ttim(27) = ttim(27) - secs(0.0) 205c 206 rk = pt5 * ( u1**2 + u2**2 + u3**2 ) 207c 208c.... get the thermodynamic properties 209c 210 if (iLSet .ne. 0)then 211 Sclr = zero 212 isc=abs(iRANS)+6 213 do n = 1, nshl 214 Sclr = Sclr + shape(:,n) * ycl(:,n,isc) 215 enddo 216 endif 217 218 ithm = 7 219 call getthm (pres, T, Sclr, 220 & rk, rho, ei, 221 & h, tmp, cv, 222 & cp, alfap, betaT, 223 & tmp, tmp) 224c 225 ttim(27) = ttim(27) + secs(0.0) 226c 227c ........Get the material properties 228c 229 call getDiff (T, cp, rho, ycl, 230 & rmu, rlm, rlm2mu, con, shape, 231 & xmudmi, xl) 232 233c.... ---------------------> Element Metrics <----------------------- 234c 235 ttim(26) = ttim(26) - secs(0.0) 236c 237 call e3metric( xl, shgl, dxidx, 238 & shg, WdetJ) 239c 240c 241 ttim(26) = ttim(26) + secs(0.0) 242c 243c.... ---------------------> Global Gradients <----------------------- 244c 245 ttim(13) = ttim(13) - secs(0.0) 246 247 g1yi = zero 248 g2yi = zero 249 g3yi = zero 250c 251 ttim(13) = ttim(13) + secs(0.0) 252 ttim(7) = ttim(7) - secs(0.0) 253c 254c.... compute the global gradient of Y-variables 255c 256c 257c Y_{,x_i}=SUM_{a=1}^nshl (N_{a,x_i}(int) Ya) 258c 259 if(nshl.eq.4) then 260 g1yi(:,1) = g1yi(:,1) + shg(:,1,1) * yl(:,1,1) 261 & + shg(:,2,1) * yl(:,2,1) 262 & + shg(:,3,1) * yl(:,3,1) 263 & + shg(:,4,1) * yl(:,4,1) 264c 265 g1yi(:,2) = g1yi(:,2) + shg(:,1,1) * yl(:,1,2) 266 & + shg(:,2,1) * yl(:,2,2) 267 & + shg(:,3,1) * yl(:,3,2) 268 & + shg(:,4,1) * yl(:,4,2) 269c 270 g1yi(:,3) = g1yi(:,3) + shg(:,1,1) * yl(:,1,3) 271 & + shg(:,2,1) * yl(:,2,3) 272 & + shg(:,3,1) * yl(:,3,3) 273 & + shg(:,4,1) * yl(:,4,3) 274c 275 g1yi(:,4) = g1yi(:,4) + shg(:,1,1) * yl(:,1,4) 276 & + shg(:,2,1) * yl(:,2,4) 277 & + shg(:,3,1) * yl(:,3,4) 278 & + shg(:,4,1) * yl(:,4,4) 279c 280 g1yi(:,5) = g1yi(:,5) + shg(:,1,1) * yl(:,1,5) 281 & + shg(:,2,1) * yl(:,2,5) 282 & + shg(:,3,1) * yl(:,3,5) 283 & + shg(:,4,1) * yl(:,4,5) 284c 285 g2yi(:,1) = g2yi(:,1) + shg(:,1,2) * yl(:,1,1) 286 & + shg(:,2,2) * yl(:,2,1) 287 & + shg(:,3,2) * yl(:,3,1) 288 & + shg(:,4,2) * yl(:,4,1) 289c 290 g2yi(:,2) = g2yi(:,2) + shg(:,1,2) * yl(:,1,2) 291 & + shg(:,2,2) * yl(:,2,2) 292 & + shg(:,3,2) * yl(:,3,2) 293 & + shg(:,4,2) * yl(:,4,2) 294c 295 g2yi(:,3) = g2yi(:,3) + shg(:,1,2) * yl(:,1,3) 296 & + shg(:,2,2) * yl(:,2,3) 297 & + shg(:,3,2) * yl(:,3,3) 298 & + shg(:,4,2) * yl(:,4,3) 299c 300 g2yi(:,4) = g2yi(:,4) + shg(:,1,2) * yl(:,1,4) 301 & + shg(:,2,2) * yl(:,2,4) 302 & + shg(:,3,2) * yl(:,3,4) 303 & + shg(:,4,2) * yl(:,4,4) 304c 305 g2yi(:,5) = g2yi(:,5) + shg(:,1,2) * yl(:,1,5) 306 & + shg(:,2,2) * yl(:,2,5) 307 & + shg(:,3,2) * yl(:,3,5) 308 & + shg(:,4,2) * yl(:,4,5) 309c 310 g3yi(:,1) = g3yi(:,1) + shg(:,1,3) * yl(:,1,1) 311 & + shg(:,2,3) * yl(:,2,1) 312 & + shg(:,3,3) * yl(:,3,1) 313 & + shg(:,4,3) * yl(:,4,1) 314c 315 g3yi(:,2) = g3yi(:,2) + shg(:,1,3) * yl(:,1,2) 316 & + shg(:,2,3) * yl(:,2,2) 317 & + shg(:,3,3) * yl(:,3,2) 318 & + shg(:,4,3) * yl(:,4,2) 319c 320 g3yi(:,3) = g3yi(:,3) + shg(:,1,3) * yl(:,1,3) 321 & + shg(:,2,3) * yl(:,2,3) 322 & + shg(:,3,3) * yl(:,3,3) 323 & + shg(:,4,3) * yl(:,4,3) 324c 325 g3yi(:,4) = g3yi(:,4) + shg(:,1,3) * yl(:,1,4) 326 & + shg(:,2,3) * yl(:,2,4) 327 & + shg(:,3,3) * yl(:,3,4) 328 & + shg(:,4,3) * yl(:,4,4) 329c 330 g3yi(:,5) = g3yi(:,5) + shg(:,1,3) * yl(:,1,5) 331 & + shg(:,2,3) * yl(:,2,5) 332 & + shg(:,3,3) * yl(:,3,5) 333 & + shg(:,4,3) * yl(:,4,5) 334c 335 else 336 do n = 1, nshl 337 g1yi(:,1) = g1yi(:,1) + shg(:,n,1) * yl(:,n,1) 338 g1yi(:,2) = g1yi(:,2) + shg(:,n,1) * yl(:,n,2) 339 g1yi(:,3) = g1yi(:,3) + shg(:,n,1) * yl(:,n,3) 340 g1yi(:,4) = g1yi(:,4) + shg(:,n,1) * yl(:,n,4) 341 g1yi(:,5) = g1yi(:,5) + shg(:,n,1) * yl(:,n,5) 342c 343 g2yi(:,1) = g2yi(:,1) + shg(:,n,2) * yl(:,n,1) 344 g2yi(:,2) = g2yi(:,2) + shg(:,n,2) * yl(:,n,2) 345 g2yi(:,3) = g2yi(:,3) + shg(:,n,2) * yl(:,n,3) 346 g2yi(:,4) = g2yi(:,4) + shg(:,n,2) * yl(:,n,4) 347 g2yi(:,5) = g2yi(:,5) + shg(:,n,2) * yl(:,n,5) 348c 349 g3yi(:,1) = g3yi(:,1) + shg(:,n,3) * yl(:,n,1) 350 g3yi(:,2) = g3yi(:,2) + shg(:,n,3) * yl(:,n,2) 351 g3yi(:,3) = g3yi(:,3) + shg(:,n,3) * yl(:,n,3) 352 g3yi(:,4) = g3yi(:,4) + shg(:,n,3) * yl(:,n,4) 353 g3yi(:,5) = g3yi(:,5) + shg(:,n,3) * yl(:,n,5) 354c 355 enddo 356 endif 357 358 359c 360c 361 divqi = zero 362 idflow = 0 363 if (((idiff >= 1) .or. isurf==1).and. 364 & (ires == 3 .or. ires==1)) then 365 ttim(32) = ttim(32) - tmr() 366c 367c CLASS please ignore all terms related to qi until after you 368c understand EVERYTHING ELSE IN THE CODE. You may run with 369c idiff=0 for the whole class and everything will be ok never 370c having understood this part. (leave it for later). 371c 372c.... compute divergence of diffusive flux vector, qi,i 373c 374 if(idiff >= 1) then 375 idflow = idflow + 4 376 do n=1, nshl 377 378 divqi(:,1) = divqi(:,1) + shg(:,n,1)*ql(:,n,1 ) 379 & + shg(:,n,2)*ql(:,n,5 ) 380 & + shg(:,n,3)*ql(:,n,9 ) 381 382 divqi(:,2) = divqi(:,2) + shg(:,n,1)*ql(:,n,2 ) 383 & + shg(:,n,2)*ql(:,n,6 ) 384 & + shg(:,n,3)*ql(:,n,10) 385 386 divqi(:,3) = divqi(:,3) + shg(:,n,1)*ql(:,n,3 ) 387 & + shg(:,n,2)*ql(:,n,7 ) 388 & + shg(:,n,3)*ql(:,n,11) 389 390 divqi(:,4) = divqi(:,4) + shg(:,n,1)*ql(:,n,4 ) 391 & + shg(:,n,2)*ql(:,n,8 ) 392 & + shg(:,n,3)*ql(:,n,12) 393 394 enddo 395 endif !end of idiff 396c 397 if (isurf .eq. 1) then 398c .... divergence of normal calculation 399 do n=1, nshl 400 divqi(:,idflow+1) = divqi(:,idflow+1) 401 & + shg(:,n,1)*ql(:,n,idflx-2) 402 & + shg(:,n,2)*ql(:,n,idflx-1) 403 & + shg(:,n,3)*ql(:,n,idflx) 404 enddo 405c .... initialization of some variables 406 Sclr = zero 407 gradh= zero 408 gyti = zero 409 sforce=zero 410 do i = 1, npro 411 do n = 1, nshl 412 Sclr(i) = Sclr(i) + shape(i,n) * ycl(i,n,6) !scalar 413c 414c .... compute the global gradient of Scalar variable 415c 416 gyti(i,1) = gyti(i,1) + shg(i,n,1) * ycl(i,n,6) 417 gyti(i,2) = gyti(i,2) + shg(i,n,2) * ycl(i,n,6) 418 gyti(i,3) = gyti(i,3) + shg(i,n,3) * ycl(i,n,6) 419c 420 enddo 421 422 if (abs (sclr(i)) .le. epsilon_ls) then 423 gradh(i,1) = 0.5/epsilon_ls * (1 424 & + cos(pi*Sclr(i)/epsilon_ls)) * gyti(i,1) 425 gradh(i,2) = 0.5/epsilon_ls * (1 426 & + cos(pi*Sclr(i)/epsilon_ls)) * gyti(i,2) 427 gradh(i,3) = 0.5/epsilon_ls * (1 428 & + cos(pi*Sclr(i)/epsilon_ls)) * gyti(i,3) 429 endif 430 enddo !end of the loop over npro 431c 432c .... surface tension force calculation 433c .. divide by density now as it gets multiplied in e3res.f, as surface 434c tension force is already in the form of force per unit volume 435c 436 437 weber(:)= Bo 438 sforce(:,1) = -(1.0/weber(:)) * divqi(:,idflow+1) !x-direction 439 & *gradh(:,1)/rho(:) 440 sforce(:,2) = -(1.0/weber(:)) * divqi(:,idflow+1) !y-direction 441 & *gradh(:,2)/rho(:) 442 sforce(:,3) = -(1.0/weber(:)) * divqi(:,idflow+1) !z-direction 443 & *gradh(:,3)/rho(:) 444 445 endif ! end of the surface tension force calculation 446 447 ttim(32) = ttim(32) + secs(0.0) 448 449 endif ! diffusive flux computation 450c 451c.... return 452c 453 ttim(20) = ttim(20) + secs(0.0) 454c 455c.... -----------------------> dist. kin energy at int. point <-------------- 456c 457c 458c if (ires .ne. 2 .and. iter.eq.1) then !only do at beginning of step 459c 460c calc exact velocity for a channel at quadrature points. 461c 462c dkei=0.0 463c 464c first guess would be this but it is wrong since it measures the 465c error outside of FEM space as well. This would be correct if we 466c wanted to measure error but for disturbance we would like to get 467c zero if the solution was nodally exact (i.e no disturbance added to 468c the base flow which is nodally exact). 469c 470c do n = 1, nshl 471c dkei = dkei + shape(:,n) * xl(:,n,2) ! this is just the y coord 472c enddo 473c dkei = (1.0-dkei*dkei) ! u_exact with u_cl=1 474c 475c 476c correct way 477c 478c do n = 1, nshl 479c dkei = dkei + shape(:,n) * (1.0-xl(:,n,2)**2) !u_ex^~ (in FEM space) 480c enddo 481c dkei = (u1-dkei)**2 +u2**2 ! u'^2+v'^2 482c dkei = dkei*WdetJ ! mult function*W*det of jacobian to 483c get this quadrature point contribution 484c dke = dke+sum(dkei) ! we move the sum over elements inside of the 485c sum over quadrature to save memory (we want 486c a scalar only) 487c endif 488 return 489 end 490 subroutine e3ivarSclr (ycl, acl, acti, 491 & shape, shgl, xl, 492 & T, cp, 493 & dxidx, Sclr, 494 & WdetJ, vort, gVnrm, 495 & g1yti, g2yti, g3yti, 496 & rho, rmu, con, 497 & rk, u1, u2, 498 & u3, shg, dwl, 499 & dist2w) 500c 501c---------------------------------------------------------------------- 502c 503c This routine computes the variables at integration point. 504c 505c input: 506c ycl (npro,nshl,ndof) : primitive variables 507c actl (npro,nshl) : time-deriv of ytl 508c dwl (npro,nshl) : distances to wall 509c shape (npro,nshl) : element shape-functions 510c shgl (npro,nsd,nshl) : element local-grad-shape-functions 511c xl (npro,nenl,nsd) : nodal coordinates at current step 512c 513c output: 514c acti (npro) : time-deriv of Scalar variable 515c Sclr (npro) : Scalar variable at intpt 516c dist2w (npro) : distance to nearest wall at intpt 517c WdetJ (npro) : weighted Jacobian at intpt 518c vort (npro) : vorticity at intpt 519c gVnrm (npro) : gradV norm at intpt 520c g1yti (npro) : grad-Sclr in direction 1 at intpt 521c g2yti (npro) : grad-Sclr in direction 2 at intpt 522c g3yti (npro) : grad-Sclr in direction 3 at intpt 523c rho (npro) : density at intpt 524c rmu (npro) : molecular viscosity 525c con (npro) : conductivity 526c rk (npro) : kinetic energy 527c u1 (npro) : x1-velocity component at intpt 528c u2 (npro) : x2-velocity component at intpt 529c u3 (npro) : x3-velocity component at intpt 530c shg (npro,nshl,nsd) : element global grad-shape-functions at intpt 531c 532c also used: 533c pres (npro) : pressure at intpt 534c T (npro) : temperature at intpt 535c ei (npro) : internal energy at intpt 536c h (npro) : enthalpy at intpt 537c alfap (npro) : expansivity at intpt 538c betaT (npro) : isothermal compressibility at intpt 539c cp (npro) : specific heat at constant pressure at intpt 540c dxidx (npro,nsd,nsd) : inverse of deformation gradient at intpt 541c divqi (npro,nflow-1) : divergence of diffusive flux 542c 543c 544c Zdenek Johan, Summer 1990. (Modified from e2ivar.f) 545c Zdenek Johan, Winter 1991. (Fortran 90) 546c Kenneth Jansen, Winter 1997. Primitive Variables 547c---------------------------------------------------------------------- 548c 549 include "common.h" 550c 551c passed arrays 552c 553 dimension ycl(npro,nshl,ndof), 554 & acl(npro,nshl,ndof), acti(npro), 555 & shape(npro,nshl), 556 & shgl(npro,nsd,nshl), xl(npro,nenl,nsd), 557 & dui(npro,nflow), 558 & g1yi(npro,nflow), 559 & g2yi(npro,nflow), g3yi(npro,nflow), 560 & shg(npro,nshl,nsd), dxidx(npro,nsd,nsd), 561 & WdetJ(npro), 562 & rho(npro), pres(npro), 563 & T(npro), ei(npro), 564 & h(npro), alfap(npro), 565 & betaT(npro), cp(npro), 566 & rk(npro), 567 & u1(npro), u2(npro), 568 & u3(npro), divqi(npro,nflow-1), 569 & ql(npro,nshl,(nflow-1)*nsd),Sclr(npro), 570 & dwl(npro,nenl), 571 & dist2w(npro), 572 & vort(npro), gVnrm(npro), 573 & rmu(npro), con(npro), 574 & g1yti(npro), 575 & g2yti(npro), g3yti(npro) 576c 577 578 dimension tmp(npro), dxdxi(npro,nsd,nsd) 579c 580 ttim(20) = ttim(20) - tmr() 581c 582c.... -------------> Primitive variables at int. point <-------------- 583c 584c.... compute primitive variables 585c 586 ttim(11) = ttim(11) - tmr() 587 588 pres = zero 589 u1 = zero 590 u2 = zero 591 u3 = zero 592 T = zero 593 Sclr = zero 594 dist2w = zero 595c 596 id = isclr + 5 597 do n = 1, nshl 598c 599c y(intp)=SUM_{a=1}^nshl (N_a(intp) Ya) 600c 601 pres = pres + shape(:,n) * ycl( :,n,1) 602 u1 = u1 + shape(:,n) * ycl( :,n,2) 603 u2 = u2 + shape(:,n) * ycl( :,n,3) 604 u3 = u3 + shape(:,n) * ycl( :,n,4) 605 T = T + shape(:,n) * ycl( :,n,5) 606 Sclr = Sclr + shape(:,n) * ycl(:,n,id) 607 if (iRANS.lt.0) then 608 dist2w = dist2w + shape(:,n) * dwl(:,n) 609 endif 610 enddo 611c 612 ttim(11) = ttim(11) + tmr() 613c 614c.... -----------------------> accel. at intp. point <---------------------- 615c 616 617 if (ires .ne. 2) then 618 ttim(12) = ttim(12) - tmr() 619c 620c.... compute time-derivative of Scalar variable 621c 622 acti = zero 623 do n = 1, nshl 624 acti(:) = acti(:) + shape(:,n) * acl(:,n,id) 625 enddo 626c 627! flops = flops + 10*nshl*npro 628 ttim(12) = ttim(12) + tmr() 629 endif !ires .ne. 2 630c 631c.... -----------------> Thermodynamic Properties <------------------- 632c 633c.... compute the kinetic energy 634c 635 ttim(27) = ttim(27) - tmr() 636c 637 rk = pt5 * ( u1**2 + u2**2 + u3**2 ) 638c 639c.... get the thermodynamic and material properties 640c 641 ithm = 7 642 call getthm (pres, T, Sclr, 643 & rk, rho, tmp, 644 & tmp, tmp, tmp, 645 & cp, tmp, tmp, 646 & tmp, tmp) 647c 648 if (iconvsclr.eq.2) rho=one 649c 650 call getDiffSclr(T, cp, rmu, 651 & tmp, tmp, con, rho, Sclr) 652 653 ttim(27) = ttim(27) + tmr() 654 655 656c 657c.... ---------------------> Element Metrics <----------------------- 658c 659 call e3metric( xl, shgl, dxidx, 660 & shg, WdetJ) 661 662 663 664c.... ---------------------> Global Gradients <----------------------- 665c 666 ttim(13) = ttim(13) - tmr() 667 668 g1yi = zero 669 g2yi = zero 670 g3yi = zero 671c 672c.... compute the global gradient of Y-variables 673c 674c 675c Y_{,x_i}=SUM_{a=1}^nshl (N_{a,x_i}(intp) Ya) 676c 677 do n = 1, nshl 678c g1yi(:,1) = g1yi(:,1) + shg(:,n,1) * ycl(:,n,1) 679 g1yi(:,2) = g1yi(:,2) + shg(:,n,1) * ycl(:,n,2) 680 g1yi(:,3) = g1yi(:,3) + shg(:,n,1) * ycl(:,n,3) 681 g1yi(:,4) = g1yi(:,4) + shg(:,n,1) * ycl(:,n,4) 682c g1yi(:,5) = g1yi(:,5) + shg(:,n,1) * ycl(:,n,5) 683c 684c g2yi(:,1) = g2yi(:,1) + shg(:,n,2) * ycl(:,n,1) 685 g2yi(:,2) = g2yi(:,2) + shg(:,n,2) * ycl(:,n,2) 686 g2yi(:,3) = g2yi(:,3) + shg(:,n,2) * ycl(:,n,3) 687 g2yi(:,4) = g2yi(:,4) + shg(:,n,2) * ycl(:,n,4) 688c g2yi(:,5) = g2yi(:,5) + shg(:,n,2) * ycl(:,n,5) 689c 690c g3yi(:,1) = g3yi(:,1) + shg(:,n,3) * ycl(:,n,1) 691 g3yi(:,2) = g3yi(:,2) + shg(:,n,3) * ycl(:,n,2) 692 g3yi(:,3) = g3yi(:,3) + shg(:,n,3) * ycl(:,n,3) 693 g3yi(:,4) = g3yi(:,4) + shg(:,n,3) * ycl(:,n,4) 694c g3yi(:,5) = g3yi(:,5) + shg(:,n,3) * ycl(:,n,5) 695c 696 enddo 697 698 699 700 g1yti = zero 701 g2yti = zero 702 g3yti = zero 703c 704c.... compute the global gradient of Scalar variable 705c 706c 707c Y_{,x_i}=SUM_{a=1}^nshl (N_{a,x_i}(intp) Ya) 708c 709 do n = 1, nshl 710 g1yti(:) = g1yti(:) + shg(:,n,1) * ycl(:,n,id) 711 g2yti(:) = g2yti(:) + shg(:,n,2) * ycl(:,n,id) 712 g3yti(:) = g3yti(:) + shg(:,n,3) * ycl(:,n,id) 713c 714 enddo 715c ********************************************************** 716c 717c.... compute vorticity at this integration point 718c 719 vort = sqrt( 720 & (g2yi(:,4)-g3yi(:,3))**2 721 & +(g3yi(:,2)-g1yi(:,4))**2 722 & +(g1yi(:,3)-g2yi(:,2))**2 723 & ) 724 if(iles.lt.0) then 725 gVnrm = sqrt(g1yi(:,2)**2+g2yi(:,2)**2+g3yi(:,2)**2 726 & +g1yi(:,3)**2+g2yi(:,3)**2+g3yi(:,3)**2 727 & +g1yi(:,4)**2+g2yi(:,4)**2+g3yi(:,4)**2) 728 endif ! safe to not define since use is only in same condtnl 729 730c *********************************************************** 731 732 ttim(7) = ttim(7) + tmr() 733 734c 735c.... return 736c 737 ttim(20) = ttim(20) + tmr() 738 return 739 end 740 741 742