Sources Esope | Cast3M

ausmw
C AUSMW     SOURCE    CB215821  16/04/21    21:15:10     8920
c----------------------------------------------------------------------
c INPUT:
c
c      rog : left density
C
c      uxg : left x-component of speed
c
c      uyg : left y-component of speed
c
c      pg  : left pressure
c
c      rod : right density
c
c      uxd : right x-component of speed
c
c      uyd : right y-component of speed
c
c      pd  : right pressure
c
c      ga  : gamma of the gas
c
c      cnx : x-component of the surface normal
c
c      cny : y-component of the surface normal
c
c OUTPUT:
c
c      dfrun : derivative of the numerical estimation of (rho u_n u_n +
c              p) (flux-function for conservative variable (rho u_n))
c              with respect to the primitive variables
c
c----------------------------------------------------------------------
       SUBROUTINE AUSMW(ROG,UXG,UYG,PG,ROD,UXD,UYD,PD,
     &                  ga,cnx,cny,dfrun)
      IMPLICIT INTEGER(I-N)
       real*8 rog,uxg,uyg,pg
       real*8 rod,uxd,uyd,pd
       real*8 cnx,cny
       real*8 jl(4,4),jr(4,4),dfrun(4)
       real*8 alpha,beta
       real*8 ga, gm1
       real*8 n_x,n_y,t_x,t_y
       real*8 un_l, un_r
       real*8 ml,mr,Mplus,Mmin
       real*8 rold_l,uold_l,vold_l,pold_l,eold_l
       real*8 rold_r,uold_r,vold_r,pold_r,eold_r
       real*8 Pplus,Pmin
       real*8 temp_l,temp_r,brac_l,brac_r
       real*8 aleft, arigh, am
       real*8 damr_l,damr_r,damu_l,damu_r
       real*8 damv_l,damv_r,damp_l,damp_r
       real*8 dmlr_l,dmlr_r,dmlu_l,dmlu_r
       real*8 dmlv_l,dmlv_r,dmlp_l,dmlp_r
       real*8 dmrr_l,dmrr_r,dmru_l,dmru_r
       real*8 dmrv_l,dmrv_r,dmrp_l,dmrp_r
       real*8 dPpr_l,dPpr_r,dPpu_l,dPpu_r
       real*8 dPpv_l,dPpv_r,dPpp_l,dPpp_r
       real*8 dPmr_l,dPmr_r,dPmu_l,dPmu_r
       real*8 dPmv_l,dPmv_r,dPmp_l,dPmp_r
       real*8 dpir_l,dpir_r,dpiu_l,dpiu_r
       real*8 dpiv_l,dpiv_r,dpip_l,dpip_r
       real*8 tcoef,bcoef
       integer i
       parameter(alpha = 0.1875D0, beta = 0.125D0)
c-----------------------
       gm1=ga-1.0d0
c-----------------------
       n_x=cnx
       n_y=cny
       t_x=-cny
       t_y=cnx
c----------------------------
c-----------------------
       rold_l=rog
       uold_l=uxg
       vold_l=uyg
       pold_l=pg
c-----------------------
       rold_r=rod
       uold_r=uxd
       vold_r=uyd
       pold_r=pd
c-----------------------
       eold_l=(uold_l*uold_l+vold_l*vold_l)/2.0d0
       eold_l=eold_l+pold_l/(gm1*rold_l)
       eold_r=(uold_r*uold_r+vold_r*vold_r)/2.0d0
       eold_r=eold_r+pold_r/(gm1*rold_r)
c-------------------------------------------------------------------
c Computation of the speed of sound and its derivatives
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
         damp_r=ga/(4.0d0*arigh*rold_r)
c-----------------------
         damr_l=-aleft/(4.0d0*rold_l)
         damu_l=0.0d0
         damv_l=0.0d0
         damp_l=ga/(4.0d0*aleft*rold_l)
c----------------------------------------------------------------------
c  computing numerical Mach number and its derivatives
c----------------------------------------------------------------------
       un_l=uold_l*n_x+vold_l*n_y
       un_r=uold_r*n_x+vold_r*n_y
c----------------------------------------
       ml=un_l/am
       mr=un_r/am
c-------------------------------------
c  Mplus and Mmin are calligraphic
c  lettes M+ and M- from the paper
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
c to both sets of primitive variables:
c 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--------
       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---------------------------------------------------------------
       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------------
         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-----------
         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------------
         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-----------
         dPmp_l=0.0d0
         dPmp_r=0.0d0
c-----------
       endif
c---------------------------------------------------------------------
c computing pmid -- p_{1/2} and its derivatives
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
       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
       dpip_r=dPpp_r*pold_l+Pmin+dPmp_r*pold_r
c---------------------------------------------------------------------
c  !!!!!!!!!!!!!!!!!!  JACOBIAN  !!!!!!!!!!!!!!!!!!!!!!!!!!!!
c---------------------------------------------------------------------
       jl(1,1)=0.0D0
       jl(1,2)=0.0D0
       jl(1,3)=0.0D0
       jl(1,4)=0.0D0
c------------------------------------
       jr(1,1)=0.0D0
       jr(1,2)=0.0D0
       jr(1,3)=0.0D0
       jr(1,4)=0.0D0
c------------------------------------
ccccc    f222222222222  -------------
c------------------------------------
c---------------------------------------------------------
c--------------------------------------------
       jl(2,1)=n_x*dpir_l
c-------------------
       jl(2,2)=n_x*dpiu_l
c-------------------
       jl(2,3)=n_x*dpiv_l
c-------------------
       jl(2,4)=n_x*dpip_l
c-------------------------------------------------------------
       jr(2,1)=n_x*dpir_r
c-------------------
       jr(2,2)=n_x*dpiu_r
c-------------------
       jr(2,3)=n_x*dpiv_r
c-------------------
       jr(2,4)=n_x*dpip_r
c-------------------------------------------------------------
c------------    f33333333333333333333   ---------------------
c-------------------------------------------------------------
       jl(3,1)=n_y*dpir_l
c-------------------
       jl(3,2)=n_y*dpiu_l
c-------------------
       jl(3,3)=n_y*dpiv_l
c-------------------
       jl(3,4)=n_y*dpip_l
c-------------------------------------------------------------
       jr(3,1)=n_y*dpir_r
c-------------------
       jr(3,2)=n_y*dpiu_r
c-------------------
       jr(3,3)=n_y*dpiv_r
c-------------------
       jr(3,4)=n_y*dpip_r
c-------------------------------------------------------------
c  ------  f44444444444444444444444444444 ---------
c-------------------------------------------------------------
       jl(4,1)=0.0D0
c---------------------
       jl(4,2)=0.0D0
c---------------------
       jl(4,3)=0.0D0
c---------------------
       jl(4,4)=0.0D0
c----------------------------------------------------------
c----------------------------------------------------------
       jr(4,1)=0.0D0
c---------------------
       jr(4,2)=0.0D0
c---------------------
       jr(4,3)=0.0D0
c---------------------
       jr(4,4)=0.0D0
c----------------------------------------------------------
       tcoef=n_x*t_y+t_x*n_y
       bcoef=n_x*t_y-t_x*n_y
c-----------------------------
c-----------------------------------
       do 11 i=1,4
         jl(i,1)=jl(i,1)+jr(i,1)
         jl(i,2)=jl(i,2)+jr(i,2)*(-tcoef/bcoef)+
     &            jr(i,3)*2.0d0*n_x*t_x/bcoef
         jl(i,3)=jl(i,3)+jr(i,2)*(-2.0d0*n_y*t_y/bcoef)+
     &            jr(i,3)*tcoef/bcoef
         jl(i,4)=jl(i,4)+jr(i,4)
 11    continue
c--------------------------------------
       do 12 i=1,4
        dfrun(i)=jl(2,i)*n_x+jl(3,i)*n_y
 12    continue
c----------------------------------------------------------------------
       return
       end