1 subroutine e3ql (ycl, 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 ycl (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,nshape,(nflow-1)*nsd) : element RHS diffusion residual 19c 20c---------------------------------------------------------------------- 21c 22 include "common.h" 23c 24 dimension ycl(npro,nshl,ndof), 25 & shp(nshl,ngauss), 26 & 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,nflow), g2yi(npro,nflow), 33 & g3yi(npro,nflow), shg(npro,nshl,nsd), 34 & dxidx(npro,nsd,nsd), WdetJ(npro), 35 & T(npro), cp(npro), 36 & u1(npro), u2(npro), 37 & u3(npro), rmu(npro), 38 & rlm(npro), rlm2mu(npro), 39 & con(npro), rminv(npro,nshl,nshl), 40 & qrl(npro,nshl,(nflow-1)*nsd) 41c 42 dimension qdi(npro,nsd*(nflow-1)), shape(npro,nshl), 43 & shdrv(npro,nsd,nshl), indx(nshl), 44 & rmass(npro,nshl,nshl), rho(npro) 45 46 real*8 tmp(npro) 47c 48c.... loop through the integration points 49c 50 rminv = zero 51 rmass = zero 52 qrl = zero 53 54 do intp = 1, ngauss 55 if (Qwt(lcsyst,intp) .eq. zero) cycle ! precaution 56c 57c.... create a matrix of shape functions (and derivatives) for each 58c element at this quadrature point. These arrays will contain 59c the correct signs for the hierarchic basis 60c 61 call getshp(shp, shgl, sgn, 62 & shape, shdrv) 63c 64c 65c.... initialize 66c 67 qdi = zero 68c 69c 70c.... calculate the integration variables necessary for the 71c formation of q 72c 73 74 call e3qvar (ycl, shape, shdrv, 75 & rho, xl, g1yi, 76 & g2yi, g3yi, shg, 77 & dxidx, WdetJ, T, 78 & cp, u1, u2, 79 & u3 ) 80c 81c 82c.... compute diffusive flux vector at this integration point 83c 84c 85c.... get material properties 86c 87 call getDiff (T, cp, rho, ycl, 88 & rmu, rlm, rlm2mu, con, shape, 89 & xmudmi, xl) 90c 91c.... compute diffusive fluxes 92c 93c.... diffusive flux in x1-direction 94c 95 qdi(:,1) = rlm2mu * g1yi(:,2) 96 & + rlm * g2yi(:,3) 97 & + rlm * g3yi(:,4) 98 qdi(:,2) = rmu * g1yi(:,3) 99 & + rmu * g2yi(:,2) 100 qdi(:,3) = rmu * g1yi(:,4) 101 & + rmu * g3yi(:,2) 102 qdi(:,4) = rlm2mu * u1 * g1yi(:,2) + rmu * u2 * g1yi(:,3) 103 & + rmu * u3 * g1yi(:,4) 104 & + rmu * u2 * g2yi(:,2) + rlm * u1 * g2yi(:,3) 105 & + rmu * u3 * g3yi(:,2) + rlm * u1 * g3yi(:,4) 106 & + con * g1yi(:,5) 107 108c 109c.... diffusive flux in x2-direction 110c 111 qdi(:,5) = rmu * g1yi(:,3) 112 & + rmu * g2yi(:,2) 113 qdi(:,6) = rlm * g1yi(:,2) 114 & + rlm2mu * g2yi(:,3) 115 & + rlm * g3yi(:,4) 116 qdi(:,7) = rmu * g2yi(:,4) 117 & + rmu * g3yi(:,3) 118 qdi(:,8) = rlm * u2 * g1yi(:,2) + rmu * u1 * g1yi(:,3) 119 & + rmu * u1 * g2yi(:,2) + rlm2mu * u2 * g2yi(:,3) 120 & + rmu * u3 * g2yi(:,4) 121 & + rmu * u3 * g3yi(:,3) + rlm * u2 * g3yi(:,4) 122 & + con * g2yi(:,5) 123c 124c.... diffusive flux in x3-direction 125c 126 qdi(:,9 ) = rmu * g1yi(:,4) 127 & + rmu * g3yi(:,2) 128 qdi(:,10) = rmu * g2yi(:,4) 129 & + rmu * g3yi(:,3) 130 qdi(:,11) = rlm * g1yi(:,2) 131 & + rlm * g2yi(:,3) 132 & + rlm2mu * g3yi(:,4) 133 qdi(:,12) = rlm * u3 * g1yi(:,2) + rmu * u1 * g1yi(:,4) 134 & + rlm * u3 * g2yi(:,3) + rmu * u2 * g2yi(:,4) 135 & + rmu * u1 * g3yi(:,2) + rmu * u2 * g3yi(:,3) 136 & + rlm2mu * u3 * g3yi(:,4) 137 & + con * g3yi(:,5) 138 139c 140c 141c.... assemble contribution of qdi to qrl,i.e., contribution to 142c each element shape function 143c 144 tmp = Qwt(lcsyst,intp) 145 if ( lcsyst .eq. 1) then 146 tmp = tmp*(three/four) 147 end if 148 149 do i=1,nshl 150 qrl(:,i,1 ) = qrl(:,i,1 )+ shape(:,i)*tmp*qdi(:,1 ) 151 qrl(:,i,2 ) = qrl(:,i,2 )+ shape(:,i)*tmp*qdi(:,2 ) 152 qrl(:,i,3 ) = qrl(:,i,3 )+ shape(:,i)*tmp*qdi(:,3 ) 153 qrl(:,i,4 ) = qrl(:,i,4 )+ shape(:,i)*tmp*qdi(:,4 ) 154 qrl(:,i,5 ) = qrl(:,i,5 )+ shape(:,i)*tmp*qdi(:,5 ) 155 qrl(:,i,6 ) = qrl(:,i,6 )+ shape(:,i)*tmp*qdi(:,6 ) 156 qrl(:,i,7 ) = qrl(:,i,7 )+ shape(:,i)*tmp*qdi(:,7 ) 157 qrl(:,i,8 ) = qrl(:,i,8 )+ shape(:,i)*tmp*qdi(:,8 ) 158 qrl(:,i,9 ) = qrl(:,i,9 )+ shape(:,i)*tmp*qdi(:,9 ) 159 qrl(:,i,10) = qrl(:,i,10)+ shape(:,i)*tmp*qdi(:,10) 160 qrl(:,i,11) = qrl(:,i,11)+ shape(:,i)*tmp*qdi(:,11) 161 qrl(:,i,12) = qrl(:,i,12)+ shape(:,i)*tmp*qdi(:,12) 162 enddo 163c 164c.... add contribution to local mass matrix 165c 166 do i=1,nshl 167 do j=1,nshl 168 rmass(:,i,j) = rmass(:,i,j)+shape(:,i)*shape(:,j)*tmp 169 enddo 170 enddo 171c 172c.... end of the loop over integration points 173c 174 enddo 175 if ( lcsyst .eq. 1) then 176 qrl = qrl/6.d0 177 rmass = rmass/6.0 178 end if 179c 180c.... find the inverse of the local mass matrix for each element 181c 182 do iel=1,npro 183 184 do i=1,nshl ! form the identy matrix 185 do j=1,nshl 186 rminv(iel,i,j) = 0.0 187 enddo 188 rminv(iel,i,i)=1.0 189 enddo 190c 191c.... LU factor the mass matrix 192c 193 call ludcmp(rmass(iel,:,:),nshl,nshl,indx,d) 194c 195c.... back substitute with the identy matrix to find the 196c matrix inverse 197c 198 do j=1,nshl 199 call lubksb(rmass(iel,:,:),nshl,nshl, 200 & indx,rminv(iel,:,j)) 201 enddo 202 enddo 203 204c 205c.... find the modal coefficients of ql by multiplying by the inverse of 206c the local mass matrix 207c 208 do iel=1,npro 209 do j=1,12 210 ql(iel,:,j) = matmul( rminv(iel,:,:),qrl(iel,:,j) ) 211 enddo 212 enddo 213c 214c.... return 215c 216 return 217 end 218