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

        subroutine e3juel (yl,     acl,    sls,    A0,
     &			   WdetJ,  rl,     rml)
c
c----------------------------------------------------------------------
c
c This routine calculates Exactly integrated Linear Tetrahedra
c Mass term (assuming U(Y) is linear it is not).
c
c input:
c  WdetJ  (npro)                : weighted Jacobian
c
c output:
c  rl     (npro,nshl,nflow)      : residual
c  rml    (npro,nshl,nflow)      : modified residual
c
c
c  note that this routine wipes out yl by putting ul into it
c  and then (in ires=1 case ) it is used again
c
c Kenneth Jansen, Winter 1997, Primitive Variables
c----------------------------------------------------------------------
c
	include "common.h"
c
        dimension yl(npro,nshl,nflow),        acl(npro,nshl,ndof),
     &            WdetJ(npro),               A0(npro,nflow,nflow),
     &            rl(npro,nshl,nflow),        rml(npro,nshl,nflow)
c
        dimension rk(npro),                  rho(npro),
     &            ei(npro),                  tmp(npro),
     &            ub(npro,nflow),             fact(npro),
     &		  fddt(npro)

	ttim(28) = ttim(28) - secs(0.0)

c
c.... --------------------->  Time term   <--------------------
c
c.... add contribution of U to rml
c
c.... compute conservative variables
c
c
c   multiply by exact mass matrix   integral N_aN_b=(I+1)*V/20
c   where 1 is a matrix with every element=1
c
c   note that the wght has 4/3 multiplier so 3/4*20=15
c

        fact=WdetJ/(Qwt(lcsyst,intp)*15.0d0)
         fct1=almi/(gami*alfi)*Dtgl  ! factor for predictor (scalar)
        if(ires.ne.1) then
         fddt=fact*fct1
         do inod=1,nshl
	  rk = pt5 * (yl(:,inod,2)**2 + yl(:,inod,3)**2 + yl(:,inod,4)**2)
c
          ithm = 6
          call getthm (yl(:,inod,1), yl(:,inod,5), sls,
     &                 rk,           rho,          ei,
     &                 tmp,          tmp,          tmp,
     &                 tmp,          tmp,          tmp,
     &                 tmp,          tmp)
c
         yl(:,inod,1) = rho
         yl(:,inod,2) = rho * yl(:,inod,2)
         yl(:,inod,3) = rho * yl(:,inod,3)
         yl(:,inod,4) = rho * yl(:,inod,4)
         yl(:,inod,5) = rho * (ei + rk)
        enddo
	    ub(:,:)=yl(:,1,:)+yl(:,2,:)
     &             +yl(:,3,:)+yl(:,4,:)

c
c  what we have now in yl is the U_b^e
c  we want to get Resm(mass)=M^e_{ab} (U^e_b(Y+epsilon P)
c                                              - U^e_b(Y))*fact/dt,
c  since only the difference between resm's is important we
c  do not have to subtract off the unperturbed U(Y) vector
c  This term is meant to carry through the effect of a perturbation
c  in Y upon dY/dt (through the predictor into the term fact)
c
c

	  do i = 1, nshl
	    rml(:,i,1) = rml(:,i,1) + fddt * (yl(:,i,1)+ub(:,1))
	    rml(:,i,2) = rml(:,i,2) + fddt * (yl(:,i,2)+ub(:,2))
	    rml(:,i,3) = rml(:,i,3) + fddt * (yl(:,i,3)+ub(:,3))
	    rml(:,i,4) = rml(:,i,4) + fddt * (yl(:,i,4)+ub(:,4))
	    rml(:,i,5) = rml(:,i,5) + fddt * (yl(:,i,5)+ub(:,5))
	  enddo
c
          flops = flops + 35*nshl*npro
	endif

c
c.... ires = 2 or 3
c
	if ((ires .eq. 1) .or. (ires .eq. 3)) then

	    ub(:,:)=acl(:,1,:)+acl(:,2,:)
     &             +acl(:,3,:)+acl(:,4,:)

	  do i = 1, nshl
	    yl(:,i,1) = fact*(acl(:,i,1)+ub(:,1))
	    yl(:,i,2) = fact*(acl(:,i,2)+ub(:,2))
	    yl(:,i,3) = fact*(acl(:,i,3)+ub(:,3))
	    yl(:,i,4) = fact*(acl(:,i,4)+ub(:,4))
	    yl(:,i,5) = fact*(acl(:,i,5)+ub(:,5))
	  enddo
c
c  what we have now in yl is the dY_a^e/dt=M^e_{ab} Y_{b,t},  must multiply by
c  A0 to get dU^e_a/dt= dU/dY(centroid) dY_a^e/dt,  take advantage of zeros
c  in A0(Prim) with comments
c
          do i = 1, nshl
            rl(:,i,1) = rl(:,i,1)
     &     + A0(:,1,1)*yl(:,i,1)
c    &     + A0(:,1,2)*yl(:,i,2)
c    &     + A0(:,1,3)*yl(:,i,3)
c    &     + A0(:,1,4)*yl(:,i,4)
     &     + A0(:,1,5)*yl(:,i,5)
c
            rl(:,i,2) = rl(:,i,2)
     &     + A0(:,2,1)*yl(:,i,1)
     &     + A0(:,2,2)*yl(:,i,2)
c    &     + A0(:,2,3)*yl(:,i,3)
c    &     + A0(:,2,4)*yl(:,i,4)
     &     + A0(:,2,5)*yl(:,i,5)
c
            rl(:,i,3) = rl(:,i,3)
     &     + A0(:,3,1)*yl(:,i,1)
c    &     + A0(:,3,2)*yl(:,i,2)
     &     + A0(:,3,3)*yl(:,i,3)
c    &     + A0(:,3,4)*yl(:,i,4)
     &     + A0(:,3,5)*yl(:,i,5)
c
            rl(:,i,4) = rl(:,i,4)
     &     + A0(:,4,1)*yl(:,i,1)
c    &     + A0(:,4,2)*yl(:,i,2)
c    &     + A0(:,4,3)*yl(:,i,3)
     &     + A0(:,4,4)*yl(:,i,4)
     &     + A0(:,4,5)*yl(:,i,5)
c
            rl(:,i,5) = rl(:,i,5)
     &     + A0(:,5,1)*yl(:,i,1)
     &     + A0(:,5,2)*yl(:,i,2)
     &     + A0(:,5,3)*yl(:,i,3)
     &     + A0(:,5,4)*yl(:,i,4)
     &     + A0(:,5,5)*yl(:,i,5)
c
          enddo
c
          flops = flops + 45*nshl*npro
	endif

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

c$$$
c$$$        subroutine e3juelSclr (ycl,     acl,    A0t,
c$$$     &			       WdetJ,  rtl,     rmtl)
c$$$c
c$$$c----------------------------------------------------------------------
c$$$c
c$$$c This routine calculates Exactly integrated Linear Tetrahedra
c$$$c Mass term (assuming U(Y) is linear it is not).
c$$$c
c$$$c input:
c$$$c  WdetJ  (npro)                : weighted Jacobian
c$$$c
c$$$c output:
c$$$c  rtl     (npro,nshl,nflow)      : residual
c$$$c  rmtl    (npro,nshl,nflow)      : modified residual
c$$$c
c$$$c
c$$$c  note that this routine wipes out ycl by putting ul into it
c$$$c  and then (in ires=1 case ) it is used again
c$$$c
c$$$c Kenneth Jansen, Winter 1997, Primitive Variables
c$$$c----------------------------------------------------------------------
c$$$c
c$$$	include "common.h"
c$$$c
c$$$        dimension ycl(npro,nshl,ndof),    acl(npro,nshl,ndof),
c$$$     &            WdetJ(npro),           A0t(npro),
c$$$     &            rtl(npro,nshl),        rmtl(npro,nshl)
c$$$
c$$$c
c$$$        dimension rk(npro),              rho(npro),
c$$$     &            ei(npro),              tmp(npro),
c$$$     &            ubt(npro),             fact(npro),
c$$$     &		  fddt(npro)
c$$$
c$$$	ttim(28) = ttim(28) - tmr()
c$$$
c$$$c
c$$$c.... --------------------->  Time term   <--------------------
c$$$c
c$$$c.... add contribution of U to rml
c$$$c
c$$$c.... compute conservative variables
c$$$c
c$$$c
c$$$c   multiply by exact mass matrix   integral N_aN_b=(I+1)*V/20
c$$$c   where 1 is a matrix with every element=1
c$$$c
c$$$c   note that the wght has 4/3 multiplier so 3/4*20=15
c$$$c
c$$$
c$$$        fact=WdetJ/(Qwt(lcsyst,intp)*15.0d0)
c$$$         fct1=almi/(gami*alfi)*Dtgl  ! factor for predictor (scaler)
c$$$c
c$$$c
c$$$c.... ires = 2 or 3
c$$$c
c$$$
c$$$	if ((ires .eq. 1) .or. (ires .eq. 3)) then
c$$$	    ubt(:)=acl(:,1,id)+acl(:,2,id)
c$$$     &            +acl(:,3,id)+acl(:,4,id)
c$$$	  do i = 1, nshl
c$$$	    ycl(:,i,id) = fact*(acl(:,i,id)+ubt(:))
c$$$
c$$$	  enddo
c$$$c
c$$$c  what we have now in ycl is the dY_a^e/dt=M^e_{ab} Y_{b,t},  must multiply by
c$$$c  A0 to get dU^e_a/dt= dU/dY(centroid) dY_a^e/dt,  take advantage of zeros
c$$$c  in A0(Prim) with comments
c$$$c
c$$$          do i = 1, nshl
c$$$            rtl(:,i) = rtl(:,i)
c$$$     &     + A0t(:)*ycl(:,i,id)
c$$$
c$$$c
c$$$          enddo
c$$$c
c$$$          flops = flops + 45*nenl*npro
c$$$	endif
c$$$
c$$$	ttim(28) = ttim(28) + tmr()
c$$$c
c$$$c.... return
c$$$c
c$$$	return
c$$$	end