Sources Esope | Cast3M

conj3w
C CONJ3W    SOURCE    CB215821  16/04/21    21:15:56     8920
c----------------------------------------------------------------------
c GENERAL DESCRIPTION:
c   This subroutine provides the jacobians which are the derivatives
c   of the numerical flux function defined at the cell interface
c   with respect to the conservative variables of the left and right
c   cells (relative to the cell interface).
c
c EQUATIONS: 3D Euler equations of gas dynamics
c
c
c REFERENCE: JCP, 129, 364-382 (1996)
c            " A Sequel to AUSM: AUSM+ ".
c----------------------------------------------------------------------
c INPUT:
c
c      alpha   -- parameter of the AUSM+ scheme in the Pressure function;
c                 ( -3/4 &lt;= alpha &lt;= 3/16 ) (imposed as a parameter)
c
c      beta    -- parameter of the AUSM+ scheme in the Mach function;
c                 ( -1/16 &lt;= beta &lt;= 1/2 )  (imposed as a parameter)
c
c      wvec_l  -- vector of the primitive variables (rho,ux,uy,uz,p) at the
c                 left cell;
c
c      wvec_r  -- vector of the primitive variables (rho,ux,uy,uz,p) at the
c                 right cell;
c
c      nvect   -- normal vector to the interface (3 components in 3D);
c
c      tvec1   -- tangential vector to the interface;
c
c      tvec2   -- tangential vector to the interface;
c
c      ga      -- ratio of the specific heats (assumed constant)
c----------------------------------------------------------------------
c
c OUTPUT:
c
c      jtl     -- jakobian matrix 5 by 5 - derivatives of the numerical
c                 flux function with respect to the conservative variables
c                 from the left cell;
c
c      jtr     -- jakobian matrix 5 by 5 - derivatives of the numerical
c                 flux function with respect to the conservative variables
c                 from the right cell.
c----------------------------------------------------------------------
       SUBROUTINE CONJ3W(JTL,JTR,WVEC_L,WVEC_R,NVECT,TVECT1,TVECT2,
     &                   ga)
      IMPLICIT INTEGER(I-N)
       real*8 wvec_l(5),wvec_r(5)
       real*8 nvect(3),tvect1(3),tvect2(3)
       real*8 jl(5,5),jr(5,5)
       real*8 wl(5,5),wr(5,5)
       real*8 jtl(5,5),jtr(5,5)
       real*8 alpha,beta
       real*8 ga, gm1
       real*8 n_x,n_y,n_z,t1_x,t1_y,t1_z,t2_x,t2_y,t2_z
       real*8 un_l,un_r
       real*8 ml,mr,Mplus,Mmin
       real*8 am
       real*8 rold_l,uold_l,vold_l,wold_l,pold_l,eold_l
       real*8 rold_r,uold_r,vold_r,wold_r,pold_r,eold_r
       real*8 Pplus,Pmin
       real*8 temp_l,temp_r,brac_l,brac_r
       real*8 aleft, arigh
       real*8 damr_l,damr_r,damu_l,damu_r
       real*8 damv_l,damv_r,damp_l,damp_r
       real*8 damw_l,damw_r
       real*8 dmlr_l,dmlr_r,dmlu_l,dmlu_r
       real*8 dmlv_l,dmlv_r,dmlp_l,dmlp_r
       real*8 dmlw_l,dmlw_r
       real*8 dmrr_l,dmrr_r,dmru_l,dmru_r
       real*8 dmrv_l,dmrv_r,dmrp_l,dmrp_r
       real*8 dmrw_l,dmrw_r
       real*8 dPpr_l,dPpr_r,dPpu_l,dPpu_r
       real*8 dPpv_l,dPpv_r,dPpp_l,dPpp_r
       real*8 dPpw_l,dPpw_r
       real*8 dPmr_l,dPmr_r,dPmu_l,dPmu_r
       real*8 dPmv_l,dPmv_r,dPmp_l,dPmp_r
       real*8 dPmw_l,dPmw_r
       real*8 dpir_l,dpir_r,dpiu_l,dpiu_r
       real*8 dpiv_l,dpiv_r,dpip_l,dpip_r
       real*8 dpiw_l,dpiw_r
       real*8 c11,c12,c13,c21,c22,c23,c31,c32,c33,det
       real*8 zc11,zc12,zc13,zc21,zc22,zc23,zc31,zc32,zc33
       integer i,j,k
       parameter(alpha = 0.1875D0, beta = 0.125D0)
c-----------------------
       gm1=ga-1.0d0
c-----------------------
       n_x=nvect(1)
       n_y=nvect(2)
       n_z=nvect(3)
c-----------------------
       t1_x=tvect1(1)
       t1_y=tvect1(2)
       t1_z=tvect1(3)
c-----------------------
       t2_x=tvect2(1)
       t2_y=tvect2(2)
       t2_z=tvect2(3)
c--------------------------------------
       rold_l=wvec_l(1)
       uold_l=wvec_l(2)
       vold_l=wvec_l(3)
       wold_l=wvec_l(4)
       pold_l=wvec_l(5)
c-----------------------
       rold_r=wvec_r(1)
       uold_r=wvec_r(2)
       vold_r=wvec_r(3)
       wold_r=wvec_r(4)
       pold_r=wvec_r(5)
c------------------------------------------------------------------
c Computation of the specific total energy on the left and right.
c------------------------------------------------------------------
       eold_l=(uold_l*uold_l+vold_l*vold_l+wold_l*wold_l)/2.0d0
       eold_l=eold_l+pold_l/(gm1*rold_l)
       eold_r=(uold_r*uold_r+vold_r*vold_r+wold_r*wold_r)/2.0d0
       eold_r=eold_r+pold_r/(gm1*rold_r)
c-------------------------------------------------------------------
c Computation of the speed of sound and its derivatives;
c numerical speed of sound at the interface is taken as an average
c of the speeds of sounds of the neighbouring cells
c---------------------------------------------------------------------
       aleft=SQRT(ga*pold_l/rold_l)
       arigh=SQRT(ga*pold_r/rold_r)
       am=0.5d0*(aleft+arigh)
c--------------------------------------------------------------------
       damr_r=-arigh/(4.0d0*rold_r)
       damu_r=0.0d0
       damv_r=0.0d0
       damw_r=0.0d0
       damp_r=ga/(4.0d0*arigh*rold_r)
c-----------------------
       damr_l=-aleft/(4.0d0*rold_l)
       damu_l=0.0d0
       damv_l=0.0d0
       damw_l=0.0d0
       damp_l=ga/(4.0d0*aleft*rold_l)
c-------------------------------------------------------------------
c Computing numerical Mach number and its derivatives,
c see p.370, under (A1).
c-------------------------------------------------------------------
       un_l=uold_l*n_x+vold_l*n_y+wold_l*n_z
       un_r=uold_r*n_x+vold_r*n_y+wold_r*n_z
c----------------------------------------
       ml=un_l/am
       mr=un_r/am
c-------------------------------------------------------------------
c  Mplus and Mmin are calligraphic lettes M+ and M- from the paper,
c  see (19a) and (19b), p.367.
c-------------------------------------------------------------------
       if(ABS(ml) .ge. 1.0d0) then
         Mplus=(ml+ABS(ml))/2.0d0
       else
         Mplus=(ml+1.0d0)*(ml+1.0d0)/4.0d0
         Mplus=Mplus+beta*(ml*ml-1.0d0)*(ml*ml-1.0d0)
       endif
c-----------------
       if(ABS(mr) .ge. 1.0d0) then
         Mmin=(mr-ABS(mr))/2.0d0
       else
         Mmin=-(mr-1.0d0)*(mr-1.0d0)/4.0d0
         Mmin=Mmin-beta*(mr*mr-1.0d0)*(mr*mr-1.0d0)
       endif
c-------------------------------------------------------------------
c Derivatives of ml and mr with respect to both sets of primitive
c variables: left  and right.
c-------------------------------------------------------------------
       temp_l=-un_l/(am*am)
       temp_r=-un_r/(am*am)
c--------
       dmlr_l=temp_l*damr_l
       dmlr_r=temp_l*damr_r
       dmrr_l=temp_r*damr_l
       dmrr_r=temp_r*damr_r
c--------
       dmlu_l=n_x/am+temp_l*damu_l
       dmlu_r=temp_l*damu_r
       dmru_l=temp_r*damu_l
       dmru_r=n_x/am+temp_r*damu_r
c--------
       dmlv_l=n_y/am+temp_l*damv_l
       dmlv_r=temp_l*damv_r
       dmrv_l=temp_r*damv_l
       dmrv_r=n_y/am+temp_r*damv_r
c--------
       dmlw_l=n_z/am+temp_l*damw_l
       dmlw_r=temp_l*damw_r
       dmrw_l=temp_r*damw_l
       dmrw_r=n_z/am+temp_r*damw_r
c--------
       dmlp_l=temp_l*damp_l
       dmlp_r=temp_l*damp_r
       dmrp_l=temp_r*damp_l
       dmrp_r=temp_r*damp_r
c---------------------------------------------------------------
c Computing the calligraphic P+ and P- with their derivatives
c see (21a) & (21b) on p.368.
c---------------------------------------------------------------
       if(ml .ge. 1.0d0) then
         Pplus = 1.0d0
       else
         if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
           Pplus=(ml+1.0d0)*(ml+1.0d0)*(2.0d0-ml)/4.0d0
           Pplus=Pplus+alpha*ml*(ml*ml-1.0d0)*(ml*ml-1.0d0)
         else
           Pplus = 0.0d0
         endif
       endif
c---------------------------------------------------------------
       if(mr .ge. 1.0d0) then
         Pmin = 0.0d0
       else
         if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
           Pmin=(mr-1.0d0)*(mr-1.0d0)*(2.0d0+mr)/4.0d0
           Pmin=Pmin-alpha*mr*(mr*mr-1.0d0)*(mr*mr-1.0d0)
         else
           Pmin = 1.0d0
         endif
       endif
c---------------------------------------------------------------
       brac_l=(ml+1.0d0)*(2.0d0-ml)/2.0d0-(ml+1.0d0)*(ml+1.0d0)/4.0d0
       brac_l=brac_l+alpha*(ml*ml-1.0d0)*(ml*ml-1.0d0)
       brac_l=brac_l+4.0d0*alpha*ml*ml*(ml*ml-1.0d0)
c--------------
       brac_r=(mr-1.0d0)*(2.0d0+mr)/2.0d0+(mr-1.0d0)*(mr-1.0d0)/4.0d0
       brac_r=brac_r-alpha*(mr*mr-1.0d0)*(mr*mr-1.0d0)
       brac_r=brac_r-4.0d0*alpha*mr*mr*(mr*mr-1.0d0)
c---------------------------------------------------------------
       if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
         dPpr_l=brac_l*dmlr_l
         dPpr_r=brac_l*dmlr_r
c------------
         dPpu_l=brac_l*dmlu_l
         dPpu_r=brac_l*dmlu_r
c------------
         dPpv_l=brac_l*dmlv_l
         dPpv_r=brac_l*dmlv_r
c------------
         dPpw_l=brac_l*dmlw_l
         dPpw_r=brac_l*dmlw_r
c------------
         dPpp_l=brac_l*dmlp_l
         dPpp_r=brac_l*dmlp_r
c------------
       else
         dPpr_l=0.0d0
         dPpr_r=0.0d0
c-----------
         dPpu_l=0.0d0
         dPpu_r=0.0d0
c-----------
         dPpv_l=0.0d0
         dPpv_r=0.0d0
c-----------
         dPpw_l=0.0d0
         dPpw_r=0.0d0
c-----------
         dPpp_l=0.0d0
         dPpp_r=0.0d0
c-----------
       endif
c---------------------------------------------------------------
       if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
         dPmr_l=brac_r*dmrr_l
         dPmr_r=brac_r*dmrr_r
c------------
         dPmu_l=brac_r*dmru_l
         dPmu_r=brac_r*dmru_r
c------------
         dPmv_l=brac_r*dmrv_l
         dPmv_r=brac_r*dmrv_r
c------------
         dPmw_l=brac_r*dmrw_l
         dPmw_r=brac_r*dmrw_r
c------------
         dPmp_l=brac_r*dmrp_l
         dPmp_r=brac_r*dmrp_r
c------------
       else
         dPmr_l=0.0d0
         dPmr_r=0.0d0
c-----------
         dPmu_l=0.0d0
         dPmu_r=0.0d0
c-----------
         dPmv_l=0.0d0
         dPmv_r=0.0d0
c-----------
         dPmw_l=0.0d0
         dPmw_r=0.0d0
c-----------
         dPmp_l=0.0d0
         dPmp_r=0.0d0
c-----------
       endif
c---------------------------------------------------------------------
c computing pmid -- p_{1/2} and its derivatives, see (20b), p.367.
c---------------------------------------------------------------------
       dpir_l=dPpr_l*pold_l+dPmr_l*pold_r
       dpiu_l=dPpu_l*pold_l+dPmu_l*pold_r
       dpiv_l=dPpv_l*pold_l+dPmv_l*pold_r
       dpiw_l=dPpw_l*pold_l+dPmw_l*pold_r
       dpip_l=dPpp_l*pold_l+Pplus+dPmp_l*pold_r
c----------------------------
       dpir_r=dPpr_r*pold_l+dPmr_r*pold_r
       dpiu_r=dPpu_r*pold_l+dPmu_r*pold_r
       dpiv_r=dPpv_r*pold_l+dPmv_r*pold_r
       dpiw_r=dPpw_r*pold_l+dPmw_r*pold_r
       dpip_r=dPpp_r*pold_l+Pmin+dPmp_r*pold_r
c---------------------------------------------------------------------
c computing JACOBIAN as a derivative of the numerical flux function
c with respect to the primitive variables.
c Notation: jl(i,j) --- is the derivative of the i-component of the
c           flux function with respect to the j-component of the
c           vector of primitive variables of the left state.
c           jr(i,j) --- is the derivative of the i-component of the
c           flux function with respect to the j-component of the
c           vector of primitive variables of the right state.
c---------------------------------------------------------------------
       jl(1,1)=0.0D0
       jl(1,2)=0.0D0
       jl(1,3)=0.0D0
       jl(1,4)=0.0D0
       jl(1,5)=0.0D0
c------------------------------------
       jr(1,1)=0.0D0
       jr(1,2)=0.0D0
       jr(1,3)=0.0D0
       jr(1,4)=0.0D0
       jr(1,5)=0.0D0
c------------------------------------
c------------------------------------
       jl(2,1)=dpir_l*n_x
c------------------------
       jl(2,2)=dpiu_l*n_x
c------------------------
       jl(2,3)=dpiv_l*n_x
c------------------------
       jl(2,4)=dpiw_l*n_x
c------------------------
       jl(2,5)=dpip_l*n_x
c-------------------------------------------------
c-------------------------------------------------
       jl(3,1)=dpir_l*n_y
c------------------------
       jl(3,2)=dpiu_l*n_y
c------------------------
       jl(3,3)=dpiv_l*n_y
c------------------------
       jl(3,4)=dpiw_l*n_y
c------------------------
       jl(3,5)=dpip_l*n_y
c-------------------------------------------------
c-------------------------------------------------
       jl(4,1)=dpir_l*n_z
c------------------------
       jl(4,2)=dpiu_l*n_z
c------------------------
       jl(4,3)=dpiv_l*n_z
c------------------------
       jl(4,4)=dpiw_l*n_z
c------------------------
       jl(4,5)=dpip_l*n_z
c-------------------------------------------------------
c derivatives with respect to the right
c set of primitive variables
c-------------------------------------------------------
       jr(2,1)=dpir_r*n_x
c------------------------
       jr(2,2)=dpiu_r*n_x
c------------------------
       jr(2,3)=dpiv_r*n_x
c------------------------
       jr(2,4)=dpiw_r*n_x
c------------------------
       jr(2,5)=dpip_r*n_x
c-------------------------------------------------------
       jr(3,1)=dpir_r*n_y
c------------------------
       jr(3,2)=dpiu_r*n_y
c------------------------
       jr(3,3)=dpiv_r*n_y
c------------------------
       jr(3,4)=dpiw_r*n_y
c------------------------
       jr(3,5)=dpip_r*n_y
c--------------------------------------------------------
       jr(4,1)=dpir_r*n_z
c------------------------
       jr(4,2)=dpiu_r*n_z
c------------------------
       jr(4,3)=dpiv_r*n_z
c------------------------
       jr(4,4)=dpiw_r*n_z
c------------------------
       jr(4,5)=dpip_r*n_z
c-------------------------------------------------------
c-------------------------------------------------------
       jl(5,1)=0.0D0
       jl(5,2)=0.0D0
       jl(5,3)=0.0D0
       jl(5,4)=0.0D0
       jl(5,5)=0.0D0
c-------------------------------------------------
c-------------------------------------------------
       jr(5,1)=0.0D0
       jr(5,2)=0.0D0
       jr(5,3)=0.0D0
       jr(5,4)=0.0D0
       jr(5,5)=0.0D0
c-------------------------------------------------------------
c  matrix wl(i,j) represents the derivative of the i-component
c  of the vector of primitive variables of the left state with
c  respect to the j-component of the vector of the conservative
c  variables of the left state.
c
c  Here: (rho, u, v, w, p) - vector of primitive variables;
c        (rho, rho u, rh o v, rho w, rho e) - conservative variables.
c-------------------------------------------------------------
       wl(1,1)=1.0d0
       wl(1,2)=0.0d0
       wl(1,3)=0.0d0
       wl(1,4)=0.0d0
       wl(1,5)=0.0d0
c------------------------------
       wl(2,1)=-uold_l/rold_l
       wl(2,2)=1.0d0/rold_l
       wl(2,3)=0.0d0
       wl(2,4)=0.0d0
       wl(2,5)=0.0d0
c------------------------------
       wl(3,1)=-vold_l/rold_l
       wl(3,2)=0.0d0
       wl(3,3)=1.0d0/rold_l
       wl(3,4)=0.0d0
       wl(3,5)=0.0d0
c------------------------------
       wl(4,1)=-wold_l/rold_l
       wl(4,2)=0.0d0
       wl(4,3)=0.0d0
       wl(4,4)=1.0d0/rold_l
       wl(4,5)=0.0d0
c------------------------------
       wl(5,1)=gm1*(uold_l*uold_l+vold_l*vold_l+wold_l*wold_l)/2.0d0
       wl(5,2)=-uold_l*gm1
       wl(5,3)=-vold_l*gm1
       wl(5,4)=-wold_l*gm1
       wl(5,5)=gm1
c------------------------------
c------------------------------
       wr(1,1)=1.0d0
       wr(1,2)=0.0d0
       wr(1,3)=0.0d0
       wr(1,4)=0.0d0
       wr(1,5)=0.0d0
c------------------------------
       wr(2,1)=-uold_r/rold_r
       wr(2,2)=1.0d0/rold_r
       wr(2,3)=0.0d0
       wr(2,4)=0.0d0
       wr(2,5)=0.0d0
c------------------------------
       wr(3,1)=-vold_r/rold_r
       wr(3,2)=0.0d0
       wr(3,3)=1.0d0/rold_r
       wr(3,4)=0.0d0
       wr(3,5)=0.0d0
c------------------------------
       wr(4,1)=-wold_r/rold_r
       wr(4,2)=0.0d0
       wr(4,3)=0.0d0
       wr(4,4)=1.0d0/rold_r
       wr(4,5)=0.0d0
c------------------------------
       wr(5,1)=gm1*(uold_r*uold_r+vold_r*vold_r+wold_r*wold_r)/2.0d0
       wr(5,2)=-uold_r*gm1
       wr(5,3)=-vold_r*gm1
       wr(5,4)=-wold_r*gm1
       wr(5,5)=gm1
c----------------------------------------------
c----------------------------------------------
       do 1 i=1,5
        do 2 j=1,5
         jtl(i,j)=0.0d0
         jtr(i,j)=0.0d0
         do 3 k=1,5
          jtl(i,j)=jtl(i,j)+jl(i,k)*wl(k,j)
          jtr(i,j)=jtr(i,j)+jr(i,k)*wr(k,j)
 3       continue
 2      continue
 1     continue
c-----------------------------------------------
c Taking in account the dependancy of variables
c-----------------------------------------------
       c11=t1_y*t2_z - t1_z*t2_y
       c12=n_y*t2_z - n_z*t2_y
       c13=n_y*t1_z - n_z*t1_y
c-------------------------------------
       c21=t1_x*t2_z - t1_z*t2_x
       c22=n_x*t2_z - n_z*t2_x
       c23=n_x*t1_z - n_z*t1_x
c-------------------------------------
       c31=t1_x*t2_y - t1_y*t2_x
       c32=n_x*t2_y - n_y*t2_x
       c33=n_x*t1_y - n_y*t1_x
       det=n_x*c11 - n_y*c21 + n_z*c31
c----------------------------------------------------------------------
       ZC11=-NVECT(1)*C11-TVECT1(1)*C12+TVECT2(1)*C13
       ZC12=-NVECT(2)*C11-TVECT1(2)*C12+TVECT2(2)*C13
       ZC13=-NVECT(3)*C11-TVECT1(3)*C12+TVECT2(3)*C13
C---------------------------------
       ZC21=NVECT(1)*C21+TVECT1(1)*C22-TVECT2(1)*C23
       ZC22=NVECT(2)*C21+TVECT1(2)*C22-TVECT2(2)*C23
       ZC23=NVECT(3)*C21+TVECT1(3)*C22-TVECT2(3)*C23
C---------------------------------
       ZC31=-NVECT(1)*C31-TVECT1(1)*C32+TVECT2(1)*C33
       ZC32=-NVECT(2)*C31-TVECT1(2)*C32+TVECT2(2)*C33
       ZC33=-NVECT(3)*C31-TVECT1(3)*C32+TVECT2(3)*C33
c---------------------------------------------------------------------
       do 11 i=1,5
         jtl(i,1)=jtl(i,1)+jtr(i,1)
         jtl(i,2)=jtl(i,2)+jtr(i,2)*zc11/det+jtr(i,3)*zc21/det+
     &            jtr(i,4)*zc31/det
         jtl(i,3)=jtl(i,3)+jtr(i,2)*zc12/det+jtr(i,3)*zc22/det+
     &            jtr(i,4)*zc32/det
         jtl(i,4)=jtl(i,4)+jtr(i,2)*zc13/det+jtr(i,3)*zc23/det+
     &            jtr(i,4)*zc33/det
         jtl(i,5)=jtl(i,5)+jtr(i,5)
 11    continue
c----------------------------------------------------------------------
       return
       end