1 subroutine e3ql (yl, dwl, shp, shgl, 2 & xl, ql, xmudmi, 3 & sgn ) 4c 5c---------------------------------------------------------------------- 6c 7c This routine computes the local diffusive flux vector using a 8c local projection algorithm 9c 10c input: 11c yl (npro,nshl,ndof) : Y variables 12c shp (nen,ngauss) : element shape-functions 13c shgl (nsd,nen,ngauss) : element local-grad-shape-functions 14c xl (npro,nshape,nsd) : nodal coordinates at current step 15c sgn (npro,nshl) : signs for reversed shape functions 16c 17c output: 18c ql (npro,nshl,nsd*nsd) : element RHS diffusion residual 19c 20c---------------------------------------------------------------------- 21c 22 use local_mass 23 include "common.h" 24c 25 dimension yl(npro,nshl,ndof), dwl(npro,nshl), 26 & shp(nshl,ngauss), shgl(nsd,nshl,ngauss), 27 & xl(npro,nenl,nsd), sgn(npro,nshl), 28 & ql(npro,nshl,idflx), xmudmi(npro,ngauss) 29c 30c local arrays 31c 32 dimension g1yi(npro,ndof), g2yi(npro,ndof), 33 & g3yi(npro,ndof), shg(npro,nshl,nsd), 34 & dxidx(npro,nsd,nsd), WdetJ(npro), 35 & rmu(npro), 36 & rminv(npro,nshl,nshl), 37 & qrl(npro,nshl,nsd*nsd) 38c 39 dimension qdi(npro,nsd*nsd), shape(npro,nshl), 40 & shdrv(npro,nsd,nshl), indx(nshl), 41 & rmass(npro,nshl,nshl) 42 43 44 real*8 tmp(npro) 45c 46c.... loop through the integration points 47c 48 rminv = zero 49 rmass = zero 50 qrl = zero 51 52 do intp = 1, ngauss 53 54 call getshp(shp, shgl, sgn, shape, shdrv) 55 56 qdi = zero 57c 58c.... calculate the integration variables 59c 60c 61 call e3qvar (yl, shdrv, 62 & xl, g1yi, 63 & g2yi, g3yi, shg, 64 & dxidx, WdetJ ) 65 66 call getdiff(dwl, yl, shape, xmudmi, xl,rmu, tmp) 67c 68c.... diffusive flux in x1-direction 69c 70 qdi(:,1) = two * rmu * g1yi(:,2) 71 qdi(:,4) = rmu * (g1yi(:,3) + g2yi(:,2)) 72 qdi(:,7) = rmu * (g1yi(:,4) + g3yi(:,2)) 73c 74c.... diffusive flux in x2-direction 75c 76 qdi(:,2) = rmu * (g1yi(:,3) + g2yi(:,2)) 77 qdi(:,5) = two * rmu * g2yi(:,3) 78 qdi(:,8) = rmu * (g2yi(:,4) + g3yi(:,3)) 79c 80c.... diffusive flux in x3-direction 81c 82 qdi(:,3) = rmu * (g1yi(:,4) + g3yi(:,2)) 83 qdi(:,6)= rmu * (g2yi(:,4) + g3yi(:,3)) 84 qdi(:,9)= two * rmu * g3yi(:,4) 85c 86c 87c.... assemble contribution of qdi to qrl,i.e., contribution to 88c each element shape function 89c 90 tmp = Qwt(lcsyst,intp) 91 if (lcsyst .eq. 1) then 92 tmp = tmp*(three/four) 93 endif 94c 95c reconsider this when hierarchic wedges come into code WDGCHECK 96c 97 98 do i=1,nshl 99 qrl(:,i,1 ) = qrl(:,i,1 )+ shape(:,i)*tmp*qdi(:,1 ) 100 qrl(:,i,2 ) = qrl(:,i,2 )+ shape(:,i)*tmp*qdi(:,2 ) 101 qrl(:,i,3 ) = qrl(:,i,3 )+ shape(:,i)*tmp*qdi(:,3 ) 102 103 qrl(:,i,4 ) = qrl(:,i,4 )+ shape(:,i)*tmp*qdi(:,4 ) 104 qrl(:,i,5 ) = qrl(:,i,5 )+ shape(:,i)*tmp*qdi(:,5 ) 105 qrl(:,i,6 ) = qrl(:,i,6 )+ shape(:,i)*tmp*qdi(:,6 ) 106 107 qrl(:,i,7 ) = qrl(:,i,7 )+ shape(:,i)*tmp*qdi(:,7 ) 108 qrl(:,i,8 ) = qrl(:,i,8 )+ shape(:,i)*tmp*qdi(:,8 ) 109 qrl(:,i,9 ) = qrl(:,i,9 )+ shape(:,i)*tmp*qdi(:,9 ) 110 enddo 111c 112c.... add contribution to local mass matrix 113c 114 115 if (have_local_mass .eq. 0) then 116 do i=1,nshl 117 do j=1,nshl 118 rmass(:,i,j) = rmass(:,i,j)+shape(:,i)*shape(:,j)*tmp 119 enddo 120 enddo 121 endif 122c 123c.... end of the loop over integration points 124c 125 enddo 126 127c 128c.... find the inverse of the local mass matrix for each element 129 130 131 if (have_local_mass .eq. 0) then 132 allocate (lmassinv(iblock)%p(npro,nshl,nshl)) 133 134 do iel=1,npro 135 do i=1,nshl ! form the identy matrix 136 do j=1,nshl 137 lmassinv(iblock)%p(iel,i,j) = 0.0 138 enddo 139 lmassinv(iblock)%p(iel,i,i)=1.0 140 enddo 141c 142c.... LU factor the mass matrix 143c 144 call ludcmp(rmass(iel,:,:),nshl,nshl,indx,d) 145c 146c.... back substitute with the identy matrix to find the 147c matrix inverse 148c 149 do j=1,nshl 150 call lubksb(rmass(iel,:,:),nshl,nshl,indx, 151 & lmassinv(iblock)%p(iel,:,j)) 152 enddo 153 enddo 154 rminv(:,:,:) = lmassinv(iblock)%p(:,:,:) 155 else 156 rminv(:,:,:) = lmassinv(iblock)%p(:,:,:) 157 endif 158c 159c.... find the modal coefficients of ql by multiplying by the inverse of 160c the local mass matrix 161c 162 do iel=1,npro 163 do j=1,9 164c do j=1, 3*nsd 165 ql(iel,:,j) = matmul( rminv(iel,:,:),qrl(iel,:,j) ) 166 enddo 167 enddo 168c 169c.... return 170c 171 return 172 end 173 174 175 176 177 subroutine e3qlSclr (yl, dwl, shp, shgl, 178 & xl, ql, sgn ) 179c 180c---------------------------------------------------------------------- 181c 182c This routine computes the local diffusive flux vector using a 183c local projection algorithm: 184c diffus * phi,i 185c 186c---------------------------------------------------------------------- 187c 188 use local_mass 189 include "common.h" 190c 191 dimension yl(npro,nshl,ndof), dwl(npro,nshl), 192 & shp(nshl,ngauss), shgl(nsd,nshl,ngauss), 193 & xl(npro,nenl,nsd), sgn(npro,nshl), 194 & ql(npro,nshl,nsd) 195c 196c local arrays 197c 198 dimension dxidx(npro,nsd,nsd), WdetJ(npro), 199 & diffus(npro), 200 & rminv(npro,nshl,nshl), 201 & qrl(npro,nshl,nsd) 202c 203 dimension qdi(npro,nsd), shape(npro,nshl), 204 & shdrv(npro,nsd,nshl), indx(nshl), 205 & rmass(npro,nshl,nshl), gradT(npro,nsd), 206 & eviscv(npro) 207 208c 209c.... loop through the integration points 210c 211 rminv = zero 212 rmass = zero 213 qrl = zero 214 215 do intp = 1, ngauss 216 217 call getshp(shp, shgl, sgn, shape, shdrv) 218 219 qdi = zero 220c 221c.... calculate the integration variables 222c 223c 224 call e3qvarSclr (yl, shdrv, xl, 225 & gradT, dxidx, WdetJ ) 226c 227c.... call function to sort out diffusivity (at end of this file) 228c 229 call getdiffsclr(dwl,shape,yl, diffus) 230c 231c.... diffusive flux in x1-direction 232c 233 qdi(:,1) = diffus * gradT(:,1) 234 qdi(:,2) = diffus * gradT(:,2) 235 qdi(:,3) = diffus * gradT(:,3) 236 237c 238c.... assemble contribution of qdi to qrl,i.e., contribution to 239c each element shape function 240c 241 tmp = Qwt(lcsyst,intp) 242 if (lcsyst .eq. 1) then 243 tmp = tmp*(three/four) 244 endif 245 246 do i=1,nshl 247 qrl(:,i,1 ) = qrl(:,i,1 )+ shape(:,i)*tmp*qdi(:,1 ) 248 qrl(:,i,2 ) = qrl(:,i,2 )+ shape(:,i)*tmp*qdi(:,2 ) 249 qrl(:,i,3 ) = qrl(:,i,3 )+ shape(:,i)*tmp*qdi(:,3 ) 250 enddo 251c 252c.... add contribution to local mass matrix 253c 254 if (have_local_mass .eq. 0) then 255 do i=1,nshl 256 do j=1,nshl 257 rmass(:,i,j)=rmass(:,i,j)+shape(:,i)*shape(:,j)*tmp 258 enddo 259 enddo 260 endif 261 262c.... end of the loop over integration points 263c 264 enddo 265 266c 267c.... find the inverse of the local mass matrix for each element 268c 269 qrl = qrl/6.d0 270c 271c.... Assuming that lmassinv was already computed for flow equations 272c 273 rmass = rmass/6.0 274c 275c.... for cubics, it cannot be precomputed, so compute and 276c save it the first time it is needed 277c 278 if (have_local_mass .eq. 0) then 279 allocate (lmassinv(iblock)%p(npro,nshl,nshl)) 280 281 do iel=1,npro 282 do i=1,nshl ! form the identy matrix 283 do j=1,nshl 284 lmassinv(iblock)%p(iel,i,j) = 0.0 285 enddo 286 lmassinv(iblock)%p(iel,i,i)=1.0 287 enddo 288c 289c.... LU factor the mass matrix 290c 291 call ludcmp(rmass(iel,:,:),nshl,nshl,indx,d) 292c 293c.... back substitute with the identy matrix to find the 294c matrix inverse 295c 296 do j=1,nshl 297 call lubksb(rmass(iel,:,:),nshl,nshl,indx, 298 & lmassinv(iblock)%p(iel,:,j)) 299 enddo 300 enddo 301 rminv(:,:,:) = lmassinv(iblock)%p(:,:,:) 302 else 303 rminv(:,:,:) = lmassinv(iblock)%p(:,:,:) 304 endif 305c 306c.... find the modal coefficients of ql by multiplying by the inverse of 307c the local mass matrix 308c 309 do iel=1,npro 310 do j=1,nsd 311 ql(iel,:,j) = matmul( rminv(iel,:,:),qrl(iel,:,j) ) 312 enddo 313 enddo 314c 315c.... return 316c 317 return 318 end 319