1 subroutine e3visc (g1yi, g2yi, g3yi, 2 & dxidx, 3 & rmu, rlm, rlm2mu, 4 & u1, u2, u3, 5 & ri, rmi, stiff, 6 & con, rlsli, compK, T) 7c 8c---------------------------------------------------------------------- 9c 10c This routine calculates the contribution of viscous and heat fluxes 11c to both RHS and LHS. 12c 13c input: 14c g1yi (npro,nflow) : grad-y in direction 1 15c g2yi (npro,nflow) : grad-y in direction 2 16c g3yi (npro,nflow) : grad-y in direction 3 17c u1 (npro) : x1-velocity component 18c u2 (npro) : x2-velocity component 19c u3 (npro) : x3-velocity component 20c rmu (npro) : viscosity 21c rlm (npro) : Lambda 22c rlm2mu (npro) : Lambda + 2 Mu 23c con (npro) : Conductivity 24c cp (npro) : specific heat at constant pressure 25c rlsli (npro,6) : Resolved Reynold's stresses 26c output: 27c ri (npro,nflow*(nsd+1)) : partial residual 28c rmi (npro,nflow*(nsd+1)) : partial modified residual 29c stiff (npro,nsd*nflow,nsd*nflow) : stiffness matrix 30c compK (npro,10) : K_ij in (eq.134) A new ... III 31c 32c 33c Zdenek Johan, Summer 1990. (Modified from e2visc.f) 34c Zdenek Johan, Winter 1991. (Fortran 90) 35c Kenneth Jansen, Winter 1997 Primitive Variables 36c---------------------------------------------------------------------- 37c 38 include "common.h" 39c 40c passed arrays 41c 42 dimension g1yi(npro,nflow), g2yi(npro,nflow), 43 & g3yi(npro,nflow), 44 & Sclr(npro), dxidx(npro,nsd,nsd), 45 & u1(npro), u2(npro), 46 & u3(npro), rho(npro), 47 & ri(npro,nflow*(nsd+1)), rmi(npro,nflow*(nsd+1)) 48c 49c 50 dimension rmu(npro), rlm(npro), 51 & rlm2mu(npro), con(npro), 52 & stiff(npro,3*nflow,3*nflow), 53 & rlsli(npro,6), compK(npro,10), 54 & f1(npro), f2(npro), f3(npro), f4(npro), 55 & f5(npro), f6(npro), T(npro), rk(npro) 56 57 ttim(23) = ttim(23) - secs(0.0) 58c 59c... dynamic model is now being accounted for in getdiff.f 60c 61c 62c.... ---------------------> Diffusivity Matrix <------------------- 63c 64 65 if (lhs .eq. 1) then 66c 67c.... K11 68c 69 stiff(:, 2, 2) = rlm2mu 70 stiff(:, 3, 3) = rmu 71 stiff(:, 4, 4) = rmu 72 stiff(:, 5, 2) = rlm2mu * u1 73 stiff(:, 5, 3) = rmu * u2 74 stiff(:, 5, 4) = rmu * u3 75 stiff(:, 5, 5) = con 76c 77c.... K12 78c 79 stiff(:, 2, 8) = rlm 80 stiff(:, 3, 7) = rmu 81 stiff(:, 5, 7) = rmu * u2 82 stiff(:, 5, 8) = rlm * u1 83c 84c.... K13 85c 86 stiff(:, 2,14) = rlm 87 stiff(:, 4,12) = rmu 88 stiff(:, 5,12) = rmu * u3 89 stiff(:, 5,14) = rlm * u1 90 91c 92c.... K21 93c 94 stiff(:, 7, 3) = rmu 95 stiff(:, 8, 2) = rlm 96 stiff(:,10, 2) = rlm * u2 97 stiff(:,10, 3) = rmu * u1 98 99c 100c.... K22 101c 102 stiff(:, 7, 7) = rmu 103 stiff(:, 8, 8) = rlm2mu 104 stiff(:, 9, 9) = rmu 105 stiff(:,10, 7) = rmu * u1 106 stiff(:,10, 8) = rlm2mu * u2 107 stiff(:,10, 9) = rmu * u3 108 stiff(:,10,10) = con 109c 110c.... K23 111c 112 stiff(:, 8,14) = rlm 113 stiff(:, 9,13) = rmu 114 stiff(:,10,13) = rmu * u3 115 stiff(:,10,14) = rlm * u2 116c 117c.... K31 118c 119 stiff(:,12, 4) = rmu 120 stiff(:,14, 2) = rlm 121 stiff(:,15, 2) = rlm * u3 122 stiff(:,15, 4) = rmu * u1 123c 124c.... K32 125c 126 stiff(:,13, 9) = rmu 127 stiff(:,14, 8) = rlm 128 stiff(:,15, 8) = rlm * u3 129 stiff(:,15, 9) = rmu * u2 130c 131c.... K33 132c 133 stiff(:,12,12) = rmu 134 stiff(:,13,13) = rmu 135 stiff(:,14,14) = rlm2mu 136 stiff(:,15,12) = rmu * u1 137 stiff(:,15,13) = rmu * u2 138 stiff(:,15,14) = rlm2mu * u3 139 stiff(:,15,15) = con 140 141 endif 142 143 if (itau .ge. 10) then ! non-diag tau, buld compK 144 145c 146c.... calculate element factors 147c 148 f1 = dxidx(:,1,1) * dxidx(:,1,1) + 149 & dxidx(:,2,1) * dxidx(:,2,1) + 150 & dxidx(:,3,1) * dxidx(:,3,1) 151 f2 = dxidx(:,1,1) * dxidx(:,1,2) + 152 & dxidx(:,2,1) * dxidx(:,2,2) + 153 & dxidx(:,3,1) * dxidx(:,3,2) 154 f3 = dxidx(:,1,2) * dxidx(:,1,2) + 155 & dxidx(:,2,2) * dxidx(:,2,2) + 156 & dxidx(:,3,2) * dxidx(:,3,2) 157 f4 = dxidx(:,1,1) * dxidx(:,1,3) + 158 & dxidx(:,2,1) * dxidx(:,2,3) + 159 & dxidx(:,3,1) * dxidx(:,3,3) 160 f5 = dxidx(:,1,2) * dxidx(:,1,3) + 161 & dxidx(:,2,2) * dxidx(:,2,3) + 162 & dxidx(:,3,2) * dxidx(:,3,3) 163 f6 = dxidx(:,1,3) * dxidx(:,1,3) + 164 & dxidx(:,2,3) * dxidx(:,2,3) + 165 & dxidx(:,3,3) * dxidx(:,3,3) 166c 167c.... correction for tetrahedra (invariance w.r.t triangular coord) 168c 169 if (lcsyst .eq. 1) then 170 f1 = ( f1 + (dxidx(:,1,1) + 171 & dxidx(:,2,1) + dxidx(:,3,1))**2 ) * pt39 172 f2 = ( f2 + (dxidx(:,1,1) + 173 & dxidx(:,2,1) + dxidx(:,3,1)) * 174 & (dxidx(:,1,2) + 175 & dxidx(:,2,2) + dxidx(:,3,2)) ) * pt39 176 f3 = ( f3 + (dxidx(:,1,2) + 177 & dxidx(:,2,2) + dxidx(:,3,2))**2 ) * pt39 178 f4 = ( f4 + (dxidx(:,1,1) + 179 & dxidx(:,2,1) + dxidx(:,3,1)) * 180 & (dxidx(:,1,3) + 181 & dxidx(:,2,3) + dxidx(:,3,3)) ) * pt39 182 f5 = ( f5 + (dxidx(:,1,2) + 183 & dxidx(:,2,2) + dxidx(:,3,2)) * 184 & (dxidx(:,1,3) + 185 & dxidx(:,2,3) + dxidx(:,3,3)) ) * pt39 186 f6 = ( f6 + (dxidx(:,1,3) + 187 & dxidx(:,2,3) + dxidx(:,3,3))**2 ) * pt39 188c 189 ! flops = flops + 36*npro 190 endif 191c 192c.... correction for wedges 193c 194 if ((lcsyst .eq. 3) .or. (lcsyst .eq. 4)) then 195 f1 = ( dxidx(:,1,1) * dxidx(:,1,1) + 196 & dxidx(:,2,1) * dxidx(:,2,1) + 197 & ( dxidx(:,1,1) + dxidx(:,2,1) )**2 ) * pt57 198 & + dxidx(:,3,1) * dxidx(:,3,1) 199 f2 = ( dxidx(:,1,1) * dxidx(:,1,2) + 200 & dxidx(:,2,1) * dxidx(:,2,2) + 201 & ( dxidx(:,1,1) + dxidx(:,2,1) ) * 202 & ( dxidx(:,1,2) + dxidx(:,2,2) ) ) * pt57 203 & + dxidx(:,3,1) * dxidx(:,3,2) 204 f3 = ( dxidx(:,1,2) * dxidx(:,1,2) + 205 & dxidx(:,2,2) * dxidx(:,2,2) + 206 & ( dxidx(:,1,2) + dxidx(:,2,2) )**2 ) * pt57 207 & + dxidx(:,3,2) * dxidx(:,3,2) 208 f4 = ( dxidx(:,1,1) * dxidx(:,1,3) + 209 & dxidx(:,2,1) * dxidx(:,2,3) + 210 & ( dxidx(:,1,1) + dxidx(:,2,1) ) * 211 & ( dxidx(:,1,3) + dxidx(:,2,3) ) ) * pt57 212 & + dxidx(:,3,1) * dxidx(:,3,3) 213 f5 = ( dxidx(:,1,2) * dxidx(:,1,3) + 214 & dxidx(:,2,2) * dxidx(:,2,3) + 215 & ( dxidx(:,1,2) + dxidx(:,2,2) ) * 216 & ( dxidx(:,1,3) + dxidx(:,2,3) ) ) * pt57 217 & + dxidx(:,3,2) * dxidx(:,3,3) 218 f6 = ( dxidx(:,1,3) * dxidx(:,1,3) + 219 & dxidx(:,2,3) * dxidx(:,2,3) + 220 & ( dxidx(:,1,3) + dxidx(:,2,3) )**2 ) * pt57 221 & + dxidx(:,3,3) * dxidx(:,3,3) 222c 223 ! flops = flops + 51*npro 224 endif 225c 226c.... calculate compact K matrix in local parent coordinates 227c.... equation 134 in "A new ... III" only w/ K^tilde_jj. Might need 228c.... complete Kij. 229 230 compK(:, 1) = f1 * T * rlm2mu + f3 * T * rmu 231 & + f6 * T * rmu 232c 233 compK(:, 2) = f2 * T * (rlm + rmu) 234 compK(:, 3) = f1 * T * rmu + f3 * T * rlm2mu 235 & + f6 * T * rmu 236c 237 compK(:, 4) = f4 * T * (rlm + rmu) 238 compK(:, 5) = f5 * T * (rlm + rmu) 239 compK(:, 6) = f1 * T * rmu + f3 * T * rmu 240 & + f6 * T * rlm2mu 241c 242 compK(:, 7) = f1 * T * rlm2mu * u1 + f2 * T * (rlm + rmu) * u2 243 & + f3 * T * rmu * u1 + f4 * T * (rlm + rmu) * u3 244 & + f6 * T * rmu * u1 245 compK(:, 8) = f1 * T * rmu * u2 + f2 * T * (rlm + rmu) * u1 246 & + f3 * T * rlm2mu * u2 + f5 * T * (rlm + rmu) * u3 247 & + f6 * T * rmu * u2 248 compK(:, 9) = f1 * T * rmu * u3 + f3 * T * rmu * u3 249 & + f4 * T * (rlm + rmu) * u1 + f5 * T * (rlm + rmu) * u2 250 & + f6 * T * rlm2mu * u3 251 252 rk=pt5*(u1**2+u2**2+u3**2) 253 254 compK(:,10) = f1 * T * (con * T + two * rmu * rk + (rlm + 255 & rmu) * u1**2) + f2 * T * (rlm + rmu) * two * u1 * u2 256 & + f3 * T * (con * T + two * rmu * rk + (rlm + rmu) * 257 & u2**2) + f4 * T * (rlm + rmu) * two * u1 * u3 258 & + f5 * T * (rlm + rmu) * two * u2 * u3 + f6 * T * (con 259 & * T + two * rmu * rk + (rlm + rmu) * u3**2) 260c 261c.... flop count 262c 263 ! flops = flops + 86*npro 264c 265c.... end of GLS 266c 267 268 endif 269c 270c.... ---------------------------> RHS <----------------------------- 271c 272c.... compute diffusive fluxes and add them to ri and rmi 273 274c 275c.... diffusive flux in x1-direction 276c 277c rmi(:,1) = zero ! already initialized 278 rmi(:,2) = rlm2mu * g1yi(:,2) 279 & + rlm * g2yi(:,3) 280 & + rlm * g3yi(:,4) 281 & - rlsli(:,1) 282 rmi(:,3) = rmu * g1yi(:,3) 283 & + rmu * g2yi(:,2) 284 & - rlsli(:,4) 285 rmi(:,4) = rmu * g1yi(:,4) 286 & + rmu * g3yi(:,2) 287 & - rlsli(:,5) 288 rmi(:,5) = rlm2mu * u1 * g1yi(:,2) + rmu * u2 * g1yi(:,3) 289 & + rmu * u3 * g1yi(:,4) 290 & + rmu * u2 * g2yi(:,2) + rlm * u1 * g2yi(:,3) 291 & + rmu * u3 * g3yi(:,2) + rlm * u1 * g3yi(:,4) 292 & + con * g1yi(:,5) 293 294c 295 ri (:,2:5) = ri (:,2:5) + rmi(:,2:5) 296c rmi(:,2:5) = rmi(:,2:5) + qdi(:,2:5) 297c 298c! flops = flops + 74*npro 299c 300c.... diffusive flux in x2-direction 301c 302c rmi(:, 6) = zero 303 rmi(:, 7) = rmu * g1yi(:,3) 304 & + rmu * g2yi(:,2) 305 & - rlsli(:,4) 306 rmi(:, 8) = rlm * g1yi(:,2) 307 & + rlm2mu * g2yi(:,3) 308 & + rlm * g3yi(:,4) 309 & - rlsli(:,2) 310 rmi(:, 9) = rmu * g2yi(:,4) 311 & + rmu * g3yi(:,3) 312 & - rlsli(:,6) 313 rmi(:,10) = rlm * u2 * g1yi(:,2) + rmu * u1 * g1yi(:,3) 314 & + rmu * u1 * g2yi(:,2) + rlm2mu * u2 * g2yi(:,3) 315 & + rmu * u3 * g2yi(:,4) 316 & + rmu * u3 * g3yi(:,3) + rlm * u2 * g3yi(:,4) 317 & + con * g2yi(:,5) 318c 319 ri (:,7:10) = ri (:,7:10) + rmi(:,7:10) 320c rmi(:,7:10) = rmi(:,7:10) + qdi(:,2:5) 321c 322c! flops = flops + 74*npro 323c 324c.... diffusive flux in x3-direction 325c 326c rmi(:,11) = zero 327 rmi(:,12) = rmu * g1yi(:,4) 328 & + rmu * g3yi(:,2) 329 & - rlsli(:,5) 330 rmi(:,13) = rmu * g2yi(:,4) 331 & + rmu * g3yi(:,3) 332 & - rlsli(:,6) 333 rmi(:,14) = rlm * g1yi(:,2) 334 & + rlm * g2yi(:,3) 335 & + rlm2mu * g3yi(:,4) 336 & - rlsli(:,3) 337 rmi(:,15) = rlm * u3 * g1yi(:,2) + rmu * u1 * g1yi(:,4) 338 & + rlm * u3 * g2yi(:,3) + rmu * u2 * g2yi(:,4) 339 & + rmu * u1 * g3yi(:,2) + rmu * u2 * g3yi(:,3) 340 & + rlm2mu * u3 * g3yi(:,4) 341 & + con * g3yi(:,5) 342c 343 ri (:,12:15) = ri (:,12:15) + rmi(:,12:15) 344c 345c! flops = flops + 74*npro 346c 347c stiff for preconditioner has been eliminated 348c all preconditioner stuff is in e3bdg.f 349c 350 351 ttim(23) = ttim(23) + secs(0.0) 352 353c 354c.... return 355c 356 return 357 end 358c 359c 360c 361 subroutine e3viscSclr (g1yti, g2yti, g3yti, 362 & rmu, Sclr, rho, 363 & rti, rmti, stifft) 364c 365c---------------------------------------------------------------------- 366c 367c This routine calculates the contribution of viscous and heat fluxes 368c to both RHS and LHS. 369c 370c input: 371c g1yti (npro) : grad-y in direction 1 372c g2yti (npro) : grad-y in direction 2 373c g3yti (npro) : grad-y in direction 3 374c rmu (npro) : viscosity 375c Sclr (npro) : scalar variable 376c 377c output: 378c rti (npro,nsd+1) : partial residual 379c rmti (npro,nsd+1) : partial modified residual 380c stifft (npro,nsd,nsd) : stiffness matrix 381c 382c 383c 384c Zdenek Johan, Summer 1990. (Modified from e2visc.f) 385c Zdenek Johan, Winter 1991. (Fortran 90) 386c Kenneth Jansen, Winter 1997 Primitive Variables 387c---------------------------------------------------------------------- 388c 389 use turbSA ! for saSigma 390 include "common.h" 391c 392c passed arrays 393c 394 dimension g1yti(npro), g2yti(npro), 395 & g3yti(npro), rmu(npro), 396 & Sclr(npro), rho(npro), 397 & rti(npro,nsd+1), rmti(npro,nsd+1), 398 & stifft(npro,3,3) 399 400 ttim(23) = ttim(23) - tmr() 401 402 if ((ilset.ne.zero) .and. (isclr.lt.3)) return 403c 404c.... ---------------------------> RHS <----------------------------- 405c 406c.... ---------------------> Diffusivity Matrix <------------------- 407c 408c if (lhs .eq. 1) then 409 410 stifft = zero 411c 412c.... K11 413c 414 stifft(:,1,1)=rmu 415 if (iRANS .lt. 0) then 416 stifft(:,1,1)=saSigmaInv*((rmu/rho)+max(zero,Sclr)) 417!Sclr is nu_til not mu and thus nu+nu_til is the total diffusion 418 endif 419c 420c.... K22 421c 422 stifft(:,2,2)=stifft(:,1,1) 423c 424c.... K33 425c 426 stifft(:,3,3)=stifft(:,1,1) 427 428c endif 429c 430c.... ---------------------------> RHS <----------------------------- 431c 432c.... compute diffusive fluxes and add them to ri and rmi 433c 434c.... diffusive fluxes 435c 436 rmti(:,1) = stifft(:,1,1) * g1yti(:) 437 rmti(:,2) = stifft(:,2,2) * g2yti(:) 438 rmti(:,3) = stifft(:,3,3) * g3yti(:) 439 440 rti (:,:) = rti (:,:) + rmti(:,:) 441c rmi(:,2:5) = rmi(:,2:5) + qdi(:,2:5) 442c 443c! flops = flops + 74*npro 444c 445 ttim(23) = ttim(23) + tmr() 446 447c 448c.... return 449c 450 return 451 end 452