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