xref: /phasta/phSolver/compressible/e3eig1.f (revision 595995161822a203c8467e0e4a253d7bd7d6df32)

        subroutine e3eig1 (rho,    T,      cp,     gamb,   c,
     &                     u1,     u2,     u3,     a1,     a2,
     &                     a3,     eb1,
     &                     dxidx,  u,      Q)
c
c----------------------------------------------------------------------
c
c This routine performs the first step of the eigenvalue decomposition
c of the Tau matrix.
c
c input:
c  rho    (npro)           : density
c  T      (npro)           : temperature
c  cp     (npro)           : specific heat at constant pressure
c  gamb   (npro)           : gamma_bar (defined in paper by Chalot et al.)
c  c      (npro)           : speed of sound
c  u1     (npro)           : x1-velocity component
c  u2     (npro)           : x2-velocity component
c  u3     (npro)           : x3-velocity component
c  a1     (npro)           : x1-acceleration component
c  a2     (npro)           : x2-acceleration component
c  a3     (npro)           : x3-acceleration component
c  eb1    (npro)           : e1_bar (defined in paper by Chalot et al.)
c  dxidx  (npro,nsd,nsd)   : inverse of deformation gradient
c
c output:
c  u      (npro)           : fluid velocity
c  a1     (npro)           : aspect ratio factor in streamline direction
c  a2     (npro)           : aspect ratio factor in normal_1 direction
c  a3     (npro)           : aspect ratio factor in normal_2 direction
c  Q      (npro,nflow,nflow) : 1st level eigenvectors of Tau matrix
c
c
c Zdenek Johan, Summer 1990.  (Modified from e2tau.f)
c Zdenek Johan, Winter 1991.  (Fortran 90)
c----------------------------------------------------------------------
c
        include "common.h"
c
        dimension rho(npro),                 T(npro),
     &            cp(npro),                  gamb(npro),
     &            c(npro),                   u1(npro),
     &            u2(npro),                  u3(npro),
     &            a1(npro),                  a2(npro),
     &            a3(npro),
     &            eb1(npro),                 dxidx(npro,nsd,nsd),
     &            u(npro),                   Q(npro,nflow,nflow)
c
        dimension Rcs(npro,nsd,nsd),         fact(npro),
     &            temp(npro)
c
c.... compute the directional cosines (streamline direction)
c
c

        where (u .ne. zero)
           fact       = one / u
           Rcs(:,1,1) = u1 * fact
           Rcs(:,1,2) = u2 * fact
           Rcs(:,1,3) = u3 * fact
        elsewhere
           Rcs(:,1,1) = one
           Rcs(:,1,2) = zero
           Rcs(:,1,3) = zero
        endwhere
c
c.... compute the directional cosines (normal acceleration direction)
c


        fact = a1 * Rcs(:,1,1) + a2 * Rcs(:,1,2) + a3 * Rcs(:,1,3)
        a1   = a1 - fact * Rcs(:,1,1)
        a2   = a2 - fact * Rcs(:,1,2)
        a3   = a3 - fact * Rcs(:,1,3)
        fact = a1**2 + a2**2 + a3**2
c
        where (fact .gt. epsM)
c
          Rcs(:,2,1) = a1
          Rcs(:,2,2) = a2
          Rcs(:,2,3) = a3
c
        elsewhere
c
          Rcs(:,2,1) =   Rcs(:,1,2)
     &                 + sign(one, Rcs(:,1,2)*Rcs(:,1,3)) * Rcs(:,1,3)
          Rcs(:,2,2) = - Rcs(:,1,1)
     &                 - sign(one, Rcs(:,1,2)*Rcs(:,1,3)) * Rcs(:,1,3)
          Rcs(:,2,3) =   sign(one, Rcs(:,1,2)*Rcs(:,1,3)) *
     &                 ( Rcs(:,1,2) - Rcs(:,1,1) )
c
          fact = Rcs(:,2,1)**2 + Rcs(:,2,2)**2 + Rcs(:,2,3)**2
c
        endwhere
c
	fact = one / sqrt(fact)
c
        Rcs(:,2,1) = Rcs(:,2,1) * fact
        Rcs(:,2,2) = Rcs(:,2,2) * fact
        Rcs(:,2,3) = Rcs(:,2,3) * fact
c
c.... compute the directional cosines (last direction)
c
        Rcs(:,3,1) = Rcs(:,1,2) * Rcs(:,2,3) - Rcs(:,1,3) * Rcs(:,2,2)
        Rcs(:,3,2) = Rcs(:,1,3) * Rcs(:,2,1) - Rcs(:,1,1) * Rcs(:,2,3)
        Rcs(:,3,3) = Rcs(:,1,1) * Rcs(:,2,2) - Rcs(:,1,2) * Rcs(:,2,1)

c
c.... calculate the element aspect ratio factors
c
        a1 = Rcs(:,1,1) * dxidx(:,1,1) + Rcs(:,1,2) * dxidx(:,1,2) +
     &       Rcs(:,1,3) * dxidx(:,1,3)
        a2 = Rcs(:,2,1) * dxidx(:,1,1) + Rcs(:,2,2) * dxidx(:,1,2) +
     &       Rcs(:,2,3) * dxidx(:,1,3)
        a3 = Rcs(:,3,1) * dxidx(:,1,1) + Rcs(:,3,2) * dxidx(:,1,2) +
     &       Rcs(:,3,3) * dxidx(:,1,3)
        dxidx(:,1,1) = a1
        dxidx(:,1,2) = a2
        dxidx(:,1,3) = a3
c
        a1 = Rcs(:,1,1) * dxidx(:,2,1) + Rcs(:,1,2) * dxidx(:,2,2) +
     &       Rcs(:,1,3) * dxidx(:,2,3)
        a2 = Rcs(:,2,1) * dxidx(:,2,1) + Rcs(:,2,2) * dxidx(:,2,2) +
     &       Rcs(:,2,3) * dxidx(:,2,3)
        a3 = Rcs(:,3,1) * dxidx(:,2,1) + Rcs(:,3,2) * dxidx(:,2,2) +
     &       Rcs(:,3,3) * dxidx(:,2,3)
        dxidx(:,2,1) = a1
        dxidx(:,2,2) = a2
        dxidx(:,2,3) = a3
c
        a1 = Rcs(:,1,1) * dxidx(:,3,1) + Rcs(:,1,2) * dxidx(:,3,2) +
     &       Rcs(:,1,3) * dxidx(:,3,3)
        a2 = Rcs(:,2,1) * dxidx(:,3,1) + Rcs(:,2,2) * dxidx(:,3,2) +
     &       Rcs(:,2,3) * dxidx(:,3,3)
        a3 = Rcs(:,3,1) * dxidx(:,3,1) + Rcs(:,3,2) * dxidx(:,3,2) +
     &       Rcs(:,3,3) * dxidx(:,3,3)
        dxidx(:,3,1) = a1
        dxidx(:,3,2) = a2
        dxidx(:,3,3) = a3

c
c... original
c

        a1 = dxidx(:,1,1)**2 + dxidx(:,2,1)**2 + dxidx(:,3,1)**2
        a2 = dxidx(:,1,2)**2 + dxidx(:,2,2)**2 + dxidx(:,3,2)**2
        a3 = dxidx(:,1,3)**2 + dxidx(:,2,3)**2 + dxidx(:,3,3)**2

c
c... change from original (analyt., could be error)
c

cc        a1 = dxidx(:,1,1)**2 + dxidx(:,1,2)**2 + dxidx(:,1,3)**2
cc        a2 = dxidx(:,2,1)**2 + dxidx(:,2,2)**2 + dxidx(:,2,3)**2
cc        a3 = dxidx(:,3,1)**2 + dxidx(:,3,2)**2 + dxidx(:,3,3)**2

c
c.... correct for tetrahedra
c
        if (lcsyst .eq. 1) then
          a1 = ( a1 + (dxidx(:,1,1) + dxidx(:,2,1) +
     &                 dxidx(:,3,1))**2 ) * pt39
          a2 = ( a2 + (dxidx(:,1,2) + dxidx(:,2,2) +
     &                 dxidx(:,3,2))**2 ) * pt39
          a3 = ( a3 + (dxidx(:,1,3) + dxidx(:,2,3) +
     &                 dxidx(:,3,3))**2 ) * pt39
c
          flops = flops + 15*npro
        endif
c
c.... set up the 1st level Eigenvectors =  R_t*Tau^*
c
        fact     =  (one / sqrt( two * rho * T )) / c
c
        Q(:,1,5) =  fact * ((c + u) * c - eb1 * gamb)
        Q(:,2,5) = -fact * (c + u * gamb) * Rcs(:,1,1)
        Q(:,3,5) = -fact * (c + u * gamb) * Rcs(:,1,2)
        Q(:,4,5) = -fact * (c + u * gamb) * Rcs(:,1,3)
        Q(:,5,5) =  fact * gamb
c
        Q(:,1,1) =  fact * ((c - u) * c - eb1 * gamb)
        Q(:,2,1) =  fact * (c - u * gamb) * Rcs(:,1,1)
        Q(:,3,1) =  fact * (c - u * gamb) * Rcs(:,1,2)
        Q(:,4,1) =  fact * (c - u * gamb) * Rcs(:,1,3)
        Q(:,5,1) =  fact * gamb
c
        Q(:,1,3) =  zero
        Q(:,2,3) = -fact * c * sqt2 * Rcs(:,2,1) ! + in original Jo/Sh code
        Q(:,3,3) = -fact * c * sqt2 * Rcs(:,2,2) ! could be error here,
        Q(:,4,3) = -fact * c * sqt2 * Rcs(:,2,3) ! but unlikely
        Q(:,5,3) =  zero
c
        Q(:,1,2) =  zero
        Q(:,2,2) =  fact * c * sqt2 * Rcs(:,3,1)
        Q(:,3,2) =  fact * c * sqt2 * Rcs(:,3,2)
        Q(:,4,2) =  fact * c * sqt2 * Rcs(:,3,3)
        Q(:,5,2) =  zero
c
        fact     =  (one / sqrt( cp * rho )) / T
c
        Q(:,1,4) =  fact * eb1
        Q(:,2,4) =  fact * u * Rcs(:,1,1)
        Q(:,3,4) =  fact * u * Rcs(:,1,2)
        Q(:,4,4) =  fact * u * Rcs(:,1,3)
        Q(:,5,4) = -fact
c
c.... flop count
c

c        do i = 1, nflow
c           do j = 1, nflow
c              Q(:,i,j) = abs(Q(:,i,j)) !make sure eigenv. are positive
c           enddo
c        enddo

        flops = flops + 203*npro
c
c.... return
c
        return
        end


        subroutine e3eig2 (u,      c,      AR1,    AR2,    AR3,
     &                     rlam,   Q,    eigmax)
c
c----------------------------------------------------------------------
c
c  This routine diagonalizes a partially reduced matrix using the
c Jacobi transformation.  This routine assumes that the original
c 5x5 matrix is already reduced to 2x2.
c
c input:
c  u      (npro)           : fluid velocity
c  c      (npro)           : speed of sound
c  AR1    (npro)           : aspect ratio factor in streamline direction
c  AR2    (npro)           : aspect ratio factor in normal_1 direction
c  AR3    (npro)           : aspect ratio factor in normal_2 direction
c  Q      (npro,nflow,nflow) : 1st level eigenvectors of Tau matrix
c
c output:
c  rlam   (npro,nflow)      : eigenvalues
c  Q      (npro,nflow,nflow) : eigenvectors
c
c
c                                  T
c  This routine finds S that      S  Atau S = Lam  (Lam is Symm.)
c
c  Then returns       rlam <- Lam     and      Q <- Q S
c
c
c Note: Jacobi transformation is extracted (and modified) from
c       "Numerical Recipes: The Art of Scientific Computing" by
c       Press, Flannery, Teukolsky and Vetterling, pp. 342-349.
c
c
c Farzin Shakib, Spring 1989.
c Zdenek Johan,  Winter 1991.  (Fortran 90)
c----------------------------------------------------------------------
c
        include "common.h"
c
        dimension u(npro),                   c(npro),
     &            AR1(npro),                 AR2(npro),
     &            AR3(npro),                 rlam(npro,nflow),
     &            Q(npro,nflow,nflow),         eigmax(npro,nflow)
c
        dimension offd(npro),                t(npro),
     &            Rcs(npro),                 Rsn(npro)
c
c.... set the reduced eigensystem
c
        offd      = (AR2 + AR3) * pt5 * c**2
c
        rlam(:,1) = AR1 * (u + c)**2 + offd
        rlam(:,2) = AR1 * u**2       + AR3 * c**2
        rlam(:,3) = AR1 * u**2       + AR2 * c**2
        rlam(:,4) = AR1 * u**2
        rlam(:,5) = AR1 * (u - c)**2 + offd
c
c.... modify for time dependent problems
c
        ! consider time term if iremoveStabTimeTerm is set to zero
        if(iremoveStabTimeTerm.eq.0) then
           tmp  = dtsfct * four * (Dtgl * Dtgl)
           rlam(:,:) = rlam(:,:) + tmp
        endif
c
c.... compute the rotation tangent ( IEEE arithmetic if offd=0 )
c
        t = pt5 * (rlam(:,1) - rlam(:,5)) / offd
        t = sign(one, t) / ( abs(t) + sqrt(one + t**2) )
c
c.... compute cosine and sin, and rotate the eigenvalues
c
        Rcs = one / sqrt( one + t**2 )
        Rsn = t * Rcs
c
        rlam(:,1) = rlam(:,1) - t * offd
        rlam(:,5) = rlam(:,5) + t * offd
c
c.... transform the Eigenvectors (all 5 components)
c
        t        = Rcs * Q(:,1,1) - Rsn * Q(:,1,5)
        Q(:,1,5) = Rsn * Q(:,1,1) + Rcs * Q(:,1,5)
        Q(:,1,1) = t
c
        t        = Rcs * Q(:,2,1) - Rsn * Q(:,2,5)
        Q(:,2,5) = Rsn * Q(:,2,1) + Rcs * Q(:,2,5)
        Q(:,2,1) = t
c
        t        = Rcs * Q(:,3,1) - Rsn * Q(:,3,5)
        Q(:,3,5) = Rsn * Q(:,3,1) + Rcs * Q(:,3,5)
        Q(:,3,1) = t
c
        t        = Rcs * Q(:,4,1) - Rsn * Q(:,4,5)
        Q(:,4,5) = Rsn * Q(:,4,1) + Rcs * Q(:,4,5)
        Q(:,4,1) = t
c
        t        = Rcs * Q(:,5,1) - Rsn * Q(:,5,5)
        Q(:,5,5) = Rsn * Q(:,5,1) + Rcs * Q(:,5,5)
        Q(:,5,1) = t

c
c.... extract maximum eigenvalues
c
        eigmax(:,1) = max(rlam(:,1), rlam(:,2), rlam(:,3),
     &                    rlam(:,4), rlam(:,5))
        eigmax(:,2) = eigmax(:,1)
        eigmax(:,3) = eigmax(:,1)
        eigmax(:,4) = eigmax(:,1)
        eigmax(:,5) = eigmax(:,1)

c
c.... flop count
c
        flops = flops + 85*npro
c
c.... return
c
        return
        end