1 subroutine e3 (yl, acl, dwl, shp, 2 & shgl, xl, rl, ql, 3 & xKebe, xGoC, xmudmi, sgn, 4 & rerrl, rlsl) 5c 6c---------------------------------------------------------------------- 7c 8c This routine calculates the residual and tangent matrix for the 9c UBar formulation of the incompressible Navier Stokes equations. 10c 11c 12c input: e a 1..5 when we think of U^e_a and U is 5 variables 13c yl (npro,nshl,ndof) : Y variables (not U) 14c acl (npro,nshl,ndof) : Y acceleration (Y_{,t}) 15c shp (nen,ngauss) : element shape-functions N_a 16c shgl (nsd,nen,ngauss) : element local-shape-functions N_{a,xi} 17c wght (ngauss) : element weight (for quadrature) 18c xl (npro,nenl,nsd) : nodal coordinates at current step (x^e_a) 19c ql (npro,nshl,nsd*nsd) : diffusive flux vector (don't worry) 20c rlsl (npro,nshl,6) : resolved Leonard stresses 21c 22c output: 23c rl (npro,nshl,nflow) : element RHS residual (G^e_a) 24c rml (npro,nshl,nflow) : element modified residual (G^e_a tilde) 25c xKebe (npro,9,nshl,nshl) : element LHS tangent mass matrix 26c xGoC (npro,4,nshl,nshl) : element LHS tangent mass matrix 27c 28c Note: This routine will calculate the element matrices for the 29c Hulbert's generalized alpha method integrator 30c 31c Mathematics done by: Michel Mallet, Farzin Shakib (version 1) 32c Farzin Shakib (version 2) 33c 34c K. E. Jansen, Winter 1999. (advective form formulation) 35c C. H. Whiting, Winter 1999. (advective form formulation) 36c---------------------------------------------------------------------- 37c 38 include "common.h" 39c 40 dimension yl(npro,nshl,ndof), 41 & acl(npro,nshl,ndof), 42 & shp(nshl,ngauss), shgl(nsd,nshl,ngauss), 43 & xl(npro,nenl,nsd), dwl(npro,nenl), 44 & rl(npro,nshl,nflow), ql(npro,nshl,idflx) 45c 46 47 dimension xKebe(npro,9,nshl,nshl), xGoC(npro,4,nshl,nshl) 48c 49c.... local declarations 50c 51 dimension g1yi(npro,ndof), g2yi(npro,ndof), 52 & g3yi(npro,ndof), shg(npro,nshl,nsd), 53 & aci(npro,3), dxidx(npro,nsd,nsd), 54 & WdetJ(npro), rho(npro), 55 & pres(npro), u1(npro), 56 & u2(npro), u3(npro), 57 & rLui(npro,nsd), uBar(npro,nsd), 58 & xmudmi(npro,ngauss), sgn(npro,nshl), 59 & shpfun(npro,nshl), shdrv(npro,nsd,nshl), 60 & rmu(npro), tauC(npro), 61 & tauM(npro), tauBar(npro), 62 & src(npro,3) 63 64 dimension rlsl(npro,nshl,6), rlsli(npro,6) 65 66 real*8 rerrl(npro,nshl,6) 67 integer aa 68 69c 70c 71c.... local reconstruction of diffusive flux vector for quadratics 72c or greater but NOT for bflux since local mass was not mapped 73c 74 if ( idiff==2 .and. ires .eq. 1 ) then 75 call e3ql (yl, dwl, shp, shgl, 76 & xl, ql, xmudmi, 77 & sgn) 78 endif 79c 80c.... loop through the integration points 81c 82 83 do intp = 1, ngauss 84 85 if (Qwt(lcsyst,intp) .eq. zero) cycle ! precaution 86c 87c.... get the hierarchic shape functions at this int point 88c 89 call getshp(shp, shgl, sgn, 90 & shpfun, shdrv) 91c 92c.... get necessary fluid properties (including eddy viscosity) 93c 94 call getdiff(dwl, yl, shpfun, xmudmi, xl, rmu, rho) 95c 96c.... calculate the integration variables 97c 98 call e3ivar (yl, acl, shpfun, 99 & shdrv, xl, 100 & aci, g1yi, g2yi, 101 & g3yi, shg, dxidx, 102 & WdetJ, rho, pres, 103 & u1, u2, u3, 104 & ql, rLui, src, 105 & rerrl, rlsl, rlsli, 106 & dwl) 107c 108c.... compute the stabilization terms 109c 110 call e3stab (rho, u1, u2, 111 & u3, dxidx, rLui, 112 & rmu, tauC, tauM, 113 & tauBar, uBar ) 114c 115c.... compute the residual contribution at this integration point 116c 117 call e3Res ( u1, u2, u3, 118 & uBar, aci, WdetJ, 119 & g1yi, g2yi, g3yi, 120 & rLui, rmu, rho, 121 & tauC, tauM, tauBar, 122 & shpfun, shg, src, 123 & rl, pres, acl, 124 & rlsli) 125c 126c.... compute the tangent matrix contribution 127c 128 if (lhs .eq. 1) then 129 call e3LHS ( u1, u2, u3, 130 & uBar, WdetJ, rho, 131 & rLui, rmu, 132 & tauC, tauM, tauBar, 133 & shpfun, shg, xKebe, 134 & xGoC ) 135 endif 136 137c 138c.... end of integration loop 139c 140 enddo 141 142c 143c.... symmetrize C 144c 145 if (lhs .eq. 1) then 146 do ib = 1, nshl 147 do iaa = 1, ib-1 148 xGoC(:,4,iaa,ib) = xGoC(:,4,ib,iaa) 149 enddo 150 enddo 151 endif 152c 153c.... return 154c 155 return 156 end 157 158 159c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 160c################################################################### 161 162 163 subroutine e3Sclr (yl, acl, shp, 164 & shgl, xl, dwl, 165 & rl, ql, xSebe, 166 & sgn, xmudmi) 167c 168c---------------------------------------------------------------------- 169c 170c This routine calculates the residual and tangent matrix for the 171c advection - diffusion equation for scalar. 172c 173c K. E. Jansen, Winter 1999. (advective form formulation) 174c C. H. Whiting, Winter 1999. (advective form formulation) 175c---------------------------------------------------------------------- 176c 177 include "common.h" 178c 179 real*8 yl(npro,nshl,ndof), acl(npro,nshl,ndof), 180 & shp(nshl,ngauss), shgl(nsd,nshl,ngauss), 181 & xl(npro,nenl,nsd), rl(npro,nshl), 182 & ql(npro,nshl,nsd), xSebe(npro,nshl,nshl), 183 & dwl(npro,nshl) 184c 185c.... local declarations 186c 187 real*8 gradS(npro,nsd), shg(npro,nshl,nsd), 188 & Sdot(npro), Sclr(npro), 189 & dxidx(npro,nsd,nsd), WdetJ(npro), 190 & u1(npro), u2(npro), u3(npro), 191 & sgn(npro,nshl), shpfun(npro,nshl), 192 & shdrv(npro,nsd,nshl), rLS(npro), 193 & tauS(npro), diffus(npro), 194 & srcL(npro), srcR(npro), 195 & gGradS(npro,nsd), dcFct(npro), 196 & giju(npro,6) 197c 198c.... Source terms sometimes take the form (beta_i)*(phi,_i). Since 199c the convective term has (u_i)*(phi,_i), it is useful to treat 200c beta_i as a "correction" to the velocity. In calculating the 201c stabilization terms, the new "modified" velocity (u_i-beta_i) is 202c then used in place of the pure velocity for stabilization terms, 203c and the source term sneaks into the RHS and LHS. 204 real*8 uMod(npro,nsd), srcRat(npro), xmudmi(npro,ngauss) 205c 206 integer aa, b 207c 208c.... local reconstruction of diffusive flux vector 209c 210 if ( idiff==2 ) then 211 call e3qlSclr (yl, dwl, shp, shgl, xl, ql, sgn) 212 endif 213c 214c.... loop through the integration points 215c 216 do intp = 1, ngauss 217 218 if (Qwt(lcsyst,intp) .eq. zero) cycle ! precaution 219c 220c.... get the hierarchic shape functions at this int point 221c 222 call getshp(shp, shgl, sgn, 223 & shpfun, shdrv) 224c 225c.... get necessary fluid properties 226c 227 call getdiffsclr(shpfun,dwl,yl,diffus) 228c 229c.... calculate the integration variables 230c 231 call e3ivarSclr(yl, acl, shpfun, 232 & shdrv, xl, xmudmi, 233 & Sclr, Sdot, gradS, 234 & shg, dxidx, WdetJ, 235 & u1, u2, u3, 236 & ql, rLS, SrcR, 237 & SrcL, uMod, dwl, 238 & diffus, srcRat) 239 240 241c 242c.... compute the stabilization terms 243c 244 call e3StabSclr (uMod, dxidx, tauS, 245 & diffus, srcR, giju, 246 & srcRat) 247c 248c... computing the DC factor for the discontinuity capturing 249c 250 if (idcsclr(1) .ne. 0) then 251 if ((idcsclr(2).eq.1 .and. isclr.eq.1) .or. 252 & (idcsclr(2).eq.2 .and. isclr.eq.2)) then ! scalar with dc 253c 254 call e3dcSclr ( gradS, giju, gGradS, 255 & rLS, tauS, srcR, 256 & dcFct) 257 endif 258 endif !end of idcsclr 259c 260c.... compute the residual contribution at this integration point 261c 262 call e3ResSclr ( uMod, gGradS, 263 & Sclr, Sdot, gradS, 264 & WdetJ, rLS, tauS, 265 & shpfun, shg, srcR, 266 & diffus, 267 & rl ) 268c 269c.... compute the tangent matrix contribution 270c 271 if (lhs .eq. 1) then 272 call e3LHSSclr ( uMod, giju, dcFct, 273 & Sclr, Sdot, gradS, 274 & WdetJ, rLS, tauS, 275 & shpfun, shg, srcL, 276 & diffus, 277 & xSebe ) 278 279 endif 280 281c 282c.... end of integration loop 283c 284 enddo 285 286c 287c.... return 288c 289 return 290 end 291