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