1 subroutine e3conv (g1yi, g2yi, g3yi, 2 & A1, A2, A3, 3 & rho, pres, T, 4 & ei, rk, u1, 5 & u2, u3, rLyi, 6 & ri, rmi, EGmass, 7 & shg, shape, WdetJ ) 8c 9c---------------------------------------------------------------------- 10c 11c This routine calculates the contribution of Galerkin part of the 12c Convective term (Time and Euler fluxes) to both RHS and LHS. 13c 14c input: 15c g1yi (npro,nflow) : grad-y in direction 1 16c g2yi (npro,nflow) : grad-y in direction 2 17c g3yi (npro,nflow) : grad-y in direction 3 18c A1 (npro,nflow,nflow) : A-1 19c A2 (npro,nflow,nflow) : A-2 20c A3 (npro,nflow,nflow) : A-3 21c rho (npro) : density 22c pres (npro) : pressure 23c T (npro) : temperature 24c ei (npro) : internal energy 25c rk (npro) : kinetic energy 26c u1 (npro) : x1-velocity component 27c u2 (npro) : x2-velocity component 28c u3 (npro) : x3-velocity component 29c shg (npro,nshl,nsd) : global grad's of shape functions 30c shape (npro,nshl) : element shape functions 31c WdetJ (npro) : weighted Jacobian determinant 32c 33c output: 34c rLyi (npro,nflow) : least-squares residual vector 35c ri (npro,nflow*(nsd+1)) : partial residual 36c rmi (npro,nflow*(nsd+1)) : partial modified residual 37c EGmass (npro,nedof,nedof) : partial LHS tangent matrix 38c 39c 40c Zdenek Johan, Summer 1990. (Modified from e2conv.f) 41c Zdenek Johan, Winter 1991. (Fortran 90) 42c Kenneth Jansen, Winter 1997 Primitive Variables 43c---------------------------------------------------------------------- 44c 45 include "common.h" 46c 47c passed arrays 48c 49 dimension g1yi(npro,nflow), g2yi(npro,nflow), 50 & g3yi(npro,nflow), 51 & A1(npro,nflow,nflow), 52 & A2(npro,nflow,nflow), A3(npro,nflow,nflow), 53 & rho(npro), pres(npro), 54 & T(npro), ei(npro), 55 & rk(npro), u1(npro), 56 & u2(npro), u3(npro), 57 & rLyi(npro,nflow), ri(npro,nflow*(nsd+1)), 58 & rmi(npro,nflow*(nsd+1)), EGmass(npro,nedof,nedof), 59 & shg(npro,nshl,nsd), shape(npro,nshl), 60 & WdetJ(npro) 61c 62c local arrays 63c 64 dimension AiNbi(npro,nflow,nflow), fact1(npro), 65 & fact2(npro), fact3(npro) 66 67 ttim(22) = ttim(22) - secs(0.0) 68 69c 70c.... ----------------------> RHS, Euler Flux <---------------------- 71c 72 if ((ires .eq. 1) .or. (ires .eq. 3)) then 73c 74c.... calculate integrated by part contribution of Euler flux (Galerkin) 75c 76 ri(:, 1) = (- u1) * rho 77 ri(:, 2) = (- u1) * rho * u1 - pres 78 ri(:, 3) = (- u1) * rho * u2 79 ri(:, 4) = (- u1) * rho * u3 80 ri(:, 5) = (- u1) * rho * (ei + rk) - u1 * pres 81c 82 ri(:, 6) = (- u2) * rho 83 ri(:, 7) = (- u2) * rho * u1 84 ri(:, 8) = (- u2) * rho * u2 - pres 85 ri(:, 9) = (- u2) * rho * u3 86 ri(:,10) = (- u2) * rho * (ei + rk) - u2 * pres 87c 88 ri(:,11) = (- u3) * rho 89 ri(:,12) = (- u3) * rho * u1 90 ri(:,13) = (- u3) * rho * u2 91 ri(:,14) = (- u3) * rho * u3 - pres 92 ri(:,15) = (- u3) * rho * (ei + rk) - u3 * pres 93c 94 flops = flops + 28*npro 95c 96 endif 97c 98c.... calculate ( A_i Y,i ) --> rLyi Commented out zeros of A matrices 99c 100 rLyi(:,1) = 101 & A1(:,1,1) * g1yi(:,1) 102 & + A1(:,1,2) * g1yi(:,2) 103c & + A1(:,1,3) * g1yi(:,3) 104c & + A1(:,1,4) * g1yi(:,4) 105 & + A1(:,1,5) * g1yi(:,5) 106 & + A2(:,1,1) * g2yi(:,1) 107c & + A2(:,1,2) * g2yi(:,2) 108 & + A2(:,1,3) * g2yi(:,3) 109c & + A2(:,1,4) * g2yi(:,4) 110 & + A2(:,1,5) * g2yi(:,5) 111 & + A3(:,1,1) * g3yi(:,1) 112c & + A3(:,1,2) * g3yi(:,2) 113c & + A3(:,1,3) * g3yi(:,3) 114 & + A3(:,1,4) * g3yi(:,4) 115 & + A3(:,1,5) * g3yi(:,5) 116 rLyi(:,2) = 117 & A1(:,2,1) * g1yi(:,1) 118 & + A1(:,2,2) * g1yi(:,2) 119c & + A1(:,2,3) * g1yi(:,3) 120c & + A1(:,2,4) * g1yi(:,4) 121 & + A1(:,2,5) * g1yi(:,5) 122 & + A2(:,2,1) * g2yi(:,1) 123 & + A2(:,2,2) * g2yi(:,2) 124 & + A2(:,2,3) * g2yi(:,3) 125c & + A2(:,2,4) * g2yi(:,4) 126 & + A2(:,2,5) * g2yi(:,5) 127 & + A3(:,2,1) * g3yi(:,1) 128 & + A3(:,2,2) * g3yi(:,2) 129c & + A3(:,2,3) * g3yi(:,3) 130 & + A3(:,2,4) * g3yi(:,4) 131 & + A3(:,2,5) * g3yi(:,5) 132 rLyi(:,3) = 133 & A1(:,3,1) * g1yi(:,1) 134 & + A1(:,3,2) * g1yi(:,2) 135 & + A1(:,3,3) * g1yi(:,3) 136c & + A1(:,3,4) * g1yi(:,4) 137 & + A1(:,3,5) * g1yi(:,5) 138 & + A2(:,3,1) * g2yi(:,1) 139c & + A2(:,3,2) * g2yi(:,2) 140 & + A2(:,3,3) * g2yi(:,3) 141c & + A2(:,3,4) * g2yi(:,4) 142 & + A2(:,3,5) * g2yi(:,5) 143 & + A3(:,3,1) * g3yi(:,1) 144c & + A3(:,3,2) * g3yi(:,2) 145 & + A3(:,3,3) * g3yi(:,3) 146 & + A3(:,3,4) * g3yi(:,4) 147 & + A3(:,3,5) * g3yi(:,5) 148 rLyi(:,4) = 149 & A1(:,4,1) * g1yi(:,1) 150 & + A1(:,4,2) * g1yi(:,2) 151c & + A1(:,4,3) * g1yi(:,3) 152 & + A1(:,4,4) * g1yi(:,4) 153 & + A1(:,4,5) * g1yi(:,5) 154 & + A2(:,4,1) * g2yi(:,1) 155c & + A2(:,4,2) * g2yi(:,2) 156 & + A2(:,4,3) * g2yi(:,3) 157 & + A2(:,4,4) * g2yi(:,4) 158 & + A2(:,4,5) * g2yi(:,5) 159 & + A3(:,4,1) * g3yi(:,1) 160c & + A3(:,4,2) * g3yi(:,2) 161c & + A3(:,4,3) * g3yi(:,3) 162 & + A3(:,4,4) * g3yi(:,4) 163 & + A3(:,4,5) * g3yi(:,5) 164 rLyi(:,5) = 165 & A1(:,5,1) * g1yi(:,1) 166 & + A1(:,5,2) * g1yi(:,2) 167 & + A1(:,5,3) * g1yi(:,3) 168 & + A1(:,5,4) * g1yi(:,4) 169 & + A1(:,5,5) * g1yi(:,5) 170 & + A2(:,5,1) * g2yi(:,1) 171 & + A2(:,5,2) * g2yi(:,2) 172 & + A2(:,5,3) * g2yi(:,3) 173 & + A2(:,5,4) * g2yi(:,4) 174 & + A2(:,5,5) * g2yi(:,5) 175 & + A3(:,5,1) * g3yi(:,1) 176 & + A3(:,5,2) * g3yi(:,2) 177 & + A3(:,5,3) * g3yi(:,3) 178 & + A3(:,5,4) * g3yi(:,4) 179 & + A3(:,5,5) * g3yi(:,5) 180c 181c.... add contribution to rmi 182c 183 if ((ires .eq. 2) .or. (ires .eq. 3)) 184 & rmi(:,16:20) = rLyi ! modified residual uses non i.b.p form of conv. 185c 186c.... ----------------------> LHS <----------------------- 187c 188 if (lhs .eq. 1) then 189c 190c.... loop through the columns 191c 192 do j = 1, nshl 193 j0 = nflow * (j - 1) 194c 195c.... compute some useful factors 196c 197 fact1 = WdetJ * shg(:,j,1) 198 fact2 = WdetJ * shg(:,j,2) 199 fact3 = WdetJ * shg(:,j,3) 200c 201c.... first compute (A_i N_b,i) 202c 203 AiNbi(:,1,1) = 204 & fact1 * A1(:,1,1) 205 & + fact2 * A2(:,1,1) 206 & + fact3 * A3(:,1,1) 207 AiNbi(:,1,2) = 208 & fact1 * A1(:,1,2) 209 & + fact2 * A2(:,1,2) 210 & + fact3 * A3(:,1,2) 211 AiNbi(:,1,3) = 212 & fact1 * A1(:,1,3) 213 & + fact2 * A2(:,1,3) 214 & + fact3 * A3(:,1,3) 215 AiNbi(:,1,4) = 216 & fact1 * A1(:,1,4) 217 & + fact2 * A2(:,1,4) 218 & + fact3 * A3(:,1,4) 219 AiNbi(:,1,5) = 220 & fact1 * A1(:,1,5) 221 & + fact2 * A2(:,1,5) 222 & + fact3 * A3(:,1,5) 223c 224 AiNbi(:,2,1) = 225 & fact1 * A1(:,2,1) 226 & + fact2 * A2(:,2,1) 227 & + fact3 * A3(:,2,1) 228 AiNbi(:,2,2) = 229 & fact1 * A1(:,2,2) 230 & + fact2 * A2(:,2,2) 231 & + fact3 * A3(:,2,2) 232 AiNbi(:,2,3) = 233 & fact1 * A1(:,2,3) 234 & + fact2 * A2(:,2,3) 235 & + fact3 * A3(:,2,3) 236 AiNbi(:,2,4) = 237 & fact1 * A1(:,2,4) 238 & + fact2 * A2(:,2,4) 239 & + fact3 * A3(:,2,4) 240 AiNbi(:,2,5) = 241 & fact1 * A1(:,2,5) 242 & + fact2 * A2(:,2,5) 243 & + fact3 * A3(:,2,5) 244c 245 AiNbi(:,3,1) = 246 & fact1 * A1(:,3,1) 247 & + fact2 * A2(:,3,1) 248 & + fact3 * A3(:,3,1) 249 AiNbi(:,3,2) = 250 & fact1 * A1(:,3,2) 251 & + fact2 * A2(:,3,2) 252 & + fact3 * A3(:,3,2) 253 AiNbi(:,3,3) = 254 & fact1 * A1(:,3,3) 255 & + fact2 * A2(:,3,3) 256 & + fact3 * A3(:,3,3) 257 AiNbi(:,3,4) = 258 & fact1 * A1(:,3,4) 259 & + fact2 * A2(:,3,4) 260 & + fact3 * A3(:,3,4) 261 AiNbi(:,3,5) = 262 & fact1 * A1(:,3,5) 263 & + fact2 * A2(:,3,5) 264 & + fact3 * A3(:,3,5) 265c 266 AiNbi(:,4,1) = 267 & fact1 * A1(:,4,1) 268 & + fact2 * A2(:,4,1) 269 & + fact3 * A3(:,4,1) 270 AiNbi(:,4,2) = 271 & fact1 * A1(:,4,2) 272 & + fact2 * A2(:,4,2) 273 & + fact3 * A3(:,4,2) 274 AiNbi(:,4,3) = 275 & fact1 * A1(:,4,3) 276 & + fact2 * A2(:,4,3) 277 & + fact3 * A3(:,4,3) 278 AiNbi(:,4,4) = 279 & fact1 * A1(:,4,4) 280 & + fact2 * A2(:,4,4) 281 & + fact3 * A3(:,4,4) 282 AiNbi(:,4,5) = 283 & fact1 * A1(:,4,5) 284 & + fact2 * A2(:,4,5) 285 & + fact3 * A3(:,4,5) 286c 287 AiNbi(:,5,1) = 288 & fact1 * A1(:,5,1) 289 & + fact2 * A2(:,5,1) 290 & + fact3 * A3(:,5,1) 291 AiNbi(:,5,2) = 292 & fact1 * A1(:,5,2) 293 & + fact2 * A2(:,5,2) 294 & + fact3 * A3(:,5,2) 295 AiNbi(:,5,3) = 296 & fact1 * A1(:,5,3) 297 & + fact2 * A2(:,5,3) 298 & + fact3 * A3(:,5,3) 299 AiNbi(:,5,4) = 300 & fact1 * A1(:,5,4) 301 & + fact2 * A2(:,5,4) 302 & + fact3 * A3(:,5,4) 303 AiNbi(:,5,5) = 304 & fact1 * A1(:,5,5) 305 & + fact2 * A2(:,5,5) 306 & + fact3 * A3(:,5,5) 307c 308c.... now loop through the row nodes and add (N_a A_i N_b,i) to 309c the tangent matrix. 310c 311 do i = 1, nshl 312 i0 = nflow * (i - 1) 313c 314c.... loop through dof's 315c 316 do jdof = 1, nflow 317 jl = j0 + jdof 318 319 EGmass(:,i0+1,jl) = EGmass(:,i0+1,jl) + 320 & shape(:,i) * AiNbi(:,1,jdof) 321 322 EGmass(:,i0+2,jl) = EGmass(:,i0+2,jl) + 323 & shape(:,i) * AiNbi(:,2,jdof) 324 325 EGmass(:,i0+3,jl) = EGmass(:,i0+3,jl) + 326 & shape(:,i) * AiNbi(:,3,jdof) 327 328 EGmass(:,i0+4,jl) = EGmass(:,i0+4,jl) + 329 & shape(:,i) * AiNbi(:,4,jdof) 330 331 EGmass(:,i0+5,jl) = EGmass(:,i0+5,jl) + 332 & shape(:,i) * AiNbi(:,5,jdof) 333 enddo 334c 335c.... end loop on rows 336c 337 enddo 338c 339c.... end loop on columns 340c 341 enddo 342c 343c.... end of LHS tangent matrix computation 344c 345 endif 346 347 ttim(22) = ttim(22) + secs(0.0) 348c 349c.... return 350c 351 return 352 end 353c 354c 355c 356 subroutine e3convSclr (g1yti, g2yti, g3yti, 357 & A1t, A2t, A3t, 358 & rho, u1, Sclr, 359 & u2, u3, rLyti, 360 & rti, rmti, EGmasst, 361 & shg, shape, WdetJ) 362c 363c---------------------------------------------------------------------- 364c 365c This routine calculates the contribution of Galerkin part of the 366c Convective term (Time and Euler fluxes) to both RHS and LHS. 367c 368c input: 369c Sclr (npro) : Scalar variable 370c g1yti (npro) : grad-y in direction 1 371c g2yti (npro) : grad-y in direction 2 372c g3yti (npro) : grad-y in direction 3 373c A1t (npro) : A-1 374c A2t (npro) : A-2 375c A3t (npro) : A-3 376c rho (npro) : density 377c u1 (npro) : x1-velocity component 378c u2 (npro) : x2-velocity component 379c u3 (npro) : x3-velocity component 380c shg (npro,nshl,nsd) : global grad's of shape functions 381c shape (npro,nshl) : element shape functions 382c WdetJ (npro) : weighted Jacobian determinant 383c 384c output: 385c rLyti (npro) : least-squares residual vector 386c rti (npro,nsd+1) : partial residual 387c rmti (npro,nsd+1) : partial modified residual 388c EGmasst (npro,nshape,nshape): partial LHS tangent matrix 389c 390c 391c Zdenek Johan, Summer 1990. (Modified from e2conv.f) 392c Zdenek Johan, Winter 1991. (Fortran 90) 393c Kenneth Jansen, Winter 1997 Primitive Variables 394c---------------------------------------------------------------------- 395c 396 include "common.h" 397c 398c passed arrays 399c 400 dimension g1yti(npro), g2yti(npro), 401 & g3yti(npro), Sclr(npro), 402 & A1t(npro), 403 & A2t(npro), A3t(npro), 404 & rho(npro), u1(npro), 405 & u2(npro), u3(npro), 406 & rLyti(npro), rti(npro,nsd+1), 407 & rmti(npro,nsd+1), EGmasst(npro,nshape,nshape), 408 & shg(npro,nshl,nsd), shape(npro,nshl), 409 & WdetJ(npro) 410c 411c local arrays 412c 413 dimension AitNbi(npro) 414 415 ttim(22) = ttim(22) - tmr(0.0) 416c 417c.... ----------------------> RHS, Euler Flux <---------------------- 418c 419 if ((ires .eq. 1) .or. (ires .eq. 3)) then 420c 421c.... calculate integrated by part contribution of Euler flux (Galerkin) 422c 423 if (iconvsclr.eq.2) then ! convective form 424c 425 rti(:, 4) = rti(:,4) + ( u1) * g1yti(:) 426 & + ( u2) * g2yti(:) 427 & + ( u3) * g3yti(:) 428c 429 else ! conservative form 430c 431 rti(:, 1) = rti(:,1) + (- u1) * rho * Sclr 432 rti(:, 2) = rti(:,2) + (- u2) * rho * Sclr 433 rti(:, 3) = rti(:,3) + (- u3) * rho * Sclr 434c 435 endif 436 437 flops = flops + 28*npro 438 439 endif 440c 441c.... calculate ( A_i Y,i ) --> rLyi 442c 443 rLyti(:) = rLyti(:) 444 & + A1t(:) * g1yti(:) 445 & + A2t(:) * g2yti(:) 446 & + A3t(:) * g3yti(:) 447 448c 449c.... add contribution to rmi 450c 451c if ((ires .eq. 2) .or. (ires .eq. 3)) 452c & rmi(:,16:20) = rLyi ! modified residual uses non i.b.p form of conv 453c 454c.... ----------------------> LHS <----------------------- 455c 456 if (lhs .eq. 1) then 457c 458c.... loop through the columns 459c 460 do j = 1, nshl 461 462c 463c.... first compute (A_i N_b,i) 464c 465 AitNbi(:) = 466 & WdetJ * shg(:,j,1) * A1t(:) 467 & + WdetJ * shg(:,j,2) * A2t(:) 468 & + WdetJ * shg(:,j,3) * A3t(:) 469 470c 471c.... now loop through the rows and add (N_a A_i N_b,i) to 472c the tangent matrix. 473c 474 do i = 1, nshl 475 476 EGmasst(:,i,j) = EGmasst(:,i,j) + shape(:,i) * AitNbi(:) 477 478 479c 480c.... end loop on rows 481c 482 enddo 483c 484c.... end loop on columns 485c 486 enddo 487c 488c.... end of LHS tangent matrix computation 489c 490 endif 491 492 ttim(22) = ttim(22) + tmr() 493c 494c.... return 495c 496 return 497 end 498 499