Sources Esope | Cast3M

flusti
C FLUSTI    SOURCE    BECC      11/05/26    21:15:25     6981
      SUBROUTINE FLUSTI(
     &     PC_L, GAM_L,
     &     PC_R, GAM_R,
     &     Z,
     &     D_L, D_B, D_A, D_R,
     &     RHO_L, P_L, U_L, UT1_L, UT2_L,
     &     RHO_B, P_B, U_B,
     &     RHO_A, P_A, U_A,
     &     RHO_R, P_R, U_R, UT1_R, UT2_R,
     &     UINT,
     &     FLURHO, FLURU, FLURT1, FLURT2, FLURET,
     &     LOGDEB, LOGAN)
*
*************************************************************************
*
* project           :  CAST3M, EUROPLEXUS...
*
* name              :  flusti
*
* description       :  euler equations for a mixture of stiffened gases
*                      flux in the non-reactive case.
*
*                      given z = x/t, we compute the flux
*                      f(x/t) - ((x/t) * u)
*                      we compute density, pressure and velocity by
*                      averaging the intermediate states using the
*                      weve speeds.
*                      we compute the other variables using these
*                      quantities and the eos
*
* language          :  fortran 77
*
* author            :  a. beccantini den/dm2s/sfme/ltmf
*
*************************************************************************
*
* called by         :
*
*
*************************************************************************
*
***** input
*
*     pc_l, gam_l     =  properties of the gas in the left
*
*     pc_r, gam_r     =  properties of the gas in the right
*
*     z               =  x / t
*
*     d_l, d_b, d_a, d_r =
*                     = wave speeds
*
*     rho_, p_, u_
*                     = density, pressure, velocity
*
*     uint            = interface speed
*
*     logdeb          = debugging ?
*
***** output
*
*     flurho, flurhou, flurt1, flurt2, fluret =
*                      fluxes for rho, rho u, rho ut1, rhout_2, rho e_t
*
*     logan          = anomaly
*
*************************************************************************
*
*     26/11/2009      created
*     25/05/2011      evolution in CAST3M
*
*************************************************************************
*
* n.b.: all variables are declared
*
C      IMPLICIT NONE
      IMPLICIT INTEGER(I-N)
      REAL*8
     &     PC_L, GAM_L,
     &     PC_R, GAM_R,
     &     Z,
     &     D_L, D_B, D_A, D_R,
     &     RHO_L, P_L, U_L, UT1_L, UT2_L,
     &     RHO_B, P_B, U_B,
     &     RHO_A, P_A, U_A,
     &     RHO_R, P_R, U_R, UT1_R, UT2_R,
     &     UINT,
     &     FLURHO, FLURU, FLURT1, FLURT2, FLURET
*
      REAL*8 FP_L, FP_B, FP_A, FP_R
     &     , GAM_Z, PC_Z
     &     , RHO_Z, P_Z, U_Z, E_Z, UT1_Z, UT2_Z
     &     , ECIN_Z
*
*     debugging ?
*
      LOGICAL LOGDEB, LOGAN
*
******** states
*
*         l
*        ***
*          |*                              r
*          | *   b                       *********
*          |  ********                  *|
*          |  |      |  *   a          * |
*          |  |      |    *************  |
*          |  |      |    |           |  |
*          |  |      |    |           |  |
*        d_l d_b   un_b   un_a       d_a d_r
*
*c
*c     test
*c
*      z = 1.13d0
*      d_l = -1.0d0
*      d_b = 1.10d0
*      d_a = 1.11d0
*      d_r = 1.12d0
*
*     speed of the contact discontinuty u_ab = u_a = u_b
*     if there is no vacuum between a and b
*
*
***** weight functions
*
      FP_L = 0.0D0
      FP_B = 0.0D0
      FP_A = 0.0D0
      FP_R = 0.0D0
*
*     this can be simplified from a computational point of view
*     using max, min, abs...
*     but in this way it is easy to read it...
*
      IF (Z .LT. D_L) THEN
         FP_L = 1.0D0
      ELSEIF(Z .LT. D_B) THEN
*        rarefaction if d_l &lt; d_b
         FP_L = (D_B - Z) / (D_B - D_L)
         FP_B = 1.0D0 - FP_L
      ELSEIF(Z .LT. U_B)THEN
*        contact discontinuty, left
         FP_B = 1.0D0
      ELSEIF(Z .LT. U_A)THEN
*        vacuum (u_a &lt; u_b)
         FP_B = (U_A - Z) / (U_A - U_B)
         FP_A = 1.0D0 - FP_B
      ELSEIF(Z .LT. D_A)THEN
*        rs, right
         FP_A = 1.0D0
      ELSEIF(Z .LT. D_R)THEN
*        rgnl rarefaction if d_r &lt; d_a
         FP_A = (D_R - Z) / (D_R - D_A)
         FP_R = 1.0D0 - FP_A
      ELSE
*        rgnl, right
         FP_R = 1.0D0
      ENDIF
*
***** shock-shock + (e.f) intermediate solution
*     or
*     er-parameterisation, with linearisation over
*     the rarefaction for p, rho, w
*
*     nb expression of gam_z and pc_z are exact except
*        in the case of vacuum (u_b &lt; z &lt; u_a).
*        indeed
*        gam_z = gam_l if z &lt; un_b
*        gam_z = gam_r if z > un_a
*
      GAM_Z = ((FP_L + FP_B) * GAM_L) +
     &     ((FP_A + FP_R) * GAM_R)
      PC_Z = ((FP_L + FP_B) * PC_L) +
     &     ((FP_A + FP_R) * PC_R)
      UT1_Z = ((FP_L + FP_B) * UT1_L) +
     &     ((FP_A + FP_R) * UT1_R)
      UT2_Z = ((FP_L + FP_B) * UT2_L) +
     &     ((FP_A + FP_R) * UT2_R)
      RHO_Z = (FP_L * RHO_L) + (FP_B * RHO_B) +
     &     (FP_A * RHO_A) + (FP_R * RHO_R)
      P_Z = (FP_L * P_L) + (FP_B * P_B) +
     &     (FP_A * P_A) + (FP_R * P_R)
      U_Z = (FP_L * U_L) + (FP_B * U_B) +
     &     (FP_A * U_A) + (FP_R * U_R)
*
      IF (LOGDEB) THEN
         WRITE(*,*) 'FP_L, FP_B, FP_A, FP_R ', FP_L, FP_B, FP_A, FP_R
         WRITE(*,*) 'GAM_Z ', GAM_Z
         WRITE(*,*) 'PC_Z ', PC_Z
         WRITE(*,*) 'UT1_Z ', UT1_Z
         WRITE(*,*) 'UT2_Z ', UT2_Z
         WRITE(*,*) 'RHO_Z ', RHO_Z
         WRITE(*,*) 'P_Z ', P_Z
         WRITE(*,*) 'U_Z ', U_Z
      ENDIF
*
***** computation of the internal energy
*
      CALL ESTIFF(RHO_Z, P_Z, GAM_Z, PC_Z, E_Z, LOGDEB, LOGAN)
      IF (LOGAN) THEN
         WRITE(*,*) 'SUBROUTINE FLUSTI.F'
         WRITE(*,*) 'ANOMALY DETECTED'
         GOTO 9999
      ENDIF
      ECIN_Z = 0.5D0 * ((U_Z * U_Z) + (UT1_Z * UT1_Z) +
     &     (UT2_Z * UT2_Z))
*
***** interfacial flux
*     according to nkonga, comput methods appl. mech engnr 190, 2000
*     \dep{u}{x} + \dep{f(u)}{x} = 0
*     z = speed of the surface
*     f_z = f - z u
*
      UINT = U_Z - Z
      FLURHO = RHO_Z * (U_Z - Z)
*      write(*,*) 'flurho', flurho
      FLURU = (FLURHO * U_Z) + P_Z
      FLURET = (FLURHO * (E_Z + ECIN_Z))
     &     + (P_Z * U_Z)
      FLURT1 = FLURHO * UT1_Z
      FLURT2 = FLURHO * UT2_Z
*
*      write(*,*) (et_z + ef_z +
*     &     ecin_z), flurho
*      write(*,*) p_z, u_z
 9999 RETURN
      END