1 subroutine e3 (yl, ycl, acl, shp, 2 & shgl, xl, rl, rml, xmudmi, 3 & BDiagl, ql, sgn, rlsl, EGmass, 4 & rerrl, ytargetl) 5c 6c---------------------------------------------------------------------- 7c 8c This routine is the 3D element routine for the N-S equations. 9c This routine calculates the RHS residual and if requested the 10c modified residual or the LHS consistent mass matrix or the block- 11c diagonal preconditioner. 12c 13c input: 14c yl (npro,nshl,nflow) : Y variables (DOES NOT CONTAIN SCALARS) 15c ycl (npro,nshl,ndof) : Y variables at current step 16c acl (npro,nshl,ndof) : Y acceleration (Y_{,t}) 17c shp (nshl,ngauss) : element shape-functions N_a 18c shgl (nsd,nshl,ngauss) : element local-shape-functions N_{a,xi} 19c xl (npro,nenl,nsd) : nodal coordinates at current step (x^e_a) 20c ql (npro,nshl,(nflow-1)*nsd) : diffusive flux vector 21c rlsl (npro,nshl,6) : resolved Leonard stresses 22c sgn (npro,nshl) : shape function sign matrix 23c 24c output: 25c rl (npro,nshl,nflow) : element RHS residual (G^e_a) 26c rml (npro,nshl,nflow) : element modified residual (G^e_a tilde) 27c EGmass (npro,nedof,nedof) : element LHS tangent mass matrix (dG^e_a 28c dY_b ) 29c BDiagl (npro,nshl,nflow,nflow) : block-diagonal preconditioner 30c 31c 32c Note: This routine will calculate the element matrices for the 33c Hulbert's generalized alpha method integrator 34c 35c Note: nedof = nflow*nshape is the total number of degrees of freedom 36c at each element node. 37c 38c Mathematics done by: Michel Mallet, Farzin Shakib (version 1) 39c Farzin Shakib (version 2) 40c 41c 42c Zdenek Johan, Summer 1990. (Modified from e2.f) 43c Zdenek Johan, Winter 1991. (Fortran 90) 44c Kenneth Jansen, Winter 1997. (Primitive Variables) 45c Chris Whiting, Winter 1998. (LHS matrix formation) 46c---------------------------------------------------------------------- 47c 48 include "common.h" 49c 50 dimension yl(npro,nshl,nflow), ycl(npro,nshl,ndof), 51 & acl(npro,nshl,ndof), rmu(npro), 52 & shp(nshl,ngauss), rlm2mu(npro), 53 & shgl(nsd,nshl,ngauss), con(npro), 54 & xl(npro,nenl,nsd), rlm(npro), 55 & rl(npro,nshl,nflow), ql(npro,nshl,idflx), 56 & rml(npro,nshl,nflow), xmudmi(npro,ngauss), 57 & BDiagl(npro,nshl,nflow,nflow), 58 & EGmass(npro,nedof,nedof),cv(npro), 59 & ytargetl(npro,nshl,nflow) 60c 61 dimension dui(npro,ndof), aci(npro,ndof) 62c 63 dimension g1yi(npro,nflow), g2yi(npro,nflow), 64 & g3yi(npro,nflow), shg(npro,nshl,nsd), 65 & divqi(npro,nflow), tau(npro,5) 66c 67 dimension dxidx(npro,nsd,nsd), WdetJ(npro) 68c 69 dimension rho(npro), pres(npro), 70 & T(npro), ei(npro), 71 & h(npro), alfap(npro), 72 & betaT(npro), DC(npro,ngauss), 73 & cp(npro), rk(npro), 74 & u1(npro), u2(npro), 75 & u3(npro), A0DC(npro,4), 76 & Sclr(npro), dVdY(npro,15), 77 & giju(npro,6), rTLS(npro), 78 & raLS(npro), A0inv(npro,15) 79c 80 dimension A0(npro,nflow,nflow), A1(npro,nflow,nflow), 81 & A2(npro,nflow,nflow), A3(npro,nflow,nflow) 82c 83 dimension rLyi(npro,nflow), sgn(npro,nshl) 84c 85 dimension ri(npro,nflow*(nsd+1)), rmi(npro,nflow*(nsd+1)), 86 & shape(npro,nshl), shdrv(npro,nsd,nshl), 87 & stiff(npro,nsd*nflow,nsd*nflow), 88 & PTau(npro,5,5), 89 & sforce(npro,3), compK(npro,10) 90c 91 dimension x(npro,3), bcool(npro) 92 93 dimension rlsl(npro,nshl,6), rlsli(npro,6) 94 95 real*8 rerrl(npro,nshl,6) 96 ttim(6) = ttim(6) - secs(0.0) 97c 98c.... local reconstruction of diffusive flux vector 99c (note: not currently included in mfg) 100 if (idiff==2 .and. (ires==3 .or. ires==1)) then 101 call e3ql (ycl, shp, shgl, xl, ql, xmudmi, sgn) 102 endif 103c 104c.... loop through the integration points 105c 106 do intp = 1, ngauss 107c 108c.... if Det. .eq. 0, do not include this point 109c 110 if (Qwt(lcsyst,intp) .eq. zero) cycle ! precaution 111c 112c.... create a matrix of shape functions (and derivatives) for each 113c element at this quadrature point. These arrays will contain 114c the correct signs for the hierarchic basis 115c 116 call getshp(shp, shgl, sgn, 117 & shape, shdrv) 118c 119c.... initialize 120c 121 ri = zero 122 rmi = zero 123 if (lhs .eq. 1) stiff = zero 124c 125c 126c.... calculate the integration variables 127c 128 ttim(8) = ttim(8) - secs(0.0) 129 call e3ivar (yl, ycl, acl, 130 & Sclr, shape, shdrv, 131 & xl, dui, aci, 132 & g1yi, g2yi, g3yi, 133 & shg, dxidx, WdetJ, 134 & rho, pres, T, 135 & ei, h, alfap, 136 & betaT, cp, rk, 137 & u1, u2, u3, 138 & ql, divqi, sgn, 139 & rLyi, !passed as a work array 140 & rmu, rlm, rlm2mu, 141 & con, rlsl, rlsli, 142 & xmudmi, sforce, cv) 143 ttim(8) = ttim(8) + secs(0.0) 144 145c 146c.... calculate the relevant matrices 147c 148 ttim(9) = ttim(9) - secs(0.0) 149 call e3mtrx (rho, pres, T, 150 & ei, h, alfap, 151 & betaT, cp, rk, 152 & u1, u2, u3, 153 & A0, A1, 154 & A2, A3, 155 & rLyi(:,1), rLyi(:,2), rLyi(:,3), ! work arrays 156 & rLyi(:,4), rLyi(:,5), A0DC, 157 & A0inv, dVdY) 158 ttim(9) = ttim(9) + tmr() 159c 160c.... calculate the convective contribution (Galerkin) 161c 162 ttim(14) = ttim(14) - secs(0.0) 163 call e3conv (g1yi, g2yi, g3yi, 164 & A1, A2, A3, 165 & rho, pres, T, 166 & ei, rk, u1, 167 & u2, u3, rLyi, 168 & ri, rmi, EGmass, 169 & shg, shape, WdetJ) 170 ttim(14) = ttim(14) + secs(0.0) 171c 172c.... calculate the diffusion contribution 173c 174 ttim(15) = ttim(15) - secs(0.0) 175 compK = zero 176 if (Navier .eq. 1) then 177 call e3visc (g1yi, g2yi, g3yi, 178 & dxidx, 179 & rmu, rlm, rlm2mu, 180 & u1, u2, u3, 181 & ri, rmi, stiff, 182 & con, rlsli, compK, T) 183 endif 184 ttim(15) = ttim(15) + secs(0.0) 185c 186c.... calculate the body force contribution 187c 188 if(isurf .ne. 1 .and. matflg(5,1).gt.0) then 189 call e3source (ri, rmi, rlyi, 190 & rho, u1, u2, 191 & u3, pres, sforce, 192 & dui, dxidx, ytargetl, 193 & xl, shape, bcool) 194 else 195 bcool=zero 196 endif 197c 198c.... calculate the least-squares contribution 199c 200 ttim(16) = ttim(16) - secs(0.0) 201 call e3LS (A1, A2, A3, 202 & rho, rmu, cp, 203 & cv, con, T, 204 & u1, u2, u3, 205 & rLyi, dxidx, tau, 206 & ri, rmi, rk, 207 & dui, aci, A0, 208 & divqi, shape, shg, 209 & EGmass, stiff, WdetJ, 210 & giju, rTLS, raLS, 211 & A0inv, dVdY, rerrl, 212 & compK, pres, PTau) 213 ttim(16) = ttim(16) + secs(0.0) 214c 215c.... Discontinuity capturing 216c 217 if(iDC.ne.0) then 218 call e3dc (g1yi, g2yi, g3yi, 219 & A0, raLS, rTLS, 220 & giju, DC, 221 & ri, rmi, stiff, A0DC) 222 if((intp.eq.1).and.(ierrcalc.eq.1).and.(nitr.eq.iter)) then 223 do i=1,npro 224 Tmax=maxval(yl(i,:,5)) 225 Tmin=minval(yl(i,:,5)) 226 rerrl(i,:,6)=(Tmax-Tmin)/T(i) 227 enddo 228! do j=1,nshl 229! rerrl(:,j,6)=rerrl(:,j,6)+DC(:,intp) !super hack to get error indicator for shock capturing 230! enddo 231 endif 232 endif 233c 234c 235c.... calculate the time derivative (mass) contribution to RHS 236c 237 if (ngauss .eq. 1 .and. nshl .eq. 4) then ! trilinear tets 238 ttim(17) = ttim(17) - secs(0.0) 239 call e3juel (yl, acl, Sclr, A0, WdetJ, rl, rml) 240 ttim(17) = ttim(17) + secs(0.0) 241 else 242 call e3massr (aci, dui, ri, rmi, A0) 243 endif 244 245c 246c.... calculate the time (mass) contribution to the LHS 247c 248 if (lhs .eq. 1) then 249 call e3massl (bcool,shape, WdetJ, A0, EGmass) 250 endif 251c 252c.... calculate the preconditioner all at once now instead of in separate 253c routines 254c 255 if(iprec.eq.1 .and. lhs.ne.1)then 256 ttim(18) = ttim(18) - secs(0.0) 257 258 if (itau.lt.10) then 259 260 call e3bdg(shape, shg, WdetJ, 261 & A1, A2, A3, 262 & A0, bcool, tau, 263 & u1, u2, u3, 264 & BDiagl, 265 & rmu, rlm2mu, con) 266 267 else 268 269 call e3bdg_nd(shape, shg, WdetJ, 270 & A1, A2, A3, 271 & A0, bcool, PTau, 272 & u1, u2, u3, 273 & BDiagl, 274 & rmu, rlm2mu, con) 275 276 endif 277 278 ttim(18) = ttim(18) + secs(0.0) 279 endif 280c 281c 282c.... multiply flux terms by shape functions and derivatives (from weight space for RHS and 283c by both the weight space and solution space for the LHS stiffness term) 284c 285 ttim(19) = ttim(19) - secs(0.0) 286 call e3wmlt (shape, shg, WdetJ, 287 & ri, rmi, rl, 288 & rml, stiff, EGmass) 289 ttim(19) = ttim(19) + secs(0.0) 290c 291c.... end of integration loop 292c 293 enddo 294 295 ttim(6) = ttim(6) + secs(0.0) 296c 297c.... return 298c 299 return 300 end 301c 302c 303c 304c 305 subroutine e3Sclr (ycl, acl, 306 & dwl, elDwl, 307 & shp, sgn, 308 & shgl, xl, 309 & rtl, rmtl, 310 & qtl, EGmasst) 311c 312c---------------------------------------------------------------------- 313c 314c This routine is the 3D element routine for the N-S equations. 315c This routine calculates the RHS residual and if requested the 316c modified residual or the LHS consistent mass matrix or the block- 317c diagonal preconditioner. 318c 319c input: e a 1..5 when we think of U^e_a and U is 5 variables 320c ycl (npro,nshl,ndof) : Y variables 321c actl (npro,nshl) : scalar variable time derivative 322c dwl (npro,nenl) : distances to wall 323c shp (nen,ngauss) : element shape-functions N_a 324c shgl (nsd,nen,ngauss) : element local-shape-functions N_{a,xi} 325c xl (npro,nenl,nsd) : nodal coordinates at current step (x^e_a) 326c qtl (npro,nshl) : diffusive flux vector (don't worry) 327c 328c output: 329c rtl (npro,nshl) : element RHS residual (G^e_a) 330c rmtl (npro,nshl) : element modified residual (G^e_a tilde) 331c EGmasst (npro,nshape,nshape) : element LHS tangent mass matrix (dG^e_a 332c dY_b ) 333c 334c 335c Note: This routine will calculate the element matrices for the 336c Hulbert's generalized alpha method integrator 337c 338c Note: nedof = nflow*nshl is the total number of degrees of freedom 339c at each element node. 340c 341c Mathematics done by: Michel Mallet, Farzin Shakib (version 1) 342c Farzin Shakib (version 2) 343c 344c 345c Zdenek Johan, Summer 1990. (Modified from e2.f) 346c Zdenek Johan, Winter 1991. (Fortran 90) 347c Kenneth Jansen, Winter 1997. (Primitive Variables) 348c Chris Whiting, Winter 1998. (LHS matrix formation) 349c---------------------------------------------------------------------- 350c 351 include "common.h" 352c 353 dimension ycl(npro,nshl,ndof), 354 & acl(npro,nshl,ndof), 355 & dwl(npro,nenl), 356 & shp(nshl,ngauss), shgl(nsd,nshl,ngauss), 357 & xl(npro,nenl,nsd), 358 & rtl(npro,nshl), qtl(npro,nshl), 359 & rmtl(npro,nshl), Diagl(npro,nshl), 360 & EGmasst(npro,nshape,nshape), 361 & dist2w(npro), sgn(npro,nshl), 362 & vort(npro), gVnrm(npro), 363 & rmu(npro), con(npro), 364 & T(npro), cp(npro), 365 & g1yti(npro), acti(npro), 366 & g2yti(npro), g3yti(npro), 367 & Sclr(npro), srcp (npro) 368 369c 370 dimension shg(npro,nshl,nsd) 371c 372 dimension dxidx(npro,nsd,nsd), WdetJ(npro) 373c 374 dimension rho(npro), rk(npro), 375 & u1(npro), u2(npro), 376 & u3(npro) 377c 378 dimension A0t(npro), A1t(npro), 379 & A2t(npro), A3t(npro) 380c 381 dimension rLyti(npro), raLSt(npro), 382 & rTLSt(npro), giju(npro,6), 383 & DCt(npro, ngauss) 384c 385 dimension rti(npro,nsd+1), rmti(npro,nsd+1), 386 & stifft(npro,nsd,nsd), 387 & shape(npro,nshl), shdrv(npro,nsd,nshl) 388 real*8 elDwl(npro) 389c 390 ttim(6) = ttim(6) - tmr() 391c 392c.... loop through the integration points 393c 394 elDwl(:)=zero 395 do intp = 1, ngauss 396c 397c.... if Det. .eq. 0, do not include this point 398c 399 if (Qwt(lcsyst,intp) .eq. zero) cycle ! precaution 400c 401c 402c.... create a matrix of shape functions (and derivatives) for each 403c element at this quadrature point. These arrays will contain 404c the correct signs for the hierarchic basis 405c 406 call getshp(shp, shgl, sgn, 407 & shape, shdrv) 408c 409c.... initialize 410c 411 rlyti = zero 412 rti = zero 413 rmti = zero 414 srcp = zero 415 stifft = zero 416c if (lhs .eq. 1) stifft = zero 417c 418c 419c.... calculate the integration variables 420c 421 ttim(8) = ttim(8) - tmr() 422c 423 call e3ivarSclr(ycl, acl, acti, 424 & shape, shdrv, xl, 425 & T, cp, 426 & dxidx, Sclr, 427 & Wdetj, vort, gVnrm, 428 & g1yti, g2yti, g3yti, 429 & rho, rmu, con, 430 & rk, u1, u2, 431 & u3, shg, dwl, 432 & dist2w) 433 ttim(8) = ttim(8) + tmr() 434 435c 436c.... calculate the source term contribution 437c 438 if(nosource.ne.1) 439 & call e3sourceSclr (Sclr, rho, rmu, 440 & dist2w, vort, gVnrm, con, 441 & g1yti, g2yti, g3yti, 442 & rti, rLyti, srcp, 443 & ycl, shape, u1, 444 & u2, u3, xl, 445 & elDwl) 446c 447 if (ilset.eq.2 .and. isclr.eq.2) then 448 rk = pt5 * ( u1**2 + u2**2 + u3**2 ) 449 endif 450c.... calculate the relevant matrices 451c 452 ttim(9) = ttim(9) - tmr() 453 call e3mtrxSclr (rho, 454 & u1, u2, u3, 455 & A0t, A1t, 456 & A2t, A3t) 457 ttim(9) = ttim(9) + tmr() 458c 459c.... calculate the convective contribution (Galerkin) 460c 461 ttim(14) = ttim(14) - tmr() 462 call e3convSclr (g1yti, g2yti, g3yti, 463 & A1t, A2t, A3t, 464 & rho, u1, Sclr, 465 & u2, u3, rLyti, 466 & rti, rmti, EGmasst, 467 & shg, shape, WdetJ) 468 ttim(14) = ttim(14) + tmr() 469c 470c.... calculate the diffusion contribution 471c 472 ttim(15) = ttim(15) - tmr() 473 if (Navier .eq. 1) then 474 call e3viscSclr (g1yti, g2yti, g3yti, 475 & rmu, Sclr, rho, 476 & rti, rmti, stifft ) 477 endif 478 ttim(15) = ttim(15) + tmr() 479c 480 if (ilset.eq.2) srcp = zero 481 482c 483c.... calculate the least-squares contribution 484c 485 ttim(16) = ttim(16) - tmr() 486 call e3LSSclr(A1t, A2t, A3t, 487 & rho, rmu, rtLSt, 488 & u1, u2, u3, 489 & rLyti, dxidx, raLSt, 490 & rti, rk, giju, 491 & acti, A0t, 492 & shape, shg, 493 & EGmasst, stifft, WdetJ, 494 & srcp) 495 ttim(16) = ttim(16) + tmr() 496c 497c******************************DC TERMS***************************** 498 if (idcsclr(1) .ne. 0) then 499 if ((idcsclr(2).eq.1 .and. isclr.eq.1) .or. 500 & (idcsclr(2).eq.2 .and. isclr.eq.2)) then ! scalar with dc 501 call e3dcSclr(g1yti, g2yti, g3yti, 502 & A0t, raLSt, rTLSt, 503 & DCt, giju, 504 & rti, rmti, stifft) 505 endif 506 endif 507c 508c****************************************************************** 509c.... calculate the time derivative (mass) contribution to RHS 510c 511 512 call e3massrSclr (acti, rti, A0t) 513c 514c.... calculate the time (mass) contribution to the LHS 515c 516 if (lhs .eq. 1) then 517 call e3masslSclr (shape, WdetJ, A0t, EGmasst,srcp) 518 endif 519c 520 521c.... multiply flux terms by shape functions and derivatives (from weight space for RHS and 522c by both the weight space and solution space for the LHS stiffness term) 523c 524 ttim(19) = ttim(19) - tmr() 525 call e3wmltSclr (shape, shg, WdetJ, 526 & rti, 527 & rtl, stifft, EGmasst) 528 ttim(19) = ttim(19) + tmr() 529c 530c.... end of the loop 531c 532 enddo 533 qptinv=one/ngauss 534 elDwl(:)=elDwl(:)*qptinv 535 536 537 ttim(6) = ttim(6) + tmr() 538c 539c.... return 540c 541 return 542 end 543 544