riemisot.f90

### TABLE OF CONTENTS

- 1. PHY/riemisot

## PHY/riemisot [ Modules ]

[ Top ] [ Modules ]

NOM

     riemisot(rg,ug,phig,alfg,rd,ud,phid,alfd,xi,r,u,p,phi,alf)

DESCRIPTION

     Routine de resolution du probleme de Riemann isotherme
     avec toutes les routines necessaires
     Routine for solving the isothermal Riemann problem
     issu de V3D (code recherche F77 Helluy/Golay) / from V3D (F77 research code Helluy/Golay)
     la loi de pression est de la forme / the pressure law is of the form
     p=alf p1 + (1-alf) p2
     p1=p0_phy + c_phy**2 (   alf /   phi *rho-rho0A_phy)
     p2=p0_phy + c_phy**2 ((1-alf)/(1-phi)*rho-rho0W_phy)
     alf et phi sont convectes a la vitesse u / alf and phi are convected at velocity u

        ENTREES / INPUT
      rg    : masse volumique a gauche / density left
      ug    : vitesse a gauche / velocity left
      phig  : fraction de masse de fluide a gauche / mass fraction of fluid at left
      alfg  : fraction de volume de fluide a gauche / volume fraction of fluid at left
      rd    : masse volumique a droite / density right
      ud    : vitesse a droite / velocity right
      phid  : fraction de masse de fluide a droite / mass fraction of fluid at right
      alfd  : fraction de volume de fluide  a droite / volume fraction of fluid at right
      xi    : x/t   (x/t=0 pour godunov)

        SORTIES / OUTPUT
       r   : masse volumique / density
       u   : vitesse / velocity
       p   : pression / pressure
       phi : fraction de masse de fluide / mass fraction of fluid
       alf : fraction de volume de fluide / volume fraction of fluid

SOURCE

    ! Ce logiciel est regi par la licence [CeCILL-B]
    ! This software is governed by the [CeCILL-B] license
    !=========================== DEBUT DES DECLARATIONS ====================
    !.1-----  Implicit, Use
      use phy_typ
      IMPLICIT NONE
    !.2-----  Declaration
      real(kind=kind(0.d0)),intent(in)  ::  rg    !! masse volumique a gauche / left density
      real(kind=kind(0.d0)),intent(in)  ::  ug    !! vitesse a gauche / velocity left
      real(kind=kind(0.d0)),intent(in)  ::  phig  !! fraction de masse de fluide a gauche / fraction of fluid mass on left
      real(kind=kind(0.d0)),intent(in)  ::  alfg  !! fraction de volume de fluide a gauche / fraction of fluid volume at left
      real(kind=kind(0.d0)),intent(in)  ::  rd    !! masse volumique a droite / density right
      real(kind=kind(0.d0)),intent(in)  ::  ud    !! vitesse a droite / velocity right
      real(kind=kind(0.d0)),intent(in)  ::  phid  !! fraction de masse de fluide a droite / mass fraction of fluid at right
      real(kind=kind(0.d0)),intent(in)  ::  alfd  !! fraction de volume de fluide  a droite / volume fraction of fluid at right
      real(kind=kind(0.d0)),intent(in)  ::  xi    !! x/t   (x/t=0 pour godunov) / (x/t=0 for godunov)
      real(kind=kind(0.d0)),intent(out) ::  r     !! masse volumique / density
      real(kind=kind(0.d0)),intent(out) ::  u     !! vitesse / velocity
      real(kind=kind(0.d0)),intent(out) ::  p     !! pression / pressure
      real(kind=kind(0.d0)),intent(out) ::  phi   !! fraction de masse de fluide / mass fraction of fluid
      real(kind=kind(0.d0)),intent(out) ::  alf   !! fraction de volume de fluide / volume fraction of fluid

      real(kind=kind(0.d0))             :: r0g,r0d,s1m,s1p,s2m,s2p
      real(kind=kind(0.d0))             :: ugini,udini,um,ui,pmin,pming,pmind
      real(kind=kind(0.d0))             :: pg,pd,pdini,pgini,pm,cd,cg
      real(kind=kind(0.d0))             :: rdini,rgini,r1,r2
      real(kind=kind(0.d0))             :: dp,drg,drd,ff,df
      integer                           :: itermax,iter
      real(kind=kind(0.d0))             :: eps=1.d-10

    !=========================== DEBUT DU CODE EXECUTABLE ==================
    !-------  Initialisations
      IF (rd.le.0.d0) write(*,*)'rd',rd
      IF (rg.le.0.d0) write(*,*)'rg',rg

      ! vitesses du son Ã  droite et Ã  gauche
      cg=dsqrt(phig*c_phy**2+(1.d0-phig)*c_phy**2)
      cd=dsqrt(phid*c_phy**2+(1.d0-phid)*c_phy**2)

      r0g=(alfg*c_phy**2*rho0A_phy+(1.d0-alfg)*c_phy**2*rho0W_phy)/cg**2
      r0d=(alfd*c_phy**2*rho0A_phy+(1.d0-alfd)*c_phy**2*rho0W_phy)/cd**2

      pming=p0_phy-cg**2*r0g
      pmind=p0_phy-cd**2*r0d
      IF (pming.le.pmind) THEN
         pmin=pmind
      ELSE
         pmin=pming
      END IF

      drg=rg-r0g
      drd=rd-r0d

      pg=p0_phy+cg**2*(rg-r0g)
      pd=p0_phy+cd**2*(rd-r0d)

      p=0.5d0*(pg+pd)

      dp=1
      ff=1
      iter=1
      itermax=100

      eps=1.d-6

      rgini=rg
      ugini=ug
      pgini=pg
      rdini=rd
      udini=ud
      pdini=pd

      DO WHILE (dabs(dp).gt.eps.and.iter.lt.itermax)

          IF (p.lt.pmin) p=pmin+eps
          iter=iter+1
          r2=(p-p0_phy)/cd**2+r0d
          r1=(p-p0_phy)/cg**2+r0g
          if(r1.le.0.d0)write(*,*)r1,p,phig
          if(r2.le.0.d0)write(*,*)r2,p,phid
          ff=ud-ug-h(rd,r2,cd)-h(rg,r1,cg)
          df=-dhdr(rd,r2,cd)/cd**2-dhdr(rg,r1,cg)/cg**2
          dp=-ff/df
          p=p+dp

      END DO

      IF (iter.eq.itermax) THEN
          write(*,*) 'Riemiso pas de convergence',dabs(dp/p0_phy)
          write(*,*) 'rg ug phig alfg pg',rg,ug,phig,alfg,pg
          write(*,*) 'rd ud phid alfd pd',rd,ud,phid,alfd,pd
          STOP
      END IF

      ! intermediate densities and velocity
      r2=(p-p0_phy)/cd**2+r0d
      r1=(p-p0_phy)/cg**2+r0g

      um=ug+h(rg,r1,cg)
      um=ud-h(rd,r2,cd)

      ! vitesse d interface
      ui=um
      pm=p

      ! wave velocities

      ! 1-shock or 1-rarefaction
      IF (r1.gt.rg) THEN
          s1m=ug-cg*dsqrt(r1/rg)
          s1p=um-cg*dsqrt(rg/r1)
      ELSE
          s1m=ug-cg
          s1p=um-cg
      END IF

      ! 2-shock or 2-rarefaction
      IF (r2.gt.rd) THEN
          s2m=ud+cd*dsqrt(r2/rd)
          s2p=um+cd*dsqrt(rd/r2)
      ELSE
          s2m=um+cd
          s2p=ud+cd
      END IF

      ! computation of the solution at xi=x/t

      IF (xi.lt.s1m) THEN
          r=rg
          u=ug
          phi=phig
          alf=alfg
          p=pg
      ELSE IF (xi.lt.s1p) THEN
          u=xi+cg
          r=rg*exp((ug-u)/cg)
          phi=phig
          alf=alfg
          p=p0_phy+cg**2*(r-r0g)
      ELSE IF (xi.lt.um) then
          r=r1
          u=um
          phi=phig
          alf=alfg
          p=p0_phy+cg**2*(r-r0g)
      ELSE IF (xi.lt.s2m) then
          r=r2
          u=um
          phi=phid
          alf=alfd
          p=p0_phy+cd**2*(r-r0d)
      ELSE if (xi.lt.s2p) then
          u=xi-cd
          r=rd*exp((u-ud)/cd)
          phi=phid
          alf=alfd
          p=p0_phy+cd**2*(r-r0d)
      ELSE
          r=rd
          u=ud
          phi=phid
          alf=alfd
          p=p0_phy+cd**2*(r-r0d)
      END IF

    CONTAINS

    !=======================================================================
    FUNCTION h(rL,rLocal,c)
    !------------------
      use phy
      IMPLICIT NONE

      real(kind=kind(0.d0)), intent(in) :: rL,rLocal,c
      real(kind=kind(0.d0))             :: h

      IF (rLocal*rL.le.0.d0) THEN
          write(*,*)'pb function h: r=',rLocal,' rL=',rL
          STOP
      END IF
      IF (rLocal.gt.rL) THEN
          h=c*(rL-rLocal)/dsqrt(rL*rLocal)
      ELSE
          h=c*log(rL/rLocal)
      END IF

    END FUNCTION h
    !=======================================================================
    FUNCTION dhdr(rL,rLocal,c)
    !------------------
      use phy
      IMPLICIT NONE

      real(kind=kind(0.d0)), intent(in) :: rL,rLocal,c
      real(kind=kind(0.d0))             :: dhdr,t2,t7

      IF(rLocal*rL.le.0.d0) THEN
          write(*,*)'pb function dh'
          STOP
      END IF

      IF (rLocal.gt.rL) THEN
          t2 = dsqrt(rL*rLocal)
          t7 = t2*t2
          dhdr = -c/t2-c*(rL-rLocal)/t7/t2*rL/2.D0
      ELSE
          dhdr=-c/rLocal
      END IF

    END FUNCTION dhdr

    END SUBROUTINE riemisot

    !===========================   FIN DE LA ROUTINE    ====================