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 endif 223! SAM wants a threshold here so we are going to take over this little used 224! error indictor for that purpose. To revert note you will want to uncomment the original 225! form of this error indicator in e3LS.f 226 if((intp.eq.1).and.(ierrcalc.eq.1).and.(nitr.eq.iter)) then 227 do i=1,npro 228 Tmax=maxval(yl(i,:,5)) 229 Tmin=minval(yl(i,:,5)) 230 rerrl(i,:,6)=(Tmax-Tmin)/T(i) 231 enddo 232! the below was somewhat suprisingly ineffective compared to above for identifying shocks. 233! it refined on each side of the shock but left the actual shock quite coarse whereas the above 234! centered well on the shock 235! do j=1,nshl 236! rerrl(:,j,6)=rerrl(:,j,6)+DC(:,intp) !super hack to get error indicator for shock capturing 237! enddo 238 endif 239c 240c 241c.... calculate the time derivative (mass) contribution to RHS 242c 243 if (ngauss .eq. 1 .and. nshl .eq. 4) then ! trilinear tets 244 ttim(17) = ttim(17) - secs(0.0) 245 call e3juel (yl, acl, Sclr, A0, WdetJ, rl, rml) 246 ttim(17) = ttim(17) + secs(0.0) 247 else 248 call e3massr (aci, dui, ri, rmi, A0) 249 endif 250 251c 252c.... calculate the time (mass) contribution to the LHS 253c 254 if (lhs .eq. 1) then 255 call e3massl (bcool,shape, WdetJ, A0, EGmass) 256 endif 257c 258c.... calculate the preconditioner all at once now instead of in separate 259c routines 260c 261 if(iprec.eq.1 .and. lhs.ne.1)then 262 ttim(18) = ttim(18) - secs(0.0) 263 264 if (itau.lt.10) then 265 266 call e3bdg(shape, shg, WdetJ, 267 & A1, A2, A3, 268 & A0, bcool, tau, 269 & u1, u2, u3, 270 & BDiagl, 271 & rmu, rlm2mu, con) 272 273 else 274 275 call e3bdg_nd(shape, shg, WdetJ, 276 & A1, A2, A3, 277 & A0, bcool, PTau, 278 & u1, u2, u3, 279 & BDiagl, 280 & rmu, rlm2mu, con) 281 282 endif 283 284 ttim(18) = ttim(18) + secs(0.0) 285 endif 286c 287c 288c.... multiply flux terms by shape functions and derivatives (from weight space for RHS and 289c by both the weight space and solution space for the LHS stiffness term) 290c 291 ttim(19) = ttim(19) - secs(0.0) 292 call e3wmlt (shape, shg, WdetJ, 293 & ri, rmi, rl, 294 & rml, stiff, EGmass) 295 ttim(19) = ttim(19) + secs(0.0) 296c 297c.... end of integration loop 298c 299 enddo 300 301 ttim(6) = ttim(6) + secs(0.0) 302c 303c.... return 304c 305 return 306 end 307c 308c 309c 310c 311 subroutine e3Sclr (ycl, acl, 312 & dwl, elDwl, 313 & shp, sgn, 314 & shgl, xl, 315 & rtl, rmtl, 316 & qtl, EGmasst) 317c 318c---------------------------------------------------------------------- 319c 320c This routine is the 3D element routine for the N-S equations. 321c This routine calculates the RHS residual and if requested the 322c modified residual or the LHS consistent mass matrix or the block- 323c diagonal preconditioner. 324c 325c input: e a 1..5 when we think of U^e_a and U is 5 variables 326c ycl (npro,nshl,ndof) : Y variables 327c actl (npro,nshl) : scalar variable time derivative 328c dwl (npro,nenl) : distances to wall 329c shp (nen,ngauss) : element shape-functions N_a 330c shgl (nsd,nen,ngauss) : element local-shape-functions N_{a,xi} 331c xl (npro,nenl,nsd) : nodal coordinates at current step (x^e_a) 332c qtl (npro,nshl) : diffusive flux vector (don't worry) 333c 334c output: 335c rtl (npro,nshl) : element RHS residual (G^e_a) 336c rmtl (npro,nshl) : element modified residual (G^e_a tilde) 337c EGmasst (npro,nshape,nshape) : element LHS tangent mass matrix (dG^e_a 338c dY_b ) 339c 340c 341c Note: This routine will calculate the element matrices for the 342c Hulbert's generalized alpha method integrator 343c 344c Note: nedof = nflow*nshl is the total number of degrees of freedom 345c at each element node. 346c 347c Mathematics done by: Michel Mallet, Farzin Shakib (version 1) 348c Farzin Shakib (version 2) 349c 350c 351c Zdenek Johan, Summer 1990. (Modified from e2.f) 352c Zdenek Johan, Winter 1991. (Fortran 90) 353c Kenneth Jansen, Winter 1997. (Primitive Variables) 354c Chris Whiting, Winter 1998. (LHS matrix formation) 355c---------------------------------------------------------------------- 356c 357 include "common.h" 358c 359 dimension ycl(npro,nshl,ndof), 360 & acl(npro,nshl,ndof), 361 & dwl(npro,nenl), 362 & shp(nshl,ngauss), shgl(nsd,nshl,ngauss), 363 & xl(npro,nenl,nsd), 364 & rtl(npro,nshl), qtl(npro,nshl), 365 & rmtl(npro,nshl), Diagl(npro,nshl), 366 & EGmasst(npro,nshape,nshape), 367 & dist2w(npro), sgn(npro,nshl), 368 & vort(npro), gVnrm(npro), 369 & rmu(npro), con(npro), 370 & T(npro), cp(npro), 371 & g1yti(npro), acti(npro), 372 & g2yti(npro), g3yti(npro), 373 & Sclr(npro), srcp (npro) 374 375c 376 dimension shg(npro,nshl,nsd) 377c 378 dimension dxidx(npro,nsd,nsd), WdetJ(npro) 379c 380 dimension rho(npro), rk(npro), 381 & u1(npro), u2(npro), 382 & u3(npro) 383c 384 dimension A0t(npro), A1t(npro), 385 & A2t(npro), A3t(npro) 386c 387 dimension rLyti(npro), raLSt(npro), 388 & rTLSt(npro), giju(npro,6), 389 & DCt(npro, ngauss) 390c 391 dimension rti(npro,nsd+1), rmti(npro,nsd+1), 392 & stifft(npro,nsd,nsd), 393 & shape(npro,nshl), shdrv(npro,nsd,nshl) 394 real*8 elDwl(npro) 395c 396 ttim(6) = ttim(6) - tmr() 397c 398c.... loop through the integration points 399c 400 elDwl(:)=zero 401 do intp = 1, ngauss 402c 403c.... if Det. .eq. 0, do not include this point 404c 405 if (Qwt(lcsyst,intp) .eq. zero) cycle ! precaution 406c 407c 408c.... create a matrix of shape functions (and derivatives) for each 409c element at this quadrature point. These arrays will contain 410c the correct signs for the hierarchic basis 411c 412 call getshp(shp, shgl, sgn, 413 & shape, shdrv) 414c 415c.... initialize 416c 417 rlyti = zero 418 rti = zero 419 rmti = zero 420 srcp = zero 421 stifft = zero 422c if (lhs .eq. 1) stifft = zero 423c 424c 425c.... calculate the integration variables 426c 427 ttim(8) = ttim(8) - tmr() 428c 429 call e3ivarSclr(ycl, acl, acti, 430 & shape, shdrv, xl, 431 & T, cp, 432 & dxidx, Sclr, 433 & Wdetj, vort, gVnrm, 434 & g1yti, g2yti, g3yti, 435 & rho, rmu, con, 436 & rk, u1, u2, 437 & u3, shg, dwl, 438 & dist2w) 439 ttim(8) = ttim(8) + tmr() 440 441c 442c.... calculate the source term contribution 443c 444 if(nosource.ne.1) 445 & call e3sourceSclr (Sclr, rho, rmu, 446 & dist2w, vort, gVnrm, con, 447 & g1yti, g2yti, g3yti, 448 & rti, rLyti, srcp, 449 & ycl, shape, u1, 450 & u2, u3, xl, 451 & elDwl) 452c 453 if (ilset.eq.2 .and. isclr.eq.2) then 454 rk = pt5 * ( u1**2 + u2**2 + u3**2 ) 455 endif 456c.... calculate the relevant matrices 457c 458 ttim(9) = ttim(9) - tmr() 459 call e3mtrxSclr (rho, 460 & u1, u2, u3, 461 & A0t, A1t, 462 & A2t, A3t) 463 ttim(9) = ttim(9) + tmr() 464c 465c.... calculate the convective contribution (Galerkin) 466c 467 ttim(14) = ttim(14) - tmr() 468 call e3convSclr (g1yti, g2yti, g3yti, 469 & A1t, A2t, A3t, 470 & rho, u1, Sclr, 471 & u2, u3, rLyti, 472 & rti, rmti, EGmasst, 473 & shg, shape, WdetJ) 474 ttim(14) = ttim(14) + tmr() 475c 476c.... calculate the diffusion contribution 477c 478 ttim(15) = ttim(15) - tmr() 479 if (Navier .eq. 1) then 480 call e3viscSclr (g1yti, g2yti, g3yti, 481 & rmu, Sclr, rho, 482 & rti, rmti, stifft ) 483 endif 484 ttim(15) = ttim(15) + tmr() 485c 486 if (ilset.eq.2) srcp = zero 487 488c 489c.... calculate the least-squares contribution 490c 491 ttim(16) = ttim(16) - tmr() 492 call e3LSSclr(A1t, A2t, A3t, 493 & rho, rmu, rtLSt, 494 & u1, u2, u3, 495 & rLyti, dxidx, raLSt, 496 & rti, rk, giju, 497 & acti, A0t, 498 & shape, shg, 499 & EGmasst, stifft, WdetJ, 500 & srcp) 501 ttim(16) = ttim(16) + tmr() 502c 503c******************************DC TERMS***************************** 504 if (idcsclr(1) .ne. 0) then 505 if ((idcsclr(2).eq.1 .and. isclr.eq.1) .or. 506 & (idcsclr(2).eq.2 .and. isclr.eq.2)) then ! scalar with dc 507 call e3dcSclr(g1yti, g2yti, g3yti, 508 & A0t, raLSt, rTLSt, 509 & DCt, giju, 510 & rti, rmti, stifft) 511 endif 512 endif 513c 514c****************************************************************** 515c.... calculate the time derivative (mass) contribution to RHS 516c 517 518 call e3massrSclr (acti, rti, A0t) 519c 520c.... calculate the time (mass) contribution to the LHS 521c 522 if (lhs .eq. 1) then 523 call e3masslSclr (shape, WdetJ, A0t, EGmasst,srcp) 524 endif 525c 526 527c.... multiply flux terms by shape functions and derivatives (from weight space for RHS and 528c by both the weight space and solution space for the LHS stiffness term) 529c 530 ttim(19) = ttim(19) - tmr() 531 call e3wmltSclr (shape, shg, WdetJ, 532 & rti, 533 & rtl, stifft, EGmasst) 534 ttim(19) = ttim(19) + tmr() 535c 536c.... end of the loop 537c 538 enddo 539 qptinv=one/ngauss 540 elDwl(:)=elDwl(:)*qptinv 541 542 543 ttim(6) = ttim(6) + tmr() 544c 545c.... return 546c 547 return 548 end 549 550