1 subroutine e3LS (A1, A2, A3, 2 & rho, rmu, cp, 3 & cv, con, T, 4 & u1, u2, u3, 5 & rLyi, dxidx, tau, 6 & ri, rmi, rk, 7 & dui, aci, A0, 8 & divqi, shape, shg, 9 & EGmass, stiff, WdetJ, 10 & giju, rTLS, raLS, 11 & A0inv, dVdY, rerrl, 12 & compK, pres, PTau) 13c 14c---------------------------------------------------------------------- 15c 16c This routine calculates the contribution of the least-squares 17c operator to the RHS vector and LHS tangent matrix. The temporary 18c results are put in ri. 19c 20c input: 21c A1 (npro,nflow,nflow) : A_1 22c A2 (npro,nflow,nflow) : A_2 23c A3 (npro,nflow,nflow) : A_3 24c rho (npro) : density 25c T (npro) : temperature 26c cp (npro) : specific heat at constant pressure 27c u1 (npro) : x1-velocity component 28c u2 (npro) : x2-velocity component 29c u3 (npro) : x3-velocity component 30c rLyi (npro,nflow) : least-squares residual vector 31c dxidx (npro,nsd,nsd) : inverse of deformation gradient 32c tau (npro,3) : stability parameter 33c PTau (npro,5,5) : matrix of stability parameters 34c rLyi (npro,nflow) : convective portion of least-squares 35c residual vector 36c divqi (npro,nflow-1) : divergence of diffusive flux 37c shape (npro,nshl) : element shape functions 38c shg (npro,nshl,nsd) : global element shape function grads 39c WdetJ (npro) : weighted jacobian determinant 40c stiff (npro,nsd*nflow,nsd*nflow) : stiffness matrix 41c EGmass(npro,nedof,nedof) : partial mass matrix 42c compK (npro,10) : K_ij in (eq.134) A new ... III 43c 44c output: 45c ri (npro,nflow*(nsd+1)) : partial residual 46c rmi (npro,nflow*(nsd+1)) : partial modified residual 47c EGmass (npro,nedof,nedof) : partial mass matrix 48c 49c 50c Zdenek Johan, Summer 1990. (Modified from e2ls.f) 51c Zdenek Johan, Winter 1991. (Fortran 90) 52c Kenneth Jansen, Winter 1997. Prim. Var. with Diag Tau 53c---------------------------------------------------------------------- 54c 55 include "common.h" 56 57c 58c passed arrays 59c 60 dimension A1(npro,nflow,nflow), A2(npro,nflow,nflow), 61 & A3(npro,nflow,nflow), cv(npro), 62 & A0(npro,nflow,nflow), rho(npro), 63 & con(npro), cp(npro), 64 & dui(npro,nflow), aci(npro,nflow), 65 & u1(npro), u2(npro), 66 & u3(npro), rk(npro), 67 & rLyi(npro,nflow), dxidx(npro,nsd,nsd), 68 & tau(npro,5), giju(npro,6), 69 & rTLS(npro), raLS(npro), 70 & ri(npro,nflow*(nsd+1)), rmi(npro,nflow*(nsd+1)), 71 & divqi(npro,nflow-1), stiff(npro,3*nflow,3*nflow), 72 & EGmass(npro,nedof,nedof),shape(npro,nshl), 73 & shg(npro,nshl,nsd), WdetJ(npro), 74 & PTau(npro,5,5), T(npro), 75 & pres(npro) 76c 77c local arrays 78c 79 dimension rLymi(npro,nflow), Atau(npro,nflow,nflow), 80 & A1tauA0(npro,nflow,nflow), A2tauA0(npro,nflow,nflow), 81 & A3tauA0(npro,nflow,nflow), fact(npro), 82 & A0inv(npro,15), dVdY(npro,15), 83 & compK(npro,10), ac1(npro), 84 & ac2(npro), ac3(npro) 85c 86 real*8 rerrl(npro,nshl,6), tmp(npro), tmp1(npro) 87 ttim(24) = ttim(24) - secs(0.0) 88c 89c 90c last step to the least squares is adding the time term. So that we 91c only have to localize one vector for each Krylov vector the modified 92c residual is quite different from the total residual. 93c 94c 95c the modified residual 96c 97 fct1=almi/gami/alfi*dtgl 98c 99c add possibility of not including time term 100c 101c if(idiff.ne.-1) 102 103 if(ires.ne.1) rLymi = rLyi + fct1*dui 104c 105 if(ires.eq.1 .or. ires .eq. 3) then 106c rLymi = rLyi 107 108 rLyi(:,1) = rLyi(:,1) 109 & + A0(:,1,1)*aci(:,1) 110c & + A0(:,1,2)*aci(:,2) 111c & + A0(:,1,3)*aci(:,3) 112c & + A0(:,1,4)*aci(:,4) 113 & + A0(:,1,5)*aci(:,5) 114c 115 rLyi(:,2) = rLyi(:,2) 116 & + A0(:,2,1)*aci(:,1) 117 & + A0(:,2,2)*aci(:,2) 118c & + A0(:,2,3)*aci(:,3) 119c & + A0(:,2,4)*aci(:,4) 120 & + A0(:,2,5)*aci(:,5) 121c 122 rLyi(:,3) = rLyi(:,3) 123 & + A0(:,3,1)*aci(:,1) 124c & + A0(:,3,2)*aci(:,2) 125 & + A0(:,3,3)*aci(:,3) 126c & + A0(:,3,4)*aci(:,4) 127 & + A0(:,3,5)*aci(:,5) 128c 129 rLyi(:,4) = rLyi(:,4) 130 & + A0(:,4,1)*aci(:,1) 131c & + A0(:,4,2)*aci(:,2) 132c & + A0(:,4,3)*aci(:,3) 133 & + A0(:,4,4)*aci(:,4) 134 & + A0(:,4,5)*aci(:,5) 135c 136 rLyi(:,5) = rLyi(:,5) 137 & + A0(:,5,1)*aci(:,1) 138 & + A0(:,5,2)*aci(:,2) 139 & + A0(:,5,3)*aci(:,3) 140 & + A0(:,5,4)*aci(:,4) 141 & + A0(:,5,5)*aci(:,5) 142c 143 endif 144c 145c.... subtract div(q) from the least squares term 146c 147 if ((idiff >= 1).and.(ires==3 .or. ires==1)) then 148c 149 if (isurf.eq.zero) then 150 rLyi(:,2) = rLyi(:,2) - divqi(:,1) 151 rLyi(:,3) = rLyi(:,3) - divqi(:,2) 152 rLyi(:,4) = rLyi(:,4) - divqi(:,3) 153 rLyi(:,5) = rLyi(:,5) - divqi(:,4) 154 endif 155 endif 156c 157c.... -------------------> error calculation <----------------- 158c 159 if((ierrcalc.eq.1).and.(nitr.eq.iter)) then 160 do ia=1,nshl 161 tmp=shape(:,ia)*WdetJ(:) 162 tmp1=shape(:,ia)*Qwt(lcsyst,intp) 163 rerrl(:,ia,1) = rerrl(:,ia,1) + 164 & tmp1(:)*rLyi(:,1)*rLyi(:,1) 165 rerrl(:,ia,2) = rerrl(:,ia,2) + 166 & tmp1(:)*rLyi(:,2)*rLyi(:,2) 167 rerrl(:,ia,3) = rerrl(:,ia,3) + 168 & tmp1(:)*rLyi(:,3)*rLyi(:,3) 169 170 rerrl(:,ia,4) = rerrl(:,ia,4) + 171 & tmp(:)*divqi(:,1)*divqi(:,1) 172 rerrl(:,ia,5) = rerrl(:,ia,5) + 173 & tmp(:)*divqi(:,2)*divqi(:,2) 174 rerrl(:,ia,6) = rerrl(:,ia,6) + 175 & tmp(:)*divqi(:,3)*divqi(:,3) 176 enddo 177 endif 178c 179c 180c.... ---------------------------> Tau <----------------------------- 181c 182c.... calculate the tau matrix 183c 184c.... in the first incarnation we will ONLY have a diagonal tau here 185 186 if (itau .lt. 10) then ! diagonal tau 187 188 call e3tau (rho, cp, rmu, 189 & u1, u2, u3, 190 & con, dxidx, rLyi, 191 & rLymi, tau, rk, 192 & giju, rTLS, raLS, 193 & A0inv, dVdY, cv) 194 195 else 196 197c.... looks like need a non-diagonal, discribed in "A new ... III" 198c.... by Hughes and Mallet. Original work - developed diffusive (and as 199c.... well advective) correction factors for 1-D a-d equation w/ hier. b. 200 201 202 ac1(:) = aci(:,2) 203 ac2(:) = aci(:,3) 204 ac3(:) = aci(:,4) 205 206 call e3tau_nd (rho, cp, rmu, T, 207 & u1, u2, u3, 208 & ac1, ac2, ac3, 209 & con, dxidx, rLyi, 210 & rLymi, PTau, rk, 211 & giju, rTLS, raLS, 212 & cv, compK, pres, 213 & A0inv, dVdY) 214 215 endif 216 217 218 ttim(25) = ttim(25) + secs(0.0) 219c 220c Warning: to save space I have put the tau times the modified residual 221c in rLymi and the tau times the total residual back in rLyi 222c 223c 224c.... ----------------------------> RHS <---------------------------- 225c 226c.... calculate (A_i^T tau Ly) 227c 228 229 if(ires.ne.1) then 230c 231c A1 * Tau L(Y): to be hit on left with Na,x in e3wmlt 232c 233 rmi(:,1) = 234 & A1(:,1,1) * rLymi(:,1) 235 & + A1(:,1,2) * rLymi(:,2) 236c & + A1(:,1,3) * rLymi(:,3) 237c & + A1(:,1,4) * rLymi(:,4) 238 & + A1(:,1,5) * rLymi(:,5) 239 & + rmi(:,1) 240 rmi(:,2) = 241 & A1(:,2,1) * rLymi(:,1) 242 & + A1(:,2,2) * rLymi(:,2) 243c & + A1(:,2,3) * rLymi(:,3) 244c & + A1(:,2,4) * rLymi(:,4) 245 & + A1(:,2,5) * rLymi(:,5) 246 & + rmi(:,2) 247 rmi(:,3) = 248 & A1(:,3,1) * rLymi(:,1) 249 & + A1(:,3,2) * rLymi(:,2) 250 & + A1(:,3,3) * rLymi(:,3) 251c & + A1(:,3,4) * rLymi(:,4) 252 & + A1(:,3,5) * rLymi(:,5) 253 & + rmi(:,3) 254 rmi(:,4) = 255 & A1(:,4,1) * rLymi(:,1) 256 & + A1(:,4,2) * rLymi(:,2) 257c & + A1(:,4,3) * rLymi(:,3) 258 & + A1(:,4,4) * rLymi(:,4) 259 & + A1(:,4,5) * rLymi(:,5) 260 & + rmi(:,4) 261 rmi(:,5) = 262 & A1(:,5,1) * rLymi(:,1) 263 & + A1(:,5,2) * rLymi(:,2) 264 & + A1(:,5,3) * rLymi(:,3) 265 & + A1(:,5,4) * rLymi(:,4) 266 & + A1(:,5,5) * rLymi(:,5) 267 & + rmi(:,5) 268c 269c A2 * Tau L(Y), to be hit on left with Na,y 270c 271 rmi(:,6) = 272 & A2(:,1,1) * rLymi(:,1) 273c & + A2(:,1,2) * rLymi(:,2) 274 & + A2(:,1,3) * rLymi(:,3) 275c & + A2(:,1,4) * rLymi(:,4) 276 & + A2(:,1,5) * rLymi(:,5) 277 & + rmi(:,6) 278 rmi(:,7) = 279 & A2(:,2,1) * rLymi(:,1) 280 & + A2(:,2,2) * rLymi(:,2) 281 & + A2(:,2,3) * rLymi(:,3) 282c & + A2(:,2,4) * rLymi(:,4) 283 & + A2(:,2,5) * rLymi(:,5) 284 & + rmi(:,7) 285 rmi(:,8) = 286 & A2(:,3,1) * rLymi(:,1) 287c & + A2(:,3,2) * rLymi(:,2) 288 & + A2(:,3,3) * rLymi(:,3) 289c & + A2(:,3,4) * rLymi(:,4) 290 & + A2(:,3,5) * rLymi(:,5) 291 & + rmi(:,8) 292 rmi(:,9) = 293 & A2(:,4,1) * rLymi(:,1) 294c & + A2(:,4,2) * rLymi(:,2) 295 & + A2(:,4,3) * rLymi(:,3) 296 & + A2(:,4,4) * rLymi(:,4) 297 & + A2(:,4,5) * rLymi(:,5) 298 & + rmi(:,9) 299 rmi(:,10) = 300 & A2(:,5,1) * rLymi(:,1) 301 & + A2(:,5,2) * rLymi(:,2) 302 & + A2(:,5,3) * rLymi(:,3) 303 & + A2(:,5,4) * rLymi(:,4) 304 & + A2(:,5,5) * rLymi(:,5) 305 & + rmi(:,10) 306c 307c A3 * Tau L(Y) to be hit on left with Na,z 308c 309 rmi(:,11) = 310 & A3(:,1,1) * rLymi(:,1) 311c & + A3(:,1,2) * rLymi(:,2) 312c & + A3(:,1,3) * rLymi(:,3) 313 & + A3(:,1,4) * rLymi(:,4) 314 & + A3(:,1,5) * rLymi(:,5) 315 & + rmi(:,11) 316 rmi(:,12) = 317 & A3(:,2,1) * rLymi(:,1) 318 & + A3(:,2,2) * rLymi(:,2) 319c & + A3(:,2,3) * rLymi(:,3) 320 & + A3(:,2,4) * rLymi(:,4) 321 & + A3(:,2,5) * rLymi(:,5) 322 & + rmi(:,12) 323 rmi(:,13) = 324 & A3(:,3,1) * rLymi(:,1) 325c & + A3(:,3,2) * rLymi(:,2) 326 & + A3(:,3,3) * rLymi(:,3) 327 & + A3(:,3,4) * rLymi(:,4) 328 & + A3(:,3,5) * rLymi(:,5) 329 & + rmi(:,13) 330 rmi(:,14) = 331 & A3(:,4,1) * rLymi(:,1) 332c & + A3(:,4,2) * rLymi(:,2) 333c & + A3(:,4,3) * rLymi(:,3) 334 & + A3(:,4,4) * rLymi(:,4) 335 & + A3(:,4,5) * rLymi(:,5) 336 & + rmi(:,14) 337 rmi(:,15) = 338 & A3(:,5,1) * rLymi(:,1) 339 & + A3(:,5,2) * rLymi(:,2) 340 & + A3(:,5,3) * rLymi(:,3) 341 & + A3(:,5,4) * rLymi(:,4) 342 & + A3(:,5,5) * rLymi(:,5) 343 & + rmi(:,15) 344 endif !ires.ne.1 345 346c 347c same thing for the real residual 348c 349 if(ires.eq.3 .or. ires .eq. 1) then ! we need the total residual 350 ri(:,1) = 351 & A1(:,1,1) * rLyi(:,1) 352 & + A1(:,1,2) * rLyi(:,2) 353c & + A1(:,1,3) * rLyi(:,3) 354c & + A1(:,1,4) * rLyi(:,4) 355 & + A1(:,1,5) * rLyi(:,5) 356 & + ri(:,1) 357 ri(:,2) = 358 & A1(:,2,1) * rLyi(:,1) 359 & + A1(:,2,2) * rLyi(:,2) 360c & + A1(:,2,3) * rLyi(:,3) 361c & + A1(:,2,4) * rLyi(:,4) 362 & + A1(:,2,5) * rLyi(:,5) 363 & + ri(:,2) 364 ri(:,3) = 365 & A1(:,3,1) * rLyi(:,1) 366 & + A1(:,3,2) * rLyi(:,2) 367 & + A1(:,3,3) * rLyi(:,3) 368c & + A1(:,3,4) * rLyi(:,4) 369 & + A1(:,3,5) * rLyi(:,5) 370 & + ri(:,3) 371 ri(:,4) = 372 & A1(:,4,1) * rLyi(:,1) 373 & + A1(:,4,2) * rLyi(:,2) 374c & + A1(:,4,3) * rLyi(:,3) 375 & + A1(:,4,4) * rLyi(:,4) 376 & + A1(:,4,5) * rLyi(:,5) 377 & + ri(:,4) 378 ri(:,5) = 379 & A1(:,5,1) * rLyi(:,1) 380 & + A1(:,5,2) * rLyi(:,2) 381 & + A1(:,5,3) * rLyi(:,3) 382 & + A1(:,5,4) * rLyi(:,4) 383 & + A1(:,5,5) * rLyi(:,5) 384 & + ri(:,5) 385c 386 ri(:,6) = 387 & A2(:,1,1) * rLyi(:,1) 388c & + A2(:,1,2) * rLyi(:,2) 389 & + A2(:,1,3) * rLyi(:,3) 390c & + A2(:,1,4) * rLyi(:,4) 391 & + A2(:,1,5) * rLyi(:,5) 392 & + ri(:,6) 393 ri(:,7) = 394 & A2(:,2,1) * rLyi(:,1) 395 & + A2(:,2,2) * rLyi(:,2) 396 & + A2(:,2,3) * rLyi(:,3) 397c & + A2(:,2,4) * rLyi(:,4) 398 & + A2(:,2,5) * rLyi(:,5) 399 & + ri(:,7) 400 ri(:,8) = 401 & A2(:,3,1) * rLyi(:,1) 402c & + A2(:,3,2) * rLyi(:,2) 403 & + A2(:,3,3) * rLyi(:,3) 404c & + A2(:,3,4) * rLyi(:,4) 405 & + A2(:,3,5) * rLyi(:,5) 406 & + ri(:,8) 407 ri(:,9) = 408 & A2(:,4,1) * rLyi(:,1) 409c & + A2(:,4,2) * rLyi(:,2) 410 & + A2(:,4,3) * rLyi(:,3) 411 & + A2(:,4,4) * rLyi(:,4) 412 & + A2(:,4,5) * rLyi(:,5) 413 & + ri(:,9) 414 ri(:,10) = 415 & A2(:,5,1) * rLyi(:,1) 416 & + A2(:,5,2) * rLyi(:,2) 417 & + A2(:,5,3) * rLyi(:,3) 418 & + A2(:,5,4) * rLyi(:,4) 419 & + A2(:,5,5) * rLyi(:,5) 420 & + ri(:,10) 421 ri(:,11) = 422 & A3(:,1,1) * rLyi(:,1) 423c & + A3(:,1,2) * rLyi(:,2) 424c & + A3(:,1,3) * rLyi(:,3) 425 & + A3(:,1,4) * rLyi(:,4) 426 & + A3(:,1,5) * rLyi(:,5) 427 & + ri(:,11) 428 ri(:,12) = 429 & A3(:,2,1) * rLyi(:,1) 430 & + A3(:,2,2) * rLyi(:,2) 431c & + A3(:,2,3) * rLyi(:,3) 432 & + A3(:,2,4) * rLyi(:,4) 433 & + A3(:,2,5) * rLyi(:,5) 434 & + ri(:,12) 435 ri(:,13) = 436 & A3(:,3,1) * rLyi(:,1) 437c & + A3(:,3,2) * rLyi(:,2) 438 & + A3(:,3,3) * rLyi(:,3) 439 & + A3(:,3,4) * rLyi(:,4) 440 & + A3(:,3,5) * rLyi(:,5) 441 & + ri(:,13) 442 ri(:,14) = 443 & A3(:,4,1) * rLyi(:,1) 444c & + A3(:,4,2) * rLyi(:,2) 445c & + A3(:,4,3) * rLyi(:,3) 446 & + A3(:,4,4) * rLyi(:,4) 447 & + A3(:,4,5) * rLyi(:,5) 448 & + ri(:,14) 449 ri(:,15) = 450 & A3(:,5,1) * rLyi(:,1) 451 & + A3(:,5,2) * rLyi(:,2) 452 & + A3(:,5,3) * rLyi(:,3) 453 & + A3(:,5,4) * rLyi(:,4) 454 & + A3(:,5,5) * rLyi(:,5) 455 & + ri(:,15) 456c 457 endif ! for ires=3 or 1 458 459c 460c.... ----------------------------> LHS <---------------------------- 461c 462 if (lhs .eq. 1) then 463c 464c.... calculate (Atau) <-- (A_1 tau) (Recall that we are using a 465c diagonal tau here) 466c 467 468 if (itau.lt.10) then 469 470 do i = 1, nflow 471 Atau(:,i,1) = A1(:,i,1)*tau(:,1) 472 Atau(:,i,2) = A1(:,i,2)*tau(:,2) 473 Atau(:,i,3) = A1(:,i,3)*tau(:,2) 474 Atau(:,i,4) = A1(:,i,4)*tau(:,2) 475 Atau(:,i,5) = A1(:,i,5)*tau(:,3) 476 enddo 477 478 else 479 480 Atau = zero 481 do i = 1, nflow 482 do j = 1, nflow 483 do k = 1, nflow 484 Atau(:,i,j) =Atau(:,i,j) + A1(:,i,k)*PTau(:,k,j) 485 enddo 486 enddo 487 enddo 488 489 endif 490c 491c.... calculate (A_1 tau A_0) (for L.S. time term of EGmass) 492c 493 do j = 1, nflow 494 do i = 1, nflow 495 A1tauA0(:,i,j) = 496 & Atau(:,i,1)*A0(:,1,j) + 497 & Atau(:,i,2)*A0(:,2,j) + 498 & Atau(:,i,3)*A0(:,3,j) + 499 & Atau(:,i,4)*A0(:,4,j) + 500 & Atau(:,i,5)*A0(:,5,j) 501 enddo 502 enddo 503c 504c.... add (A_1 tau A_1) to stiff [1,1] 505c 506 do j = 1, nflow 507 do i = 1, nflow 508 stiff(:,i,j) = stiff(:,i,j) + ( 509 & Atau(:,i,1)*A1(:,1,j) 510 & + Atau(:,i,2)*A1(:,2,j) 511 & + Atau(:,i,3)*A1(:,3,j) 512 & + Atau(:,i,4)*A1(:,4,j) 513 & + Atau(:,i,5)*A1(:,5,j) 514 & ) 515 enddo 516 enddo 517c 518c.... add (A_1 tau A_2) to stiff [1,2] 519c 520 do j = 1, nflow 521 do i = 1, nflow 522 stiff(:,i,j+5) = stiff(:,i,j+5) + ( 523 & Atau(:,i,1)*A2(:,1,j) 524 & + Atau(:,i,2)*A2(:,2,j) 525 & + Atau(:,i,3)*A2(:,3,j) 526 & + Atau(:,i,4)*A2(:,4,j) 527 & + Atau(:,i,5)*A2(:,5,j) 528 & ) 529 enddo 530 enddo 531c 532c.... add (A_1 tau A_3) to stiff [1,3] 533c 534 do j = 1, nflow 535 do i = 1, nflow 536 stiff(:,i,j+10) = stiff(:,i,j+10) + ( 537 & Atau(:,i,1)*A3(:,1,j) 538 & + Atau(:,i,2)*A3(:,2,j) 539 & + Atau(:,i,3)*A3(:,3,j) 540 & + Atau(:,i,4)*A3(:,4,j) 541 & + Atau(:,i,5)*A3(:,5,j) 542 & ) 543 enddo 544 enddo 545c 546c.... calculate (Atau) <-- (A_2 tau) (Recall that we are using a 547c diagonal tau here) 548c 549 if (itau.lt.10) then 550 551 do i = 1, nflow 552 Atau(:,i,1) = A2(:,i,1)*tau(:,1) 553 Atau(:,i,2) = A2(:,i,2)*tau(:,2) 554 Atau(:,i,3) = A2(:,i,3)*tau(:,2) 555 Atau(:,i,4) = A2(:,i,4)*tau(:,2) 556 Atau(:,i,5) = A2(:,i,5)*tau(:,3) 557 enddo 558 559 else 560 Atau = zero 561 do i = 1, nflow 562 do j = 1, nflow 563 do k = 1, nflow 564 Atau(:,i,j) = Atau(:,i,j) + A2(:,i,k)*PTau(:,k,j) 565 enddo 566 enddo 567 enddo 568 569 endif 570c 571c.... calculate (A_2 tau A_0) (for L.S. time term of EGmass) 572c 573 do j = 1, nflow 574 do i = 1, nflow 575 A2tauA0(:,i,j) = 576 & Atau(:,i,1)*A0(:,1,j) + 577 & Atau(:,i,2)*A0(:,2,j) + 578 & Atau(:,i,3)*A0(:,3,j) + 579 & Atau(:,i,4)*A0(:,4,j) + 580 & Atau(:,i,5)*A0(:,5,j) 581 enddo 582 enddo 583c 584c.... add (A_2 tau A_1) to stiff [2,1] 585c 586 do j = 1, nflow 587 do i = 1, nflow 588 stiff(:,i+5,j) = stiff(:,i+5,j) + ( 589 & Atau(:,i,1)*A1(:,1,j) 590 & + Atau(:,i,2)*A1(:,2,j) 591 & + Atau(:,i,3)*A1(:,3,j) 592 & + Atau(:,i,4)*A1(:,4,j) 593 & + Atau(:,i,5)*A1(:,5,j) 594 & ) 595 enddo 596 enddo 597c 598c.... add (A_2 tau A_2) to stiff [2,2] 599c 600 do j = 1, nflow 601 do i = 1, nflow 602 stiff(:,i+5,j+5) = stiff(:,i+5,j+5) + ( 603 & Atau(:,i,1)*A2(:,1,j) 604 & + Atau(:,i,2)*A2(:,2,j) 605 & + Atau(:,i,3)*A2(:,3,j) 606 & + Atau(:,i,4)*A2(:,4,j) 607 & + Atau(:,i,5)*A2(:,5,j) 608 & ) 609 enddo 610 enddo 611c 612c.... add (A_2 tau A_3) to stiff [2,3] 613c 614 do j = 1, nflow 615 do i = 1, nflow 616 stiff(:,i+5,j+10) = stiff(:,i+5,j+10) + ( 617 & Atau(:,i,1)*A3(:,1,j) 618 & + Atau(:,i,2)*A3(:,2,j) 619 & + Atau(:,i,3)*A3(:,3,j) 620 & + Atau(:,i,4)*A3(:,4,j) 621 & + Atau(:,i,5)*A3(:,5,j) 622 & ) 623 enddo 624 enddo 625c 626c.... calculate (Atau) <-- (A_3 tau) (Recall that we are using a 627c diagonal tau here) 628c 629 if (itau.lt.10) then 630 631 do i = 1, nflow 632 Atau(:,i,1) = A3(:,i,1)*tau(:,1) 633 Atau(:,i,2) = A3(:,i,2)*tau(:,2) 634 Atau(:,i,3) = A3(:,i,3)*tau(:,2) 635 Atau(:,i,4) = A3(:,i,4)*tau(:,2) 636 Atau(:,i,5) = A3(:,i,5)*tau(:,3) 637 enddo 638 639 else 640 Atau = zero 641 do i = 1, nflow 642 do j = 1, nflow 643 do k = 1, nflow 644 Atau(:,i,j) = Atau(:,i,j) + A3(:,i,k)*PTau(:,k,j) 645 enddo 646 enddo 647 enddo 648 649 endif 650c 651c.... calculate (A_3 tau A_0) (for L.S. time term of EGmass) 652c 653 do j = 1, nflow 654 do i = 1, nflow 655 A3tauA0(:,i,j) = 656 & Atau(:,i,1)*A0(:,1,j) + 657 & Atau(:,i,2)*A0(:,2,j) + 658 & Atau(:,i,3)*A0(:,3,j) + 659 & Atau(:,i,4)*A0(:,4,j) + 660 & Atau(:,i,5)*A0(:,5,j) 661 enddo 662 enddo 663c 664c.... add (A_3 tau A_1) to stiff [3,1] 665c 666 do j = 1, nflow 667 do i = 1, nflow 668 stiff(:,i+10,j) = stiff(:,i+10,j) + ( 669 & Atau(:,i,1)*A1(:,1,j) 670 & + Atau(:,i,2)*A1(:,2,j) 671 & + Atau(:,i,3)*A1(:,3,j) 672 & + Atau(:,i,4)*A1(:,4,j) 673 & + Atau(:,i,5)*A1(:,5,j) 674 & ) 675 enddo 676 enddo 677c 678c.... add (A_3 tau A_2) to stiff [3,2] 679c 680 do j = 1, nflow 681 do i = 1, nflow 682 stiff(:,i+10,j+5) = stiff(:,i+10,j+5) + ( 683 & Atau(:,i,1)*A2(:,1,j) 684 & + Atau(:,i,2)*A2(:,2,j) 685 & + Atau(:,i,3)*A2(:,3,j) 686 & + Atau(:,i,4)*A2(:,4,j) 687 & + Atau(:,i,5)*A2(:,5,j) 688 & ) 689 enddo 690 enddo 691c 692c.... add (A_3 tau A_3) to stiff [3,3] 693c 694 do j = 1, nflow 695 do i = 1, nflow 696 stiff(:,i+10,j+10) = stiff(:,i+10,j+10) + ( 697 & Atau(:,i,1)*A3(:,1,j) 698 & + Atau(:,i,2)*A3(:,2,j) 699 & + Atau(:,i,3)*A3(:,3,j) 700 & + Atau(:,i,4)*A3(:,4,j) 701 & + Atau(:,i,5)*A3(:,5,j) 702 & ) 703 enddo 704 enddo 705c 706c.... add least squares time term to the LHS tangent mass matrix 707c 708c 709c.... loop through rows (nodes i) 710c 711 do i = 1, nshl 712 i0 = nflow * (i - 1) 713c 714c.... first calculate (Atau) <-- (N_a,i A_i tau A_0) 715c ( use Atau to conserve space ) 716c 717 do idof = 1, nflow 718 do jdof = 1, nflow 719 Atau(:,idof,jdof) = 720 & shg(:,i,1) * A1tauA0(:,idof,jdof) + 721 & shg(:,i,2) * A2tauA0(:,idof,jdof) + 722 & shg(:,i,3) * A3tauA0(:,idof,jdof) 723 enddo 724 enddo 725c 726c.... loop through column nodes, add (N_a,i A_i tau N_b) to EGmass 727c 728 do j = 1, nshl 729 j0 = nflow * (j - 1) 730c 731c.... compute the factors 732c 733 fact = shape(:,j) * WdetJ * almi/gami/alfi*dtgl 734c 735c.... loop through d.o.f.'s 736c 737 do idof = 1, nflow 738 il = i0 + idof 739 740 EGmass(:,il,j0+1) = EGmass(:,il,j0+1) + 741 & fact * Atau(:,idof,1) 742 EGmass(:,il,j0+2) = EGmass(:,il,j0+2) + 743 & fact * Atau(:,idof,2) 744 EGmass(:,il,j0+3) = EGmass(:,il,j0+3) + 745 & fact * Atau(:,idof,3) 746 EGmass(:,il,j0+4) = EGmass(:,il,j0+4) + 747 & fact * Atau(:,idof,4) 748 EGmass(:,il,j0+5) = EGmass(:,il,j0+5) + 749 & fact * Atau(:,idof,5) 750 enddo 751c 752c.... end loop on column nodes 753c 754 enddo 755c 756c.... end loop on row nodes 757c 758 enddo 759c 760c.... end LHS computation 761c 762 endif 763 764 ttim(24) = ttim(24) + secs(0.0) 765c 766c.... return 767c 768 return 769 end 770c 771c 772c 773 subroutine e3LSSclr (A1t, A2t, A3t, 774 & rho, rmu, rTLSt, 775 & u1, u2, u3, 776 & rLyti, dxidx, raLSt, 777 & rti, rk, giju, 778 & acti, A0t, 779 & shape, shg, 780 & EGmasst, stifft, WdetJ, 781 & srcp) 782c 783c---------------------------------------------------------------------- 784c 785c This routine calculates the contribution of the least-squares 786c operator to the RHS vector and LHS tangent matrix. The temporary 787c results are put in ri. 788c 789c input: 790c A0t (npro) : A_0 791c A1t (npro) : A_1 792c A2t (npro) : A_2 793c A3t (npro) : A_3 794c acti (npro) : time-deriv. of Sclr 795c rho (npro) : density 796c rmu (npro) : molecular viscosity 797c rk (npro) : kinetic energy 798c u1 (npro) : x1-velocity component 799c u2 (npro) : x2-velocity component 800c u3 (npro) : x3-velocity component 801c rLyti (npro) : least-squares residual vector 802c dxidx (npro,nsd,nsd) : inverse of deformation gradient 803c taut (npro) : stability parameter 804c rLyti (npro) : convective portion of least-squares 805c residual vector 806c divqti (npro,1) : divergence of diffusive flux 807c shape (npro,nshl) : element shape functions 808c shg (npro,nshl,nsd) : global element shape function grads 809c WdetJ (npro) : weighted jacobian determinant 810c stifft (npro,nsd,nsd) : stiffness matrix 811c EGmasst(npro,nshape,nshape): partial mass matrix 812c 813c output: 814c rti (npro,nsd+1) : partial residual 815c EGmasst(npro,nshape,nshape): partial mass matrix 816c 817c 818c Zdenek Johan, Summer 1990. (Modified from e2ls.f) 819c Zdenek Johan, Winter 1991. (Fortran 90) 820c Kenneth Jansen, Winter 1997. Prim. Var. with Diag Tau 821c---------------------------------------------------------------------- 822c 823 include "common.h" 824 825c 826c passed arrays 827c 828 dimension A1t(npro), A2t(npro), 829 & A3t(npro), 830 & A0t(npro), rho(npro), 831 & acti(npro), rmu(npro), 832 & u1(npro), u2(npro), 833 & u3(npro), rk(npro), 834 & rLyti(npro), dxidx(npro,nsd,nsd), 835 & taut(npro), raLSt(npro), 836 & rti(npro,nsd+1), rTLSt(npro), 837 & stifft(npro,3,3), giju(npro,6), 838 & EGmasst(npro,nshape,nshape), 839 & shape(npro,nshl), 840 & shg(npro,nshl,nsd), WdetJ(npro), 841 & srcp(npro) 842c 843c local arrays 844c 845 dimension rLymti(npro), Ataut(npro), 846 & A1tautA0(npro), A2tautA0(npro), 847 & A3tautA0(npro), fact(npro) 848 849 ttim(24) = ttim(24) - tmr() 850c 851 if(ivart.lt.2) return 852c 853c last step to the least squares is adding the time term. So that we 854c only have to localize one vector for each Krylov vector the modified 855c residual is quite different from the total residual. 856c 857c 858c the modified residual 859c 860 fct1=almi/gami/alfi*dtgl 861c 862c add possibility of not including time term 863c 864c if(idiff.ne.-1) 865c rLymti = rLyti + fct1*duti 866 867 if((ires.eq.1 .or. ires .eq. 3).and. idiff.ne.-1) then 868 869 rLyti(:) = rLyti(:) + A0t(:)*acti(:) 870 871 endif 872c 873c.... subtract div(q) from the least squares term 874c 875c if ((idiff >= 1).and.(ires==3 .or. ires==1)) then 876c rLyi(:) = rLyi(:) - divqti(:) 877c endif 878c 879c.... ---------------------------> Tau <----------------------------- 880c 881c.... calculate the tau matrix 882c 883 884c 885c.... we will use the same tau as used for momentum equations here 886c 887 ttim(25) = ttim(25) - tmr() 888 889 call e3tauSclr(rho, rmu, A0t, 890 & u1, u2, u3, 891 & dxidx, rLyti, rLymti, 892 & taut, rk, raLSt, 893 & rTLSt, giju) 894 895 ttim(25) = ttim(25) + tmr() 896c 897c Warning: to save space I have put the tau times the modified residual 898c in rLymi and the tau times the total residual back in rLyi 899c 900c 901c.... ----------------------------> RHS <---------------------------- 902c 903c.... calculate (A_i^T tau Ly) 904c 905 906c if(ires.ne.1) then 907c 908c A1 * Tau L(Y): to be hit on left with Na,x in e3wmlt 909c 910c rmti(:,1) = A1t(:) * rLymti(:) 911c 912c 913c A2 * Tau L(Y), to be hit on left with Na,y 914c 915c rmti(:,2) = A2t(:) * rLymti(:) 916c 917c 918c A3 * Tau L(Y) to be hit on left with Na,z 919c 920c rmti(:,3) = A3t(:) * rLymti(:) 921c 922c endif !ires.ne.1 923 924c 925c same thing for the real residual 926c 927 if(ires.eq.3 .or. ires .eq. 1) then ! we need the total residual 928 rti(:,1) = rti(:,1) + A1t(:) * rLyti(:) 929 930 rti(:,2) = rti(:,2) + A2t(:) * rLyti(:) 931 932 rti(:,3) = rti(:,3) + A3t(:) * rLyti(:) 933 934 endif ! for ires=3 or 1 935 936c 937c.... ----------------------------> LHS <---------------------------- 938c 939 if (lhs .eq. 1) then 940c 941c 942c.... calculate (Atau) <-- (A_1 tau) 943c 944 945 Ataut(:) = A1t(:)*taut(:) 946 947c 948c.... calculate (A_1 tau (A_0-srcp)) (for L.S. time term of EGmass) 949c 950 951 A1tautA0(:) = Ataut(:)*(a0t(:)*fct1-srcp(:)) 952 953c 954c.... add (A_1 tau A_1) to stiff [1,1] 955c 956 957 stifft(:,1,1) = stifft(:,1,1) + Ataut(:)*A1t(:) 958c stifft(:,1,1) = Ataut(:)*A1t(:) 959c 960c.... add (A_1 tau A_2) to stiff [1,2] 961c 962 963 stifft(:,1,2) = stifft(:,1,2) + Ataut(:)*A2t(:) 964c stifft(:,1,2) = Ataut(:)*A2t(:) 965c 966c.... add (A_1 tau A_3) to stiff [1,3] 967c 968 969 stifft(:,1,3) = stifft(:,1,3) + Ataut(:)*A3t(:) 970c stifft(:,1,3) = Ataut(:)*A3t(:) 971c 972c.... calculate (Atau) <-- (A_2 tau) 973c 974 975 Ataut(:) = A2t(:)*taut(:) 976 977c 978c.... calculate (A_2 tau (A_0-srcp)) (for L.S. time term of EGmass) 979c 980 981 A2tautA0(:) = Ataut(:)*(a0t(:)*fct1-srcp(:)) 982 983c 984c.... add (A_2 tau A_1) to stiff [2,1] 985c 986 987 stifft(:,2,1) = stifft(:,1,2) 988c 989c.... add (A_2 tau A_2) to stiff [2,2] 990c 991 992 stifft(:,2,2) = stifft(:,2,2) + Ataut(:)*A2t(:) 993 994c 995c.... add (A_2 tau A_3) to stiff [2,3] 996c 997 998 stifft(:,2,3) = stifft(:,2,3) + Ataut(:)*A3t(:) 999 1000c 1001c.... calculate (Atau) <-- (A_3 tau) 1002c 1003 1004 Ataut(:) = A3t(:)*taut(:) 1005 1006c 1007c.... calculate (A_3 tau (A_0-srcp)) (for L.S. time term of EGmass) 1008c 1009 1010 A3tautA0(:) = Ataut(:)*(a0t(:)*fct1-srcp(:)) 1011 1012c 1013c.... add (A_3 tau A_1) to stiff [3,1] 1014c 1015 1016 stifft(:,3,1) = stifft(:,1,3) 1017 1018c 1019c.... add (A_3 tau A_2) to stiff [3,2] 1020c 1021 1022 stifft(:,3,2) = stifft(:,2,3) 1023 1024c 1025c.... add (A_3 tau A_3) to stiff [3,3] 1026c 1027 1028 stifft(:,3,3) = stifft(:,3,3) + Ataut(:)*A3t(:) 1029 1030c 1031c.... add least squares time term to the LHS tangent mass matrix 1032c 1033c 1034c.... loop through rows (nodes i) 1035c 1036 do ia = 1, nshl 1037c 1038c.... first calculate (Atau) <-- (N_a,i A_i tau A_0) 1039c ( use Atau to conserve space ) 1040c 1041 1042 Ataut(:) = 1043 & shg(:,ia,1) * A1tautA0(:) + 1044 & shg(:,ia,2) * A2tautA0(:) + 1045 & shg(:,ia,3) * A3tautA0(:) 1046 1047c 1048c.... loop through column nodes, add (N_a,i A_i tau N_b) to EGmass 1049c 1050 do jb = 1, nshl 1051 1052 fact = shape(:,jb) * WdetJ 1053 1054 EGmasst(:,ia,jb) = EGmasst(:,ia,jb) + fact * Ataut(:) 1055 1056c 1057c.... end loop on column nodes 1058c 1059 1060 enddo 1061c 1062c.... end loop on row nodes 1063c 1064 enddo 1065c 1066c.... end LHS computation 1067c 1068 endif 1069 1070 ttim(24) = ttim(24) + tmr() 1071c 1072c.... return 1073c 1074 return 1075 end 1076 1077