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