1 subroutine e3mtrx (rho, pres, T, 2 & ei, h, alfap, 3 & betaT, cp, rk, 4 & u1, u2, u3, 5 & A0, A1, 6 & A2, A3, 7 & e2p, e3p, e4p, 8 & drdp, drdT, A0DC, 9 & A0inv, dVdY) 10c 11c---------------------------------------------------------------------- 12c 13c This routine sets up the necessary matrices at the integration point. 14c 15c input: 16c rho (npro) : density 17c pres (npro) : pressure 18c T (npro) : temperature 19c ei (npro) : internal energy 20c h (npro) : enthalpy 21c alfap (npro) : expansivity 22c betaT (npro) : isothermal compressibility 23c cp (npro) : specific heat at constant pressure 24c c (npro) : speed of sound 25c rk (npro) : kinetic energy 26c u1 (npro) : x1-velocity component 27c u2 (npro) : x2-velocity component 28c u3 (npro) : x3-velocity component 29c 30c output: 31c A0 (npro,nflow,nflow) : A0 matrix 32c A1 (npro,nflow,nflow) : A_1 matrix 33c A2 (npro,nflow,nflow) : A_2 matrix 34c A3 (npro,nflow,nflow) : A_3 matrix 35c 36c Note: the definition of the matrices can be found in 37c thesis by Hauke. 38c 39c Zdenek Johan, Summer 1990. (Modified from e2mtrx.f) 40c Zdenek Johan, Winter 1991. (Fortran 90) 41c Kenneth Jansen, Winter 1997 Primitive Variables 42c---------------------------------------------------------------------- 43c 44 include "common.h" 45c 46c 47c passed arrays 48c 49 dimension rho(npro), pres(npro), 50 & T(npro), ei(npro), 51 & h(npro), alfap(npro), 52 & betaT(npro), 53 & cp(npro), 54 & rk(npro), 55 & u1(npro), u2(npro), 56 & u3(npro), fact1(npro), 57 & A0(npro,nflow,nflow), dVdY(npro,15), 58 & A1(npro,nflow,nflow), A2(npro,nflow,nflow), 59 & A3(npro,nflow,nflow), A0DC(npro,4), 60 & A0inv(npro,15), d(npro), 61 & fact2(npro), s1(npro), 62 & e1bar(npro), e2bar(npro), 63 & e3bar(npro), e4bar(npro), 64 & e5bar(npro), c1bar(npro), 65 & c2bar(npro), cv(npro), 66 & c3bar(npro), u12(npro), 67 & u31(npro), u23(npro) 68c 69c local work arrays that are passed shared space 70c 71 dimension e2p(npro), 72 & e3p(npro), e4p(npro), 73 & drdp(npro), drdT(npro) 74 75 ttim(21) = ttim(21) - secs(0.0) 76c 77c.... initialize 78c 79 A0 = zero 80 A1 = zero 81 A2 = zero 82 A3 = zero 83c 84c.... set up the constants 85c 86c 87 drdp = rho * betaT 88 drdT = -rho * alfap 89 A0(:,5,1) = drdp * (h + rk) - alfap * T ! e1p 90c A0(:,5,1) = drdp * (ei + rk) + betaT * pres - alfap * T ! e1p 91 e2p = A0(:,5,1) + one 92 e3p = rho * ( h + rk) 93 e4p = drdT * (h + rk) + rho * cp 94c 95c 96c.... Calculate A0 97c 98 A0(:,1,1) = drdp 99c A0(:,1,2) = zero 100c A0(:,1,3) = zero 101c A0(:,1,4) = zero 102 A0(:,1,5) = drdT 103c 104 A0(:,2,1) = drdp * u1 105 A0(:,2,2) = rho 106c A0(:,2,3) = zero 107c A0(:,2,4) = zero 108 A0(:,2,5) = drdT * u1 109c 110 A0(:,3,1) = drdp * u2 111c A0(:,3,2) = zero 112 A0(:,3,3) = rho 113c A0(:,3,4) = zero 114 A0(:,3,5) = drdT * u2 115c 116 A0(:,4,1) = drdp * u3 117c A0(:,4,2) = zero 118c A0(:,4,3) = zero 119 A0(:,4,4) = rho 120 A0(:,4,5) = drdT * u3 121c 122covered above A0(:,5,1) = drdp * u1 123 A0(:,5,2) = rho * u1 124 A0(:,5,3) = rho * u2 125 A0(:,5,4) = rho * u3 126 A0(:,5,5) = e4p 127c 128 ! flops = flops + 67*npro 129c 130c.... Calculate A-tilde-1, A-tilde-2 and A-tilde-3 131c 132 A1(:,1,1) = drdp * u1 133 A1(:,1,2) = rho 134c A1(:,1,3) = zero 135c A1(:,1,4) = zero 136 A1(:,1,5) = drdT * u1 137c 138 A1(:,2,1) = drdp * u1 * u1 +1 139 A1(:,2,2) = two *rho * u1 140c A1(:,2,3) = zero 141c A1(:,2,4) = zero 142 A1(:,2,5) = drdT * u1 * u1 143c 144 A1(:,3,1) = drdp * u1 * u2 145 A1(:,3,2) = rho * u2 146 A1(:,3,3) = rho * u1 147c A1(:,3,4) = zero 148 A1(:,3,5) = drdT * u1 * u2 149c 150 A1(:,4,1) = drdp * u1 * u3 151 A1(:,4,2) = rho * u3 152c A1(:,4,3) = zero 153 A1(:,4,4) = rho * u1 154 A1(:,4,5) = drdT * u1 * u3 155c 156 A1(:,5,1) = u1 * e2p 157 A1(:,5,2) = e3p + rho * u1 * u1 158 A1(:,5,3) = rho * u1 * u2 159 A1(:,5,4) = rho * u1 * u3 160 A1(:,5,5) = u1 * e4p 161c 162 ! flops = flops + 35*npro 163c 164 A2(:,1,1) = drdp * u2 165c A2(:,1,2) = zero 166 A2(:,1,3) = rho 167c A2(:,1,4) = zero 168 A2(:,1,5) = drdT * u2 169c 170 A2(:,2,1) = drdp * u1 * u2 171 A2(:,2,2) = rho * u2 172 A2(:,2,3) = rho * u1 173c A2(:,2,4) = zero 174 A2(:,2,5) = drdT * u1 * u2 175c 176 A2(:,3,1) = drdp * u2 * u2 +1 177c A2(:,3,2) = zero 178 A2(:,3,3) = two * rho * u2 179c A2(:,3,4) = zero 180 A2(:,3,5) = drdT * u2 * u2 181c 182 A2(:,4,1) = drdp * u2 * u3 183c A2(:,4,2) = zero 184 A2(:,4,3) = rho * u3 185 A2(:,4,4) = rho * u2 186 A2(:,4,5) = drdT * u2 * u3 187c 188 A2(:,5,1) = u2 * e2p 189 A2(:,5,2) = rho * u1 * u2 190 A2(:,5,3) = e3p + rho * u2 * u2 191 A2(:,5,4) = rho * u2 * u3 192 A2(:,5,5) = u2 * e4p 193c 194 ! flops = flops + 35*npro 195c 196 A3(:,1,1) = drdp * u3 197c A3(:,1,2) = zero 198c A3(:,1,3) = zero 199 A3(:,1,4) = rho 200 A3(:,1,5) = drdT * u3 201c 202 A3(:,2,1) = drdp * u1 * u3 203 A3(:,2,2) = rho * u3 204c A3(:,2,3) = zero 205 A3(:,2,4) = rho * u1 206 A3(:,2,5) = drdT * u1 * u3 207c 208 A3(:,3,1) = drdp * u3 * u2 209c A3(:,3,2) = zero 210 A3(:,3,3) = rho * u3 211 A3(:,3,4) = rho * u2 212 A3(:,3,5) = drdT * u3 * u2 213c 214 A3(:,4,1) = drdp * u3 * u3 +1 215c A3(:,4,2) = zero 216c A3(:,4,3) = zero 217 A3(:,4,4) = two *rho * u3 218 A3(:,4,5) = drdT * u3 * u3 219c 220 A3(:,5,1) = u3 * e2p 221 A3(:,5,2) = rho * u1 * u3 222 A3(:,5,3) = rho * u2 * u3 223 A3(:,5,4) = e3p + rho * u3 * u3 224 A3(:,5,5) = u3 * e4p 225c 226 ! flops = flops + 35*npro 227 ttim(21) = ttim(21) + secs(0.0) 228 229c 230c.... return 231c 232 if (idc .ne. 0) then 233c.... for Discountinuity Capturing Term 234c 235c.... calculation of A0^DC matrix 236c 237c.... Ref P-163 of the handout 238c 239 s1 = one/(rho**2 * betaT * T) 240 cv = cp - (alfap**2 * T/rho/betaT) 241 A0DC(:,1) = (rho * betaT)**2*s1 242 A0DC(:,2) = -rho*alfap*rho*betaT*s1 243 A0DC(:,3) = rho/T 244 A0DC(:,4) = (-rho*alfap)**2 * s1 + (rho*cv/T**2) 245c 246c.... calculation of A0^tilda^inverse matrix 247c.... ref P-169 of the hand out 248c 249 fact1 = one/(rho*cv*T**2) 250 d = alfap*T/rho/betaT 251 e1bar = h - rk 252 e2bar = e1bar - d 253 e3bar = e2bar - cv * T 254 e4bar = e2bar - 2* cv * T 255 e5bar = e1bar**2 - 2*e1bar*d + 2*rk*cv*T + cp*T/rho/betaT 256 c1bar = u1**2 + cv * T 257 c2bar = u2**2 + cv * T 258 c3bar = u3**2 + cv * T 259 u12 = u1 * u2 260 u31 = u3 * u1 261 u23 = u2 * u3 262 A0inv( :,1) = e5bar*fact1 263 A0inv( :,2) = c1bar*fact1 264 A0inv( :,3) = c2bar*fact1 265 A0inv( :,4) = c3bar*fact1 266 A0inv( :,5) = 1*fact1 267 A0inv( :,6) = u1*e3bar*fact1 268 A0inv( :,7) = u2*e3bar*fact1 269 A0inv( :,8) = u3*e3bar*fact1 270 A0inv( :,9) = -e2bar*fact1 271 A0inv(:,10) = u12*fact1 272 A0inv(:,11) = u31*fact1 273 A0inv(:,12) = -u1*fact1 274 A0inv(:,13) = u23*fact1 275 A0inv(:,14) = -u2*fact1 276 A0inv(:,15) = -u3*fact1 277c 278c.....calculation of dV/dY (derivative of entropy variables w.r.to primitive 279c 280 fact1 = 1/T 281 fact2 = fact1/T 282 dVdY( :,1) = fact1/rho 283 dVdY( :,2) = -fact1*u1 284 dVdY( :,3) = fact1 285 dVdY( :,4) = -fact1*u2 286 dVdY( :,5) = zero 287 dVdY( :,6) = fact1 288 dVdY( :,7) = -fact1*u3 289 dVdY( :,8) = zero 290 dVdY( :,9) = zero 291 dVdY(:,10) = fact1 292 dVdY(:,11) = -(h-rk)*fact2 293 dVdY(:,12) = -fact2*u1 294 dVdY(:,13) = -fact2*u2 295 dVdY(:,14) = -fact2*u3 296 dVdY(:,15) = fact2 297 298 endif !end of idc.ne.0 299 300 return 301 end 302c 303c 304 subroutine e3mtrxSclr (rho, 305 & u1, u2, u3, 306 & A0t, A1t, 307 & A2t, A3t ) 308c 309c---------------------------------------------------------------------- 310c 311c This routine sets up the necessary matrices at the integration point. 312c 313c input: 314c rho (npro) : density 315c u1 (npro) : x1-velocity component 316c u2 (npro) : x2-velocity component 317c u3 (npro) : x3-velocity component 318c 319c output: 320c A0t (npro) : A_0 "matrix" 321c A1t (npro) : A_1 "matrix" 322c A2t (npro) : A_2 "matrix" 323c A3t (npro) : A_3 "matrix" 324c 325c Note: the definition of the matrices can be found in 326c thesis by Hauke. 327c 328c Zdenek Johan, Summer 1990. (Modified from e2mtrx.f) 329c Zdenek Johan, Winter 1991. (Fortran 90) 330c Kenneth Jansen, Winter 1997 Primitive Variables 331c---------------------------------------------------------------------- 332c 333 include "common.h" 334c 335c 336c passed arrays 337c 338 dimension rho(npro), 339 & u1(npro), u2(npro), 340 & u3(npro), 341 & A0t(npro), 342 & A1t(npro), A2t(npro), 343 & A3t(npro) 344c 345 if (iconvsclr.eq.2) then !convective form 346 A0t(:) = one 347 A1t(:) = u1(:) 348 A2t(:) = u2(:) 349 A3t(:) = u3(:) 350 else !conservative form 351 A0t(:) = rho(:) 352 A1t(:) = rho(:)*u1(:) 353 A2t(:) = rho(:)*u2(:) 354 A3t(:) = rho(:)*u3(:) 355 endif 356 357c 358c.... return 359c 360 return 361 end 362 363 364