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

        subroutine e3 (yl,      ycl,     acl,     shp,
     &                 shgl,    xl,      rl,      rml,    xmudmi,
     &                 BDiagl,  ql,      sgn,     rlsl,   EGmass,
     &                 rerrl,   ytargetl)
c
c----------------------------------------------------------------------
c
c This routine is the 3D element routine for the N-S equations.
c This routine calculates the RHS residual and if requested the
c modified residual or the LHS consistent mass matrix or the block-
c diagonal preconditioner.
c
c input:
c  yl     (npro,nshl,nflow)     : Y variables  (DOES NOT CONTAIN SCALARS)
c  ycl    (npro,nshl,ndof)      : Y variables at current step
c  acl    (npro,nshl,ndof)      : Y acceleration (Y_{,t})
c  shp    (nshl,ngauss)       : element shape-functions  N_a
c  shgl   (nsd,nshl,ngauss)   : element local-shape-functions N_{a,xi}
c  xl     (npro,nenl,nsd)       : nodal coordinates at current step (x^e_a)
c  ql     (npro,nshl,(nflow-1)*nsd) : diffusive flux vector
c  rlsl   (npro,nshl,6)       : resolved Leonard stresses
c  sgn    (npro,nshl)         : shape function sign matrix
c
c output:
c  rl     (npro,nshl,nflow)      : element RHS residual    (G^e_a)
c  rml    (npro,nshl,nflow)      : element modified residual  (G^e_a tilde)
c  EGmass (npro,nedof,nedof)    : element LHS tangent mass matrix (dG^e_a
c                                                                  dY_b  )
c  BDiagl (npro,nshl,nflow,nflow) : block-diagonal preconditioner
c
c
c Note: This routine will calculate the element matrices for the
c        Hulbert's generalized alpha method integrator
c
c Note: nedof = nflow*nshape is the total number of degrees of freedom
c       at each element node.
c
c Mathematics done by:  Michel Mallet, Farzin Shakib (version 1)
c                       Farzin Shakib                (version 2)
c
c
c Zdenek Johan, Summer 1990.   (Modified from e2.f)
c Zdenek Johan, Winter 1991.   (Fortran 90)
c Kenneth Jansen, Winter 1997. (Primitive Variables)
c Chris Whiting, Winter 1998.  (LHS matrix formation)
c----------------------------------------------------------------------
c
        include "common.h"
c
        dimension yl(npro,nshl,nflow),     ycl(npro,nshl,ndof),
     &            acl(npro,nshl,ndof),     rmu(npro),
     &            shp(nshl,ngauss),        rlm2mu(npro),
     &            shgl(nsd,nshl,ngauss),   con(npro),
     &            xl(npro,nenl,nsd),       rlm(npro),
     &            rl(npro,nshl,nflow),     ql(npro,nshl,idflx),
     &            rml(npro,nshl,nflow),    xmudmi(npro,ngauss),
     &            BDiagl(npro,nshl,nflow,nflow),
     &            EGmass(npro,nedof,nedof),cv(npro),
     &            ytargetl(npro,nshl,nflow)
c
        dimension dui(npro,ndof),            aci(npro,ndof)
c
        dimension g1yi(npro,nflow),          g2yi(npro,nflow),
     &            g3yi(npro,nflow),          shg(npro,nshl,nsd),
     &            divqi(npro,nflow),       tau(npro,5)
c
        dimension dxidx(npro,nsd,nsd),       WdetJ(npro)
c
        dimension rho(npro),                 pres(npro),
     &            T(npro),                   ei(npro),
     &            h(npro),                   alfap(npro),
     &            betaT(npro),               DC(npro,ngauss),
     &            cp(npro),                  rk(npro),
     &            u1(npro),                  u2(npro),
     &            u3(npro),                  A0DC(npro,4),
     &            Sclr(npro),                dVdY(npro,15),
     &            giju(npro,6),              rTLS(npro),
     &            raLS(npro),                A0inv(npro,15)
c
        dimension A0(npro,nflow,nflow),      A1(npro,nflow,nflow),
     &            A2(npro,nflow,nflow),      A3(npro,nflow,nflow)
c
        dimension rLyi(npro,nflow),          sgn(npro,nshl)
c
        dimension ri(npro,nflow*(nsd+1)),    rmi(npro,nflow*(nsd+1)),
     &            shape(npro,nshl),          shdrv(npro,nsd,nshl),
     &            stiff(npro,nsd*nflow,nsd*nflow),
     &            PTau(npro,5,5),
     &            sforce(npro,3),      compK(npro,10)
c
        dimension x(npro,3),              bcool(npro)

        dimension rlsl(npro,nshl,6),      rlsli(npro,6)

        real*8    rerrl(npro,nshl,6)
        ttim(6) = ttim(6) - secs(0.0)
c
c.... local reconstruction of diffusive flux vector
c (note: not currently included in mfg)
        if (idiff==2 .and. (ires==3 .or. ires==1)) then
           call e3ql (ycl, shp, shgl, xl, ql, xmudmi, sgn)
        endif
c
c.... loop through the integration points
c
        do intp = 1, ngauss
c
c.... if Det. .eq. 0, do not include this point
c
        if (Qwt(lcsyst,intp) .eq. zero) cycle          ! precaution
c
c.... create a matrix of shape functions (and derivatives) for each
c     element at this quadrature point. These arrays will contain
c     the correct signs for the hierarchic basis
c
        call getshp(shp,          shgl,      sgn,
     &              shape,        shdrv)
c
c.... initialize
c
        ri  = zero
        rmi = zero
        if (lhs .eq. 1) stiff = zero
c
c
c.... calculate the integration variables
c
        ttim(8) = ttim(8) - secs(0.0)
        call e3ivar (yl,              ycl,             acl,
     &               Sclr,            shape,           shdrv,
     &               xl,              dui,             aci,
     &               g1yi,            g2yi,            g3yi,
     &               shg,             dxidx,           WdetJ,
     &               rho,             pres,            T,
     &               ei,              h,               alfap,
     &               betaT,           cp,              rk,
     &               u1,              u2,              u3,
     &               ql,              divqi,           sgn,
     &               rLyi,  !passed as a work array
     &               rmu,             rlm,             rlm2mu,
     &               con,             rlsl,            rlsli,
     &               xmudmi,          sforce,          cv)
        ttim(8) = ttim(8) + secs(0.0)

c
c.... calculate the relevant matrices
c
        ttim(9) = ttim(9) - secs(0.0)
        call e3mtrx (rho,             pres,           T,
     &               ei,              h,              alfap,
     &               betaT,           cp,             rk,
     &               u1,              u2,             u3,
     &               A0,              A1,
     &               A2,              A3,
     &               rLyi(:,1),       rLyi(:,2),      rLyi(:,3),  ! work arrays
     &               rLyi(:,4),       rLyi(:,5),      A0DC,
     &               A0inv,           dVdY)
        ttim(9) = ttim(9) + tmr()
c
c.... calculate the convective contribution (Galerkin)
c
        ttim(14) = ttim(14) - secs(0.0)
        call e3conv (g1yi,            g2yi,            g3yi,
     &               A1,              A2,              A3,
     &               rho,             pres,            T,
     &               ei,              rk,              u1,
     &               u2,              u3,              rLyi,
     &               ri,              rmi,             EGmass,
     &               shg,             shape,           WdetJ)
        ttim(14) = ttim(14) + secs(0.0)
c
c.... calculate the diffusion contribution
c
        ttim(15) = ttim(15) - secs(0.0)
        compK = zero
        if (Navier .eq. 1) then
        call e3visc (g1yi,            g2yi,            g3yi,
     &               dxidx,
     &               rmu,             rlm,             rlm2mu,
     &               u1,              u2,              u3,
     &               ri,              rmi,             stiff,
     &               con, rlsli,     compK, T)
        endif
        ttim(15) = ttim(15) + secs(0.0)
c
c.... calculate the body force contribution
c
        if(isurf .ne. 1 .and. matflg(5,1).gt.0) then
        call e3source (ri,            rmi,           rlyi,
     &                 rho,           u1,            u2,
     &                 u3,            pres,          sforce,
     &                 dui,           dxidx,         ytargetl,
     &                 xl,            shape,         bcool)
        else
           bcool=zero
        endif
c
c.... calculate the least-squares contribution
c
        ttim(16) = ttim(16) - secs(0.0)
        call 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)
        ttim(16) = ttim(16) + secs(0.0)
c
c....  Discontinuity capturing
c
        if(iDC.ne.0) then
          call e3dc  (g1yi,          g2yi,          g3yi,
     &                A0,            raLS,          rTLS,
     &                giju,          DC,
     &                ri,            rmi,           stiff, A0DC)
        endif
! SAM wants a threshold here so we are going to take over this little used
! error indictor for that purpose.  To revert note you will want to uncomment the original
! form of this error indicator in e3LS.f
        if((intp.eq.1).and.(ierrcalc.eq.1).and.(nitr.eq.iter))  then
          do i=1,npro
             Tmax=maxval(yl(i,:,5))
             Tmin=minval(yl(i,:,5))
             rerrl(i,:,6)=(Tmax-Tmin)/T(i)
          enddo
! the below was somewhat suprisingly ineffective compared to above for identifying shocks.
! it  refined on each side of the shock but left the actual shock quite coarse whereas the above
! centered well on the shock
!          do j=1,nshl
!            rerrl(:,j,6)=rerrl(:,j,6)+DC(:,intp)  !super hack to get error indicator for shock capturing
!          enddo
        endif
c
c
c.... calculate the time derivative (mass) contribution to RHS
c
        if (ngauss .eq. 1 .and. nshl .eq. 4) then    ! trilinear tets
           ttim(17) = ttim(17) - secs(0.0)
           call e3juel (yl, acl, Sclr, A0, WdetJ, rl, rml)
           ttim(17) = ttim(17) + secs(0.0)
        else
           call e3massr (aci, dui, ri,  rmi, A0)
        endif

c
c.... calculate the time (mass) contribution to the LHS
c
        if (lhs .eq. 1) then
           call e3massl (bcool,shape,  WdetJ,   A0,  EGmass)
        endif
c
c....  calculate the preconditioner all at once now instead of in separate
c      routines
c
       if(iprec.eq.1 .and. lhs.ne.1)then
          ttim(18) = ttim(18) - secs(0.0)

          if (itau.lt.10) then

             call e3bdg(shape,       shg,             WdetJ,
     &		  A1,	       A2,	        A3,
     &		  A0,          bcool,           tau,
     &            u1,          u2,              u3,
     &            BDiagl,
     &            rmu,         rlm2mu,          con)

          else

             call e3bdg_nd(shape,       shg,             WdetJ,
     &		  A1,	       A2,	        A3,
     &		  A0,          bcool,           PTau,
     &            u1,          u2,              u3,
     &            BDiagl,
     &            rmu,         rlm2mu,          con)

          endif

       ttim(18) = ttim(18) + secs(0.0)
       endif
c
c
c.... multiply flux terms by shape functions and derivatives (from weight space for RHS and
c     by both the weight space and solution space for the LHS stiffness term)
c
       ttim(19) = ttim(19) - secs(0.0)
       call e3wmlt (shape,         shg,       WdetJ,
     &              ri,            rmi,       rl,
     &              rml,           stiff,     EGmass)
       ttim(19) = ttim(19) + secs(0.0)
c
c.... end of integration loop
c
      enddo

      ttim(6) = ttim(6) + secs(0.0)
c
c.... return
c
      return
      end
c
c
c
c
       subroutine e3Sclr (ycl,          acl,
     &                    dwl,         elDwl,
     &                    shp,         sgn,
     &                    shgl,        xl,
     &                    rtl,         rmtl,
     &                    qtl,         EGmasst)
c
c----------------------------------------------------------------------
c
c This routine is the 3D element routine for the N-S equations.
c This routine calculates the RHS residual and if requested the
c modified residual or the LHS consistent mass matrix or the block-
c diagonal preconditioner.
c
c input:    e    a   1..5   when we think of U^e_a  and U is 5 variables
c  ycl      (npro,nshl,ndof)     : Y variables
c  actl    (npro,nshl)          : scalar variable time derivative
c  dwl     (npro,nenl)          : distances to wall
c  shp     (nen,ngauss)          : element shape-functions  N_a
c  shgl    (nsd,nen,ngauss)      : element local-shape-functions N_{a,xi}
c  xl      (npro,nenl,nsd)      : nodal coordinates at current step (x^e_a)
c  qtl     (npro,nshl)          : diffusive flux vector (don't worry)
c
c output:
c  rtl     (npro,nshl)          : element RHS residual    (G^e_a)
c  rmtl    (npro,nshl)          : element modified residual  (G^e_a tilde)
c  EGmasst (npro,nshape,nshape) : element LHS tangent mass matrix (dG^e_a
c                                                                  dY_b  )
c
c
c Note: This routine will calculate the element matrices for the
c        Hulbert's generalized alpha method integrator
c
c Note: nedof = nflow*nshl is the total number of degrees of freedom
c       at each element node.
c
c Mathematics done by:  Michel Mallet, Farzin Shakib (version 1)
c                       Farzin Shakib                (version 2)
c
c
c Zdenek Johan, Summer 1990.   (Modified from e2.f)
c Zdenek Johan, Winter 1991.   (Fortran 90)
c Kenneth Jansen, Winter 1997. (Primitive Variables)
c Chris Whiting, Winter 1998.  (LHS matrix formation)
c----------------------------------------------------------------------
c
        include "common.h"
c
        dimension ycl(npro,nshl,ndof),
     &            acl(npro,nshl,ndof),
     &            dwl(npro,nenl),
     &            shp(nshl,ngauss),         shgl(nsd,nshl,ngauss),
     &            xl(npro,nenl,nsd),
     &            rtl(npro,nshl),           qtl(npro,nshl),
     &            rmtl(npro,nshl),          Diagl(npro,nshl),
     &            EGmasst(npro,nshape,nshape),
     &            dist2w(npro),             sgn(npro,nshl),
     &            vort(npro),               gVnrm(npro),
     &            rmu(npro),                con(npro),
     &            T(npro),                  cp(npro),
     &            g1yti(npro),              acti(npro),
     &            g2yti(npro),              g3yti(npro),
     &            Sclr(npro),               srcp (npro)

c
        dimension shg(npro,nshl,nsd)
c
        dimension dxidx(npro,nsd,nsd),     WdetJ(npro)
c
        dimension rho(npro),               rk(npro),
     &            u1(npro),                u2(npro),
     &            u3(npro)
c
        dimension A0t(npro),               A1t(npro),
     &            A2t(npro),               A3t(npro)
c
        dimension rLyti(npro),             raLSt(npro),
     &            rTLSt(npro),             giju(npro,6),
     &            DCt(npro, ngauss)
c
        dimension rti(npro,nsd+1),         rmti(npro,nsd+1),
     &            stifft(npro,nsd,nsd),
     &            shape(npro,nshl),        shdrv(npro,nsd,nshl)
        real*8    elDwl(npro)
c
        ttim(6) = ttim(6) - tmr()
c
c.... loop through the integration points
c
        elDwl(:)=zero
        do intp = 1, ngauss
c
c.... if Det. .eq. 0, do not include this point
c
        if (Qwt(lcsyst,intp) .eq. zero) cycle          ! precaution
c
c
c.... create a matrix of shape functions (and derivatives) for each
c     element at this quadrature point. These arrays will contain
c     the correct signs for the hierarchic basis
c
        call getshp(shp,          shgl,      sgn,
     &              shape,        shdrv)
c
c.... initialize
c
        rlyti = zero
        rti   = zero
        rmti  = zero
        srcp  = zero
        stifft = zero
c        if (lhs .eq. 1) stifft = zero
c
c
c.... calculate the integration variables
c
        ttim(8) = ttim(8) - tmr()
c
        call e3ivarSclr(ycl,              acl,        acti,
     &                  shape,           shdrv,      xl,
     &                  T,               cp,
     &                  dxidx,           Sclr,
     &                  Wdetj,           vort,       gVnrm,
     &                  g1yti,           g2yti,      g3yti,
     &                  rho,             rmu,        con,
     &                  rk,              u1,         u2,
     &                  u3,              shg,        dwl,
     &                  dist2w)
        ttim(8) = ttim(8) + tmr()

c
c.... calculate the source term contribution
c
        if(nosource.ne.1)
     &  call e3sourceSclr (Sclr,    rho,    rmu,
     &                     dist2w,  vort,   gVnrm,   con,
     &                     g1yti,  g2yti,   g3yti,
     &                     rti,    rLyti,   srcp,
     &                     ycl,     shape,   u1,
     &                     u2,     u3,	     xl,
     &                     elDwl)
c
         if (ilset.eq.2 .and. isclr.eq.2) then
          rk = pt5 * ( u1**2 + u2**2 + u3**2 )
         endif
c.... calculate the relevant matrices
c
        ttim(9) = ttim(9) - tmr()
        call e3mtrxSclr (rho,
     &                   u1,            u2,         u3,
     &                   A0t,           A1t,
     &                   A2t,           A3t)
        ttim(9) = ttim(9) + tmr()
c
c.... calculate the convective contribution (Galerkin)
c
        ttim(14) = ttim(14) - tmr()
        call e3convSclr (g1yti,        g2yti,       g3yti,
     &                   A1t,          A2t,         A3t,
     &                   rho,          u1,          Sclr,
     &                   u2,           u3,          rLyti,
     &                   rti,          rmti,        EGmasst,
     &                   shg,          shape,       WdetJ)
        ttim(14) = ttim(14) + tmr()
c
c.... calculate the diffusion contribution
c
        ttim(15) = ttim(15) - tmr()
        if (Navier .eq. 1) then
        call e3viscSclr (g1yti,        g2yti,         g3yti,
     &                   rmu,          Sclr,          rho,
     &                   rti,          rmti,          stifft )
        endif
         ttim(15) = ttim(15) + tmr()
c
         if (ilset.eq.2)  srcp = zero

c
c.... calculate the least-squares contribution
c
        ttim(16) = ttim(16) - tmr()
        call e3LSSclr(A1t,             A2t,             A3t,
     &                rho,             rmu,             rtLSt,
     &                u1,              u2,              u3,
     &                rLyti,           dxidx,           raLSt,
     &                rti,             rk,              giju,
     &                acti,            A0t,
     &                shape,           shg,
     &                EGmasst,         stifft,          WdetJ,
     &                srcp)
        ttim(16) = ttim(16) + tmr()
c
c******************************DC TERMS*****************************
        if (idcsclr(1) .ne. 0) then
           if ((idcsclr(2).eq.1 .and. isclr.eq.1) .or.
     &         (idcsclr(2).eq.2 .and. isclr.eq.2)) then   ! scalar with dc
              call e3dcSclr(g1yti, g2yti,          g3yti,
     &             A0t,            raLSt,          rTLSt,
     &             DCt,            giju,
     &             rti,            rmti,           stifft)
           endif
        endif
c
c******************************************************************
c.... calculate the time derivative (mass) contribution to RHS
c

           call e3massrSclr (acti, rti, A0t)
c
c.... calculate the time (mass) contribution to the LHS
c
         if (lhs .eq. 1) then
            call e3masslSclr (shape,  WdetJ,   A0t,  EGmasst,srcp)
         endif
c

c.... multiply flux terms by shape functions and derivatives (from weight space for RHS and
c     by both the weight space and solution space for the LHS stiffness term)
c
       ttim(19) = ttim(19) - tmr()
       call e3wmltSclr (shape,           shg,             WdetJ,
     &                  rti,
     &                  rtl,             stifft,          EGmasst)
       ttim(19) = ttim(19) + tmr()
c
c.... end of the loop
c
      enddo
	qptinv=one/ngauss
	elDwl(:)=elDwl(:)*qptinv


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