1 subroutine e3tau (rho, cp, rmu, 2 & u1, u2, u3, 3 & con, dxidx, rLyi, 4 & rLymi, tau, rk, 5 & giju, rTLS, raLS, 6 & A0inv, dVdY, cv) 7c 8c---------------------------------------------------------------------- 9c 10c This routine computes the diagonal Tau for least-squares operator. 11c We receive the regular residual L operator and a 12c modified residual L operator, calculate tau, and return values for 13c tau and tau times these operators to maintain the format of past 14c ENSA codes 15c 16c input: 17c rho (npro) : density 18c T (npro) : temperature 19c cp (npro) : specific heat at constant pressure 20c u1 (npro) : x1-velocity component 21c u2 (npro) : x2-velocity component 22c u3 (npro) : x3-velocity component 23c dxidx (npro,nsd,nsd) : inverse of deformation gradient 24c rLyi (npro,nflow) : least-squares residual vector 25c rLymi (npro,nflow) : modified least-squares residual vector 26c 27c output: 28c rLyi (npro,nflow) : least-squares residual vector 29c rLymi (npro,nflow) : modified least-squares residual vector 30c tau (npro,3) : 3 taus 31c 32c 33c Zdenek Johan, Summer 1990. (Modified from e2tau.f) 34c Zdenek Johan, Winter 1991. (Fortran 90) 35c---------------------------------------------------------------------- 36c 37 include "common.h" 38c 39 dimension rho(npro), con(npro), 40 & cp(npro), u1(npro), 41 & u2(npro), u3(npro), 42 & dxidx(npro,nsd,nsd), rk(npro), 43 & tau(npro,5), rLyi(npro,nflow), 44 & rLymi(npro,nflow), dVdY(npro,15), 45 & rTLS(npro), raLS(npro), 46 & rLyitemp(npro,nflow), detgijI(npro) 47c 48 dimension rmu(npro), cv(npro), 49 & gijd(npro,6), uh1(npro), 50 & fact(npro), h2o2u(npro), giju(npro,6), 51 & A0inv(npro,15),gijdu(npro,6) 52c 53 call e3gijd( dxidx, gijd ) ! both diagonal taus need this 54 if(itau.eq.1) then ! tau proposed by for nearly incompressible 55c flows by Guillermo Hauke 56c 57c next form the diffusive length scale |u| h_1 = 2 ( ui (gijd)-1 uj)^{1/2} 58c 59c dividing factor for the inverse of gijd 60c 61 fact = gijd(:,1) * gijd(:,3) * gijd(:,6) 62 & - gijd(:,1) * gijd(:,5) * gijd(:,5) 63 & - gijd(:,3) * gijd(:,4) * gijd(:,4) 64 & - gijd(:,6) * gijd(:,2) * gijd(:,2) 65 & + gijd(:,2) * gijd(:,4) * gijd(:,5) * two 66c 67 uh1= u1*u1*(gijd(:,3)*gijd(:,6)-gijd(:,5)*gijd(:,5)) 68 & + u2*u2*(gijd(:,1)*gijd(:,6)-gijd(:,4)*gijd(:,4)) 69 & + u3*u3*(gijd(:,1)*gijd(:,3)-gijd(:,2)*gijd(:,2)) 70 & + two *(u1*u2*(gijd(:,4)*gijd(:,5)-gijd(:,2)*gijd(:,6)) 71 & + u1*u3*(gijd(:,2)*gijd(:,5)-gijd(:,4)*gijd(:,3)) 72 & + u2*u3*(gijd(:,4)*gijd(:,2)-gijd(:,1)*gijd(:,5))) 73c 74c at this point we have (u h1)^2 * inverse coefficient^2 / 4 75c 76c ^ fact 77c 78 uh1= two * sqrt(uh1/fact) 79c 80c next form the advective length scale |u|/h_2 = 2 ( ui (gijd) uj)^{1/2} 81c 82 h2o2u = u1*u1*gijd(:,1) 83 & + u2*u2*gijd(:,3) 84 & + u3*u3*gijd(:,6) 85 & +(u1*u2*gijd(:,2) 86 & + u1*u3*gijd(:,4) 87 & + u2*u3*gijd(:,5))*two + 1.0e-15 !FIX FOR INVALID MESHES 88c 89c at this point we have (2 u / h_2)^2 90c 91 92c call tnanqe(h2o2u,1,"riaconv ") 93 94 h2o2u = one / sqrt(h2o2u) ! this flips it over leaves it h_2/(2u) 95c 96c.... diffusive corrections 97 98c 99c.... cell Reynold number 100c 101 fact=pt5*rho*uh1/rmu 102 103c 104c.... continuity tau 105c 106 tau(:,1)=pt5*uh1*min(one,fact)*taucfct 107c 108c... momentum tau 109c 110 if (iremoveStabTimeTerm.eq.0 ) then 111 dts=one/(Dtgl*dtsfct) 112 tau(:,2)=min(dts,h2o2u) 113 endif 114 tau(:,2)=tau(:,2)/rho 115c 116c.... energy tau cv=cp/gamma assumed 117c 118c tau(:,3)=gamma*tau(:,2)/cp 119 tau(:,3)=tau(:,2)/cv 120c 121c.... diffusive corrections 122c 123 if (ipord == 1) then 124 celt = pt66 125 else if (ipord == 2) then 126 celt = pt33 127c celt = pt33*0.04762 128 else if (ipord == 3) then 129 celt = pt33 !.02 just a guess... 130 else if (ipord >= 4) then 131 celt = .008 ! yet another guess... 132 endif 133c 134c fact=h2o2u*h2o2u*rk*pt66/rmu 135 fact=h2o2u*h2o2u*rk*celt/rmu 136c 137 tau(:,2)=min(tau(:,2),fact) 138 tau(:,3)=min(tau(:,3),fact*rmu/con)*temper 139c 140 else if(itau.eq.0) then ! tau proposed by Farzin and Shakib 141c 142c... momentum tau 143c 144 145c 146c... higher order element diffusive correction 147c 148 if (ipord == 1) then 149 fff = 36.0d0 150 else if (ipord == 2) then 151 fff = 60.0d0 152c fff = 36.0d0 153 else if (ipord == 3) then 154 fff = 128.0d0 155c fff = 144.0d0 156 endif 157 if (iremoveStabTimeTerm.eq.1 ) then 158 dts = zero 159 else 160 dts = dtsfct*Dtgl 161 endif 162 tau(:,2)=rho**2*((two*dts)**2 163 & + (u1*(u1*gijd(:,1) + two*(u2*gijd(:,2)+u3*gijd(:,4))) 164 & + u2*(u2*gijd(:,3) + two*u3*gijd(:,5)) 165 & + u3*u3*gijd(:,6))) 166 & +fff*rmu**2*(gijd(:,1)**2 + gijd(:,3)**2 + gijd(:,6)**2 + 167 & two*(gijd(:,2)**2 + gijd(:,4)**2 + gijd(:,5)**2)) 168 fact=sqrt(tau(:,2)) 169 tau(:,1)=pt125*fact/(rho*(gijd(:,1)+gijd(:,3)+gijd(:,6)))*taucfct 170 tau(:,2)=one/fact 171c 172c.... energy tau cv=cp/gamma assumed 173c 174 tau(:,3)=tau(:,2)/cv*temper 175c 176 endif 177c 178c... finally multiply this tau matrix times the 179c two residual vectors 180c 181c ... calculate (tau Ly) --> (rLyi) 182c ... storing rLyi for the DC term 183 if(iDC .ne. 0) rLyitemp=rLyi 184 185 if(ires.eq.3 .or. ires .eq. 1) then 186 rLyi(:,1) = tau(:,1) * rLyi(:,1) 187 rLyi(:,2) = tau(:,2) * rLyi(:,2) 188 rLyi(:,3) = tau(:,2) * rLyi(:,3) 189 rLyi(:,4) = tau(:,2) * rLyi(:,4) 190 rLyi(:,5) = tau(:,3) * rLyi(:,5) 191 endif 192c 193 if(iDC .ne. 0) then 194c..... calculation of rTLS & raLS for DC term 195c 196c..... calculation of (Ly - S).tau^tilda*(Ly - S) 197c 198 rTLS = rLYItemp(:,1) * (rLyi(:,1)*dVdY(:,1) 199 & + dVdY(:,2)*rLyi(:,2) + dVdY(:,4)*rLyi(:,3) 200 & + rLyi(:,4)*dVdY(:,7) + dVdY(:,11)*rLyi(:,5)) 201 & + rLyitemp(:,2) * (rLyi(:,2)*dVdY(:,3) 202 & + rLyi(:,3)*dVdY(:,5) + dVdY(:,8)*rLyi(:,4) 203 & + rLyi(:,5)*dVdY(:,12)) 204 & + rLyitemp(:,3) * (rLyi(:,3)*dVdY(:,6) 205 & + dVdY(:,9)*rLyi(:,4) + dVdY(:,13)*rLyi(:,5)) 206 & + rLyitemp(:,4) * (rLyi(:,4)*dVdY(:,10) 207 & + dVdY(:,14)*rLyi(:,5)) 208 & + rLyitemp(:,5) * (dVdY(:,15)*rLyi(:,5)) 209 210c 211c...... calculation of (Ly - S).A0inv*(Ly - S) 212c 213 raLS = two*rLyitemp(:,4)*rLyitemp(:,5)*A0inv(:,15) 214 & + two*rLyitemp(:,3)*rLyitemp(:,5)*A0inv(:,14) 215 & + two*rLyitemp(:,1)*rLyitemp(:,2)*A0inv( :,6) 216 & + two*rLyitemp(:,2)*rLyitemp(:,3)*A0inv(:,10) 217 & + two*rLyitemp(:,2)*rLyitemp(:,4)*A0inv(:,11) 218 & + two*rLyitemp(:,1)*rLyitemp(:,3)*A0inv( :,7) 219 & + two*rLyitemp(:,3)*rLyitemp(:,4)*A0inv(:,13) 220 & + two*rLyitemp(:,2)*rLyitemp(:,5)*A0inv(:,12) 221 & + two*rLyitemp(:,1)*rLyitemp(:,4)*A0inv( :,8) 222 & + two*rLyitemp(:,1)*rLyitemp(:,5)*A0inv( :,9) 223 & + rLyitemp(:,1)**2*A0inv(:,1) 224 & + rLyitemp(:,2)**2*A0inv(:,2) 225 & + rLyitemp(:,3)**2*A0inv(:,3) 226 & + rLyitemp(:,4)**2*A0inv(:,4) 227 & + rLyitemp(:,5)**2*A0inv(:,5) 228c 229c.....****************calculation of giju for DC term*************** 230c 231c.... for the notation of different numbering 232c 233 gijdu(:,1)=gijd(:,1) 234 gijdu(:,2)=gijd(:,3) 235 gijdu(:,3)=gijd(:,6) 236 gijdu(:,4)=gijd(:,2) 237 gijdu(:,5)=gijd(:,4) 238 gijdu(:,6)=gijd(:,5) 239c 240c 241 detgijI = one/( 242 & gijdu(:,1) * gijdu(:,2) * gijdu(:,3) 243 & - gijdu(:,1) * gijdu(:,6) * gijdu(:,6) 244 & - gijdu(:,4) * gijdu(:,4) * gijdu(:,3) 245 & + gijdu(:,4) * gijdu(:,5) * gijdu(:,6) * two 246 & - gijdu(:,5) * gijdu(:,5) * gijdu(:,2) 247 & ) 248 giju(:,1) = detgijI * (gijdu(:,2)*gijdu(:,3) 249 & - gijdu(:,6)**2) 250 giju(:,2) = detgijI * (gijdu(:,1)*gijdu(:,3) 251 & - gijdu(:,5)**2) 252 giju(:,3) = detgijI * (gijdu(:,1)*gijdu(:,2) 253 & - gijdu(:,4)**2) 254 giju(:,4) = detgijI * (gijdu(:,5)*gijdu(:,6) 255 & - gijdu(:,4)*gijdu(:,3) ) 256 giju(:,5) = detgijI * (gijdu(:,4)*gijdu(:,6) 257 & - gijdu(:,5)*gijdu(:,2) ) 258 giju(:,6) = detgijI * (gijdu(:,4)*gijdu(:,5) 259 & - gijdu(:,1)*gijdu(:,6) ) 260c 261 endif ! end of iDC.ne.0 262c 263c.... calculate (tau Lym) --> (rLymi) 264c 265 if(ires.ne.1 ) then 266 rLymi(:,1) = tau(:,1) * rLymi(:,1) 267 rLymi(:,2) = tau(:,2) * rLymi(:,2) 268 rLymi(:,3) = tau(:,2) * rLymi(:,3) 269 rLymi(:,4) = tau(:,2) * rLymi(:,4) 270 rLymi(:,5) = tau(:,3) * rLymi(:,5) 271 endif 272c 273c INCORRECT NOW ! flops = flops + 255*npro 274c 275c 276c.... return 277c 278 return 279 end 280c 281c 282 283 284 subroutine e3tau_nd (rho, cp, rmu, T, 285 & u1, u2, u3, 286 & a1, a2, a3, 287 & con, dxidx, rLyi, 288 & rLymi, Tau, rk, 289 & giju, rTLS, raLS, 290 & cv, compK, pres, 291 & A0inv, dVdY) 292c 293c---------------------------------------------------------------------- 294c 295c This routine computes the diagonal Tau for least-squares operator. 296c We receive the regular residual L operator and a 297c modified residual L operator, calculate tau, and return values for 298c tau and tau times these operators to maintain the format of past 299c ENSA codes 300c 301c input: 302c rho (npro) : density 303c T (npro) : temperature 304c cp (npro) : specific heat at constant pressure 305c u1 (npro) : x1-velocity component 306c u2 (npro) : x2-velocity component 307c u3 (npro) : x3-velocity component 308c dxidx (npro,nsd,nsd) : inverse of deformation gradient 309c rLyi (npro,nflow) : least-squares residual vector 310c rLymi (npro,nflow) : modified least-squares residual vector 311c 312c output: 313c rLyi (npro,nflow) : least-squares residual vector 314c rLymi (npro,nflow) : modified least-squares residual vector 315c Mtau (npro,5,5) : Matrix of Stabilized Parameters 316c 317c 318c Zdenek Johan, Summer 1990. (Modified from e2tau.f) 319c Zdenek Johan, Winter 1991. (Fortran 90) 320c---------------------------------------------------------------------- 321c 322 include "common.h" 323c 324 dimension rho(npro), con(npro), 325 & cp(npro), u1(npro), 326 & u2(npro), u3(npro), 327 & dxidx(npro,nsd,nsd), rk(npro), 328 & rLyi(npro,nflow), 329 & rLymi(npro,nflow), dVdY(npro,15), 330 & rTLS(npro), raLS(npro), 331 & rLyitemp(npro,nflow), detgijI(npro) 332c 333 dimension rmu(npro), cv(npro), 334 & gijd(npro,6), 335 & fact(npro), giju(npro,6), 336 & A0inv(npro,15),gijdu(npro,6), compK(npro,10), 337 & a1(npro), a2(npro), a3(npro), 338 & T(npro), pres(npro), alphap(npro), 339 & betaT(npro), gamb(npro), c(npro), 340 & u(npro), eb1(npro), Q(npro,5,5), 341 & rlam(npro,5), eigmax(npro,5), Pe(npro), 342 & Ak(npro), B(npro), D(npro), E(npro), 343 & STau(npro,15), Tau(npro,nflow,nflow) 344 345 346c... build some necessary quantities at integration point: 347 348c alfap = one / T 349c betaT = one / pres 350 gamb = gamma1 351 c = sqrt( (gamma * Rgas) * T ) 352 u = sqrt(u1**2+u2**2+u3**2) 353 eb1 = cp*T - pt5*(u1**2+u2**2+u3**2) 354 355c... get eigenvectors for diagonalizing Tau_adv (e.v) before final rotations 356c... Note: the ordering of eigenvectors corresponds to the following 357c... ordering of the eigenvalues in the 1-st Euler Jacobian rotated to 358c... the streamline coordinates 359 360c |u+c | 361c | u | 362c | u | 363c | u | ! for origins of this see Beam, Warming, Hyett 364c | u-c| ! Mathematics of Computation vol. 29 No. 132 p. 1037 365 366c... different ordering assumed in Shakib/Johan entropy vars. code 367 368 369 370 call e3eig1 (rho, T, cp, 371 & gamb, c, u1, 372 & u2, u3, a1, 373 & a2, a3, eb1, 374 & dxidx, u, Q) 375 376c... Perform a Jacobi rotation on the Lambda matrix as well as adjust 377c... the eigenvectors. Tau still remains in entropy variables. 378 379 380 381 call e3eig2 (u, c, a1, 382 & a2, a3, rlam, 383 & Q, eigmax) 384 385c 386c.... invert the eigenvalues 387c 388 389 390 where (rlam .gt. ((epsM**2) * eigmax)) 391 rlam = one / sqrt(rlam) 392 elsewhere 393 rlam = zero 394 endwhere 395 396c 397c.... flop count 398c 399 ! flops = flops + 65*npro 400 401c.... reduce the diffusion contribution 402c 403 404 if (Navier .eq. 1) then 405c 406c.... calculate sigma 407c 408 409 do i = 1, nflow ! diff. corr for every mode of Tau 410 411 c(:) = pt33 * ( 412 & Q(:,2,i) * ( compK(:, 1) * Q(:,2,i) + compK(:, 2) * Q(:,3,i) 413 & + compK(:, 4) * Q(:,4,i) + compK(:, 7) * Q(:,5,i) ) 414 & + Q(:,3,i) * ( compK(:, 2) * Q(:,2,i) + compK(:, 3) * Q(:,3,i) 415 & + compK(:, 5) * Q(:,4,i) + compK(:, 8) * Q(:,5,i) ) 416 & + Q(:,4,i) * ( compK(:, 4) * Q(:,2,i) + compK(:, 5) * Q(:,3,i) 417 & + compK(:, 6) * Q(:,4,i) + compK(:, 9) * Q(:,5,i) ) 418 & + Q(:,5,i) * ( compK(:, 7) * Q(:,2,i) + compK(:, 8) * Q(:,3,i) 419 & + compK(:, 9) * Q(:,4,i) + compK(:,10) * Q(:,5,i) ) 420 & ) 421 422c... get Peclet Number 423 424 425 Pe(:) = one / (c(:)*rlam(:,i)) 426 427 428c 429c... branch out into different polynomial orders 430c 431 432 433 if (ipord == 1) then ! linears - l = l^a*(coth(Pe)-1/Pe) 434 ! approx. l = l^a*min(1/3*Pe,1) 435 rlam(:,i) = rlam(:,i)*min(pt33*Pe(:),one) 436 437 endif 438 439 if (ipord == 2) then ! quads - Codina, CMAME' 92 440 ! approx. l = l^a*min(1/6*Pe,1/2) 441 rlam(:,i) = rlam(:,i)*min(pt33*pt5*Pe(:),pt5) 442 443 endif 444 445 if (ipord == 3) then ! cubics - Recent Work 446 ! l = l^a*1/Pe*b 447 ! b is a real root of the 448 ! following polynomial: 449 ! b^3+(-Pe*coth(Pe)b^2+(6/15*Pe^2-Pe*coth(Pe)+1)b+ 450 ! +(-1/15*Pe^3*coth(Pe)+6/15*Pe^2-Pe*coth(Pe)+1) = 0 451 ! proceed setting up newton iteration w/ initial 452 ! guess coming from the quad estimate. 453 454 Ak(:)=1.0 ! get parameters for iteration 455 456 where(Pe.lt.5.0) ! approx. to hyp. cothangent 457 alphap = exp(Pe) 458 betaT = exp(-Pe) 459 c = (alphap+betaT)/(alphap-betaT) 460 elsewhere 461 c = one 462 endwhere 463 464 B= -Pe*c + Ak 465 D= 0.4 * (Pe**2) + B 466 E=-0.066666667 * (Pe**3) * c +D 467 468 ! initial guess, use betaT 469 betaT(:) = Pe(:)*min(pt33*pt5*Pe(:),pt5) 470 471 ! iteration - 3 iterations sufficient 472 do j = 1, 3 473 474 betaT=betaT-(Ak*betaT**3+B*betaT**2+D*betaT+E)/(3*Ak 475 & *betaT**2+2*B*betaT+D) 476 enddo 477 478 rlam(:,i) = rlam(:,i)*(1/Pe(:))*betaT(:) 479 480 endif ! done w/ different polynomial orders 481 482 enddo ! done with loop over dof's 483 484 endif ! done with diffusive correction 485 486 487c 488c.... form Tau (symmetric - entropy variables - then convert) 489c 490 STau(:, 1) = rlam(:,1) * Q(:,1,1) * Q(:,1,1) + 491 & rlam(:,2) * Q(:,1,2) * Q(:,1,2) + 492 & rlam(:,3) * Q(:,1,3) * Q(:,1,3) + 493 & rlam(:,4) * Q(:,1,4) * Q(:,1,4) + 494 & rlam(:,5) * Q(:,1,5) * Q(:,1,5) 495c 496 STau(:, 2) = rlam(:,1) * Q(:,1,1) * Q(:,2,1) + 497 & rlam(:,2) * Q(:,1,2) * Q(:,2,2) + 498 & rlam(:,3) * Q(:,1,3) * Q(:,2,3) + 499 & rlam(:,4) * Q(:,1,4) * Q(:,2,4) + 500 & rlam(:,5) * Q(:,1,5) * Q(:,2,5) 501c 502 STau(:, 3) = rlam(:,1) * Q(:,2,1) * Q(:,2,1) + 503 & rlam(:,2) * Q(:,2,2) * Q(:,2,2) + 504 & rlam(:,3) * Q(:,2,3) * Q(:,2,3) + 505 & rlam(:,4) * Q(:,2,4) * Q(:,2,4) + 506 & rlam(:,5) * Q(:,2,5) * Q(:,2,5) 507c 508 STau(:, 4) = rlam(:,1) * Q(:,1,1) * Q(:,3,1) + 509 & rlam(:,2) * Q(:,1,2) * Q(:,3,2) + 510 & rlam(:,3) * Q(:,1,3) * Q(:,3,3) + 511 & rlam(:,4) * Q(:,1,4) * Q(:,3,4) + 512 & rlam(:,5) * Q(:,1,5) * Q(:,3,5) 513c 514 STau(:, 5) = rlam(:,1) * Q(:,2,1) * Q(:,3,1) + 515 & rlam(:,2) * Q(:,2,2) * Q(:,3,2) + 516 & rlam(:,3) * Q(:,2,3) * Q(:,3,3) + 517 & rlam(:,4) * Q(:,2,4) * Q(:,3,4) + 518 & rlam(:,5) * Q(:,2,5) * Q(:,3,5) 519c 520 STau(:, 6) = rlam(:,1) * Q(:,3,1) * Q(:,3,1) + 521 & rlam(:,2) * Q(:,3,2) * Q(:,3,2) + 522 & rlam(:,3) * Q(:,3,3) * Q(:,3,3) + 523 & rlam(:,4) * Q(:,3,4) * Q(:,3,4) + 524 & rlam(:,5) * Q(:,3,5) * Q(:,3,5) 525c 526 STau(:, 7) = rlam(:,1) * Q(:,1,1) * Q(:,4,1) + 527 & rlam(:,2) * Q(:,1,2) * Q(:,4,2) + 528 & rlam(:,3) * Q(:,1,3) * Q(:,4,3) + 529 & rlam(:,4) * Q(:,1,4) * Q(:,4,4) + 530 & rlam(:,5) * Q(:,1,5) * Q(:,4,5) 531c 532 STau(:, 8) = rlam(:,1) * Q(:,2,1) * Q(:,4,1) + 533 & rlam(:,2) * Q(:,2,2) * Q(:,4,2) + 534 & rlam(:,3) * Q(:,2,3) * Q(:,4,3) + 535 & rlam(:,4) * Q(:,2,4) * Q(:,4,4) + 536 & rlam(:,5) * Q(:,2,5) * Q(:,4,5) 537c 538 STau(:, 9) = rlam(:,1) * Q(:,3,1) * Q(:,4,1) + 539 & rlam(:,2) * Q(:,3,2) * Q(:,4,2) + 540 & rlam(:,3) * Q(:,3,3) * Q(:,4,3) + 541 & rlam(:,4) * Q(:,3,4) * Q(:,4,4) + 542 & rlam(:,5) * Q(:,3,5) * Q(:,4,5) 543c 544 STau(:,10) = rlam(:,1) * Q(:,4,1) * Q(:,4,1) + 545 & rlam(:,2) * Q(:,4,2) * Q(:,4,2) + 546 & rlam(:,3) * Q(:,4,3) * Q(:,4,3) + 547 & rlam(:,4) * Q(:,4,4) * Q(:,4,4) + 548 & rlam(:,5) * Q(:,4,5) * Q(:,4,5) 549c 550 STau(:,11) = rlam(:,1) * Q(:,1,1) * Q(:,5,1) + 551 & rlam(:,2) * Q(:,1,2) * Q(:,5,2) + 552 & rlam(:,3) * Q(:,1,3) * Q(:,5,3) + 553 & rlam(:,4) * Q(:,1,4) * Q(:,5,4) + 554 & rlam(:,5) * Q(:,1,5) * Q(:,5,5) 555c 556 STau(:,12) = rlam(:,1) * Q(:,2,1) * Q(:,5,1) + 557 & rlam(:,2) * Q(:,2,2) * Q(:,5,2) + 558 & rlam(:,3) * Q(:,2,3) * Q(:,5,3) + 559 & rlam(:,4) * Q(:,2,4) * Q(:,5,4) + 560 & rlam(:,5) * Q(:,2,5) * Q(:,5,5) 561c 562 STau(:,13) = rlam(:,1) * Q(:,3,1) * Q(:,5,1) + 563 & rlam(:,2) * Q(:,3,2) * Q(:,5,2) + 564 & rlam(:,3) * Q(:,3,3) * Q(:,5,3) + 565 & rlam(:,4) * Q(:,3,4) * Q(:,5,4) + 566 & rlam(:,5) * Q(:,3,5) * Q(:,5,5) 567c 568 STau(:,14) = rlam(:,1) * Q(:,4,1) * Q(:,5,1) + 569 & rlam(:,2) * Q(:,4,2) * Q(:,5,2) + 570 & rlam(:,3) * Q(:,4,3) * Q(:,5,3) + 571 & rlam(:,4) * Q(:,4,4) * Q(:,5,4) + 572 & rlam(:,5) * Q(:,4,5) * Q(:,5,5) 573c 574 STau(:,15) = rlam(:,1) * Q(:,5,1) * Q(:,5,1) + 575 & rlam(:,2) * Q(:,5,2) * Q(:,5,2) + 576 & rlam(:,3) * Q(:,5,3) * Q(:,5,3) + 577 & rlam(:,4) * Q(:,5,4) * Q(:,5,4) + 578 & rlam(:,5) * Q(:,5,5) * Q(:,5,5) 579 580 581c 582c... Form Primitive Variable Tau as [dY/dV]*tilde{Tau}, 583c... see G.Hauke's thesis appendix for [dY/dV] matrix 584c 585 betaT = cp*T + pt5*(u1**2+u2**2+u3**2) !reuse betaT 586 587 Tau(:,1,1) = rho*T*(STau(:,1)+u1*STau(:,2)+ 588 & u2*STau(:,4)+u3*STau(:,7)+betaT*STau(:,11)) 589 Tau(:,1,2) = rho*T*(STau(:,2)+u1*STau(:,3)+ 590 & u2*STau(:,5)+u3*STau(:,8)+betaT*STau(:,12)) 591 Tau(:,1,3) = rho*T*(STau(:,4)+u1*STau(:,5)+ 592 & u2*STau(:,6)+u3*STau(:,9)+betaT*STau(:,13)) 593 Tau(:,1,4) = rho*T*(STau(:,7)+u1*STau(:,8)+ 594 & u2*STau(:,9)+u3*STau(:,10)+betaT*STau(:,14)) 595 Tau(:,1,5) = rho*T*(STau(:,11)+u1*STau(:,12)+ 596 & u2*STau(:,13)+u3*STau(:,14)+betaT*STau(:,15)) 597 598 599 Tau(:,2,1) = T(:)*(STau(:,2) + u1(:)*STau(:,11)) 600 Tau(:,2,2) = T(:)*(STau(:,3) + u1(:)*STau(:,12)) 601 Tau(:,2,3) = T(:)*(STau(:,5) + u1(:)*STau(:,13)) 602 Tau(:,2,4) = T(:)*(STau(:,8) + u1(:)*STau(:,14)) 603 Tau(:,2,5) = T(:)*(STau(:,12) + u1(:)*STau(:,15)) 604 605 Tau(:,3,1) = T(:)*(STau(:,4) + u2(:)*STau(:,11)) 606 Tau(:,3,2) = T(:)*(STau(:,5) + u2(:)*STau(:,12)) 607 Tau(:,3,3) = T(:)*(STau(:,6) + u2(:)*STau(:,13)) 608 Tau(:,3,4) = T(:)*(STau(:,9) + u2(:)*STau(:,14)) 609 Tau(:,3,5) = T(:)*(STau(:,13) + u2(:)*STau(:,15)) 610 611 Tau(:,4,1) = T(:)*(STau(:,7) + u3(:)*STau(:,11)) 612 Tau(:,4,2) = T(:)*(STau(:,8) + u3(:)*STau(:,12)) 613 Tau(:,4,3) = T(:)*(STau(:,9) + u3(:)*STau(:,13)) 614 Tau(:,4,4) = T(:)*(STau(:,10) + u3(:)*STau(:,14)) 615 Tau(:,4,5) = T(:)*(STau(:,14) + u3(:)*STau(:,15)) 616 617 betaT = T**2 618 619 Tau(:,5,1) = betaT(:)*STau(:,11) 620 Tau(:,5,2) = betaT(:)*STau(:,12) 621 Tau(:,5,3) = betaT(:)*STau(:,13) 622 Tau(:,5,4) = betaT(:)*STau(:,14) 623 Tau(:,5,5) = betaT(:)*STau(:,15) 624 625 626c 627c... done with conversion to pressure primitive variables 628c... now need to interface with the rest of the computations 629c 630 631c... finally multiply this tau matrix times the 632c two residual vectors 633c 634c ... calculate (tau Ly) --> (rLyi) 635c ... storing rLyi for the DC term 636 637 638 if(iDC .ne. 0) rLyitemp=rLyi 639 640 if(ires.eq.3 .or. ires .eq. 1) then 641 eigmax = rLyi !reuse 642 rLyi = zero 643 do i = 1, nflow 644 do j = 1, nflow 645 rLyi(:,i) = rLyi(:,i) + Tau(:,i,j) * eigmax(:,j) 646 enddo 647 enddo 648 endif 649 650 651 if(iDC .ne. 0) then 652c.....calculation of rTLS & raLS for DC term 653c 654c.....calculation of (Ly - S).tau^tilda*(Ly - S) 655c 656 rTLS = rLYItemp(:,1) * (rLyi(:,1)*dVdY(:,1) 657 & + dVdY(:,2)*rLyi(:,2) + dVdY(:,4)*rLyi(:,3) 658 & + rLyi(:,4)*dVdY(:,7) + dVdY(:,11)*rLyi(:,5)) 659 & + rLyitemp(:,2) * (rLyi(:,2)*dVdY(:,3) 660 & + rLyi(:,3)*dVdY(:,5) + dVdY(:,8)*rLyi(:,4) 661 & + rLyi(:,5)*dVdY(:,12)) 662 & + rLyitemp(:,3) * (rLyi(:,3)*dVdY(:,6) 663 & + dVdY(:,9)*rLyi(:,4) + dVdY(:,13)*rLyi(:,5)) 664 & + rLyitemp(:,4) * (rLyi(:,4)*dVdY(:,10) 665 & + dVdY(:,14)*rLyi(:,5)) 666 & + rLyitemp(:,5) * (dVdY(:,15)*rLyi(:,5)) 667 668c 669c...... calculation of (Ly - S).A0inv*(Ly - S) 670c 671 raLS = two*rLyitemp(:,4)*rLyitemp(:,5)*A0inv(:,15) 672 & + two*rLyitemp(:,3)*rLyitemp(:,5)*A0inv(:,14) 673 & + two*rLyitemp(:,1)*rLyitemp(:,2)*A0inv( :,6) 674 & + two*rLyitemp(:,2)*rLyitemp(:,3)*A0inv(:,10) 675 & + two*rLyitemp(:,2)*rLyitemp(:,4)*A0inv(:,11) 676 & + two*rLyitemp(:,1)*rLyitemp(:,3)*A0inv( :,7) 677 & + two*rLyitemp(:,3)*rLyitemp(:,4)*A0inv(:,13) 678 & + two*rLyitemp(:,2)*rLyitemp(:,5)*A0inv(:,12) 679 & + two*rLyitemp(:,1)*rLyitemp(:,4)*A0inv( :,8) 680 & + two*rLyitemp(:,1)*rLyitemp(:,5)*A0inv( :,9) 681 & + rLyitemp(:,1)**2*A0inv(:,1) 682 & + rLyitemp(:,2)**2*A0inv(:,2) 683 & + rLyitemp(:,3)**2*A0inv(:,3) 684 & + rLyitemp(:,4)**2*A0inv(:,4) 685 & + rLyitemp(:,5)**2*A0inv(:,5) 686c 687c.....****************calculation of giju for DC term*************** 688c 689c.... for the notation of different numbering 690c 691 call e3gijd( dxidx, gijd ) 692 693 gijdu(:,1)=gijd(:,1) 694 gijdu(:,2)=gijd(:,3) 695 gijdu(:,3)=gijd(:,6) 696 gijdu(:,4)=gijd(:,2) 697 gijdu(:,5)=gijd(:,4) 698 gijdu(:,6)=gijd(:,5) 699c 700c 701 detgijI = one/( 702 & gijdu(:,1) * gijdu(:,2) * gijdu(:,3) 703 & - gijdu(:,1) * gijdu(:,6) * gijdu(:,6) 704 & - gijdu(:,4) * gijdu(:,4) * gijdu(:,3) 705 & + gijdu(:,4) * gijdu(:,5) * gijdu(:,6) * two 706 & - gijdu(:,5) * gijdu(:,5) * gijdu(:,2) 707 & ) 708 giju(:,1) = detgijI * (gijdu(:,2)*gijdu(:,3) 709 & - gijdu(:,6)**2) 710 giju(:,2) = detgijI * (gijdu(:,1)*gijdu(:,3) 711 & - gijdu(:,5)**2) 712 giju(:,3) = detgijI * (gijdu(:,1)*gijdu(:,2) 713 & - gijdu(:,4)**2) 714 giju(:,4) = detgijI * (gijdu(:,5)*gijdu(:,6) 715 & - gijdu(:,4)*gijdu(:,3) ) 716 giju(:,5) = detgijI * (gijdu(:,4)*gijdu(:,6) 717 & - gijdu(:,5)*gijdu(:,2) ) 718 giju(:,6) = detgijI * (gijdu(:,4)*gijdu(:,5) 719 & - gijdu(:,1)*gijdu(:,6) ) 720c 721 endif ! end of iDC.ne.0 722c 723c.... calculate (tau Lym) --> (rLymi) 724c 725 if(ires.ne.1 ) then 726 eigmax = rLymi 727 rLymi = zero 728 do i = 1, nflow 729 do j = 1, nflow 730 rLymi(:,i) = rLymi(:,i) + Tau(:,i,j) * eigmax(:,j) 731 enddo 732 enddo 733 endif 734c 735c INCORRECT NOW ! flops = flops + 255*npro 736c 737c 738c.... return 739c 740 return 741 end 742c 743 744 745 subroutine e3tau_nd_modal (rho, cp, rmu, T, 746 & u1, u2, u3, 747 & a1, a2, a3, 748 & con, dxidx, rLyi, 749 & rLymi, Tau, rk, 750 & giju, rTLS, raLS, 751 & cv, compK, pres, 752 & A0inv, dVdY) 753c 754c---------------------------------------------------------------------- 755c 756c This routine computes the diagonal Tau for least-squares operator. 757c 758c We receive the regular residual L operator and a 759c modified residual L operator, calculate tau, and return values for 760c tau and tau times these operators to maintain the format of past 761c ENSA codes 762c 763c input: 764c rho (npro) : density 765c T (npro) : temperature 766c cp (npro) : specific heat at constant pressure 767c u1 (npro) : x1-velocity component 768c u2 (npro) : x2-velocity component 769c u3 (npro) : x3-velocity component 770c dxidx (npro,nsd,nsd) : inverse of deformation gradient 771c rLyi (npro,nflow) : least-squares residual vector 772c rLymi (npro,nflow) : modified least-squares residual vector 773c 774c output: 775c rLyi (npro,nflow) : least-squares residual vector 776c rLymi (npro,nflow) : modified least-squares residual vector 777c Mtau (npro,5,5) : Matrix of Stabilized Parameters 778c 779c 780c Zdenek Johan, Summer 1990. (Modified from e2tau.f) 781c Zdenek Johan, Winter 1991. (Fortran 90) 782c---------------------------------------------------------------------- 783c 784 include "common.h" 785c 786 dimension rho(npro), con(npro), 787 & cp(npro), u1(npro), 788 & u2(npro), u3(npro), 789 & dxidx(npro,nsd,nsd), rk(npro), 790 & rLyi(npro,nflow,ipord), 791 & rLymi(npro,nflow,ipord), dVdY(npro,15), 792 & rTLS(npro), raLS(npro), 793 & rLyitemp(npro,nflow), detgijI(npro) 794c 795 dimension rmu(npro), cv(npro), 796 & gijd(npro,6), 797 & fact(npro), giju(npro,6), 798 & A0inv(npro,15),gijdu(npro,6), compK(npro,10), 799 & a1(npro), a2(npro), a3(npro), 800 & T(npro), pres(npro), alphap(npro), 801 & betaT(npro), gamb(npro), c(npro), 802 & u(npro), eb1(npro), Q(npro,5,5), 803 & rlam(npro,5), eigmax(npro,5), Pe(npro), 804 & STau(npro,15, ipord), Tau(npro,nflow,nflow,ipord), 805 & rlamtmp(npro,5,ipord) 806 807 808c... build some necessary quantities at integration point: 809 810c alfap = one / T 811c betaT = one / pres 812 gamb = gamma1 813 c = sqrt( (gamma * Rgas) * T ) 814 u = sqrt(u1**2+u2**2+u3**2) 815 eb1 = cp*T - pt5*(u1**2+u2**2+u3**2) 816 817c... get eigenvectors for diagonalizing Tau_adv (e.v) before final rotations 818c... Note: the ordering of eigenvectors corresponds to the following 819c... ordering of the eigenvalues in the 1-st Euler Jacobian rotated to 820c... the streamline coordinates 821 822c |u+c | 823c | u | 824c | u | 825c | u | ! for origins of this see Beam, Warming, Hyett 826c | u-c| ! Mathematics of Computation vol. 29 No. 132 p. 1037 827 828c... different ordering assumed in Shakib/Johan entropy vars. code 829 830 831 832 call e3eig1 (rho, T, cp, 833 & gamb, c, u1, 834 & u2, u3, a1, 835 & a2, a3, eb1, 836 & dxidx, u, Q) 837 838c... Perform a Jacobi rotation on the Lambda matrix as well as adjust 839c... the eigenvectors. Tau still remains in entropy variables. 840 841 842 843 call e3eig2 (u, c, a1, 844 & a2, a3, rlam, 845 & Q, eigmax) 846 847c 848c.... invert the eigenvalues 849c 850 851 852 where (rlam .gt. ((epsM**2) * eigmax)) 853 rlam = one / sqrt(rlam) 854 elsewhere 855 rlam = zero 856 endwhere 857 858 do i = 1, ipord 859 rlamtmp(:,:,i) = rlam(:,:) 860 enddo 861 862c 863c.... flop count 864c 865 ! flops = flops + 65*npro 866 867c.... reduce the diffusion contribution 868c 869 870 if (Navier .eq. 1) then 871c 872c.... calculate sigma 873c 874 875 do i = 1, nflow ! diff. corr for every mode of Tau 876 877 c(:) = pt33 * ( 878 & Q(:,2,i) * ( compK(:, 1) * Q(:,2,i) + compK(:, 2) * Q(:,3,i) 879 & + compK(:, 4) * Q(:,4,i) + compK(:, 7) * Q(:,5,i) ) 880 & + Q(:,3,i) * ( compK(:, 2) * Q(:,2,i) + compK(:, 3) * Q(:,3,i) 881 & + compK(:, 5) * Q(:,4,i) + compK(:, 8) * Q(:,5,i) ) 882 & + Q(:,4,i) * ( compK(:, 4) * Q(:,2,i) + compK(:, 5) * Q(:,3,i) 883 & + compK(:, 6) * Q(:,4,i) + compK(:, 9) * Q(:,5,i) ) 884 & + Q(:,5,i) * ( compK(:, 7) * Q(:,2,i) + compK(:, 8) * Q(:,3,i) 885 & + compK(:, 9) * Q(:,4,i) + compK(:,10) * Q(:,5,i) ) 886 & ) 887 888c... get Peclet Number 889 890 891 Pe(:) = one / (c(:)*rlam(:,i)) 892 893 894c 895c... branch out into different polynomial orders 896c 897 898 899 if (ipord == 1) then ! linears - l = l^a*(coth(Pe)-1/Pe) 900 ! approx. l = l^a*min(1/3*Pe,1) 901 rlamtmp(:,i,1) = rlamtmp(:,i,1)*min(pt33*Pe(:),one) 902 903 endif 904 905 if (ipord == 2) then 906 907 rlamtmp(:,i,1) = rlamtmp(:,i,1)*min((1.0/15.0)*Pe(:),pt33) 908 rlamtmp(:,i,2) = rlamtmp(:,i,2)*min((1.0/12.0)*Pe(:),pt5) 909 910 endif 911 912 if (ipord == 3) then ! cubics - Recent Work 913 914 do ii = 1, npro 915 if (Pe(ii).lt.3.0) then 916 rlamtmp(ii,i,1) = rlamtmp(ii,i,1)* 917 & 0.00054*Pe(ii)**2 918 endif 919 920 if ((Pe(ii).ge.3).and.(Pe(ii).lt.17.20)) then 921 rlamtmp(ii,i,1) = rlamtmp(ii,i,1)*(0.0114*Pe(ii) 922 & -0.0294) 923 endif 924 925 if (Pe(ii).ge.17.20) then 926 rlamtmp(ii,i,1) = rlamtmp(ii,i,1)*(1.0/6.0) 927 endif 928 929 enddo 930 931 rlamtmp(:,i,2) = rlamtmp(:,i,2)*min((1.0/45.0)*Pe(:) 932 & ,0.2d0) 933 rlamtmp(:,i,3) = rlamtmp(:,i,3)*min((1.0/25.0)*Pe(:) 934 & ,pt33) 935 936 937 endif ! done w/ different polynomial orders 938 939 enddo ! done with loop over dof's 940 941 endif ! done with diffusive correction 942 943 944c 945c.... form Tau (symmetric - entropy variables - then convert) 946c 947 do i = 1, ipord 948 949 STau(:, 1, i) = rlamtmp(:,1,i) * Q(:,1,1) * Q(:,1,1) + 950 & rlamtmp(:,2,i) * Q(:,1,2) * Q(:,1,2) + 951 & rlamtmp(:,3,i) * Q(:,1,3) * Q(:,1,3) + 952 & rlamtmp(:,4,i) * Q(:,1,4) * Q(:,1,4) + 953 & rlamtmp(:,5,i) * Q(:,1,5) * Q(:,1,5) 954c 955 STau(:, 2, i) = rlamtmp(:,1,i) * Q(:,1,1) * Q(:,2,1) + 956 & rlamtmp(:,2,i) * Q(:,1,2) * Q(:,2,2) + 957 & rlamtmp(:,3,i) * Q(:,1,3) * Q(:,2,3) + 958 & rlamtmp(:,4,i) * Q(:,1,4) * Q(:,2,4) + 959 & rlamtmp(:,5,i) * Q(:,1,5) * Q(:,2,5) 960c 961 STau(:, 3, i) = rlamtmp(:,1,i) * Q(:,2,1) * Q(:,2,1) + 962 & rlamtmp(:,2,i) * Q(:,2,2) * Q(:,2,2) + 963 & rlamtmp(:,3,i) * Q(:,2,3) * Q(:,2,3) + 964 & rlamtmp(:,4,i) * Q(:,2,4) * Q(:,2,4) + 965 & rlamtmp(:,5,i) * Q(:,2,5) * Q(:,2,5) 966c 967 STau(:, 4, i) = rlamtmp(:,1,i) * Q(:,1,1) * Q(:,3,1) + 968 & rlamtmp(:,2,i) * Q(:,1,2) * Q(:,3,2) + 969 & rlamtmp(:,3,i) * Q(:,1,3) * Q(:,3,3) + 970 & rlamtmp(:,4,i) * Q(:,1,4) * Q(:,3,4) + 971 & rlamtmp(:,5,i) * Q(:,1,5) * Q(:,3,5) 972c 973 STau(:, 5, i) = rlamtmp(:,1,i) * Q(:,2,1) * Q(:,3,1) + 974 & rlamtmp(:,2,i) * Q(:,2,2) * Q(:,3,2) + 975 & rlamtmp(:,3,i) * Q(:,2,3) * Q(:,3,3) + 976 & rlamtmp(:,4,i) * Q(:,2,4) * Q(:,3,4) + 977 & rlamtmp(:,5,i) * Q(:,2,5) * Q(:,3,5) 978c 979 STau(:, 6, i) = rlamtmp(:,1,i) * Q(:,3,1) * Q(:,3,1) + 980 & rlamtmp(:,2,i) * Q(:,3,2) * Q(:,3,2) + 981 & rlamtmp(:,3,i) * Q(:,3,3) * Q(:,3,3) + 982 & rlamtmp(:,4,i) * Q(:,3,4) * Q(:,3,4) + 983 & rlamtmp(:,5,i) * Q(:,3,5) * Q(:,3,5) 984c 985 STau(:, 7, i) = rlamtmp(:,1,i) * Q(:,1,1) * Q(:,4,1) + 986 & rlamtmp(:,2,i) * Q(:,1,2) * Q(:,4,2) + 987 & rlamtmp(:,3,i) * Q(:,1,3) * Q(:,4,3) + 988 & rlamtmp(:,4,i) * Q(:,1,4) * Q(:,4,4) + 989 & rlamtmp(:,5,i) * Q(:,1,5) * Q(:,4,5) 990c 991 STau(:, 8, i) = rlamtmp(:,1,i) * Q(:,2,1) * Q(:,4,1) + 992 & rlamtmp(:,2,i) * Q(:,2,2) * Q(:,4,2) + 993 & rlamtmp(:,3,i) * Q(:,2,3) * Q(:,4,3) + 994 & rlamtmp(:,4,i) * Q(:,2,4) * Q(:,4,4) + 995 & rlamtmp(:,5,i) * Q(:,2,5) * Q(:,4,5) 996c 997 STau(:, 9, i) = rlamtmp(:,1,i) * Q(:,3,1) * Q(:,4,1) + 998 & rlamtmp(:,2,i) * Q(:,3,2) * Q(:,4,2) + 999 & rlamtmp(:,3,i) * Q(:,3,3) * Q(:,4,3) + 1000 & rlamtmp(:,4,i) * Q(:,3,4) * Q(:,4,4) + 1001 & rlamtmp(:,5,i) * Q(:,3,5) * Q(:,4,5) 1002c 1003 STau(:,10, i) = rlamtmp(:,1,i) * Q(:,4,1) * Q(:,4,1) + 1004 & rlamtmp(:,2,i) * Q(:,4,2) * Q(:,4,2) + 1005 & rlamtmp(:,3,i) * Q(:,4,3) * Q(:,4,3) + 1006 & rlamtmp(:,4,i) * Q(:,4,4) * Q(:,4,4) + 1007 & rlamtmp(:,5,i) * Q(:,4,5) * Q(:,4,5) 1008c 1009 STau(:,11, i) = rlamtmp(:,1,i) * Q(:,1,1) * Q(:,5,1) + 1010 & rlamtmp(:,2,i) * Q(:,1,2) * Q(:,5,2) + 1011 & rlamtmp(:,3,i) * Q(:,1,3) * Q(:,5,3) + 1012 & rlamtmp(:,4,i) * Q(:,1,4) * Q(:,5,4) + 1013 & rlamtmp(:,5,i) * Q(:,1,5) * Q(:,5,5) 1014c 1015 STau(:,12, i) = rlamtmp(:,1,i) * Q(:,2,1) * Q(:,5,1) + 1016 & rlamtmp(:,2,i) * Q(:,2,2) * Q(:,5,2) + 1017 & rlamtmp(:,3,i) * Q(:,2,3) * Q(:,5,3) + 1018 & rlamtmp(:,4,i) * Q(:,2,4) * Q(:,5,4) + 1019 & rlamtmp(:,5,i) * Q(:,2,5) * Q(:,5,5) 1020c 1021 STau(:,13, i) = rlamtmp(:,1,i) * Q(:,3,1) * Q(:,5,1) + 1022 & rlamtmp(:,2,i) * Q(:,3,2) * Q(:,5,2) + 1023 & rlamtmp(:,3,i) * Q(:,3,3) * Q(:,5,3) + 1024 & rlamtmp(:,4,i) * Q(:,3,4) * Q(:,5,4) + 1025 & rlamtmp(:,5,i) * Q(:,3,5) * Q(:,5,5) 1026c 1027 STau(:,14, i) = rlamtmp(:,1,i) * Q(:,4,1) * Q(:,5,1) + 1028 & rlamtmp(:,2,i) * Q(:,4,2) * Q(:,5,2) + 1029 & rlamtmp(:,3,i) * Q(:,4,3) * Q(:,5,3) + 1030 & rlamtmp(:,4,i) * Q(:,4,4) * Q(:,5,4) + 1031 & rlamtmp(:,5,i) * Q(:,4,5) * Q(:,5,5) 1032c 1033 STau(:,15, i) = rlamtmp(:,1,i) * Q(:,5,1) * Q(:,5,1) + 1034 & rlamtmp(:,2,i) * Q(:,5,2) * Q(:,5,2) + 1035 & rlamtmp(:,3,i) * Q(:,5,3) * Q(:,5,3) + 1036 & rlamtmp(:,4,i) * Q(:,5,4) * Q(:,5,4) + 1037 & rlamtmp(:,5,i) * Q(:,5,5) * Q(:,5,5) 1038 1039 enddo 1040 1041c 1042c... Form Primitive Variable Tau as [dY/dV]*tilde{Tau}, 1043c... see G.Hauke's thesis appendix for [dY/dV] matrix 1044c 1045 do k = 1, ipord 1046 1047 betaT = cp*T + pt5*(u1**2+u2**2+u3**2) !reuse betaT 1048 1049 Tau(:,1,1,k) = rho*T*(STau(:,1,k)+u1*STau(:,2,k)+ 1050 & u2*STau(:,4,k)+u3*STau(:,7,k)+betaT*STau(:,11,k)) 1051 Tau(:,1,2,k) = rho*T*(STau(:,2,k)+u1*STau(:,3,k)+ 1052 & u2*STau(:,5,k)+u3*STau(:,8,k)+betaT*STau(:,12,k)) 1053 Tau(:,1,3,k) = rho*T*(STau(:,4,k)+u1*STau(:,5,k)+ 1054 & u2*STau(:,6,k)+u3*STau(:,9,k)+betaT*STau(:,13,k)) 1055 Tau(:,1,4,k) = rho*T*(STau(:,7,k)+u1*STau(:,8,k)+ 1056 & u2*STau(:,9,k)+u3*STau(:,10,k)+betaT*STau(:,14,k)) 1057 Tau(:,1,5,k) = rho*T*(STau(:,11,k)+u1*STau(:,12,k)+ 1058 & u2*STau(:,13,k)+u3*STau(:,14,k)+betaT*STau(:,15,k)) 1059 1060 1061 Tau(:,2,1,k) = T(:)*(STau(:,2,k) + u1(:)*STau(:,11,k)) 1062 Tau(:,2,2,k) = T(:)*(STau(:,3,k) + u1(:)*STau(:,12,k)) 1063 Tau(:,2,3,k) = T(:)*(STau(:,5,k) + u1(:)*STau(:,13,k)) 1064 Tau(:,2,4,k) = T(:)*(STau(:,8,k) + u1(:)*STau(:,14,k)) 1065 Tau(:,2,5,k) = T(:)*(STau(:,12,k) + u1(:)*STau(:,15,k)) 1066 1067 Tau(:,3,1,k) = T(:)*(STau(:,4,k) + u2(:)*STau(:,11,k)) 1068 Tau(:,3,2,k) = T(:)*(STau(:,5,k) + u2(:)*STau(:,12,k)) 1069 Tau(:,3,3,k) = T(:)*(STau(:,6,k) + u2(:)*STau(:,13,k)) 1070 Tau(:,3,4,k) = T(:)*(STau(:,9,k) + u2(:)*STau(:,14,k)) 1071 Tau(:,3,5,k) = T(:)*(STau(:,13,k) + u2(:)*STau(:,15,k)) 1072 1073 Tau(:,4,1,k) = T(:)*(STau(:,7,k) + u3(:)*STau(:,11,k)) 1074 Tau(:,4,2,k) = T(:)*(STau(:,8,k) + u3(:)*STau(:,12,k)) 1075 Tau(:,4,3,k) = T(:)*(STau(:,9,k) + u3(:)*STau(:,13,k)) 1076 Tau(:,4,4,k) = T(:)*(STau(:,10,k) + u3(:)*STau(:,14,k)) 1077 Tau(:,4,5,k) = T(:)*(STau(:,14,k) + u3(:)*STau(:,15,k)) 1078 1079 betaT = T**2 1080 1081 Tau(:,5,1,k) = betaT(:)*STau(:,11,k) 1082 Tau(:,5,2,k) = betaT(:)*STau(:,12,k) 1083 Tau(:,5,3,k) = betaT(:)*STau(:,13,k) 1084 Tau(:,5,4,k) = betaT(:)*STau(:,14,k) 1085 Tau(:,5,5,k) = betaT(:)*STau(:,15,k) 1086 1087 enddo 1088 1089c 1090c... done with conversion to pressure primitive variables 1091c... now need to interface with the rest of the computations 1092c 1093 1094c... finally multiply this tau matrix times the 1095c two residual vectors 1096c 1097c ... calculate (tau Ly) --> (rLyi) 1098c ... storing rLyi for the DC term 1099 1100 1101 if(iDC .ne. 0) rLyitemp(:,:)=rLyi(:,:,1) 1102 1103 if(ires.eq.3 .or. ires .eq. 1) then 1104 eigmax(:,:) = rLyi(:,:,1) !reuse 1105 rLyi = zero 1106 do k = 1, ipord 1107 do i = 1, nflow 1108 do j = 1, nflow 1109 rLyi(:,i,k) = rLyi(:,i,k)+Tau(:,i,j,k)*eigmax(:,j) 1110 enddo 1111 enddo 1112 enddo 1113 endif 1114 1115 1116 if(iDC .ne. 0) then 1117c.....calculation of rTLS & raLS for DC term 1118c 1119c.....calculation of (Ly - S).tau^tilda*(Ly - S) 1120c 1121 rTLS = rLYItemp(:,1) * (rLyi(:,1,1)*dVdY(:,1) 1122 & + dVdY(:,2)*rLyi(:,2,1) + dVdY(:,4)*rLyi(:,3,1) 1123 & + rLyi(:,4,1)*dVdY(:,7) + dVdY(:,11)*rLyi(:,5,1)) 1124 & + rLyitemp(:,2) * (rLyi(:,2,1)*dVdY(:,3) 1125 & + rLyi(:,3,1)*dVdY(:,5) + dVdY(:,8)*rLyi(:,4,1) 1126 & + rLyi(:,5,1)*dVdY(:,12)) 1127 & + rLyitemp(:,3) * (rLyi(:,3,1)*dVdY(:,6) 1128 & + dVdY(:,9)*rLyi(:,4,1) + dVdY(:,13)*rLyi(:,5,1)) 1129 & + rLyitemp(:,4) * (rLyi(:,4,1)*dVdY(:,10) 1130 & + dVdY(:,14)*rLyi(:,5,1)) 1131 & + rLyitemp(:,5) * (dVdY(:,15)*rLyi(:,5,1)) 1132 1133c 1134c...... calculation of (Ly - S).A0inv*(Ly - S) 1135c 1136 raLS = two*rLyitemp(:,4)*rLyitemp(:,5)*A0inv(:,15) 1137 & + two*rLyitemp(:,3)*rLyitemp(:,5)*A0inv(:,14) 1138 & + two*rLyitemp(:,1)*rLyitemp(:,2)*A0inv( :,6) 1139 & + two*rLyitemp(:,2)*rLyitemp(:,3)*A0inv(:,10) 1140 & + two*rLyitemp(:,2)*rLyitemp(:,4)*A0inv(:,11) 1141 & + two*rLyitemp(:,1)*rLyitemp(:,3)*A0inv( :,7) 1142 & + two*rLyitemp(:,3)*rLyitemp(:,4)*A0inv(:,13) 1143 & + two*rLyitemp(:,2)*rLyitemp(:,5)*A0inv(:,12) 1144 & + two*rLyitemp(:,1)*rLyitemp(:,4)*A0inv( :,8) 1145 & + two*rLyitemp(:,1)*rLyitemp(:,5)*A0inv( :,9) 1146 & + rLyitemp(:,1)**2*A0inv(:,1) 1147 & + rLyitemp(:,2)**2*A0inv(:,2) 1148 & + rLyitemp(:,3)**2*A0inv(:,3) 1149 & + rLyitemp(:,4)**2*A0inv(:,4) 1150 & + rLyitemp(:,5)**2*A0inv(:,5) 1151c 1152c.....****************calculation of giju for DC term*************** 1153c 1154c.... for the notation of different numbering 1155c 1156 gijdu(:,1)=gijd(:,1) 1157 gijdu(:,2)=gijd(:,3) 1158 gijdu(:,3)=gijd(:,6) 1159 gijdu(:,4)=gijd(:,2) 1160 gijdu(:,5)=gijd(:,4) 1161 gijdu(:,6)=gijd(:,5) 1162c 1163c 1164 detgijI = one/( 1165 & gijdu(:,1) * gijdu(:,2) * gijdu(:,3) 1166 & - gijdu(:,1) * gijdu(:,6) * gijdu(:,6) 1167 & - gijdu(:,4) * gijdu(:,4) * gijdu(:,3) 1168 & + gijdu(:,4) * gijdu(:,5) * gijdu(:,6) * two 1169 & - gijdu(:,5) * gijdu(:,5) * gijdu(:,2) 1170 & ) 1171 giju(:,1) = detgijI * (gijdu(:,2)*gijdu(:,3) 1172 & - gijdu(:,6)**2) 1173 giju(:,2) = detgijI * (gijdu(:,1)*gijdu(:,3) 1174 & - gijdu(:,5)**2) 1175 giju(:,3) = detgijI * (gijdu(:,1)*gijdu(:,2) 1176 & - gijdu(:,4)**2) 1177 giju(:,4) = detgijI * (gijdu(:,5)*gijdu(:,6) 1178 & - gijdu(:,4)*gijdu(:,3) ) 1179 giju(:,5) = detgijI * (gijdu(:,4)*gijdu(:,6) 1180 & - gijdu(:,5)*gijdu(:,2) ) 1181 giju(:,6) = detgijI * (gijdu(:,4)*gijdu(:,5) 1182 & - gijdu(:,1)*gijdu(:,6) ) 1183c 1184 endif ! end of iDC.ne.0 1185c 1186c.... calculate (tau Lym) --> (rLymi) 1187c 1188 if(ires.ne.1 ) then 1189 eigmax(:,:) = rLymi(:,:,1) 1190 rLymi = zero 1191 do k = 1, ipord 1192 do i = 1, nflow 1193 do j = 1, nflow 1194 rLymi(:,i,k) = rLymi(:,i,k)+Tau(:,i,j,k)*eigmax(:,j) 1195 enddo 1196 enddo 1197 enddo 1198 endif 1199c 1200c INCORRECT NOW ! flops = flops + 255*npro 1201c 1202c 1203c.... return 1204c 1205 return 1206 end 1207c 1208 1209 1210 1211 subroutine e3tauSclr(rho, rmu, A0t, 1212 & u1, u2, u3, 1213 & dxidx, rLyti, rLymti, 1214 & taut, rk, raLSt, 1215 & rTLSt, giju) 1216c 1217c---------------------------------------------------------------------- 1218c 1219c This routine computes the diagonal Tau for least-squares operator. 1220c This Tau is the one proposed for nearly incompressible flows by 1221c Franca and Frey. We receive the regular residual L operator and a 1222c modified residual L operator, calculate tau, and return values for 1223c tau and tau times these operators to maintain the format of past 1224c ENSA codes 1225c 1226c input: 1227c rho (npro) : density 1228c T (npro) : temperature 1229c cp (npro) : specific heat at constant pressure 1230c u1 (npro) : x1-velocity component 1231c u2 (npro) : x2-velocity component 1232c u3 (npro) : x3-velocity component 1233c dxidx (npro,nsd,nsd) : inverse of deformation gradient 1234c rLyti (npro) : least-squares residual vector 1235c rLymti (npro) : modified least-squares residual vector 1236c 1237c output: 1238c rLyti (npro,nflow) : least-squares residual vector 1239c rLymti (npro,nflow) : modified least-squares residual vector 1240c tau (npro,3) : 3 taus 1241c 1242c 1243c Zdenek Johan, Summer 1990. (Modified from e2tau.f) 1244c Zdenek Johan, Winter 1991. (Fortran 90) 1245c---------------------------------------------------------------------- 1246c 1247 use turbSA 1248 include "common.h" 1249c 1250 dimension rho(npro), T(npro), 1251 & cp(npro), u1(npro), 1252 & u2(npro), u3(npro), 1253 & dxidx(npro,nsd,nsd), rk(npro), 1254 & taut(npro), rLyti(npro), 1255 & rLymti(npro) 1256c 1257 dimension rmu(npro), A0t(npro), 1258 & gijd(npro,6), uh1(npro), 1259 & fact(npro), h2o2u(npro), 1260 & rlytitemp(npro), raLSt(npro), 1261 & rTLSt(npro), gijdu(npro,6), 1262 & giju(npro,6), detgijI(npro) 1263c 1264c 1265 call e3gijd( dxidx, gijd ) 1266 1267c 1268c next form the diffusive length scale |u| h_1 = 2 ( ui (gijd)-1 uj)^{1/2} 1269c 1270c dividing factor for the inverse of gijd 1271c 1272 fact = gijd(:,1) * gijd(:,3) * gijd(:,6) 1273 & - gijd(:,1) * gijd(:,5) * gijd(:,5) 1274 & - gijd(:,3) * gijd(:,4) * gijd(:,4) 1275 & - gijd(:,6) * gijd(:,2) * gijd(:,2) 1276 & + gijd(:,2) * gijd(:,4) * gijd(:,5) * two 1277c 1278 uh1= u1*u1*(gijd(:,3)*gijd(:,6)-gijd(:,5)*gijd(:,5)) 1279 & + u2*u2*(gijd(:,1)*gijd(:,6)-gijd(:,4)*gijd(:,4)) 1280 & + u3*u3*(gijd(:,1)*gijd(:,3)-gijd(:,2)*gijd(:,2)) 1281 & + two *(u1*u2*(gijd(:,4)*gijd(:,5)-gijd(:,2)*gijd(:,6)) 1282 & + u1*u3*(gijd(:,2)*gijd(:,5)-gijd(:,4)*gijd(:,3)) 1283 & + u2*u3*(gijd(:,4)*gijd(:,2)-gijd(:,1)*gijd(:,5))) 1284c 1285c at this point we have (u h1)^2 * inverse coefficient^2 / 4 1286c 1287c ^ fact 1288c 1289 1290 uh1= two * sqrt(uh1/fact) 1291 1292c 1293c next form the advective length scale |u|/h_2 = 2 ( ui (gijd) uj)^{1/2} 1294c 1295 h2o2u = u1*u1*gijd(:,1) 1296 & + u2*u2*gijd(:,3) 1297 & + u3*u3*gijd(:,6) 1298 & +(u1*u2*gijd(:,2) 1299 & + u1*u3*gijd(:,4) 1300 & + u2*u3*gijd(:,5))*two + 1.0e-15 !FIX FOR INVALID MESHES 1301c 1302c at this point we have (2 u / h_2)^2 1303c 1304 1305c call tnanqe(h2o2u,1,"riaconv ") 1306 1307 h2o2u = one / sqrt(h2o2u) ! this flips it over leaves it h_2/(2u) 1308c 1309c... momentum tau 1310c 1311c 1312c... rmu will now hold the total (molecular plus eddy) viscosity... 1313 dts=one/(Dtgl*dtsfct) 1314 if(iremoveStabTimeTerm.gt.0) dts = dts*100000 ! remove time term from scalar 1315! Duct code had this dts=1.0e16 1316 taut(:)=min(dts,h2o2u) 1317 taut(:)=taut(:)/rho 1318 taut(:)=min(taut(:),h2o2u*h2o2u*rk*pt66*saSigma/rmu) 1319c 1320c... finally multiply this tau matrix times the 1321c two residual vectors 1322c 1323c.... calculate (tau Lyt) --> (rLyti) 1324c 1325c.... storing rLyi for the DC term 1326 rLytitemp=rLyti 1327 1328 if(ires.eq.3 .or. ires .eq. 1) then 1329 rLyti(:) = taut(:) * rLyti(:) 1330 1331 endif 1332 if(iDCSclr(1) .ne. 0) then 1333c..... calculation of rTLS & raLS for DC term 1334c..... calculation of (Ly - S).tau^tilda*(Ly - S) 1335c 1336 rTLSt = rLYtItemp(:)*rLyti(:) 1337c 1338c...... calculation of (Ly - S).A0inv*(Ly - S) 1339c 1340 raLSt = rLYtItemp(:) * rLYtItemp(:) 1341c.....*****************calculation of giju for DC term****************** 1342c 1343c.... for the notation of different numbering 1344c 1345 gijdu(:,1)=gijd(:,1) 1346 gijdu(:,2)=gijd(:,3) 1347 gijdu(:,3)=gijd(:,6) 1348 gijdu(:,4)=gijd(:,2) 1349 gijdu(:,5)=gijd(:,4) 1350 gijdu(:,6)=gijd(:,5) 1351c 1352c we are going to need this in the dc factor later so we calculate it. 1353c 1354 detgijI = one/( 1355 & gijdu(:,1) * gijdu(:,2) * gijdu(:,3) 1356 & - gijdu(:,1) * gijdu(:,6) * gijdu(:,6) 1357 & - gijdu(:,4) * gijdu(:,4) * gijdu(:,3) 1358 & + gijdu(:,4) * gijdu(:,5) * gijdu(:,6) * two 1359 & - gijdu(:,5) * gijdu(:,5) * gijdu(:,2) 1360 & ) 1361 giju(:,1) = detgijI * (gijdu(:,2)*gijdu(:,3) 1362 & - gijdu(:,6)**2) 1363 giju(:,2) = detgijI * (gijdu(:,1)*gijdu(:,3) 1364 & - gijdu(:,5)**2) 1365 giju(:,3) = detgijI * (gijdu(:,1)*gijdu(:,2) 1366 & - gijdu(:,4)**2) 1367 giju(:,4) = detgijI * (gijdu(:,5)*gijdu(:,6) 1368 & - gijdu(:,4)*gijdu(:,3) ) 1369 giju(:,5) = detgijI * (gijdu(:,4)*gijdu(:,6) 1370 & - gijdu(:,5)*gijdu(:,2) ) 1371 giju(:,6) = detgijI * (gijdu(:,4)*gijdu(:,5) 1372 & - gijdu(:,1)*gijdu(:,6) ) 1373c 1374 endif ! end of iDCSclr(1).ne.0 1375c 1376c.... calculate (tau Lym) --> (rLymi) 1377c 1378c if(ires.ne.1 ) then 1379c rLymi(:,1) = tau(:,1) * rLymi(:,1) 1380c rLymi(:,2) = tau(:,2) * rLymi(:,2) 1381c rLymi(:,3) = tau(:,2) * rLymi(:,3) 1382c rLymi(:,4) = tau(:,2) * rLymi(:,4) 1383c rLymi(:,5) = tau(:,3) * rLymi(:,5) 1384c endif 1385c 1386c INCORRECT NOW ! flops = flops + 255*npro 1387c 1388c 1389c.... return 1390c 1391 return 1392 end 1393 1394c----------------------------------------------------------------------- 1395c get the metric tensor g_{ij}=xi_{k,i} xi_{k,j}. 1396c----------------------------------------------------------------------- 1397 subroutine e3gijd( dxidx, gijd ) 1398 1399 include "common.h" 1400 1401 real*8 dxidx(npro,nsd,nsd), gijd(npro,6), 1402 & tmp1(npro), tmp2(npro), 1403 & tmp3(npro) 1404c form metric tensor g_{ij}=xi_{k,i} xi_{k,j}. It is a symmetric 1405c tensor so we only form 6 components and use symmetric matrix numbering. 1406c (d for down since giju=[gijd]^{-1}) 1407c (Note FARZIN and others use numbering of 1,2,3 being diagonal 456 off) 1408 if (lcsyst .ge. 2) then 1409 1410 gijd(:,1) = dxidx(:,1,1) * dxidx(:,1,1) 1411 & + dxidx(:,2,1) * dxidx(:,2,1) 1412 & + dxidx(:,3,1) * dxidx(:,3,1) 1413c 1414 gijd(:,2) = dxidx(:,1,1) * dxidx(:,1,2) 1415 & + dxidx(:,2,1) * dxidx(:,2,2) 1416 & + dxidx(:,3,1) * dxidx(:,3,2) 1417c 1418 gijd(:,3) = dxidx(:,1,2) * dxidx(:,1,2) 1419 & + dxidx(:,2,2) * dxidx(:,2,2) 1420 & + dxidx(:,3,2) * dxidx(:,3,2) 1421c 1422 gijd(:,4) = dxidx(:,1,1) * dxidx(:,1,3) 1423 & + dxidx(:,2,1) * dxidx(:,2,3) 1424 & + dxidx(:,3,1) * dxidx(:,3,3) 1425c 1426 gijd(:,5) = dxidx(:,1,2) * dxidx(:,1,3) 1427 & + dxidx(:,2,2) * dxidx(:,2,3) 1428 & + dxidx(:,3,2) * dxidx(:,3,3) 1429c 1430 gijd(:,6) = dxidx(:,1,3) * dxidx(:,1,3) 1431 & + dxidx(:,2,3) * dxidx(:,2,3) 1432 & + dxidx(:,3,3) * dxidx(:,3,3) 1433c 1434 else if (lcsyst .eq. 1) then 1435c 1436c There is an invariance problem with tets 1437c It is fixed by the following modifications to gijd 1438c 1439 1440 c1 = 1.259921049894873D+00 1441 c2 = 6.299605249474365D-01 1442c 1443 tmp1(:) = c1 * dxidx(:,1,1) + c2 *(dxidx(:,2,1)+dxidx(:,3,1)) 1444 tmp2(:) = c1 * dxidx(:,2,1) + c2 *(dxidx(:,1,1)+dxidx(:,3,1)) 1445 tmp3(:) = c1 * dxidx(:,3,1) + c2 *(dxidx(:,1,1)+dxidx(:,2,1)) 1446 gijd(:,1) = dxidx(:,1,1) * tmp1 1447 1 + dxidx(:,2,1) * tmp2 1448 2 + dxidx(:,3,1) * tmp3 1449c 1450 tmp1(:) = c1 * dxidx(:,1,2) + c2 *(dxidx(:,2,2)+dxidx(:,3,2)) 1451 tmp2(:) = c1 * dxidx(:,2,2) + c2 *(dxidx(:,1,2)+dxidx(:,3,2)) 1452 tmp3(:) = c1 * dxidx(:,3,2) + c2 *(dxidx(:,1,2)+dxidx(:,2,2)) 1453 gijd(:,2) = dxidx(:,1,1) * tmp1 1454 1 + dxidx(:,2,1) * tmp2 1455 2 + dxidx(:,3,1) * tmp3 1456c 1457 gijd(:,3) = dxidx(:,1,2) * tmp1 1458 1 + dxidx(:,2,2) * tmp2 1459 2 + dxidx(:,3,2) * tmp3 1460c 1461 tmp1(:) = c1 * dxidx(:,1,3) + c2 *(dxidx(:,2,3)+dxidx(:,3,3)) 1462 tmp2(:) = c1 * dxidx(:,2,3) + c2 *(dxidx(:,1,3)+dxidx(:,3,3)) 1463 tmp3(:) = c1 * dxidx(:,3,3) + c2 *(dxidx(:,1,3)+dxidx(:,2,3)) 1464 gijd(:,4) = dxidx(:,1,1) * tmp1 1465 1 + dxidx(:,2,1) * tmp2 1466 2 + dxidx(:,3,1) * tmp3 1467c 1468 gijd(:,5) = dxidx(:,1,2) * tmp1 1469 1 + dxidx(:,2,2) * tmp2 1470 2 + dxidx(:,3,2) * tmp3 1471c 1472 gijd(:,6) = dxidx(:,1,3) * tmp1 1473 1 + dxidx(:,2,3) * tmp2 1474 2 + dxidx(:,3,3) * tmp3 1475c 1476 else 1477c This is just the hex copied from above. I have 1478c to find my notes on invariance factors for wedges 1479c 1480 1481 gijd(:,1) = dxidx(:,1,1) * dxidx(:,1,1) 1482 & + dxidx(:,2,1) * dxidx(:,2,1) 1483 & + dxidx(:,3,1) * dxidx(:,3,1) 1484c 1485 gijd(:,2) = dxidx(:,1,1) * dxidx(:,1,2) 1486 & + dxidx(:,2,1) * dxidx(:,2,2) 1487 & + dxidx(:,3,1) * dxidx(:,3,2) 1488c 1489 gijd(:,3) = dxidx(:,1,2) * dxidx(:,1,2) 1490 & + dxidx(:,2,2) * dxidx(:,2,2) 1491 & + dxidx(:,3,2) * dxidx(:,3,2) 1492c 1493 gijd(:,4) = dxidx(:,1,1) * dxidx(:,1,3) 1494 & + dxidx(:,2,1) * dxidx(:,2,3) 1495 & + dxidx(:,3,1) * dxidx(:,3,3) 1496c 1497 gijd(:,5) = dxidx(:,1,2) * dxidx(:,1,3) 1498 & + dxidx(:,2,2) * dxidx(:,2,3) 1499 & + dxidx(:,3,2) * dxidx(:,3,3) 1500c 1501 gijd(:,6) = dxidx(:,1,3) * dxidx(:,1,3) 1502 & + dxidx(:,2,3) * dxidx(:,2,3) 1503 & + dxidx(:,3,3) * dxidx(:,3,3) 1504 endif 1505c 1506 return 1507 end 1508 1509