1 subroutine e3juel (yl, acl, sls, A0, 2 & WdetJ, rl, rml) 3c 4c---------------------------------------------------------------------- 5c 6c This routine calculates Exactly integrated Linear Tetrahedra 7c Mass term (assuming U(Y) is linear it is not). 8c 9c input: 10c WdetJ (npro) : weighted Jacobian 11c 12c output: 13c rl (npro,nshl,nflow) : residual 14c rml (npro,nshl,nflow) : modified residual 15c 16c 17c note that this routine wipes out yl by putting ul into it 18c and then (in ires=1 case ) it is used again 19c 20c Kenneth Jansen, Winter 1997, Primitive Variables 21c---------------------------------------------------------------------- 22c 23 include "common.h" 24c 25 dimension yl(npro,nshl,nflow), acl(npro,nshl,ndof), 26 & WdetJ(npro), A0(npro,nflow,nflow), 27 & rl(npro,nshl,nflow), rml(npro,nshl,nflow) 28c 29 dimension rk(npro), rho(npro), 30 & ei(npro), tmp(npro), 31 & ub(npro,nflow), fact(npro), 32 & fddt(npro) 33 34 ttim(28) = ttim(28) - secs(0.0) 35 36c 37c.... ---------------------> Time term <-------------------- 38c 39c.... add contribution of U to rml 40c 41c.... compute conservative variables 42c 43c 44c multiply by exact mass matrix integral N_aN_b=(I+1)*V/20 45c where 1 is a matrix with every element=1 46c 47c note that the wght has 4/3 multiplier so 3/4*20=15 48c 49 50 fact=WdetJ/(Qwt(lcsyst,intp)*15.0d0) 51 fct1=almi/(gami*alfi)*Dtgl ! factor for predictor (scalar) 52 if(ires.ne.1) then 53 fddt=fact*fct1 54 do inod=1,nshl 55 rk = pt5 * (yl(:,inod,2)**2 + yl(:,inod,3)**2 + yl(:,inod,4)**2) 56c 57 ithm = 6 58 call getthm (yl(:,inod,1), yl(:,inod,5), sls, 59 & rk, rho, ei, 60 & tmp, tmp, tmp, 61 & tmp, tmp, tmp, 62 & tmp, tmp) 63c 64 yl(:,inod,1) = rho 65 yl(:,inod,2) = rho * yl(:,inod,2) 66 yl(:,inod,3) = rho * yl(:,inod,3) 67 yl(:,inod,4) = rho * yl(:,inod,4) 68 yl(:,inod,5) = rho * (ei + rk) 69 enddo 70 ub(:,:)=yl(:,1,:)+yl(:,2,:) 71 & +yl(:,3,:)+yl(:,4,:) 72 73c 74c what we have now in yl is the U_b^e 75c we want to get Resm(mass)=M^e_{ab} (U^e_b(Y+epsilon P) 76c - U^e_b(Y))*fact/dt, 77c since only the difference between resm's is important we 78c do not have to subtract off the unperturbed U(Y) vector 79c This term is meant to carry through the effect of a perturbation 80c in Y upon dY/dt (through the predictor into the term fact) 81c 82c 83 84 do i = 1, nshl 85 rml(:,i,1) = rml(:,i,1) + fddt * (yl(:,i,1)+ub(:,1)) 86 rml(:,i,2) = rml(:,i,2) + fddt * (yl(:,i,2)+ub(:,2)) 87 rml(:,i,3) = rml(:,i,3) + fddt * (yl(:,i,3)+ub(:,3)) 88 rml(:,i,4) = rml(:,i,4) + fddt * (yl(:,i,4)+ub(:,4)) 89 rml(:,i,5) = rml(:,i,5) + fddt * (yl(:,i,5)+ub(:,5)) 90 enddo 91c 92! flops = flops + 35*nshl*npro 93 endif 94 95c 96c.... ires = 2 or 3 97c 98 if ((ires .eq. 1) .or. (ires .eq. 3)) then 99 100 ub(:,:)=acl(:,1,:)+acl(:,2,:) 101 & +acl(:,3,:)+acl(:,4,:) 102 103 do i = 1, nshl 104 yl(:,i,1) = fact*(acl(:,i,1)+ub(:,1)) 105 yl(:,i,2) = fact*(acl(:,i,2)+ub(:,2)) 106 yl(:,i,3) = fact*(acl(:,i,3)+ub(:,3)) 107 yl(:,i,4) = fact*(acl(:,i,4)+ub(:,4)) 108 yl(:,i,5) = fact*(acl(:,i,5)+ub(:,5)) 109 enddo 110c 111c what we have now in yl is the dY_a^e/dt=M^e_{ab} Y_{b,t}, must multiply by 112c A0 to get dU^e_a/dt= dU/dY(centroid) dY_a^e/dt, take advantage of zeros 113c in A0(Prim) with comments 114c 115 do i = 1, nshl 116 rl(:,i,1) = rl(:,i,1) 117 & + A0(:,1,1)*yl(:,i,1) 118c & + A0(:,1,2)*yl(:,i,2) 119c & + A0(:,1,3)*yl(:,i,3) 120c & + A0(:,1,4)*yl(:,i,4) 121 & + A0(:,1,5)*yl(:,i,5) 122c 123 rl(:,i,2) = rl(:,i,2) 124 & + A0(:,2,1)*yl(:,i,1) 125 & + A0(:,2,2)*yl(:,i,2) 126c & + A0(:,2,3)*yl(:,i,3) 127c & + A0(:,2,4)*yl(:,i,4) 128 & + A0(:,2,5)*yl(:,i,5) 129c 130 rl(:,i,3) = rl(:,i,3) 131 & + A0(:,3,1)*yl(:,i,1) 132c & + A0(:,3,2)*yl(:,i,2) 133 & + A0(:,3,3)*yl(:,i,3) 134c & + A0(:,3,4)*yl(:,i,4) 135 & + A0(:,3,5)*yl(:,i,5) 136c 137 rl(:,i,4) = rl(:,i,4) 138 & + A0(:,4,1)*yl(:,i,1) 139c & + A0(:,4,2)*yl(:,i,2) 140c & + A0(:,4,3)*yl(:,i,3) 141 & + A0(:,4,4)*yl(:,i,4) 142 & + A0(:,4,5)*yl(:,i,5) 143c 144 rl(:,i,5) = rl(:,i,5) 145 & + A0(:,5,1)*yl(:,i,1) 146 & + A0(:,5,2)*yl(:,i,2) 147 & + A0(:,5,3)*yl(:,i,3) 148 & + A0(:,5,4)*yl(:,i,4) 149 & + A0(:,5,5)*yl(:,i,5) 150c 151 enddo 152c 153! flops = flops + 45*nshl*npro 154 endif 155 156 ttim(28) = ttim(28) + tmr() 157c 158c.... return 159c 160 return 161 end 162 163c$$$ 164c$$$ subroutine e3juelSclr (ycl, acl, A0t, 165c$$$ & WdetJ, rtl, rmtl) 166c$$$c 167c$$$c---------------------------------------------------------------------- 168c$$$c 169c$$$c This routine calculates Exactly integrated Linear Tetrahedra 170c$$$c Mass term (assuming U(Y) is linear it is not). 171c$$$c 172c$$$c input: 173c$$$c WdetJ (npro) : weighted Jacobian 174c$$$c 175c$$$c output: 176c$$$c rtl (npro,nshl,nflow) : residual 177c$$$c rmtl (npro,nshl,nflow) : modified residual 178c$$$c 179c$$$c 180c$$$c note that this routine wipes out ycl by putting ul into it 181c$$$c and then (in ires=1 case ) it is used again 182c$$$c 183c$$$c Kenneth Jansen, Winter 1997, Primitive Variables 184c$$$c---------------------------------------------------------------------- 185c$$$c 186c$$$ include "common.h" 187c$$$c 188c$$$ dimension ycl(npro,nshl,ndof), acl(npro,nshl,ndof), 189c$$$ & WdetJ(npro), A0t(npro), 190c$$$ & rtl(npro,nshl), rmtl(npro,nshl) 191c$$$ 192c$$$c 193c$$$ dimension rk(npro), rho(npro), 194c$$$ & ei(npro), tmp(npro), 195c$$$ & ubt(npro), fact(npro), 196c$$$ & fddt(npro) 197c$$$ 198c$$$ ttim(28) = ttim(28) - tmr() 199c$$$ 200c$$$c 201c$$$c.... ---------------------> Time term <-------------------- 202c$$$c 203c$$$c.... add contribution of U to rml 204c$$$c 205c$$$c.... compute conservative variables 206c$$$c 207c$$$c 208c$$$c multiply by exact mass matrix integral N_aN_b=(I+1)*V/20 209c$$$c where 1 is a matrix with every element=1 210c$$$c 211c$$$c note that the wght has 4/3 multiplier so 3/4*20=15 212c$$$c 213c$$$ 214c$$$ fact=WdetJ/(Qwt(lcsyst,intp)*15.0d0) 215c$$$ fct1=almi/(gami*alfi)*Dtgl ! factor for predictor (scaler) 216c$$$c 217c$$$c 218c$$$c.... ires = 2 or 3 219c$$$c 220c$$$ 221c$$$ if ((ires .eq. 1) .or. (ires .eq. 3)) then 222c$$$ ubt(:)=acl(:,1,id)+acl(:,2,id) 223c$$$ & +acl(:,3,id)+acl(:,4,id) 224c$$$ do i = 1, nshl 225c$$$ ycl(:,i,id) = fact*(acl(:,i,id)+ubt(:)) 226c$$$ 227c$$$ enddo 228c$$$c 229c$$$c what we have now in ycl is the dY_a^e/dt=M^e_{ab} Y_{b,t}, must multiply by 230c$$$c A0 to get dU^e_a/dt= dU/dY(centroid) dY_a^e/dt, take advantage of zeros 231c$$$c in A0(Prim) with comments 232c$$$c 233c$$$ do i = 1, nshl 234c$$$ rtl(:,i) = rtl(:,i) 235c$$$ & + A0t(:)*ycl(:,i,id) 236c$$$ 237c$$$c 238c$$$ enddo 239c$$$c 240c$$$ ! flops = flops + 45*nenl*npro 241c$$$ endif 242c$$$ 243c$$$ ttim(28) = ttim(28) + tmr() 244c$$$c 245c$$$c.... return 246c$$$c 247c$$$ return 248c$$$ end 249