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, 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 g1yti (npro) : grad-Sclr in direction 1 at intpt 520c g2yti (npro) : grad-Sclr in direction 2 at intpt 521c g3yti (npro) : grad-Sclr in direction 3 at intpt 522c rho (npro) : density at intpt 523c rmu (npro) : molecular viscosity 524c con (npro) : conductivity 525c rk (npro) : kinetic energy 526c u1 (npro) : x1-velocity component at intpt 527c u2 (npro) : x2-velocity component at intpt 528c u3 (npro) : x3-velocity component at intpt 529c shg (npro,nshl,nsd) : element global grad-shape-functions at intpt 530c 531c also used: 532c pres (npro) : pressure at intpt 533c T (npro) : temperature at intpt 534c ei (npro) : internal energy at intpt 535c h (npro) : enthalpy at intpt 536c alfap (npro) : expansivity at intpt 537c betaT (npro) : isothermal compressibility at intpt 538c cp (npro) : specific heat at constant pressure at intpt 539c dxidx (npro,nsd,nsd) : inverse of deformation gradient at intpt 540c divqi (npro,nflow-1) : divergence of diffusive flux 541c 542c 543c Zdenek Johan, Summer 1990. (Modified from e2ivar.f) 544c Zdenek Johan, Winter 1991. (Fortran 90) 545c Kenneth Jansen, Winter 1997. Primitive Variables 546c---------------------------------------------------------------------- 547c 548 include "common.h" 549c 550c passed arrays 551c 552 dimension ycl(npro,nshl,ndof), 553 & acl(npro,nshl,ndof), acti(npro), 554 & shape(npro,nshl), 555 & shgl(npro,nsd,nshl), xl(npro,nenl,nsd), 556 & dui(npro,nflow), 557 & g1yi(npro,nflow), 558 & g2yi(npro,nflow), g3yi(npro,nflow), 559 & shg(npro,nshl,nsd), dxidx(npro,nsd,nsd), 560 & WdetJ(npro), 561 & rho(npro), pres(npro), 562 & T(npro), ei(npro), 563 & h(npro), alfap(npro), 564 & betaT(npro), cp(npro), 565 & rk(npro), 566 & u1(npro), u2(npro), 567 & u3(npro), divqi(npro,nflow-1), 568 & ql(npro,nshl,(nflow-1)*nsd),Sclr(npro), 569 & dwl(npro,nenl), 570 & dist2w(npro), 571 & vort(npro), 572 & rmu(npro), con(npro), 573 & g1yti(npro), 574 & g2yti(npro), g3yti(npro) 575c 576 577 dimension tmp(npro), dxdxi(npro,nsd,nsd) 578c 579 ttim(20) = ttim(20) - tmr() 580c 581c.... -------------> Primitive variables at int. point <-------------- 582c 583c.... compute primitive variables 584c 585 ttim(11) = ttim(11) - tmr() 586 587 pres = zero 588 u1 = zero 589 u2 = zero 590 u3 = zero 591 T = zero 592 Sclr = zero 593 dist2w = zero 594c 595 id = isclr + 5 596 do n = 1, nshl 597c 598c y(intp)=SUM_{a=1}^nshl (N_a(intp) Ya) 599c 600 pres = pres + shape(:,n) * ycl( :,n,1) 601 u1 = u1 + shape(:,n) * ycl( :,n,2) 602 u2 = u2 + shape(:,n) * ycl( :,n,3) 603 u3 = u3 + shape(:,n) * ycl( :,n,4) 604 T = T + shape(:,n) * ycl( :,n,5) 605 Sclr = Sclr + shape(:,n) * ycl(:,n,id) 606 if (iRANS.lt.0) then 607 dist2w = dist2w + shape(:,n) * dwl(:,n) 608 endif 609 enddo 610c 611 ttim(11) = ttim(11) + tmr() 612c 613c.... -----------------------> accel. at intp. point <---------------------- 614c 615 616 if (ires .ne. 2) then 617 ttim(12) = ttim(12) - tmr() 618c 619c.... compute time-derivative of Scalar variable 620c 621 acti = zero 622 do n = 1, nshl 623 acti(:) = acti(:) + shape(:,n) * acl(:,n,id) 624 enddo 625c 626 flops = flops + 10*nshl*npro 627 ttim(12) = ttim(12) + tmr() 628 endif !ires .ne. 2 629c 630c.... -----------------> Thermodynamic Properties <------------------- 631c 632c.... compute the kinetic energy 633c 634 ttim(27) = ttim(27) - tmr() 635c 636 rk = pt5 * ( u1**2 + u2**2 + u3**2 ) 637c 638c.... get the thermodynamic and material properties 639c 640 ithm = 7 641 call getthm (pres, T, Sclr, 642 & rk, rho, tmp, 643 & tmp, tmp, tmp, 644 & cp, tmp, tmp, 645 & tmp, tmp) 646c 647 if (iconvsclr.eq.2) rho=one 648c 649 call getDiffSclr(T, cp, rmu, 650 & tmp, tmp, con, rho, Sclr) 651 652 ttim(27) = ttim(27) + tmr() 653 654 655c 656c.... ---------------------> Element Metrics <----------------------- 657c 658 call e3metric( xl, shgl, dxidx, 659 & shg, WdetJ) 660 661 662 663c.... ---------------------> Global Gradients <----------------------- 664c 665 ttim(13) = ttim(13) - tmr() 666 667 g1yi = zero 668 g2yi = zero 669 g3yi = zero 670c 671c.... compute the global gradient of Y-variables 672c 673c 674c Y_{,x_i}=SUM_{a=1}^nshl (N_{a,x_i}(intp) Ya) 675c 676 do n = 1, nshl 677c g1yi(:,1) = g1yi(:,1) + shg(:,n,1) * ycl(:,n,1) 678c g1yi(:,2) = g1yi(:,2) + shg(:,n,1) * ycl(:,n,2) 679 g1yi(:,3) = g1yi(:,3) + shg(:,n,1) * ycl(:,n,3) 680 g1yi(:,4) = g1yi(:,4) + shg(:,n,1) * ycl(:,n,4) 681c g1yi(:,5) = g1yi(:,5) + shg(:,n,1) * ycl(:,n,5) 682c 683c g2yi(:,1) = g2yi(:,1) + shg(:,n,2) * ycl(:,n,1) 684 g2yi(:,2) = g2yi(:,2) + shg(:,n,2) * ycl(:,n,2) 685c g2yi(:,3) = g2yi(:,3) + shg(:,n,2) * ycl(:,n,3) 686 g2yi(:,4) = g2yi(:,4) + shg(:,n,2) * ycl(:,n,4) 687c g2yi(:,5) = g2yi(:,5) + shg(:,n,2) * ycl(:,n,5) 688c 689c g3yi(:,1) = g3yi(:,1) + shg(:,n,3) * ycl(:,n,1) 690 g3yi(:,2) = g3yi(:,2) + shg(:,n,3) * ycl(:,n,2) 691 g3yi(:,3) = g3yi(:,3) + shg(:,n,3) * ycl(:,n,3) 692c g3yi(:,4) = g3yi(:,4) + shg(:,n,3) * ycl(:,n,4) 693c g3yi(:,5) = g3yi(:,5) + shg(:,n,3) * ycl(:,n,5) 694c 695 enddo 696 697 698 699 g1yti = zero 700 g2yti = zero 701 g3yti = zero 702c 703c.... compute the global gradient of Scalar variable 704c 705c 706c Y_{,x_i}=SUM_{a=1}^nshl (N_{a,x_i}(intp) Ya) 707c 708 do n = 1, nshl 709 g1yti(:) = g1yti(:) + shg(:,n,1) * ycl(:,n,id) 710 g2yti(:) = g2yti(:) + shg(:,n,2) * ycl(:,n,id) 711 g3yti(:) = g3yti(:) + shg(:,n,3) * ycl(:,n,id) 712c 713 enddo 714c ********************************************************** 715c 716c.... compute vorticity at this integration point 717c 718 vort = sqrt( 719 & (g2yi(:,4)-g3yi(:,3))**2 720 & +(g3yi(:,2)-g1yi(:,4))**2 721 & +(g1yi(:,3)-g2yi(:,2))**2 722 & ) 723c *********************************************************** 724 725 ttim(7) = ttim(7) + tmr() 726 727c 728c.... return 729c 730 ttim(20) = ttim(20) + tmr() 731 return 732 end 733 734 735