xref: /phasta/phSolver/incompressible/e3res.f (revision 595995161822a203c8467e0e4a253d7bd7d6df32)

      subroutine e3Res ( u1,        u2,         u3,
     &                   uBar,      aci,        WdetJ,
     &                   g1yi,      g2yi,       g3yi,
     &                   rLui,      rmu,        rho,
     &                   tauC,      tauM,       tauBar,
     &                   shpfun,    shg,        src,
     &                   rl,        pres,       acl,
     &                   rlsli)
c------------------------------------------------------------------------
c
c  This routine computes the residual vector at an
c  integration point.
c
c  input:
c     u1(npro)                  : x1-velocity
c     u2(npro)                  : x2-velocity
c     u3(npro)                  : x3-velocity
c     uBar(npro,3)              : u - tauM * Li
c     aci(npro,3)               : acceleration
c     rlsli(npro,6)             : resolved Leonard stresses
c     WdetJ(npro)               : weighted jacobian determinant
c     g1yi(npro,ndof)              : x1-gradient of variables
c     g2yi(npro,ndof)              : x2-gradient of variables
c     g3yi(npro,ndof)              : x3-gradient of variables
c     rLui(npro,3)              : total residual of NS equations
c     rmu(npro)                 : fluid viscosity
c     rho(npro)                 : density
c     tauC(npro)                : continuity tau
c     tauM(npro)                : momentum tau
c     tauBar(npro)              : additional tau
c     shpfun(npro,nshl)         : element shape functions
c     shg(npro,nshl,nsd)        : global grad of element shape functions
c     src(npro,nsd)             : body force term
c
c  output:
c     rl(npro,nshl,nflow)
c
c------------------------------------------------------------------------
      include "common.h"

      dimension u1(npro),         u2(npro),         u3(npro),
     &          uBar(npro,nsd),   aci(npro,nsd),    WdetJ(npro),
     &          g1yi(npro,nflow), g2yi(npro,nflow), g3yi(npro,nflow),
     &          rLui(npro,nsd),   rmu(npro),        rho(npro),
     &          tauC(npro),       tauM(npro),       tauBar(npro),
     &          shpfun(npro,nshl),shg(npro,nshl,nsd), src(npro,nsd),
     &          pres(npro)

      dimension rl(npro,nshl,nflow),
     &          acl(npro,nshl,ndof),
     &          rlsli(npro,6)
c
c.... local declarations
c
      real*8    tmp1(npro),   tmp2(npro),          tmp3(npro),
     &          tmp(npro),    rGNa(npro,nsd,nsd),  rNa(npro,nsd),
     &          locmass(npro,nshl),omega(3)

      integer aa
c
c.... initialize multipliers for Na and Na_{,i}
c
      rNa  = zero
      rGNa = zero
c
c.... compute the Na multiplier
c
      tmps     = one-flmpr  ! consistant mass factor

c
c no density yet...it comes later
c
      rNa(:,1) = aci(:,1)  * tmps
     &         - src(:,1)
      rNa(:,2) = aci(:,2)  * tmps
     &         - src(:,2)
      rNa(:,3) = aci(:,3)  * tmps
     &         - src(:,3)

c
c.... rotating frame terms if needed
c

      if(matflg(6,1).eq.1) then ! rotation

         omega(1)=datmat(1,6,1)
         omega(2)=datmat(2,6,1)
         omega(3)=datmat(3,6,1)
c
c no density yet...it comes later
c

         rNa(:,1) = rNa(:,1) + (omega(2)-omega(3)) * tauM * rLui(:,1)
         rNa(:,2) = rNa(:,2) + (omega(3)-omega(1)) * tauM * rLui(:,2)
         rNa(:,3) = rNa(:,3) + (omega(1)-omega(2)) * tauM * rLui(:,3)
      endif


c
c.... compute the Na,i multiplier
c
      tmp  = -pres + tauC * (g1yi(:,2) + g2yi(:,3) + g3yi(:,4))
      tmp1 =  rmu * ( g2yi(:,2) + g1yi(:,3) )
      tmp2 =  rmu * ( g3yi(:,3) + g2yi(:,4) )
      tmp3 =  rmu * ( g1yi(:,4) + g3yi(:,2) )


      if(iconvflow.eq.2) then  ! advective form (NO IBP either)
c
c no density yet...it comes later
c
         rNa(:,1) = rNa(:,1)
     &            + ubar(:,1) * g1yi(:,2)
     &            + ubar(:,2) * g2yi(:,2)
     &            + ubar(:,3) * g3yi(:,2)
         rNa(:,2) = rNa(:,2)
     &            + ubar(:,1) * g1yi(:,3)
     &            + ubar(:,2) * g2yi(:,3)
     &            + ubar(:,3) * g3yi(:,3)
         rNa(:,3) = rNa(:,3)
     &            + ubar(:,1) * g1yi(:,4)
     &            + ubar(:,2) * g2yi(:,4)
     &            + ubar(:,3) * g3yi(:,4)

         rGNa(:,1,1) = two * rmu * g1yi(:,2) + tmp
         rGNa(:,1,2) = tmp1
         rGNa(:,1,3) = tmp3
         rGNa(:,2,1) = tmp1
         rGNa(:,2,2) = two * rmu * g2yi(:,3) + tmp
         rGNa(:,2,3) = tmp2
         rGNa(:,3,1) = tmp3
         rGNa(:,3,2) = tmp2
         rGNa(:,3,3) = two * rmu * g3yi(:,4) + tmp
      else   ! conservative form (with IBP)

c                                            IBP conservative convection
c                                                      |||||
c                                                      vvvvv
         rGNa(:,1,1) = two * rmu * g1yi(:,2) + tmp - u1(:)*u1(:)*rho(:)
         rGNa(:,1,2) = tmp1                        - u1(:)*u2(:)*rho(:)
         rGNa(:,1,3) = tmp3                        - u1(:)*u3(:)*rho(:)
         rGNa(:,2,1) = tmp1                        - u1(:)*u2(:)*rho(:)
         rGNa(:,2,2) = two * rmu * g2yi(:,3) + tmp - u2(:)*u2(:)*rho(:)
         rGNa(:,2,3) = tmp2                        - u3(:)*u2(:)*rho(:)
         rGNa(:,3,1) = tmp3                        - u1(:)*u3(:)*rho(:)
         rGNa(:,3,2) = tmp2                        - u3(:)*u2(:)*rho(:)
         rGNa(:,3,3) = two * rmu * g3yi(:,4) + tmp - u3(:)*u3(:)*rho(:)
      endif
      if((iLES.gt.10).and.(iLES.lt.20)) then    ! bard
         rGNa(:,1,1) = rGNa(:,1,1) - rlsli(:,1)*rho(:)
         rGNa(:,1,2) = rGNa(:,1,2) - rlsli(:,4)*rho(:)
         rGNa(:,1,3) = rGNa(:,1,3) - rlsli(:,5)*rho(:)
         rGNa(:,2,1) = rGNa(:,2,1) - rlsli(:,4)*rho(:)
         rGNa(:,2,2) = rGNa(:,2,2) - rlsli(:,2)*rho(:)
         rGNa(:,2,3) = rGNa(:,2,3) - rlsli(:,6)*rho(:)
         rGNa(:,3,1) = rGNa(:,3,1) - rlsli(:,5)*rho(:)
         rGNa(:,3,2) = rGNa(:,3,2) - rlsli(:,6)*rho(:)
         rGNa(:,3,3) = rGNa(:,3,3) - rlsli(:,3)*rho(:)
      endif

      tmp1        = tauM * rLui(:,1)
      tmp2        = tauM * rLui(:,2)
      tmp3        = tauM * rLui(:,3)

      rGNa(:,1,1) = rGNa(:,1,1) + tmp1 * u1
      rGNa(:,1,2) = rGNa(:,1,2) + tmp1 * u2
      rGNa(:,1,3) = rGNa(:,1,3) + tmp1 * u3
      rGNa(:,2,1) = rGNa(:,2,1) + tmp2 * u1
      rGNa(:,2,2) = rGNa(:,2,2) + tmp2 * u2
      rGNa(:,2,3) = rGNa(:,2,3) + tmp2 * u3
      rGNa(:,3,1) = rGNa(:,3,1) + tmp3 * u1
      rGNa(:,3,2) = rGNa(:,3,2) + tmp3 * u2
      rGNa(:,3,3) = rGNa(:,3,3) + tmp3 * u3

      if(iconvflow.eq.1) then
c
c... get the u_j w_{i,i} term in there to match A_j^T w_{i,j} tau L_i
c    to match the SUPG of incompressible limit
c
         rGNa(:,1,1) = rGNa(:,1,1) + tmp1 * u1
         rGNa(:,1,2) = rGNa(:,1,2) + tmp2 * u1
         rGNa(:,1,3) = rGNa(:,1,3) + tmp3 * u1
         rGNa(:,2,1) = rGNa(:,2,1) + tmp1 * u2
         rGNa(:,2,2) = rGNa(:,2,2) + tmp2 * u2
         rGNa(:,2,3) = rGNa(:,2,3) + tmp3 * u2
         rGNa(:,3,1) = rGNa(:,3,1) + tmp1 * u3
         rGNa(:,3,2) = rGNa(:,3,2) + tmp2 * u3
         rGNa(:,3,3) = rGNa(:,3,3) + tmp3 * u3
      endif

      if(iconvflow.eq.2) then  ! advective form has a taubar term to restore con
         tmp1 = tauBar
     &     * ( rLui(:,1) * g1yi(:,2)
     &       + rLui(:,2) * g2yi(:,2)
     &       + rLui(:,3) * g3yi(:,2) )
         tmp2 = tauBar
     &     * ( rLui(:,1) * g1yi(:,3)
     &       + rLui(:,2) * g2yi(:,3)
     &       + rLui(:,3) * g3yi(:,3) )
         tmp3 = tauBar
     &     * ( rLui(:,1) * g1yi(:,4)
     &       + rLui(:,2) * g2yi(:,4)
     &       + rLui(:,3) * g3yi(:,4) )

         rGNa(:,1,1) = rGNa(:,1,1) + tmp1 * rLui(:,1)
         rGNa(:,1,2) = rGNa(:,1,2) + tmp1 * rLui(:,2)
         rGNa(:,1,3) = rGNa(:,1,3) + tmp1 * rLui(:,3)
         rGNa(:,2,1) = rGNa(:,2,1) + tmp2 * rLui(:,1)
         rGNa(:,2,2) = rGNa(:,2,2) + tmp2 * rLui(:,2)
         rGNa(:,2,3) = rGNa(:,2,3) + tmp2 * rLui(:,3)
         rGNa(:,3,1) = rGNa(:,3,1) + tmp3 * rLui(:,1)
         rGNa(:,3,2) = rGNa(:,3,2) + tmp3 * rLui(:,2)
         rGNa(:,3,3) = rGNa(:,3,3) + tmp3 * rLui(:,3)
      endif   ! end of advective form
c
c.... everything that gets multiplied by rNa was supposed
c     to have density multiplying it.  Do it now.

      rNa(:,1) = rNa(:,1) * rho
      rNa(:,2) = rNa(:,2) * rho
      rNa(:,3) = rNa(:,3) * rho

c
c
c.... multiply the residual pieces by the weight space
c
      do aa = 1,nshl
c
c.... continuity
c
         rl(:,aa,4) = rl(:,aa,4) + WdetJ
     &              * ( shg(:,aa,1) * uBar(:,1)
     &                + shg(:,aa,2) * uBar(:,2)
     &                + shg(:,aa,3) * uBar(:,3) )
c
c.... momentum
c
         rl(:,aa,1) = rl(:,aa,1) - WdetJ
     &              * ( shpfun(:,aa) * rNa(:,1)
     &                + shg(:,aa,1) * rGNa(:,1,1)
     &                + shg(:,aa,2) * rGNa(:,1,2)
     &                + shg(:,aa,3) * rGNa(:,1,3) )
         rl(:,aa,2) = rl(:,aa,2) - WdetJ
     &              * ( shpfun(:,aa) * rNa(:,2)
     &                + shg(:,aa,1) * rGNa(:,2,1)
     &                + shg(:,aa,2) * rGNa(:,2,2)
     &                + shg(:,aa,3) * rGNa(:,2,3) )
         rl(:,aa,3) = rl(:,aa,3) - WdetJ
     &              * ( shpfun(:,aa) * rNa(:,3)
     &                + shg(:,aa,1) * rGNa(:,3,1)
     &                + shg(:,aa,2) * rGNa(:,3,2)
     &                + shg(:,aa,3) * rGNa(:,3,3) )

      enddo
c
c.... return
c
      return
      end



c------------------------------------------------------------------------
c
c     calculate the residual for the advection-diffusion equation
c
c------------------------------------------------------------------------
      subroutine e3ResSclr ( uMod,              gGradS,
     &                       Sclr,		Sdot,	gradS,
     &                       WdetJ,		rLS,	tauS,
     &                       shpfun,            shg,    src,
     &                       diffus,
     &                       rl )
c
      include "common.h"

      real*8    uMod(npro,nsd),   gGradS(npro, nsd),
     &          Sclr(npro),       Sdot(npro),	gradS(npro,nsd),
     &          WdetJ(npro),      rLS(npro),	rho(npro),
     &          tauS(npro),       shpfun(npro,nshl), src(npro),
     &          shg(npro,nshl,3), rl(npro,nshl)

      real*8    diffus(npro)
c
c.... local declarations
c
      real*8    rGNa(npro,nsd),   rNa(npro),  rcp(npro), tmp(npro)

      integer   aa
c
c.... initialize multipliers for Na and Na_{,i}
c
      rNa  = zero
      rGNa = zero
c
c.... Na multiplier
c
      tmps     = one-flmpr  ! consistant mass factor
      rcp = one ! rho * cp


         rNa = rcp*(tmps*Sdot + uMod(:,1) * gradS(:,1)
     &                         + uMod(:,2) * gradS(:,2)
     &                         + uMod(:,3) * gradS(:,3) )
     &        - src



      tmp = rcp * tauS * (rLS -src)
c
c.... Na,i multiplier
c
      rGNa(:,1) = diffus * gradS(:,1) + uMod(:,1) * tmp
      rGNa(:,2) = diffus * gradS(:,2) + uMod(:,2) * tmp
      rGNa(:,3) = diffus * gradS(:,3) + uMod(:,3) * tmp
c
      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
c
c.... add the contribution of DC to residual
c
            rGNa(:,1) = rGNa(:,1) + gGradS(:,1) ! gGradS is
            rGNa(:,2) = rGNa(:,2) + gGradS(:,2) ! g^{ij}*Y_{j}*dcFct
            rGNa(:,3) = rGNa(:,3) + gGradS(:,3) ! calculated in e3dc.f
c
         endif
      endif                     ! end of idcsclr
c
c.... multiply the residual pieces by the weight space
c
      do aa = 1,nshl
c
         rl(:,aa) = rl(:,aa)	- WdetJ
     &                        * ( shpfun(:,aa) * rNa(:)
     &                        + shg(:,aa,1) * rGNa(:,1)
     &                        + shg(:,aa,2) * rGNa(:,2)
     &                        + shg(:,aa,3) * rGNa(:,3) )

      enddo
c
c.... return
c
      return
      end



c----------------------------------------------------------------------
c     Calculate the strong PDE residual.
c----------------------------------------------------------------------
      subroutine e3resStrongPDE(
     &     aci,  u1,   u2,   u3,   Temp, rho,  xx,
     &           g1yi, g2yi, g3yi,
     &     rLui, src, divqi)
      include "common.h"
c     INPUTS
      double precision, intent(in), dimension(npro,nsd) ::
     &     aci, xx
      double precision, intent(in), dimension(npro,nflow) ::
     &     g1yi, g2yi, g3yi
      double precision, intent(in), dimension(npro) ::
     &     u1, u2, u3, Temp, rho
c     OUTPUTS
      double precision, intent(out), dimension(npro,nsd) ::
     &     rLui, src
c     LOCALS
      double precision, dimension(npro) ::
     &     divu
      double precision, dimension(npro,nsd) ::
     &     divqi
      double precision, dimension(nsd) ::
     &     omega
c.... compute source term
      src = zero
      if(matflg(5,1) .ge. 1) then
c        body force contribution to src
         bfx      = datmat(1,5,1) ! Boussinesq, g*alfap
         bfy      = datmat(2,5,1)
         bfz      = datmat(3,5,1)
         select case ( matflg(5,1) )
         case ( 1 )             ! standard linear body force
            src(:,1) = bfx
            src(:,2) = bfy
            src(:,3) = bfz
         case ( 2 )             ! boussinesq body force
            Tref = datmat(2,2,1)
            src(:,1) = bfx * (Temp(:)-Tref)
            src(:,2) = bfy * (Temp(:)-Tref)
            src(:,3) = bfz * (Temp(:)-Tref)
         case ( 3 )             ! user specified f(x,y,z)
            call e3source(xx, src)
         end select
      endif
c
      if(matflg(6,1).eq.1) then
c        coriolis force contribution to src
         omega(1)=datmat(1,6,1)
         omega(2)=datmat(2,6,1)
         omega(3)=datmat(3,6,1)
c  note that we calculate f as if it contains the usual source
c  plus the Coriolis and the centrifugal forces taken to the rhs (sign change)
c  as long as we are doing SUPG with no accounting for these terms in the
c  LHS this is the only change (which will find its way to the RHS momentum
c  equation (both Galerkin and SUPG parts)).
c
c  uncomment later if you want rotation always about z axis
c                 orig_src - om x om x r       - two om x u
c
c$$$          src(:,1)=src(:,1)+omega(3)*omega(3)*xx(:,1)+two*omega(3)*u2
c$$$          src(:,2)=src(:,2)+omega(3)*omega(3)*xx(:,2)-two*omega(3)*u1
c
c more general for testing
c
         src(:,1)=src(:,1)
     &        -omega(2)*(omega(1)*xx(:,2)-omega(2)*xx(:,1))
     &        -omega(3)*(omega(1)*xx(:,3)-omega(3)*xx(:,1))
     &        -two*(omega(2)*u3-omega(3)*u2)
         src(:,2)=src(:,2)
     &        -omega(1)*(omega(2)*xx(:,1)-omega(1)*xx(:,2))
     &        -omega(3)*(omega(2)*xx(:,3)-omega(3)*xx(:,2))
     &        -two*(omega(3)*u1-omega(1)*u3)
         src(:,3)=src(:,3)
     &        -omega(1)*(omega(3)*xx(:,1)-omega(1)*xx(:,3))
     &        -omega(2)*(omega(3)*xx(:,2)-omega(2)*xx(:,3))
     &        -two*(omega(1)*u2-omega(2)*u1)
      endif
c     calculate momentum residual
      rLui(:,1) =(aci(:,1) + u1 * g1yi(:,2)
     &     + u2 * g2yi(:,2)
     &     + u3 * g3yi(:,2) - src(:,1) ) * rho
     &     + g1yi(:,1)
     &        - divqi(:,1)
      rLui(:,2) =(aci(:,2) + u1 * g1yi(:,3)
     &     + u2 * g2yi(:,3)
     &     + u3 * g3yi(:,3) - src(:,2) ) * rho
     &     + g2yi(:,1)
     &        - divqi(:,2)
      rLui(:,3) =(aci(:,3) + u1 * g1yi(:,4)
     &     + u2 * g2yi(:,4)
     &     + u3 * g3yi(:,4) - src(:,3) ) * rho
     &     + g3yi(:,1)
     &        - divqi(:,3)
      if(iconvflow.eq.1) then
         divu(:)  = (g1yi(:,2) + g2yi(:,3) + g3yi(:,4))*rho
         rLui(:,1)=rlui(:,1)+u1(:)*divu(:)
         rLui(:,2)=rlui(:,2)+u2(:)*divu(:)
         rLui(:,3)=rlui(:,3)+u3(:)*divu(:)
      endif
c
      return
      end subroutine e3resStrongPDE