xref: /phasta/phSolver/compressible/e3ls.f (revision 0d39a63aa7ba54dac51aef7a29581103537187d1)

        subroutine e3LS   (A1,        A2,          A3,
     &                     rho,       rmu,         cp,
     &                     cv,        con,         T,
     &                     u1,        u2,          u3,
     &                     rLyi,      dxidx,       tau,
     &                     ri,        rmi,         rk,
     &                     dui,       aci,         A0,
     &                     divqi,     shape,       shg,
     &                     EGmass,    stiff,       WdetJ,
     &                     giju,      rTLS,        raLS,
     &                     A0inv,     dVdY,        rerrl,
     &                     compK,     pres,        PTau)
c
c----------------------------------------------------------------------
c
c This routine calculates the contribution of the least-squares
c operator to the RHS vector and LHS tangent matrix. The temporary
c results are put in ri.
c
c input:
c  A1    (npro,nflow,nflow)     : A_1
c  A2    (npro,nflow,nflow)     : A_2
c  A3    (npro,nflow,nflow)     : A_3
c  rho   (npro)               : density
c  T     (npro)               : temperature
c  cp    (npro)               : specific heat at constant pressure
c  u1    (npro)               : x1-velocity component
c  u2    (npro)               : x2-velocity component
c  u3    (npro)               : x3-velocity component
c  rLyi  (npro,nflow)          : least-squares residual vector
c  dxidx (npro,nsd,nsd)       : inverse of deformation gradient
c  tau   (npro,3)             : stability parameter
c  PTau  (npro,5,5)           : matrix of stability parameters
c  rLyi  (npro,nflow)          : convective portion of least-squares
c                               residual vector
c  divqi (npro,nflow-1)        : divergence of diffusive flux
c  shape (npro,nshl)        : element shape functions
c  shg   (npro,nshl,nsd)    : global element shape function grads
c  WdetJ (npro)               : weighted jacobian determinant
c  stiff (npro,nsd*nflow,nsd*nflow) : stiffness matrix
c  EGmass(npro,nedof,nedof)   : partial mass matrix
c  compK (npro,10)             : K_ij in (eq.134) A new ... III
c
c output:
c  ri     (npro,nflow*(nsd+1)) : partial residual
c  rmi    (npro,nflow*(nsd+1)) : partial modified residual
c  EGmass (npro,nedof,nedof)  : partial mass matrix
c
c
c Zdenek Johan, Summer 1990. (Modified from e2ls.f)
c Zdenek Johan, Winter 1991. (Fortran 90)
c Kenneth Jansen, Winter 1997. Prim. Var. with Diag Tau
c----------------------------------------------------------------------
c
        include "common.h"

c
c  passed arrays
c
        dimension A1(npro,nflow,nflow),    A2(npro,nflow,nflow),
     &            A3(npro,nflow,nflow),    cv(npro),
     &            A0(npro,nflow,nflow),    rho(npro),
     &            con(npro),               cp(npro),
     &            dui(npro,nflow),         aci(npro,nflow),
     &            u1(npro),                u2(npro),
     &            u3(npro),                rk(npro),
     &            rLyi(npro,nflow),        dxidx(npro,nsd,nsd),
     &            tau(npro,5),             giju(npro,6),
     &            rTLS(npro),              raLS(npro),
     &            ri(npro,nflow*(nsd+1)),  rmi(npro,nflow*(nsd+1)),
     &            divqi(npro,nflow-1),     stiff(npro,3*nflow,3*nflow),
     &            EGmass(npro,nedof,nedof),shape(npro,nshl),
     &            shg(npro,nshl,nsd),      WdetJ(npro),
     &            PTau(npro,5,5),          T(npro),
     &            pres(npro)
c
c local arrays
c
        dimension rLymi(npro,nflow),         Atau(npro,nflow,nflow),
     &            A1tauA0(npro,nflow,nflow), A2tauA0(npro,nflow,nflow),
     &            A3tauA0(npro,nflow,nflow), fact(npro),
     &            A0inv(npro,15),            dVdY(npro,15),
     &            compK(npro,10),            ac1(npro),
     &            ac2(npro),                 ac3(npro)
c
        real*8    rerrl(npro,nshl,6), tmp(npro), tmp1(npro)
        ttim(24) = ttim(24) - secs(0.0)
c
c
c last step to the least squares is adding the time term.  So that we
c only have to localize one vector for each Krylov vector the modified
c residual is quite different from the total residual.
c
c
c the modified residual
c
       fct1=almi/gami/alfi*dtgl
c
c  add possibility of not including time term
c
c       if(idiff.ne.-1)

       if(ires.ne.1) rLymi = rLyi + fct1*dui
c
       if(ires.eq.1 .or. ires .eq. 3) then
c       rLymi = rLyi

        rLyi(:,1) = rLyi(:,1)
     &            + A0(:,1,1)*aci(:,1)
c    &            + A0(:,1,2)*aci(:,2)
c    &            + A0(:,1,3)*aci(:,3)
c    &            + A0(:,1,4)*aci(:,4)
     &            + A0(:,1,5)*aci(:,5)
c
        rLyi(:,2) = rLyi(:,2)
     &            + A0(:,2,1)*aci(:,1)
     &            + A0(:,2,2)*aci(:,2)
c    &            + A0(:,2,3)*aci(:,3)
c    &            + A0(:,2,4)*aci(:,4)
     &            + A0(:,2,5)*aci(:,5)
c
        rLyi(:,3) = rLyi(:,3)
     &            + A0(:,3,1)*aci(:,1)
c    &            + A0(:,3,2)*aci(:,2)
     &            + A0(:,3,3)*aci(:,3)
c    &            + A0(:,3,4)*aci(:,4)
     &            + A0(:,3,5)*aci(:,5)
c
        rLyi(:,4) = rLyi(:,4)
     &            + A0(:,4,1)*aci(:,1)
c    &            + A0(:,4,2)*aci(:,2)
c    &            + A0(:,4,3)*aci(:,3)
     &            + A0(:,4,4)*aci(:,4)
     &            + A0(:,4,5)*aci(:,5)
c
        rLyi(:,5) = rLyi(:,5)
     &            + A0(:,5,1)*aci(:,1)
     &            + A0(:,5,2)*aci(:,2)
     &            + A0(:,5,3)*aci(:,3)
     &            + A0(:,5,4)*aci(:,4)
     &            + A0(:,5,5)*aci(:,5)
c
      endif
c
c.... subtract div(q) from the least squares term
c
      if ((idiff >= 1).and.(ires==3 .or. ires==1)) then
c
      if (isurf.eq.zero) then
         rLyi(:,2) = rLyi(:,2) - divqi(:,1)
         rLyi(:,3) = rLyi(:,3) - divqi(:,2)
         rLyi(:,4) = rLyi(:,4) - divqi(:,3)
         rLyi(:,5) = rLyi(:,5) - divqi(:,4)
      endif
      endif
c
c.... -------------------> error calculation  <-----------------
c
       if((ierrcalc.eq.1).and.(nitr.eq.iter)) then
          do ia=1,nshl
             tmp=shape(:,ia)*WdetJ(:)
             tmp1=shape(:,ia)*Qwt(lcsyst,intp)
             rerrl(:,ia,1) = rerrl(:,ia,1) +
     &                       tmp1(:)*rLyi(:,1)*rLyi(:,1)
             rerrl(:,ia,2) = rerrl(:,ia,2) +
     &                       tmp1(:)*rLyi(:,2)*rLyi(:,2)
             rerrl(:,ia,3) = rerrl(:,ia,3) +
     &                       tmp1(:)*rLyi(:,3)*rLyi(:,3)

             rerrl(:,ia,4) = rerrl(:,ia,4) +
     &                       tmp(:)*divqi(:,1)*divqi(:,1)
             rerrl(:,ia,5) = rerrl(:,ia,5) +
     &                       tmp(:)*divqi(:,2)*divqi(:,2)
! SAM wants a threshold here so we are going to take over this little used
! error indictor for that purpose
! commented for ShockError             rerrl(:,ia,6) = rerrl(:,ia,6) +
! commented for ShockError     &                       tmp(:)*divqi(:,3)*divqi(:,3)
          enddo
       endif
c
c
c.... --------------------------->  Tau  <-----------------------------
c
c.... calculate the tau matrix
c
c.... in the first incarnation we will ONLY have a diagonal tau here

       if (itau .lt. 10) then    ! diagonal tau

          call e3tau  (rho,             cp,		rmu,
     &         u1,              u2,             u3,
     &         con,             dxidx,          rLyi,
     &         rLymi,           tau,            rk,
     &         giju,            rTLS,           raLS,
     &         A0inv,           dVdY,           cv)

       else

c.... looks like need a non-diagonal, discribed in "A new ... III"
c.... by Hughes and Mallet. Original work - developed diffusive (and as
c.... well advective) correction factors for 1-D a-d equation w/ hier. b.


          ac1(:) = aci(:,2)
          ac2(:) = aci(:,3)
          ac3(:) = aci(:,4)

          call e3tau_nd  (rho,       cp,  rmu,   T,
     &         u1,              u2,             u3,
     &         ac1,             ac2,             ac3,
     &         con,             dxidx,          rLyi,
     &         rLymi,           PTau,           rk,
     &         giju,            rTLS,           raLS,
     &         cv,              compK,          pres,
     &         A0inv,           dVdY)

       endif


        ttim(25) = ttim(25) + secs(0.0)
c
c Warning:  to save space I have put the tau times the modified residual
c           in rLymi and the tau times the total residual back in rLyi
c
c
c.... ---------------------------->  RHS  <----------------------------
c
c.... calculate (A_i^T tau Ly)
c

       if(ires.ne.1) then
c
c  A1 * Tau L(Y):  to be hit on left with Na,x in e3wmlt
c
        rmi(:,1) =
     &               A1(:,1,1) * rLymi(:,1)
     &             + A1(:,1,2) * rLymi(:,2)
c    &             + A1(:,1,3) * rLymi(:,3)
c    &             + A1(:,1,4) * rLymi(:,4)
     &             + A1(:,1,5) * rLymi(:,5)
     &             + rmi(:,1)
        rmi(:,2) =
     &               A1(:,2,1) * rLymi(:,1)
     &             + A1(:,2,2) * rLymi(:,2)
c    &             + A1(:,2,3) * rLymi(:,3)
c    &             + A1(:,2,4) * rLymi(:,4)
     &             + A1(:,2,5) * rLymi(:,5)
     &             + rmi(:,2)
        rmi(:,3) =
     &               A1(:,3,1) * rLymi(:,1)
     &             + A1(:,3,2) * rLymi(:,2)
     &             + A1(:,3,3) * rLymi(:,3)
c    &             + A1(:,3,4) * rLymi(:,4)
     &             + A1(:,3,5) * rLymi(:,5)
     &             + rmi(:,3)
        rmi(:,4) =
     &               A1(:,4,1) * rLymi(:,1)
     &             + A1(:,4,2) * rLymi(:,2)
c    &             + A1(:,4,3) * rLymi(:,3)
     &             + A1(:,4,4) * rLymi(:,4)
     &             + A1(:,4,5) * rLymi(:,5)
     &             + rmi(:,4)
        rmi(:,5) =
     &               A1(:,5,1) * rLymi(:,1)
     &             + A1(:,5,2) * rLymi(:,2)
     &             + A1(:,5,3) * rLymi(:,3)
     &             + A1(:,5,4) * rLymi(:,4)
     &             + A1(:,5,5) * rLymi(:,5)
     &             + rmi(:,5)
c
c  A2 * Tau L(Y),  to be hit on left with Na,y
c
        rmi(:,6) =
     &               A2(:,1,1) * rLymi(:,1)
c    &             + A2(:,1,2) * rLymi(:,2)
     &             + A2(:,1,3) * rLymi(:,3)
c    &             + A2(:,1,4) * rLymi(:,4)
     &             + A2(:,1,5) * rLymi(:,5)
     &             + rmi(:,6)
        rmi(:,7) =
     &               A2(:,2,1) * rLymi(:,1)
     &             + A2(:,2,2) * rLymi(:,2)
     &             + A2(:,2,3) * rLymi(:,3)
c    &             + A2(:,2,4) * rLymi(:,4)
     &             + A2(:,2,5) * rLymi(:,5)
     &             + rmi(:,7)
        rmi(:,8) =
     &               A2(:,3,1) * rLymi(:,1)
c    &             + A2(:,3,2) * rLymi(:,2)
     &             + A2(:,3,3) * rLymi(:,3)
c    &             + A2(:,3,4) * rLymi(:,4)
     &             + A2(:,3,5) * rLymi(:,5)
     &             + rmi(:,8)
        rmi(:,9) =
     &               A2(:,4,1) * rLymi(:,1)
c    &             + A2(:,4,2) * rLymi(:,2)
     &             + A2(:,4,3) * rLymi(:,3)
     &             + A2(:,4,4) * rLymi(:,4)
     &             + A2(:,4,5) * rLymi(:,5)
     &             + rmi(:,9)
        rmi(:,10) =
     &               A2(:,5,1) * rLymi(:,1)
     &             + A2(:,5,2) * rLymi(:,2)
     &             + A2(:,5,3) * rLymi(:,3)
     &             + A2(:,5,4) * rLymi(:,4)
     &             + A2(:,5,5) * rLymi(:,5)
     &             + rmi(:,10)
c
c  A3 * Tau L(Y)  to be hit on left with Na,z
c
        rmi(:,11) =
     &               A3(:,1,1) * rLymi(:,1)
c    &             + A3(:,1,2) * rLymi(:,2)
c    &             + A3(:,1,3) * rLymi(:,3)
     &             + A3(:,1,4) * rLymi(:,4)
     &             + A3(:,1,5) * rLymi(:,5)
     &             + rmi(:,11)
        rmi(:,12) =
     &               A3(:,2,1) * rLymi(:,1)
     &             + A3(:,2,2) * rLymi(:,2)
c    &             + A3(:,2,3) * rLymi(:,3)
     &             + A3(:,2,4) * rLymi(:,4)
     &             + A3(:,2,5) * rLymi(:,5)
     &             + rmi(:,12)
        rmi(:,13) =
     &               A3(:,3,1) * rLymi(:,1)
c    &             + A3(:,3,2) * rLymi(:,2)
     &             + A3(:,3,3) * rLymi(:,3)
     &             + A3(:,3,4) * rLymi(:,4)
     &             + A3(:,3,5) * rLymi(:,5)
     &             + rmi(:,13)
        rmi(:,14) =
     &               A3(:,4,1) * rLymi(:,1)
c    &             + A3(:,4,2) * rLymi(:,2)
c    &             + A3(:,4,3) * rLymi(:,3)
     &             + A3(:,4,4) * rLymi(:,4)
     &             + A3(:,4,5) * rLymi(:,5)
     &             + rmi(:,14)
        rmi(:,15) =
     &               A3(:,5,1) * rLymi(:,1)
     &             + A3(:,5,2) * rLymi(:,2)
     &             + A3(:,5,3) * rLymi(:,3)
     &             + A3(:,5,4) * rLymi(:,4)
     &             + A3(:,5,5) * rLymi(:,5)
     &             + rmi(:,15)
      endif  !ires.ne.1

c
c  same thing for the real residual
c
       if(ires.eq.3 .or. ires .eq. 1) then  ! we need the total residual
        ri(:,1) =
     &               A1(:,1,1) * rLyi(:,1)
     &             + A1(:,1,2) * rLyi(:,2)
c    &             + A1(:,1,3) * rLyi(:,3)
c    &             + A1(:,1,4) * rLyi(:,4)
     &             + A1(:,1,5) * rLyi(:,5)
     &             + ri(:,1)
        ri(:,2) =
     &               A1(:,2,1) * rLyi(:,1)
     &             + A1(:,2,2) * rLyi(:,2)
c    &             + A1(:,2,3) * rLyi(:,3)
c    &             + A1(:,2,4) * rLyi(:,4)
     &             + A1(:,2,5) * rLyi(:,5)
     &             + ri(:,2)
        ri(:,3) =
     &               A1(:,3,1) * rLyi(:,1)
     &             + A1(:,3,2) * rLyi(:,2)
     &             + A1(:,3,3) * rLyi(:,3)
c    &             + A1(:,3,4) * rLyi(:,4)
     &             + A1(:,3,5) * rLyi(:,5)
     &             + ri(:,3)
        ri(:,4) =
     &               A1(:,4,1) * rLyi(:,1)
     &             + A1(:,4,2) * rLyi(:,2)
c    &             + A1(:,4,3) * rLyi(:,3)
     &             + A1(:,4,4) * rLyi(:,4)
     &             + A1(:,4,5) * rLyi(:,5)
     &             + ri(:,4)
        ri(:,5) =
     &               A1(:,5,1) * rLyi(:,1)
     &             + A1(:,5,2) * rLyi(:,2)
     &             + A1(:,5,3) * rLyi(:,3)
     &             + A1(:,5,4) * rLyi(:,4)
     &             + A1(:,5,5) * rLyi(:,5)
     &             + ri(:,5)
c
        ri(:,6) =
     &               A2(:,1,1) * rLyi(:,1)
c    &             + A2(:,1,2) * rLyi(:,2)
     &             + A2(:,1,3) * rLyi(:,3)
c    &             + A2(:,1,4) * rLyi(:,4)
     &             + A2(:,1,5) * rLyi(:,5)
     &             + ri(:,6)
        ri(:,7) =
     &               A2(:,2,1) * rLyi(:,1)
     &             + A2(:,2,2) * rLyi(:,2)
     &             + A2(:,2,3) * rLyi(:,3)
c    &             + A2(:,2,4) * rLyi(:,4)
     &             + A2(:,2,5) * rLyi(:,5)
     &             + ri(:,7)
        ri(:,8) =
     &               A2(:,3,1) * rLyi(:,1)
c    &             + A2(:,3,2) * rLyi(:,2)
     &             + A2(:,3,3) * rLyi(:,3)
c    &             + A2(:,3,4) * rLyi(:,4)
     &             + A2(:,3,5) * rLyi(:,5)
     &             + ri(:,8)
        ri(:,9) =
     &               A2(:,4,1) * rLyi(:,1)
c    &             + A2(:,4,2) * rLyi(:,2)
     &             + A2(:,4,3) * rLyi(:,3)
     &             + A2(:,4,4) * rLyi(:,4)
     &             + A2(:,4,5) * rLyi(:,5)
     &             + ri(:,9)
        ri(:,10) =
     &               A2(:,5,1) * rLyi(:,1)
     &             + A2(:,5,2) * rLyi(:,2)
     &             + A2(:,5,3) * rLyi(:,3)
     &             + A2(:,5,4) * rLyi(:,4)
     &             + A2(:,5,5) * rLyi(:,5)
     &             + ri(:,10)
        ri(:,11) =
     &               A3(:,1,1) * rLyi(:,1)
c    &             + A3(:,1,2) * rLyi(:,2)
c    &             + A3(:,1,3) * rLyi(:,3)
     &             + A3(:,1,4) * rLyi(:,4)
     &             + A3(:,1,5) * rLyi(:,5)
     &             + ri(:,11)
        ri(:,12) =
     &               A3(:,2,1) * rLyi(:,1)
     &             + A3(:,2,2) * rLyi(:,2)
c    &             + A3(:,2,3) * rLyi(:,3)
     &             + A3(:,2,4) * rLyi(:,4)
     &             + A3(:,2,5) * rLyi(:,5)
     &             + ri(:,12)
        ri(:,13) =
     &               A3(:,3,1) * rLyi(:,1)
c    &             + A3(:,3,2) * rLyi(:,2)
     &             + A3(:,3,3) * rLyi(:,3)
     &             + A3(:,3,4) * rLyi(:,4)
     &             + A3(:,3,5) * rLyi(:,5)
     &             + ri(:,13)
        ri(:,14) =
     &               A3(:,4,1) * rLyi(:,1)
c    &             + A3(:,4,2) * rLyi(:,2)
c    &             + A3(:,4,3) * rLyi(:,3)
     &             + A3(:,4,4) * rLyi(:,4)
     &             + A3(:,4,5) * rLyi(:,5)
     &             + ri(:,14)
        ri(:,15) =
     &               A3(:,5,1) * rLyi(:,1)
     &             + A3(:,5,2) * rLyi(:,2)
     &             + A3(:,5,3) * rLyi(:,3)
     &             + A3(:,5,4) * rLyi(:,4)
     &             + A3(:,5,5) * rLyi(:,5)
     &             + ri(:,15)
c
       endif ! for ires=3 or 1

c
c.... ---------------------------->  LHS  <----------------------------
c
       if (lhs .eq. 1) then
c
c.... calculate (Atau) <-- (A_1 tau) (Recall that we are using a
c                                     diagonal tau here)
c

          if (itau.lt.10) then

             do i = 1, nflow
                Atau(:,i,1) = A1(:,i,1)*tau(:,1)
                Atau(:,i,2) = A1(:,i,2)*tau(:,2)
                Atau(:,i,3) = A1(:,i,3)*tau(:,2)
                Atau(:,i,4) = A1(:,i,4)*tau(:,2)
                Atau(:,i,5) = A1(:,i,5)*tau(:,3)
             enddo

          else

             Atau = zero
             do i = 1, nflow
                do j = 1, nflow
                   do k = 1, nflow
                      Atau(:,i,j) =Atau(:,i,j) + A1(:,i,k)*PTau(:,k,j)
                   enddo
                enddo
             enddo

          endif
c
c.... calculate (A_1 tau A_0) (for L.S. time term of EGmass)
c
       do j = 1, nflow
          do i = 1, nflow
             A1tauA0(:,i,j) =
     &            Atau(:,i,1)*A0(:,1,j) +
     &            Atau(:,i,2)*A0(:,2,j) +
     &            Atau(:,i,3)*A0(:,3,j) +
     &            Atau(:,i,4)*A0(:,4,j) +
     &            Atau(:,i,5)*A0(:,5,j)
          enddo
       enddo
c
c.... add (A_1 tau A_1) to stiff [1,1]
c
       do j = 1, nflow
          do i = 1, nflow
             stiff(:,i,j) = stiff(:,i,j) + (
     &              Atau(:,i,1)*A1(:,1,j)
     &            + Atau(:,i,2)*A1(:,2,j)
     &            + Atau(:,i,3)*A1(:,3,j)
     &            + Atau(:,i,4)*A1(:,4,j)
     &            + Atau(:,i,5)*A1(:,5,j)
     &            )
          enddo
       enddo
c
c.... add (A_1 tau A_2) to stiff [1,2]
c
       do j = 1, nflow
          do i = 1, nflow
             stiff(:,i,j+5) = stiff(:,i,j+5) + (
     &              Atau(:,i,1)*A2(:,1,j)
     &            + Atau(:,i,2)*A2(:,2,j)
     &            + Atau(:,i,3)*A2(:,3,j)
     &            + Atau(:,i,4)*A2(:,4,j)
     &            + Atau(:,i,5)*A2(:,5,j)
     &            )
          enddo
       enddo
c
c.... add (A_1 tau A_3) to stiff [1,3]
c
       do j = 1, nflow
          do i = 1, nflow
             stiff(:,i,j+10) = stiff(:,i,j+10) + (
     &              Atau(:,i,1)*A3(:,1,j)
     &            + Atau(:,i,2)*A3(:,2,j)
     &            + Atau(:,i,3)*A3(:,3,j)
     &            + Atau(:,i,4)*A3(:,4,j)
     &            + Atau(:,i,5)*A3(:,5,j)
     &            )
          enddo
       enddo
c
c.... calculate (Atau) <-- (A_2 tau) (Recall that we are using a
c                                     diagonal tau here)
c
       if (itau.lt.10) then

          do i = 1, nflow
             Atau(:,i,1) = A2(:,i,1)*tau(:,1)
             Atau(:,i,2) = A2(:,i,2)*tau(:,2)
             Atau(:,i,3) = A2(:,i,3)*tau(:,2)
             Atau(:,i,4) = A2(:,i,4)*tau(:,2)
             Atau(:,i,5) = A2(:,i,5)*tau(:,3)
          enddo

       else
          Atau = zero
          do i = 1, nflow
             do j = 1, nflow
                do k = 1, nflow
                   Atau(:,i,j) = Atau(:,i,j) + A2(:,i,k)*PTau(:,k,j)
                enddo
             enddo
          enddo

       endif
c
c.... calculate (A_2 tau A_0) (for L.S. time term of EGmass)
c
       do j = 1, nflow
          do i = 1, nflow
             A2tauA0(:,i,j) =
     &            Atau(:,i,1)*A0(:,1,j) +
     &            Atau(:,i,2)*A0(:,2,j) +
     &            Atau(:,i,3)*A0(:,3,j) +
     &            Atau(:,i,4)*A0(:,4,j) +
     &            Atau(:,i,5)*A0(:,5,j)
          enddo
       enddo
c
c.... add (A_2 tau A_1) to stiff [2,1]
c
       do j = 1, nflow
          do i = 1, nflow
             stiff(:,i+5,j) = stiff(:,i+5,j) + (
     &              Atau(:,i,1)*A1(:,1,j)
     &            + Atau(:,i,2)*A1(:,2,j)
     &            + Atau(:,i,3)*A1(:,3,j)
     &            + Atau(:,i,4)*A1(:,4,j)
     &            + Atau(:,i,5)*A1(:,5,j)
     &            )
          enddo
       enddo
c
c.... add (A_2 tau A_2) to stiff [2,2]
c
       do j = 1, nflow
          do i = 1, nflow
             stiff(:,i+5,j+5) = stiff(:,i+5,j+5) + (
     &              Atau(:,i,1)*A2(:,1,j)
     &            + Atau(:,i,2)*A2(:,2,j)
     &            + Atau(:,i,3)*A2(:,3,j)
     &            + Atau(:,i,4)*A2(:,4,j)
     &            + Atau(:,i,5)*A2(:,5,j)
     &            )
          enddo
       enddo
c
c.... add (A_2 tau A_3) to stiff [2,3]
c
       do j = 1, nflow
          do i = 1, nflow
             stiff(:,i+5,j+10) = stiff(:,i+5,j+10) + (
     &              Atau(:,i,1)*A3(:,1,j)
     &            + Atau(:,i,2)*A3(:,2,j)
     &            + Atau(:,i,3)*A3(:,3,j)
     &            + Atau(:,i,4)*A3(:,4,j)
     &            + Atau(:,i,5)*A3(:,5,j)
     &            )
          enddo
       enddo
c
c.... calculate (Atau) <-- (A_3 tau) (Recall that we are using a
c                                     diagonal tau here)
c
       if (itau.lt.10) then

          do i = 1, nflow
             Atau(:,i,1) = A3(:,i,1)*tau(:,1)
             Atau(:,i,2) = A3(:,i,2)*tau(:,2)
             Atau(:,i,3) = A3(:,i,3)*tau(:,2)
             Atau(:,i,4) = A3(:,i,4)*tau(:,2)
             Atau(:,i,5) = A3(:,i,5)*tau(:,3)
          enddo

       else
          Atau = zero
          do i = 1, nflow
             do j = 1, nflow
                do k = 1, nflow
                   Atau(:,i,j) = Atau(:,i,j) + A3(:,i,k)*PTau(:,k,j)
                enddo
             enddo
          enddo

       endif
c
c.... calculate (A_3 tau A_0) (for L.S. time term of EGmass)
c
       do j = 1, nflow
          do i = 1, nflow
             A3tauA0(:,i,j) =
     &            Atau(:,i,1)*A0(:,1,j) +
     &            Atau(:,i,2)*A0(:,2,j) +
     &            Atau(:,i,3)*A0(:,3,j) +
     &            Atau(:,i,4)*A0(:,4,j) +
     &            Atau(:,i,5)*A0(:,5,j)
          enddo
       enddo
c
c.... add (A_3 tau A_1) to stiff [3,1]
c
       do j = 1, nflow
          do i = 1, nflow
             stiff(:,i+10,j) = stiff(:,i+10,j) + (
     &              Atau(:,i,1)*A1(:,1,j)
     &            + Atau(:,i,2)*A1(:,2,j)
     &            + Atau(:,i,3)*A1(:,3,j)
     &            + Atau(:,i,4)*A1(:,4,j)
     &            + Atau(:,i,5)*A1(:,5,j)
     &            )
          enddo
       enddo
c
c.... add (A_3 tau A_2) to stiff [3,2]
c
       do j = 1, nflow
          do i = 1, nflow
             stiff(:,i+10,j+5) = stiff(:,i+10,j+5) + (
     &              Atau(:,i,1)*A2(:,1,j)
     &            + Atau(:,i,2)*A2(:,2,j)
     &            + Atau(:,i,3)*A2(:,3,j)
     &            + Atau(:,i,4)*A2(:,4,j)
     &            + Atau(:,i,5)*A2(:,5,j)
     &            )
          enddo
       enddo
c
c.... add (A_3 tau A_3) to stiff [3,3]
c
       do j = 1, nflow
          do i = 1, nflow
             stiff(:,i+10,j+10) = stiff(:,i+10,j+10) + (
     &              Atau(:,i,1)*A3(:,1,j)
     &            + Atau(:,i,2)*A3(:,2,j)
     &            + Atau(:,i,3)*A3(:,3,j)
     &            + Atau(:,i,4)*A3(:,4,j)
     &            + Atau(:,i,5)*A3(:,5,j)
     &            )
          enddo
       enddo
c
c.... add least squares time term to the LHS tangent mass matrix
c
c
c.... loop through rows (nodes i)
c
       do i = 1, nshl
          i0 = nflow * (i - 1)
c
c.... first calculate (Atau) <-- (N_a,i A_i tau A_0)
c     ( use Atau to conserve space )
c
          do idof = 1, nflow
             do jdof = 1, nflow
                Atau(:,idof,jdof) =
     &               shg(:,i,1) * A1tauA0(:,idof,jdof) +
     &               shg(:,i,2) * A2tauA0(:,idof,jdof) +
     &               shg(:,i,3) * A3tauA0(:,idof,jdof)
             enddo
          enddo
c
c.... loop through column nodes, add (N_a,i A_i tau N_b) to EGmass
c
          do j = 1, nshl
             j0 = nflow * (j - 1)
c
c.... compute the factors
c
            fact = shape(:,j) * WdetJ * almi/gami/alfi*dtgl
c
c.... loop through d.o.f.'s
c
            do idof = 1, nflow
               il = i0 + idof

               EGmass(:,il,j0+1) = EGmass(:,il,j0+1) +
     &                             fact * Atau(:,idof,1)
               EGmass(:,il,j0+2) = EGmass(:,il,j0+2) +
     &                             fact * Atau(:,idof,2)
               EGmass(:,il,j0+3) = EGmass(:,il,j0+3) +
     &                             fact * Atau(:,idof,3)
               EGmass(:,il,j0+4) = EGmass(:,il,j0+4) +
     &                             fact * Atau(:,idof,4)
               EGmass(:,il,j0+5) = EGmass(:,il,j0+5) +
     &                             fact * Atau(:,idof,5)
            enddo
c
c.... end loop on column nodes
c
         enddo
c
c.... end loop on row nodes
c
       enddo
c
c.... end LHS computation
c
       endif

       ttim(24) = ttim(24) + secs(0.0)
c
c.... return
c
        return
        end
c
c
c
        subroutine e3LSSclr  (A1t,     A2t,     A3t,
     &                        rho,     rmu,     rTLSt,
     &                        u1,      u2,      u3,
     &                        rLyti,   dxidx,   raLSt,
     &                        rti,     rk,      giju,
     &                        acti,    A0t,
     &                        shape,   shg,
     &                        EGmasst, stifft,  WdetJ,
     &                        srcp)
c
c----------------------------------------------------------------------
c
c This routine calculates the contribution of the least-squares
c operator to the RHS vector and LHS tangent matrix. The temporary
c results are put in ri.
c
c input:
c  A0t    (npro)              : A_0
c  A1t    (npro)              : A_1
c  A2t    (npro)              : A_2
c  A3t    (npro)              : A_3
c  acti   (npro)              : time-deriv. of Sclr
c  rho    (npro)              : density
c  rmu    (npro)              : molecular viscosity
c  rk     (npro)              : kinetic energy
c  u1     (npro)              : x1-velocity component
c  u2     (npro)              : x2-velocity component
c  u3     (npro)              : x3-velocity component
c  rLyti  (npro)              : least-squares residual vector
c  dxidx  (npro,nsd,nsd)      : inverse of deformation gradient
c  taut   (npro)              : stability parameter
c  rLyti  (npro)              : convective portion of least-squares
c                               residual vector
c  divqti (npro,1)            : divergence of diffusive flux
c  shape  (npro,nshl)         : element shape functions
c  shg    (npro,nshl,nsd)     : global element shape function grads
c  WdetJ  (npro)              : weighted jacobian determinant
c  stifft (npro,nsd,nsd)      : stiffness matrix
c  EGmasst(npro,nshape,nshape): partial mass matrix
c
c output:
c  rti    (npro,nsd+1)        : partial residual
c  EGmasst(npro,nshape,nshape): partial mass matrix
c
c
c Zdenek Johan, Summer 1990. (Modified from e2ls.f)
c Zdenek Johan, Winter 1991. (Fortran 90)
c Kenneth Jansen, Winter 1997. Prim. Var. with Diag Tau
c----------------------------------------------------------------------
c
        include "common.h"

c
c  passed arrays
c
        dimension A1t(npro),                 A2t(npro),
     &            A3t(npro),
     &            A0t(npro),                 rho(npro),
     &            acti(npro),                rmu(npro),
     &            u1(npro),                  u2(npro),
     &            u3(npro),                  rk(npro),
     &            rLyti(npro),               dxidx(npro,nsd,nsd),
     &            taut(npro),                raLSt(npro),
     &            rti(npro,nsd+1),           rTLSt(npro),
     &            stifft(npro,3,3),          giju(npro,6),
     &            EGmasst(npro,nshape,nshape),
     &            shape(npro,nshl),
     &            shg(npro,nshl,nsd),        WdetJ(npro),
     &            srcp(npro)
c
c local arrays
c
        dimension rLymti(npro),          Ataut(npro),
     &            A1tautA0(npro),        A2tautA0(npro),
     &            A3tautA0(npro),        fact(npro)

        ttim(24) = ttim(24) - tmr()
c
       if(ivart.lt.2) return
c
c last step to the least squares is adding the time term.  So that we
c only have to localize one vector for each Krylov vector the modified
c residual is quite different from the total residual.
c
c
c the modified residual
c
       fct1=almi/gami/alfi*dtgl
c
c  add possibility of not including time term
c
c       if(idiff.ne.-1)
c       rLymti = rLyti + fct1*duti

       if((ires.eq.1 .or. ires .eq. 3).and. idiff.ne.-1) then

        rLyti(:) = rLyti(:) + A0t(:)*acti(:)

      endif
c
c.... subtract div(q) from the least squares term
c
c      if ((idiff >= 1).and.(ires==3 .or. ires==1)) then
c         rLyi(:) = rLyi(:) - divqti(:)
c      endif
c
c.... --------------------------->  Tau  <-----------------------------
c
c.... calculate the tau matrix
c

c
c.... we will use the same tau as used for momentum equations here
c
        ttim(25) = ttim(25) - tmr()

          call e3tauSclr(rho,         rmu,        A0t,
     &                   u1,          u2,         u3,
     &                   dxidx,       rLyti,      rLymti,
     &                   taut,        rk,         raLSt,
     &                   rTLSt,       giju)

        ttim(25) = ttim(25) + tmr()
c
c Warning:  to save space I have put the tau times the modified residual
c           in rLymi and the tau times the total residual back in rLyi
c
c
c.... ---------------------------->  RHS  <----------------------------
c
c.... calculate (A_i^T tau Ly)
c

c      if(ires.ne.1) then
c
c  A1 * Tau L(Y):  to be hit on left with Na,x in e3wmlt
c
c        rmti(:,1) =   A1t(:) * rLymti(:)
c
c
c  A2 * Tau L(Y),  to be hit on left with Na,y
c
c        rmti(:,2) = A2t(:) * rLymti(:)
c
c
c  A3 * Tau L(Y)  to be hit on left with Na,z
c
c        rmti(:,3) = A3t(:) * rLymti(:)
c
c      endif  !ires.ne.1

c
c  same thing for the real residual
c
       if(ires.eq.3 .or. ires .eq. 1) then  ! we need the total residual
        rti(:,1) = rti(:,1) + A1t(:) * rLyti(:)

        rti(:,2) = rti(:,2) + A2t(:) * rLyti(:)

        rti(:,3) = rti(:,3) + A3t(:) * rLyti(:)

       endif ! for ires=3 or 1

c
c.... ---------------------------->  LHS  <----------------------------
c
       if (lhs .eq. 1) then
c
c
c.... calculate (Atau) <-- (A_1 tau)
c

          Ataut(:) = A1t(:)*taut(:)

c
c.... calculate (A_1 tau (A_0-srcp)) (for L.S. time term of EGmass)
c

             A1tautA0(:) = Ataut(:)*(a0t(:)*fct1-srcp(:))

c
c.... add (A_1 tau A_1) to stiff [1,1]
c

             stifft(:,1,1) = stifft(:,1,1) + Ataut(:)*A1t(:)
c            stifft(:,1,1) = Ataut(:)*A1t(:)
c
c.... add (A_1 tau A_2) to stiff [1,2]
c

             stifft(:,1,2) = stifft(:,1,2) + Ataut(:)*A2t(:)
c            stifft(:,1,2) =  Ataut(:)*A2t(:)
c
c.... add (A_1 tau A_3) to stiff [1,3]
c

             stifft(:,1,3) = stifft(:,1,3) + Ataut(:)*A3t(:)
c            stifft(:,1,3) =  Ataut(:)*A3t(:)
c
c.... calculate (Atau) <-- (A_2 tau)
c

          Ataut(:) = A2t(:)*taut(:)

c
c.... calculate (A_2 tau (A_0-srcp)) (for L.S. time term of EGmass)
c

             A2tautA0(:) = Ataut(:)*(a0t(:)*fct1-srcp(:))

c
c.... add (A_2 tau A_1) to stiff [2,1]
c

             stifft(:,2,1) = stifft(:,1,2)
c
c.... add (A_2 tau A_2) to stiff [2,2]
c

             stifft(:,2,2) = stifft(:,2,2) + Ataut(:)*A2t(:)

c
c.... add (A_2 tau A_3) to stiff [2,3]
c

             stifft(:,2,3) = stifft(:,2,3) + Ataut(:)*A3t(:)

c
c.... calculate (Atau) <-- (A_3 tau)
c

          Ataut(:) = A3t(:)*taut(:)

c
c.... calculate (A_3 tau (A_0-srcp)) (for L.S. time term of EGmass)
c

             A3tautA0(:) = Ataut(:)*(a0t(:)*fct1-srcp(:))

c
c.... add (A_3 tau A_1) to stiff [3,1]
c

             stifft(:,3,1) = stifft(:,1,3)

c
c.... add (A_3 tau A_2) to stiff [3,2]
c

             stifft(:,3,2) = stifft(:,2,3)

c
c.... add (A_3 tau A_3) to stiff [3,3]
c

             stifft(:,3,3) = stifft(:,3,3) + Ataut(:)*A3t(:)

c
c.... add least squares time term to the LHS tangent mass matrix
c
c
c.... loop through rows (nodes i)
c
       do ia = 1, nshl
c
c.... first calculate (Atau) <-- (N_a,i A_i tau A_0)
c     ( use Atau to conserve space )
c

                Ataut(:) =
     &               shg(:,ia,1) * A1tautA0(:) +
     &               shg(:,ia,2) * A2tautA0(:) +
     &               shg(:,ia,3) * A3tautA0(:)

c
c.... loop through column nodes, add (N_a,i A_i tau N_b) to EGmass
c
          do jb = 1, nshl

            fact = shape(:,jb) * WdetJ

            EGmasst(:,ia,jb) = EGmasst(:,ia,jb) + fact * Ataut(:)

c
c.... end loop on column nodes
c

          enddo
c
c.... end loop on row nodes
c
       enddo
c
c.... end LHS computation
c
       endif

       ttim(24) = ttim(24) + tmr()
c
c.... return
c
        return
        end