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 ) 8! 9!---------------------------------------------------------------------- 10! 11! This routine calculates the contribution of Galerkin part of the 12! Convective term (Time and Euler fluxes) to both RHS and LHS. 13! 14! input: 15! g1yi (npro,nflow) : grad-y in direction 1 16! g2yi (npro,nflow) : grad-y in direction 2 17! g3yi (npro,nflow) : grad-y in direction 3 18! A1 (npro,nflow,nflow) : A-1 19! A2 (npro,nflow,nflow) : A-2 20! A3 (npro,nflow,nflow) : A-3 21! rho (npro) : density 22! pres (npro) : pressure 23! T (npro) : temperature 24! ei (npro) : internal energy 25! rk (npro) : kinetic energy 26! u1 (npro) : x1-velocity component 27! u2 (npro) : x2-velocity component 28! u3 (npro) : x3-velocity component 29! shg (npro,nshl,nsd) : global grad's of shape functions 30! shape (npro,nshl) : element shape functions 31! WdetJ (npro) : weighted Jacobian determinant 32! 33! output: 34! rLyi (npro,nflow) : least-squares residual vector 35! ri (npro,nflow*(nsd+1)) : partial residual 36! rmi (npro,nflow*(nsd+1)) : partial modified residual 37! EGmass (npro,nedof,nedof) : partial LHS tangent matrix 38! 39! 40! Zdenek Johan, Summer 1990. (Modified from e2conv.f) 41! Zdenek Johan, Winter 1991. (Fortran 90) 42! Kenneth Jansen, Winter 1997 Primitive Variables 43!---------------------------------------------------------------------- 44! 45 include "common.h" 46! 47! passed arrays 48! 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) 61! 62! local arrays 63! 64 dimension AiNbi(npro,nflow,nflow), fact1(npro), 65 & fact2(npro), fact3(npro) 66 67 ttim(22) = ttim(22) - secs(0.0) 68 69! 70!.... ----------------------> RHS, Euler Flux <---------------------- 71! 72 if ((ires .eq. 1) .or. (ires .eq. 3)) then 73! 74!.... calculate integrated by part contribution of Euler flux (Galerkin) 75! 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 81! 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 87! 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 93! 94! flops = flops + 28*npro 95! 96 endif 97! 98!.... calculate ( A_i Y,i ) --> rLyi Commented out zeros of A matrices 99! 100 rLyi(:,1) = 101 & A1(:,1,1) * g1yi(:,1) 102 & + A1(:,1,2) * g1yi(:,2) 103! & + A1(:,1,3) * g1yi(:,3) 104! & + A1(:,1,4) * g1yi(:,4) 105 & + A1(:,1,5) * g1yi(:,5) 106 & + A2(:,1,1) * g2yi(:,1) 107! & + A2(:,1,2) * g2yi(:,2) 108 & + A2(:,1,3) * g2yi(:,3) 109! & + A2(:,1,4) * g2yi(:,4) 110 & + A2(:,1,5) * g2yi(:,5) 111 & + A3(:,1,1) * g3yi(:,1) 112! & + A3(:,1,2) * g3yi(:,2) 113! & + 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) 119! & + A1(:,2,3) * g1yi(:,3) 120! & + 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) 125! & + A2(:,2,4) * g2yi(:,4) 126 & + A2(:,2,5) * g2yi(:,5) 127 & + A3(:,2,1) * g3yi(:,1) 128 & + A3(:,2,2) * g3yi(:,2) 129! & + 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) 136! & + A1(:,3,4) * g1yi(:,4) 137 & + A1(:,3,5) * g1yi(:,5) 138 & + A2(:,3,1) * g2yi(:,1) 139! & + A2(:,3,2) * g2yi(:,2) 140 & + A2(:,3,3) * g2yi(:,3) 141! & + A2(:,3,4) * g2yi(:,4) 142 & + A2(:,3,5) * g2yi(:,5) 143 & + A3(:,3,1) * g3yi(:,1) 144! & + 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) 151! & + A1(:,4,3) * g1yi(:,3) 152 & + A1(:,4,4) * g1yi(:,4) 153 & + A1(:,4,5) * g1yi(:,5) 154 & + A2(:,4,1) * g2yi(:,1) 155! & + 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) 160! & + A3(:,4,2) * g3yi(:,2) 161! & + 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) 180! 181!.... add contribution to rmi 182! 183 if ((ires .eq. 2) .or. (ires .eq. 3)) 184 & rmi(:,16:20) = rLyi ! modified residual uses non i.b.p form of conv. 185! 186!.... ----------------------> LHS <----------------------- 187! 188 if (lhs .eq. 1) then 189! 190!.... loop through the columns 191! 192 do j = 1, nshl 193 j0 = nflow * (j - 1) 194! 195!.... compute some useful factors 196! 197 fact1 = WdetJ * shg(:,j,1) 198 fact2 = WdetJ * shg(:,j,2) 199 fact3 = WdetJ * shg(:,j,3) 200! 201!.... first compute (A_i N_b,i) 202! 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) 223 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) 244 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) 265 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) 286 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) 307! 308!.... now loop through the row nodes and add (N_a A_i N_b,i) to 309! the tangent matrix. 310! 311 do i = 1, nshl 312 i0 = nflow * (i - 1) 313! 314!.... loop through dof's 315! 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 334! 335!.... end loop on rows 336! 337 enddo 338! 339!.... end loop on columns 340! 341 enddo 342! 343!.... end of LHS tangent matrix computation 344! 345 endif 346 347 ttim(22) = ttim(22) + secs(0.0) 348! 349!.... return 350! 351 return 352 end 353! 354! 355! 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) 362! 363!---------------------------------------------------------------------- 364! 365! This routine calculates the contribution of Galerkin part of the 366! Convective term (Time and Euler fluxes) to both RHS and LHS. 367! 368! input: 369! Sclr (npro) : Scalar variable 370! g1yti (npro) : grad-y in direction 1 371! g2yti (npro) : grad-y in direction 2 372! g3yti (npro) : grad-y in direction 3 373! A1t (npro) : A-1 374! A2t (npro) : A-2 375! A3t (npro) : A-3 376! rho (npro) : density 377! u1 (npro) : x1-velocity component 378! u2 (npro) : x2-velocity component 379! u3 (npro) : x3-velocity component 380! shg (npro,nshl,nsd) : global grad's of shape functions 381! shape (npro,nshl) : element shape functions 382! WdetJ (npro) : weighted Jacobian determinant 383! 384! output: 385! rLyti (npro) : least-squares residual vector 386! rti (npro,nsd+1) : partial residual 387! rmti (npro,nsd+1) : partial modified residual 388! EGmasst (npro,nshape,nshape): partial LHS tangent matrix 389! 390! 391! Zdenek Johan, Summer 1990. (Modified from e2conv.f) 392! Zdenek Johan, Winter 1991. (Fortran 90) 393! Kenneth Jansen, Winter 1997 Primitive Variables 394!---------------------------------------------------------------------- 395! 396 include "common.h" 397! 398! passed arrays 399! 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) 410! 411! local arrays 412! 413 dimension AitNbi(npro) 414 415 ttim(22) = ttim(22) - tmr() !tmr(0.0) 416! 417!.... ----------------------> RHS, Euler Flux <---------------------- 418! 419 if ((ires .eq. 1) .or. (ires .eq. 3)) then 420! 421!.... calculate integrated by part contribution of Euler flux (Galerkin) 422! 423 if (iconvsclr.eq.2) then ! convective form 424! 425 rti(:, 4) = rti(:,4) + ( u1) * g1yti(:) 426 & + ( u2) * g2yti(:) 427 & + ( u3) * g3yti(:) 428! 429 else ! conservative form 430! 431 rti(:, 1) = rti(:,1) + (- u1) * rho * Sclr 432 rti(:, 2) = rti(:,2) + (- u2) * rho * Sclr 433 rti(:, 3) = rti(:,3) + (- u3) * rho * Sclr 434! 435 endif 436 437! flops = flops + 28*npro 438 439 endif 440! 441!.... calculate ( A_i Y,i ) --> rLyi 442! 443 rLyti(:) = rLyti(:) 444 & + A1t(:) * g1yti(:) 445 & + A2t(:) * g2yti(:) 446 & + A3t(:) * g3yti(:) 447 448! 449!.... add contribution to rmi 450! 451! if ((ires .eq. 2) .or. (ires .eq. 3)) 452! & rmi(:,16:20) = rLyi ! modified residual uses non i.b.p form of conv 453! 454!.... ----------------------> LHS <----------------------- 455! 456 if (lhs .eq. 1) then 457! 458!.... loop through the columns 459! 460 do j = 1, nshl 461 462! 463!.... first compute (A_i N_b,i) 464! 465 AitNbi(:) = 466 & WdetJ * shg(:,j,1) * A1t(:) 467 & + WdetJ * shg(:,j,2) * A2t(:) 468 & + WdetJ * shg(:,j,3) * A3t(:) 469 470! 471!.... now loop through the rows and add (N_a A_i N_b,i) to 472! the tangent matrix. 473! 474 do i = 1, nshl 475 476 EGmasst(:,i,j) = EGmasst(:,i,j) + shape(:,i) * AitNbi(:) 477 478 479! 480!.... end loop on rows 481! 482 enddo 483! 484!.... end loop on columns 485! 486 enddo 487! 488!.... end of LHS tangent matrix computation 489! 490 endif 491 492 ttim(22) = ttim(22) + tmr() 493! 494!.... return 495! 496 return 497 end 498 499