xref: /phasta/phSolver/compressible/e3conv.f (revision 513954ef803c300cddba2bb96b4a5dac0b93489a)

        subroutine e3conv (g1yi,   g2yi,   g3yi,
     &                     A1,     A2,     A3,
     &                     rho,    pres,   T,
     &                     ei,     rk,     u1,
     &                     u2,     u3,     rLyi,
     &                     ri,     rmi,    EGmass,
     &                     shg,    shape,  WdetJ )
!
!----------------------------------------------------------------------
!
! This routine calculates the contribution of Galerkin part of the
! Convective term (Time and Euler fluxes) to both RHS and LHS.
!
! input:
!  g1yi   (npro,nflow)    : grad-y in direction 1
!  g2yi   (npro,nflow)    : grad-y in direction 2
!  g3yi   (npro,nflow)    : grad-y in direction 3
!  A1    (npro,nflow,nflow)  : A-1
!  A2    (npro,nflow,nflow)  : A-2
!  A3    (npro,nflow,nflow)  : A-3
!  rho    (npro)         : density
!  pres   (npro)         : pressure
!  T      (npro)         : temperature
!  ei     (npro)         : internal energy
!  rk     (npro)         : kinetic energy
!  u1     (npro)         : x1-velocity component
!  u2     (npro)         : x2-velocity component
!  u3     (npro)         : x3-velocity component
!  shg    (npro,nshl,nsd) : global grad's of shape functions
!  shape  (npro,nshl)  : element shape functions
!  WdetJ  (npro)         : weighted Jacobian determinant
!
! output:
!  rLyi   (npro,nflow)           : least-squares residual vector
!  ri     (npro,nflow*(nsd+1))   : partial residual
!  rmi    (npro,nflow*(nsd+1))   : partial modified residual
!  EGmass (npro,nedof,nedof)    : partial LHS tangent matrix
!
!
! Zdenek Johan, Summer 1990. (Modified from e2conv.f)
! Zdenek Johan, Winter 1991. (Fortran 90)
! Kenneth Jansen, Winter 1997 Primitive Variables
!----------------------------------------------------------------------
!
        include "common.h"
!
!  passed arrays
!
        dimension g1yi(npro,nflow),           g2yi(npro,nflow),
     &            g3yi(npro,nflow),
     &            A1(npro,nflow,nflow),
     &            A2(npro,nflow,nflow),       A3(npro,nflow,nflow),
     &            rho(npro),                pres(npro),
     &            T(npro),                  ei(npro),
     &            rk(npro),                 u1(npro),
     &            u2(npro),                 u3(npro),
     &            rLyi(npro,nflow),          ri(npro,nflow*(nsd+1)),
     &            rmi(npro,nflow*(nsd+1)),   EGmass(npro,nedof,nedof),
     &            shg(npro,nshl,nsd),       shape(npro,nshl),
     &            WdetJ(npro)
!
!  local arrays
!
        dimension AiNbi(npro,nflow,nflow),     fact1(npro),
     &            fact2(npro),               fact3(npro)

        ttim(22) = ttim(22) - secs(0.0)

!
!.... ---------------------->  RHS, Euler Flux  <----------------------
!
        if ((ires .eq. 1) .or. (ires .eq. 3)) then
!
!.... calculate integrated by part contribution of Euler flux (Galerkin)
!
          ri(:, 1) = (- u1) * rho
          ri(:, 2) = (- u1) * rho * u1 - pres
          ri(:, 3) = (- u1) * rho * u2
          ri(:, 4) = (- u1) * rho * u3
          ri(:, 5) = (- u1) * rho * (ei + rk) - u1 * pres
!
          ri(:, 6) = (- u2) * rho
          ri(:, 7) = (- u2) * rho * u1
          ri(:, 8) = (- u2) * rho * u2 - pres
          ri(:, 9) = (- u2) * rho * u3
          ri(:,10) = (- u2) * rho * (ei + rk) - u2 * pres
!
          ri(:,11) = (- u3) * rho
          ri(:,12) = (- u3) * rho * u1
          ri(:,13) = (- u3) * rho * u2
          ri(:,14) = (- u3) * rho * u3 - pres
          ri(:,15) = (- u3) * rho * (ei + rk) - u3 * pres
!
          flops = flops + 28*npro
!
        endif
!
!.... calculate ( A_i Y,i ) --> rLyi   Commented out zeros of A matrices
!
        rLyi(:,1) =
     &              A1(:,1,1) * g1yi(:,1)
     &            + A1(:,1,2) * g1yi(:,2)
!    &            + A1(:,1,3) * g1yi(:,3)
!    &            + A1(:,1,4) * g1yi(:,4)
     &            + A1(:,1,5) * g1yi(:,5)
     &            + A2(:,1,1) * g2yi(:,1)
!    &            + A2(:,1,2) * g2yi(:,2)
     &            + A2(:,1,3) * g2yi(:,3)
!    &            + A2(:,1,4) * g2yi(:,4)
     &            + A2(:,1,5) * g2yi(:,5)
     &            + A3(:,1,1) * g3yi(:,1)
!    &            + A3(:,1,2) * g3yi(:,2)
!    &            + A3(:,1,3) * g3yi(:,3)
     &            + A3(:,1,4) * g3yi(:,4)
     &            + A3(:,1,5) * g3yi(:,5)
        rLyi(:,2) =
     &              A1(:,2,1) * g1yi(:,1)
     &            + A1(:,2,2) * g1yi(:,2)
!    &            + A1(:,2,3) * g1yi(:,3)
!    &            + A1(:,2,4) * g1yi(:,4)
     &            + A1(:,2,5) * g1yi(:,5)
     &            + A2(:,2,1) * g2yi(:,1)
     &            + A2(:,2,2) * g2yi(:,2)
     &            + A2(:,2,3) * g2yi(:,3)
!    &            + A2(:,2,4) * g2yi(:,4)
     &            + A2(:,2,5) * g2yi(:,5)
     &            + A3(:,2,1) * g3yi(:,1)
     &            + A3(:,2,2) * g3yi(:,2)
!    &            + A3(:,2,3) * g3yi(:,3)
     &            + A3(:,2,4) * g3yi(:,4)
     &            + A3(:,2,5) * g3yi(:,5)
        rLyi(:,3) =
     &              A1(:,3,1) * g1yi(:,1)
     &            + A1(:,3,2) * g1yi(:,2)
     &            + A1(:,3,3) * g1yi(:,3)
!    &            + A1(:,3,4) * g1yi(:,4)
     &            + A1(:,3,5) * g1yi(:,5)
     &            + A2(:,3,1) * g2yi(:,1)
!    &            + A2(:,3,2) * g2yi(:,2)
     &            + A2(:,3,3) * g2yi(:,3)
!    &            + A2(:,3,4) * g2yi(:,4)
     &            + A2(:,3,5) * g2yi(:,5)
     &            + A3(:,3,1) * g3yi(:,1)
!    &            + A3(:,3,2) * g3yi(:,2)
     &            + A3(:,3,3) * g3yi(:,3)
     &            + A3(:,3,4) * g3yi(:,4)
     &            + A3(:,3,5) * g3yi(:,5)
        rLyi(:,4) =
     &              A1(:,4,1) * g1yi(:,1)
     &            + A1(:,4,2) * g1yi(:,2)
!    &            + A1(:,4,3) * g1yi(:,3)
     &            + A1(:,4,4) * g1yi(:,4)
     &            + A1(:,4,5) * g1yi(:,5)
     &            + A2(:,4,1) * g2yi(:,1)
!    &            + A2(:,4,2) * g2yi(:,2)
     &            + A2(:,4,3) * g2yi(:,3)
     &            + A2(:,4,4) * g2yi(:,4)
     &            + A2(:,4,5) * g2yi(:,5)
     &            + A3(:,4,1) * g3yi(:,1)
!    &            + A3(:,4,2) * g3yi(:,2)
!    &            + A3(:,4,3) * g3yi(:,3)
     &            + A3(:,4,4) * g3yi(:,4)
     &            + A3(:,4,5) * g3yi(:,5)
        rLyi(:,5) =
     &              A1(:,5,1) * g1yi(:,1)
     &            + A1(:,5,2) * g1yi(:,2)
     &            + A1(:,5,3) * g1yi(:,3)
     &            + A1(:,5,4) * g1yi(:,4)
     &            + A1(:,5,5) * g1yi(:,5)
     &            + A2(:,5,1) * g2yi(:,1)
     &            + A2(:,5,2) * g2yi(:,2)
     &            + A2(:,5,3) * g2yi(:,3)
     &            + A2(:,5,4) * g2yi(:,4)
     &            + A2(:,5,5) * g2yi(:,5)
     &            + A3(:,5,1) * g3yi(:,1)
     &            + A3(:,5,2) * g3yi(:,2)
     &            + A3(:,5,3) * g3yi(:,3)
     &            + A3(:,5,4) * g3yi(:,4)
     &            + A3(:,5,5) * g3yi(:,5)
!
!.... add contribution to rmi
!
        if ((ires .eq. 2) .or. (ires .eq. 3))
     &    rmi(:,16:20) = rLyi  ! modified residual uses non i.b.p form of conv.
!
!.... ---------------------->  LHS   <-----------------------
!
        if (lhs .eq. 1) then
!
!.... loop through the columns
!
        do j = 1, nshl
           j0 = nflow * (j - 1)
!
!.... compute some useful factors
!
           fact1 = WdetJ * shg(:,j,1)
           fact2 = WdetJ * shg(:,j,2)
           fact3 = WdetJ * shg(:,j,3)
!
!.... first compute (A_i N_b,i)
!
           AiNbi(:,1,1) =
     &                    fact1 * A1(:,1,1)
     &                  + fact2 * A2(:,1,1)
     &                  + fact3 * A3(:,1,1)
           AiNbi(:,1,2) =
     &                    fact1 * A1(:,1,2)
     &                  + fact2 * A2(:,1,2)
     &                  + fact3 * A3(:,1,2)
           AiNbi(:,1,3) =
     &                    fact1 * A1(:,1,3)
     &                  + fact2 * A2(:,1,3)
     &                  + fact3 * A3(:,1,3)
           AiNbi(:,1,4) =
     &                    fact1 * A1(:,1,4)
     &                  + fact2 * A2(:,1,4)
     &                  + fact3 * A3(:,1,4)
           AiNbi(:,1,5) =
     &                    fact1 * A1(:,1,5)
     &                  + fact2 * A2(:,1,5)
     &                  + fact3 * A3(:,1,5)

           AiNbi(:,2,1) =
     &                    fact1 * A1(:,2,1)
     &                  + fact2 * A2(:,2,1)
     &                  + fact3 * A3(:,2,1)
           AiNbi(:,2,2) =
     &                    fact1 * A1(:,2,2)
     &                  + fact2 * A2(:,2,2)
     &                  + fact3 * A3(:,2,2)
           AiNbi(:,2,3) =
     &                    fact1 * A1(:,2,3)
     &                  + fact2 * A2(:,2,3)
     &                  + fact3 * A3(:,2,3)
           AiNbi(:,2,4) =
     &                    fact1 * A1(:,2,4)
     &                  + fact2 * A2(:,2,4)
     &                  + fact3 * A3(:,2,4)
           AiNbi(:,2,5) =
     &                    fact1 * A1(:,2,5)
     &                  + fact2 * A2(:,2,5)
     &                  + fact3 * A3(:,2,5)

           AiNbi(:,3,1) =
     &                    fact1 * A1(:,3,1)
     &                  + fact2 * A2(:,3,1)
     &                  + fact3 * A3(:,3,1)
           AiNbi(:,3,2) =
     &                    fact1 * A1(:,3,2)
     &                  + fact2 * A2(:,3,2)
     &                  + fact3 * A3(:,3,2)
           AiNbi(:,3,3) =
     &                    fact1 * A1(:,3,3)
     &                  + fact2 * A2(:,3,3)
     &                  + fact3 * A3(:,3,3)
           AiNbi(:,3,4) =
     &                    fact1 * A1(:,3,4)
     &                  + fact2 * A2(:,3,4)
     &                  + fact3 * A3(:,3,4)
           AiNbi(:,3,5) =
     &                    fact1 * A1(:,3,5)
     &                  + fact2 * A2(:,3,5)
     &                  + fact3 * A3(:,3,5)

           AiNbi(:,4,1) =
     &                    fact1 * A1(:,4,1)
     &                  + fact2 * A2(:,4,1)
     &                  + fact3 * A3(:,4,1)
           AiNbi(:,4,2) =
     &                    fact1 * A1(:,4,2)
     &                  + fact2 * A2(:,4,2)
     &                  + fact3 * A3(:,4,2)
           AiNbi(:,4,3) =
     &                    fact1 * A1(:,4,3)
     &                  + fact2 * A2(:,4,3)
     &                  + fact3 * A3(:,4,3)
           AiNbi(:,4,4) =
     &                    fact1 * A1(:,4,4)
     &                  + fact2 * A2(:,4,4)
     &                  + fact3 * A3(:,4,4)
           AiNbi(:,4,5) =
     &                    fact1 * A1(:,4,5)
     &                  + fact2 * A2(:,4,5)
     &                  + fact3 * A3(:,4,5)

           AiNbi(:,5,1) =
     &                    fact1 * A1(:,5,1)
     &                  + fact2 * A2(:,5,1)
     &                  + fact3 * A3(:,5,1)
           AiNbi(:,5,2) =
     &                    fact1 * A1(:,5,2)
     &                  + fact2 * A2(:,5,2)
     &                  + fact3 * A3(:,5,2)
           AiNbi(:,5,3) =
     &                    fact1 * A1(:,5,3)
     &                  + fact2 * A2(:,5,3)
     &                  + fact3 * A3(:,5,3)
           AiNbi(:,5,4) =
     &                    fact1 * A1(:,5,4)
     &                  + fact2 * A2(:,5,4)
     &                  + fact3 * A3(:,5,4)
           AiNbi(:,5,5) =
     &                    fact1 * A1(:,5,5)
     &                  + fact2 * A2(:,5,5)
     &                  + fact3 * A3(:,5,5)
!
!.... now loop through the row nodes and add (N_a A_i N_b,i) to
!     the tangent matrix.
!
           do i = 1, nshl
              i0 = nflow * (i - 1)
!
!.... loop through dof's
!
              do jdof = 1, nflow
                 jl = j0 + jdof

                 EGmass(:,i0+1,jl) = EGmass(:,i0+1,jl) +
     &                               shape(:,i) * AiNbi(:,1,jdof)

                 EGmass(:,i0+2,jl) = EGmass(:,i0+2,jl) +
     &                               shape(:,i) * AiNbi(:,2,jdof)

                 EGmass(:,i0+3,jl) = EGmass(:,i0+3,jl) +
     &                               shape(:,i) * AiNbi(:,3,jdof)

                 EGmass(:,i0+4,jl) = EGmass(:,i0+4,jl) +
     &                               shape(:,i) * AiNbi(:,4,jdof)

                 EGmass(:,i0+5,jl) = EGmass(:,i0+5,jl) +
     &                               shape(:,i) * AiNbi(:,5,jdof)
              enddo
!
!.... end loop on rows
!
           enddo
!
!.... end loop on columns
!
        enddo
!
!.... end of LHS tangent matrix computation
!
        endif

        ttim(22) = ttim(22) + secs(0.0)
!
!.... return
!
        return
        end
!
!
!
        subroutine e3convSclr (g1yti,   g2yti,   g3yti,
     &                         A1t,     A2t,     A3t,
     &                         rho,     u1,      Sclr,
     &                         u2,      u3,      rLyti,
     &                         rti,     rmti,    EGmasst,
     &                         shg,     shape,   WdetJ)
!
!----------------------------------------------------------------------
!
! This routine calculates the contribution of Galerkin part of the
! Convective term (Time and Euler fluxes) to both RHS and LHS.
!
! input:
!  Sclr   (npro)          : Scalar variable
!  g1yti  (npro)          : grad-y in direction 1
!  g2yti  (npro)          : grad-y in direction 2
!  g3yti  (npro)          : grad-y in direction 3
!  A1t    (npro)          : A-1
!  A2t    (npro)          : A-2
!  A3t    (npro)          : A-3
!  rho    (npro)          : density
!  u1     (npro)          : x1-velocity component
!  u2     (npro)          : x2-velocity component
!  u3     (npro)          : x3-velocity component
!  shg    (npro,nshl,nsd) : global grad's of shape functions
!  shape  (npro,nshl)     : element shape functions
!  WdetJ  (npro)          : weighted Jacobian determinant
!
! output:
!  rLyti   (npro)         : least-squares residual vector
!  rti     (npro,nsd+1)   : partial residual
!  rmti    (npro,nsd+1)   : partial modified residual
!  EGmasst (npro,nshape,nshape): partial LHS tangent matrix
!
!
! Zdenek Johan, Summer 1990. (Modified from e2conv.f)
! Zdenek Johan, Winter 1991. (Fortran 90)
! Kenneth Jansen, Winter 1997 Primitive Variables
!----------------------------------------------------------------------
!
        include "common.h"
!
!  passed arrays
!
        dimension g1yti(npro),            g2yti(npro),
     &            g3yti(npro),            Sclr(npro),
     &            A1t(npro),
     &            A2t(npro),              A3t(npro),
     &            rho(npro),              u1(npro),
     &            u2(npro),               u3(npro),
     &            rLyti(npro),            rti(npro,nsd+1),
     &            rmti(npro,nsd+1),       EGmasst(npro,nshape,nshape),
     &            shg(npro,nshl,nsd),     shape(npro,nshl),
     &            WdetJ(npro)
!
!  local arrays
!
        dimension AitNbi(npro)

        ttim(22) = ttim(22) - tmr() !tmr(0.0)
!
!.... ---------------------->  RHS, Euler Flux  <----------------------
!
        if ((ires .eq. 1) .or. (ires .eq. 3)) then
!
!.... calculate integrated by part contribution of Euler flux (Galerkin)
!
           if (iconvsclr.eq.2) then ! convective form
!
              rti(:, 4) = rti(:,4) + ( u1) * g1yti(:)
     &                             + ( u2) * g2yti(:)
     &                             + ( u3) * g3yti(:)
!
           else                 ! conservative form
!
              rti(:, 1) = rti(:,1) + (- u1) * rho * Sclr
              rti(:, 2) = rti(:,2) + (- u2) * rho * Sclr
              rti(:, 3) = rti(:,3) + (- u3) * rho * Sclr
!
           endif

           flops = flops + 28*npro

        endif
!
!.... calculate ( A_i Y,i ) --> rLyi
!
        rLyti(:) = rLyti(:)
     &            + A1t(:) * g1yti(:)
     &            + A2t(:) * g2yti(:)
     &            + A3t(:) * g3yti(:)

!
!.... add contribution to rmi
!
!        if ((ires .eq. 2) .or. (ires .eq. 3))
!     &    rmi(:,16:20) = rLyi  ! modified residual uses non i.b.p form of conv
!
!.... ---------------------->  LHS   <-----------------------
!
        if (lhs .eq. 1) then
!
!.... loop through the columns
!
        do j = 1, nshl

!
!.... first compute (A_i N_b,i)
!
           AitNbi(:) =
     &                    WdetJ * shg(:,j,1) * A1t(:)
     &                  + WdetJ * shg(:,j,2) * A2t(:)
     &                  + WdetJ * shg(:,j,3) * A3t(:)

!
!.... now loop through the rows and add (N_a A_i N_b,i) to
!     the tangent matrix.
!
           do i = 1, nshl

             EGmasst(:,i,j) = EGmasst(:,i,j) +  shape(:,i) * AitNbi(:)


!
!.... end loop on rows
!
           enddo
!
!.... end loop on columns
!
        enddo
!
!.... end of LHS tangent matrix computation
!
        endif

        ttim(22) = ttim(22) + tmr()
!
!.... return
!
        return
        end