Esope source files

   1 : C GMAT2D    SOURCE    CHAT      05/01/13    00:18:34     5004           
   2 :        SUBROUTINE GMAT2D(GMT,WVECT,GA,V_INF,DIAM)
   3 : C
   4 : C      Voir kon171.eso
   5 : C
   6 : c-----------------------------------------------------
   7 : c  EQUATIONS:  2D Euler Equations of Gas Dynamics
   8 : c
   9 : c  REFERENCE:  AIAA Journal, Sept. 1998
  10 : c              "Low-Diffusion Flux-Splitting Methods 
  11 : c               for Flows at All Speeds" 
  12 : c-----------------------------------------------------
  13 : c  This subroutine provides the preconditioning matrix
  14 : c  which is the product of the two matrices:
  15 : c     a)  the 'original' preconditioning matrix
  16 : c          given in the Reference and
  17 : c     b)  the jacobian matrix which is
  18 : c  
  19 : c                  dW/dU
  20 : c        with W=(p,u,v,T) (T-temperature) and
  21 : c             U=(rho,rho*u,rho*v,rho*e).
  22 : c-----------------------------------------------------
  23 : c  INPUT:
  24 : c         
  25 : c        wvect   --  vector of the primitive variables
  26 : c                    (rho, u, v, p);
  27 : c
  28 : c        ga      --  ratio of the specific heats;
  29 : c
  30 : c        v_inf   --  reference velocity
  31 : c
  32 : c        diam    --  diameter of the cell
  33 : c
  34 : c  OUTPUT:
  35 : c     
  36 : c        gmt     -- preconditioning matrix 4 by 4.
  37 : c-----------------------------------------------------
  38 :       IMPLICIT INTEGER(I-N)
  39 :        real*8 wvect(4),gmt(4,4)
  40 :        real*8 ga,gm1,a
  41 :        real*8 rho,p,u,v, v_inf
  42 :        real*8 rho_p,ur,theta
  43 :        real*8 epsil,h,qq,m,t1,t2
  44 :        real*8 diam,lambda,mref,mref2,mach,mach2,dtaum1
  45 :        parameter(epsil=1.0d0) 
  46 : c--------------------------------------------------
  47 : c  Input
  48 : c--------------------------------------------------
  49 :        rho=wvect(1)
  50 :        u=wvect(2)
  51 :        v=wvect(3)
  52 :        p=wvect(4)
  53 : c--------------------------------------------------
  54 :        gm1=ga-1.0d0
  55 :        a=sqrt(ga*p/rho)
  56 : c--------------------------------------------------
  57 :        rho_p=ga/(a*a)     
  58 :        h=0.5d0*(u*u+v*v)+ga*p/(gm1*rho)
  59 : c--------------------------------------------------
  60 : c  Reference velocity
  61 : c--------------------------------------------------
  62 :        qq=sqrt(u*u+v*v)
  63 :        if(qq .lt. (epsil*v_inf)) then
  64 :           ur = epsil*v_inf
  65 :        else
  66 :           ur = qq
  67 :        endif
  68 :        if(ur .gt. a)then
  69 :           ur=a
  70 :        endif
  71 : c----------------------------------------------------
  72 :        theta=1.0d0/(ur*ur) - (1.0d0-ga)*rho_p/ga
  73 : c----------------------------------------------------       
  74 :        m=(1.0d0-ga)*(theta-rho_p)
  75 :        t1=1.0d0+m*(h-qq*qq)+(theta-rho_p)*ga/rho_p
  76 :        t2=h*(theta-rho_p)*((1.0d0-ga)*(h-qq*qq)+ga/rho_p)
  77 : c----- We compute u+a------------
  78 :        mref=ur/a
  79 :        mref2=mref*mref
  80 :        mach=qq/a
  81 :        mach2=mach*mach
  82 :        lambda=(1.0D0-mref2)**2
  83 :        lambda=(lambda*mach2)+(4.0D0*mref2)
  84 :        lambda=sqrt(lambda)*a
  85 :        lambda=((1.0D0+mref2)*qq)+lambda
  86 :        lambda=lambda*0.5D0
  87 :        dtaum1=lambda/diam
  88 : c--------------------------------
  89 :        gmt(1,1)=t1 * dtaum1
  90 :        gmt(1,2)=u*m * dtaum1
  91 :        gmt(1,3)=v*m * dtaum1
  92 :        gmt(1,4)=-1.0D0 * m * dtaum1
  93 : c-------------------------------
  94 :        gmt(2,1)=u*(t1-1.0d0) * dtaum1
  95 :        gmt(2,2)=(1.0d0+u*u*m) * dtaum1
  96 :        gmt(2,3)=u*v*m * dtaum1
  97 :        gmt(2,4)=-1.0D0 * m*u * dtaum1
  98 : c-------------------------------
  99 :        gmt(3,1)=v*(t1-1.0d0) * dtaum1
 100 :        gmt(3,2)=u*v*m * dtaum1
 101 :        gmt(3,3)=(1.0d0+v*v*m) * dtaum1
 102 :        gmt(3,4)=-1.0D0 * m*v * dtaum1
 103 : c-------------------------------
 104 :        gmt(4,1)=t2 * dtaum1
 105 :        gmt(4,2)=u*h*m * dtaum1
 106 :        gmt(4,3)=v*h*m * dtaum1
 107 :        gmt(4,4)=(-m*h+1.0d0) * dtaum1
 108 : c-------------------------------
 109 :        return
 110 :        end
 111 : 
 112 : 
 113 : 
 114 : 
 115 : 
 116 : 
 117 :  
 118 :