1 subroutine e3eig1 (rho, T, cp, gamb, c, 2 & u1, u2, u3, a1, a2, 3 & a3, eb1, 4 & dxidx, u, Q) 5c 6c---------------------------------------------------------------------- 7c 8c This routine performs the first step of the eigenvalue decomposition 9c of the Tau matrix. 10c 11c input: 12c rho (npro) : density 13c T (npro) : temperature 14c cp (npro) : specific heat at constant pressure 15c gamb (npro) : gamma_bar (defined in paper by Chalot et al.) 16c c (npro) : speed of sound 17c u1 (npro) : x1-velocity component 18c u2 (npro) : x2-velocity component 19c u3 (npro) : x3-velocity component 20c a1 (npro) : x1-acceleration component 21c a2 (npro) : x2-acceleration component 22c a3 (npro) : x3-acceleration component 23c eb1 (npro) : e1_bar (defined in paper by Chalot et al.) 24c dxidx (npro,nsd,nsd) : inverse of deformation gradient 25c 26c output: 27c u (npro) : fluid velocity 28c a1 (npro) : aspect ratio factor in streamline direction 29c a2 (npro) : aspect ratio factor in normal_1 direction 30c a3 (npro) : aspect ratio factor in normal_2 direction 31c Q (npro,nflow,nflow) : 1st level eigenvectors of Tau matrix 32c 33c 34c Zdenek Johan, Summer 1990. (Modified from e2tau.f) 35c Zdenek Johan, Winter 1991. (Fortran 90) 36c---------------------------------------------------------------------- 37c 38 include "common.h" 39c 40 dimension rho(npro), T(npro), 41 & cp(npro), gamb(npro), 42 & c(npro), u1(npro), 43 & u2(npro), u3(npro), 44 & a1(npro), a2(npro), 45 & a3(npro), 46 & eb1(npro), dxidx(npro,nsd,nsd), 47 & u(npro), Q(npro,nflow,nflow) 48c 49 dimension Rcs(npro,nsd,nsd), fact(npro), 50 & temp(npro) 51c 52c.... compute the directional cosines (streamline direction) 53c 54c 55 56 where (u .ne. zero) 57 fact = one / u 58 Rcs(:,1,1) = u1 * fact 59 Rcs(:,1,2) = u2 * fact 60 Rcs(:,1,3) = u3 * fact 61 elsewhere 62 Rcs(:,1,1) = one 63 Rcs(:,1,2) = zero 64 Rcs(:,1,3) = zero 65 endwhere 66c 67c.... compute the directional cosines (normal acceleration direction) 68c 69 70 71 fact = a1 * Rcs(:,1,1) + a2 * Rcs(:,1,2) + a3 * Rcs(:,1,3) 72 a1 = a1 - fact * Rcs(:,1,1) 73 a2 = a2 - fact * Rcs(:,1,2) 74 a3 = a3 - fact * Rcs(:,1,3) 75 fact = a1**2 + a2**2 + a3**2 76c 77 where (fact .gt. epsM) 78c 79 Rcs(:,2,1) = a1 80 Rcs(:,2,2) = a2 81 Rcs(:,2,3) = a3 82c 83 elsewhere 84c 85 Rcs(:,2,1) = Rcs(:,1,2) 86 & + sign(one, Rcs(:,1,2)*Rcs(:,1,3)) * Rcs(:,1,3) 87 Rcs(:,2,2) = - Rcs(:,1,1) 88 & - sign(one, Rcs(:,1,2)*Rcs(:,1,3)) * Rcs(:,1,3) 89 Rcs(:,2,3) = sign(one, Rcs(:,1,2)*Rcs(:,1,3)) * 90 & ( Rcs(:,1,2) - Rcs(:,1,1) ) 91c 92 fact = Rcs(:,2,1)**2 + Rcs(:,2,2)**2 + Rcs(:,2,3)**2 93c 94 endwhere 95c 96 fact = one / sqrt(fact) 97c 98 Rcs(:,2,1) = Rcs(:,2,1) * fact 99 Rcs(:,2,2) = Rcs(:,2,2) * fact 100 Rcs(:,2,3) = Rcs(:,2,3) * fact 101c 102c.... compute the directional cosines (last direction) 103c 104 Rcs(:,3,1) = Rcs(:,1,2) * Rcs(:,2,3) - Rcs(:,1,3) * Rcs(:,2,2) 105 Rcs(:,3,2) = Rcs(:,1,3) * Rcs(:,2,1) - Rcs(:,1,1) * Rcs(:,2,3) 106 Rcs(:,3,3) = Rcs(:,1,1) * Rcs(:,2,2) - Rcs(:,1,2) * Rcs(:,2,1) 107 108c 109c.... calculate the element aspect ratio factors 110c 111 a1 = Rcs(:,1,1) * dxidx(:,1,1) + Rcs(:,1,2) * dxidx(:,1,2) + 112 & Rcs(:,1,3) * dxidx(:,1,3) 113 a2 = Rcs(:,2,1) * dxidx(:,1,1) + Rcs(:,2,2) * dxidx(:,1,2) + 114 & Rcs(:,2,3) * dxidx(:,1,3) 115 a3 = Rcs(:,3,1) * dxidx(:,1,1) + Rcs(:,3,2) * dxidx(:,1,2) + 116 & Rcs(:,3,3) * dxidx(:,1,3) 117 dxidx(:,1,1) = a1 118 dxidx(:,1,2) = a2 119 dxidx(:,1,3) = a3 120c 121 a1 = Rcs(:,1,1) * dxidx(:,2,1) + Rcs(:,1,2) * dxidx(:,2,2) + 122 & Rcs(:,1,3) * dxidx(:,2,3) 123 a2 = Rcs(:,2,1) * dxidx(:,2,1) + Rcs(:,2,2) * dxidx(:,2,2) + 124 & Rcs(:,2,3) * dxidx(:,2,3) 125 a3 = Rcs(:,3,1) * dxidx(:,2,1) + Rcs(:,3,2) * dxidx(:,2,2) + 126 & Rcs(:,3,3) * dxidx(:,2,3) 127 dxidx(:,2,1) = a1 128 dxidx(:,2,2) = a2 129 dxidx(:,2,3) = a3 130c 131 a1 = Rcs(:,1,1) * dxidx(:,3,1) + Rcs(:,1,2) * dxidx(:,3,2) + 132 & Rcs(:,1,3) * dxidx(:,3,3) 133 a2 = Rcs(:,2,1) * dxidx(:,3,1) + Rcs(:,2,2) * dxidx(:,3,2) + 134 & Rcs(:,2,3) * dxidx(:,3,3) 135 a3 = Rcs(:,3,1) * dxidx(:,3,1) + Rcs(:,3,2) * dxidx(:,3,2) + 136 & Rcs(:,3,3) * dxidx(:,3,3) 137 dxidx(:,3,1) = a1 138 dxidx(:,3,2) = a2 139 dxidx(:,3,3) = a3 140 141c 142c... original 143c 144 145 a1 = dxidx(:,1,1)**2 + dxidx(:,2,1)**2 + dxidx(:,3,1)**2 146 a2 = dxidx(:,1,2)**2 + dxidx(:,2,2)**2 + dxidx(:,3,2)**2 147 a3 = dxidx(:,1,3)**2 + dxidx(:,2,3)**2 + dxidx(:,3,3)**2 148 149c 150c... change from original (analyt., could be error) 151c 152 153cc a1 = dxidx(:,1,1)**2 + dxidx(:,1,2)**2 + dxidx(:,1,3)**2 154cc a2 = dxidx(:,2,1)**2 + dxidx(:,2,2)**2 + dxidx(:,2,3)**2 155cc a3 = dxidx(:,3,1)**2 + dxidx(:,3,2)**2 + dxidx(:,3,3)**2 156 157c 158c.... correct for tetrahedra 159c 160 if (lcsyst .eq. 1) then 161 a1 = ( a1 + (dxidx(:,1,1) + dxidx(:,2,1) + 162 & dxidx(:,3,1))**2 ) * pt39 163 a2 = ( a2 + (dxidx(:,1,2) + dxidx(:,2,2) + 164 & dxidx(:,3,2))**2 ) * pt39 165 a3 = ( a3 + (dxidx(:,1,3) + dxidx(:,2,3) + 166 & dxidx(:,3,3))**2 ) * pt39 167c 168! flops = flops + 15*npro 169 endif 170c 171c.... set up the 1st level Eigenvectors = R_t*Tau^* 172c 173 fact = (one / sqrt( two * rho * T )) / c 174c 175 Q(:,1,5) = fact * ((c + u) * c - eb1 * gamb) 176 Q(:,2,5) = -fact * (c + u * gamb) * Rcs(:,1,1) 177 Q(:,3,5) = -fact * (c + u * gamb) * Rcs(:,1,2) 178 Q(:,4,5) = -fact * (c + u * gamb) * Rcs(:,1,3) 179 Q(:,5,5) = fact * gamb 180c 181 Q(:,1,1) = fact * ((c - u) * c - eb1 * gamb) 182 Q(:,2,1) = fact * (c - u * gamb) * Rcs(:,1,1) 183 Q(:,3,1) = fact * (c - u * gamb) * Rcs(:,1,2) 184 Q(:,4,1) = fact * (c - u * gamb) * Rcs(:,1,3) 185 Q(:,5,1) = fact * gamb 186c 187 Q(:,1,3) = zero 188 Q(:,2,3) = -fact * c * sqt2 * Rcs(:,2,1) ! + in original Jo/Sh code 189 Q(:,3,3) = -fact * c * sqt2 * Rcs(:,2,2) ! could be error here, 190 Q(:,4,3) = -fact * c * sqt2 * Rcs(:,2,3) ! but unlikely 191 Q(:,5,3) = zero 192c 193 Q(:,1,2) = zero 194 Q(:,2,2) = fact * c * sqt2 * Rcs(:,3,1) 195 Q(:,3,2) = fact * c * sqt2 * Rcs(:,3,2) 196 Q(:,4,2) = fact * c * sqt2 * Rcs(:,3,3) 197 Q(:,5,2) = zero 198c 199 fact = (one / sqrt( cp * rho )) / T 200c 201 Q(:,1,4) = fact * eb1 202 Q(:,2,4) = fact * u * Rcs(:,1,1) 203 Q(:,3,4) = fact * u * Rcs(:,1,2) 204 Q(:,4,4) = fact * u * Rcs(:,1,3) 205 Q(:,5,4) = -fact 206c 207c.... flop count 208c 209 210c do i = 1, nflow 211c do j = 1, nflow 212c Q(:,i,j) = abs(Q(:,i,j)) !make sure eigenv. are positive 213c enddo 214c enddo 215 216 ! flops = flops + 203*npro 217c 218c.... return 219c 220 return 221 end 222 223 224 subroutine e3eig2 (u, c, AR1, AR2, AR3, 225 & rlam, Q, eigmax) 226c 227c---------------------------------------------------------------------- 228c 229c This routine diagonalizes a partially reduced matrix using the 230c Jacobi transformation. This routine assumes that the original 231c 5x5 matrix is already reduced to 2x2. 232c 233c input: 234c u (npro) : fluid velocity 235c c (npro) : speed of sound 236c AR1 (npro) : aspect ratio factor in streamline direction 237c AR2 (npro) : aspect ratio factor in normal_1 direction 238c AR3 (npro) : aspect ratio factor in normal_2 direction 239c Q (npro,nflow,nflow) : 1st level eigenvectors of Tau matrix 240c 241c output: 242c rlam (npro,nflow) : eigenvalues 243c Q (npro,nflow,nflow) : eigenvectors 244c 245c 246c T 247c This routine finds S that S Atau S = Lam (Lam is Symm.) 248c 249c Then returns rlam <- Lam and Q <- Q S 250c 251c 252c Note: Jacobi transformation is extracted (and modified) from 253c "Numerical Recipes: The Art of Scientific Computing" by 254c Press, Flannery, Teukolsky and Vetterling, pp. 342-349. 255c 256c 257c Farzin Shakib, Spring 1989. 258c Zdenek Johan, Winter 1991. (Fortran 90) 259c---------------------------------------------------------------------- 260c 261 include "common.h" 262c 263 dimension u(npro), c(npro), 264 & AR1(npro), AR2(npro), 265 & AR3(npro), rlam(npro,nflow), 266 & Q(npro,nflow,nflow), eigmax(npro,nflow) 267c 268 dimension offd(npro), t(npro), 269 & Rcs(npro), Rsn(npro) 270c 271c.... set the reduced eigensystem 272c 273 offd = (AR2 + AR3) * pt5 * c**2 274c 275 rlam(:,1) = AR1 * (u + c)**2 + offd 276 rlam(:,2) = AR1 * u**2 + AR3 * c**2 277 rlam(:,3) = AR1 * u**2 + AR2 * c**2 278 rlam(:,4) = AR1 * u**2 279 rlam(:,5) = AR1 * (u - c)**2 + offd 280c 281c.... modify for time dependent problems 282c 283! consider time term if iremoveStabTimeTerm is set to zero 284 if(iremoveStabTimeTerm.eq.0) then 285 tmp = dtsfct * four * (Dtgl * Dtgl) 286 rlam(:,:) = rlam(:,:) + tmp 287 endif 288c 289c.... compute the rotation tangent ( IEEE arithmetic if offd=0 ) 290c 291 t = pt5 * (rlam(:,1) - rlam(:,5)) / offd 292 t = sign(one, t) / ( abs(t) + sqrt(one + t**2) ) 293c 294c.... compute cosine and sin, and rotate the eigenvalues 295c 296 Rcs = one / sqrt( one + t**2 ) 297 Rsn = t * Rcs 298c 299 rlam(:,1) = rlam(:,1) - t * offd 300 rlam(:,5) = rlam(:,5) + t * offd 301c 302c.... transform the Eigenvectors (all 5 components) 303c 304 t = Rcs * Q(:,1,1) - Rsn * Q(:,1,5) 305 Q(:,1,5) = Rsn * Q(:,1,1) + Rcs * Q(:,1,5) 306 Q(:,1,1) = t 307c 308 t = Rcs * Q(:,2,1) - Rsn * Q(:,2,5) 309 Q(:,2,5) = Rsn * Q(:,2,1) + Rcs * Q(:,2,5) 310 Q(:,2,1) = t 311c 312 t = Rcs * Q(:,3,1) - Rsn * Q(:,3,5) 313 Q(:,3,5) = Rsn * Q(:,3,1) + Rcs * Q(:,3,5) 314 Q(:,3,1) = t 315c 316 t = Rcs * Q(:,4,1) - Rsn * Q(:,4,5) 317 Q(:,4,5) = Rsn * Q(:,4,1) + Rcs * Q(:,4,5) 318 Q(:,4,1) = t 319c 320 t = Rcs * Q(:,5,1) - Rsn * Q(:,5,5) 321 Q(:,5,5) = Rsn * Q(:,5,1) + Rcs * Q(:,5,5) 322 Q(:,5,1) = t 323 324c 325c.... extract maximum eigenvalues 326c 327 eigmax(:,1) = max(rlam(:,1), rlam(:,2), rlam(:,3), 328 & rlam(:,4), rlam(:,5)) 329 eigmax(:,2) = eigmax(:,1) 330 eigmax(:,3) = eigmax(:,1) 331 eigmax(:,4) = eigmax(:,1) 332 eigmax(:,5) = eigmax(:,1) 333 334c 335c.... flop count 336c 337 ! flops = flops + 85*npro 338c 339c.... return 340c 341 return 342 end 343