fluxnum.f90

### TABLE OF CONTENTS

- 1. PHY/fluxnum

## PHY/fluxnum [ Modules ]

[ Top ] [ Modules ]

NOM

     fluxnum(wL,wR,vn,flux,itypflux)

DESCRIPTION

     Calcul du flux numerique
     Calculating the numerical flux
     cas 0: Rusanov
     cas 1 : Godunov
     Cas 2 : Flux centre
     Cas 3 : VF-Roe

        ENTREES / INPUT
     wL   : objet de type "state" Etat de la partie gauche / left state
     wR   : objet de type "state" Etat de la partie droite / right state
     vn   : objet de type "vect_3", normale Ã  la face / normal vector to the face
     itypflux   : type de flux / type of flux

        SORTIES / OUTPUT
     flux : objet de type "vect_nvar", flux numerique / numerical flux

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
      use phy ,only:  riemisot, svriemann !compflux, vpmax, vitson, energie, prim2bal, riemann_vfroe,phy_typflux
      IMPLICIT NONE
    !.2-----  Declaration
      integer                     :: itypflux !! type de flux / type of flux
      type(state), intent(in)     :: wL , wR  !! Etat gauche et droit / left and right state
      type(vect_3),intent(in)     :: vn  !! Vecteur normal a la face / Normal vector to the face
      type(vect_nvar),intent(out) :: flux !! Flux

      real(kind=kind(0.d0))       :: unL,unR,rL,rR,pL,pR
      real(kind=kind(0.d0))       :: hL,hR,h,zL,zR,hL_reco,hR_reco
      real(kind=kind(0.d0))       :: psi,psiR,psiL,r,un,p,xi,alf
      type(vect_3)                :: uL,uR,u
    !=========================== DEBUT DU CODE EXECUTABLE ==================
    !-------  Initialisations
      flux%v=0.d0
      u%v=0.d0
      uL%v=0.d0
      uR%v=0.d0

      SELECT CASE (itypflux)
      !---------------------------------------------------
          !~CASE (0) ! flux de HLL ou Rusanov
          !~    CALL compflux(wL,vn,fluxL)
          !~    CALL compflux(wR,vn,fluxR)
          !~    SELECT CASE (imodel_phy)
          !~        CASE DEFAULT
          !~            vv=vn%v(1)**2+vn%v(2)**2+vn%v(3)**2-1.d0
          !~            IF (dabs(vv).gt.1.d-10) THEN
          !~               write(*,*) vv
          !~               STOP
          !~            END IF
          !~            vp=0.d0    ! pour plus tard: maillage mobile...
          !~            uL%v(1)=wL%vprim%v(2)
          !~            uL%v(2)=wL%vprim%v(3)
          !~            uL%v(3)=wL%vprim%v(4)
          !~            CALL vitson(wL,cL)
          !~            unL=uL%v(1)*vn%v(1)+uL%v(2)*vn%v(2)+uL%v(3)*vn%v(3)
          !~            uR%v(1)=wR%vprim%v(2)
          !~            uR%v(2)=wR%vprim%v(3)
          !~            uR%v(3)=wR%vprim%v(4)
          !~            CALL vitson(wR,cR)
          !~            unR=uR%v(1)*vn%v(1)+uR%v(2)*vn%v(2)+uR%v(3)*vn%v(3)
          !~            smin=unL-vp-cL
          !~            IF (smin.gt.unR-vp-cR) smin=unR-vp-cR
          !~            smax=unL-vp+cL
          !~            IF (smax.lt.unR-vp+cR) smax=unR-vp+cR
          !~            vmax=dabs(smin)
          !~            IF (vmax.lt.dabs(smax)) vmax=dabs(smax)
          !~    END SELECT
          !~    ! pour rusanov, dÃ©commenter ces lignes : indispensable pour le bifluide !!!!!
          !~    smax=vmax
          !~    smin=-vmax
          !~    IF (smin.ge.0.d0) THEN
          !~       flux%v=fluxL%v
          !~    ELSE IF (smax.le.0.d0) THEN
          !~       flux%v=fluxR%v
          !~    ELSE
          !~       flux%v=(smax*fluxL%v-smin*fluxR%v+smin*smax*(wR%vbal%v-wL%vbal%v))/(smax-smin)
          !~    END IF
      !---------------------------------------------------
          CASE (1) ! flux de godunov

              ! calculs communs
              xi=0.d0

              SELECT CASE (imodel_phy)
    !---------------------------------------------------
                 CASE (2) ! bifluide isotherme non conservatif
                  rL=wL%vprim%v(1)
                  uL%v(1)=wL%vprim%v(2)
                  uL%v(2)=wL%vprim%v(3)
                  uL%v(3)=wL%vprim%v(4)
                  pL=wL%vprim%v(5)
                  unL=uL%v(1)*vn%v(1)+uL%v(2)*vn%v(2)+uL%v(3)*vn%v(3)

                  rR=wR%vprim%v(1)
                  uR%v(1)=wR%vprim%v(2)
                  uR%v(2)=wR%vprim%v(3)
                  uR%v(3)=wR%vprim%v(4)
                  pR=wR%vprim%v(5)
                  unR=uR%v(1)*vn%v(1)+uR%v(2)*vn%v(2)+uR%v(3)*vn%v(3)
                    psiL=wL%vbal%v(5)
                    psiR=wR%vbal%v(5)
                    !if(wL%chi.eq.0.d0.and.wR%chi.ne.0.d0)then
                    !  psiR=psiL
                    !  rR  =rL
                    !endif
                    !if(wR%chi.eq.0.d0.and.wL%chi.ne.0.d0)then
                    !  psiL=psiR
                    !  rL  =rR
                    !endif
                     CALL riemisot(rL,unL,psiL,psiL,rR,unR,psiR,psiR,xi,r,un,p,psi,alf)
                     u%v=un*vn%v
                     IF (un.gt.0.d0) THEN
                        u%v=u%v+uL%v-unL*vn%v
                     ELSE
                        u%v=u%v+uR%v-unR*vn%v
                     END IF
                     flux%v(:)=0.d0
                     flux%v(1)=r*un
                     flux%v(2)=r*u%v(1)*un+p*vn%v(1)
                     flux%v(3)=r*u%v(2)*un+p*vn%v(2)
                     flux%v(4)=r*u%v(3)*un+p*vn%v(3)
                     flux%v(6)=(r*un*un/2+r*(1.d0+log(r))*c_phy*c_phy)*un  ! flux d entropie
                     vi_phy = un   ! stockage de la vitesse du contact pour le schema nc
                     flux%v(5)=0.d0

    !---------------------------------------------------
                 CASE (1) ! SaintVenant
                     hL=wL%vprim%v(1)
                     uL%v(1)=wL%vprim%v(2)
                     uL%v(2)=wL%vprim%v(3)
                     unL=uL%v(1)*vn%v(1)+uL%v(2)*vn%v(2)
                     zL=wL%vprim%v(5)

                     hR=wR%vprim%v(1)
                     uR%v(1)=wR%vprim%v(2)
                     uR%v(2)=wR%vprim%v(3)
                     unR=uR%v(1)*vn%v(1)+uR%v(2)*vn%v(2)
                     zR=wR%vprim%v(5)

    !Reconstruction Hydrostatique
                     IF(hL.LE.hzero_phy)hL=0.d0
                     IF(hR.LE.hzero_phy)hR=0.d0
                     hL_reco=max(0.d0,hL+zL-max(zL,zR))
                     hR_reco=max(0.d0,hR+zR-max(zL,zR))
    !secu!
                     IF(hL_reco.LE.hzero_phy)THEN
                       hL_reco=0.d0
                       unL=0.d0
                     ENDIF
                     IF(hR_reco.LE.hzero_phy)THEN
                       hR_reco=0.d0
                       unR=0.d0
                     ENDIF

                     xi=0.d0
                     call svriemann(hL_reco,unL,hR_reco,unR,xi,h,un)
    !secu!
                     if(h.LE.hzero_phy)then
                      h=0.d0
                      un=0.d0
                     endif

    !projection u_riemann sur normale a la face
                     u%v=un*vn%v
    !what is it? : construction du vecteur vitesse etoile a partir de la solution
    !scalaire du probleme de riemann et des normales aux faces
                     if (un.gt.0.d0) then
                      u%v(1)=u%v(1)+uL%v(1)-unL*vn%v(1)
                      u%v(2)=u%v(2)+uL%v(2)-unL*vn%v(2)
                     else
                      u%v(1)=u%v(1)+uR%v(1)-unR*vn%v(1)
                      u%v(2)=u%v(2)+uR%v(2)-unR*vn%v(2)
                     end if

                     flux%v(:)=0.d0
                     flux%v(1)=h*un
                     flux%v(2)=h*u%v(1)*un+gpes_phy*vn%v(1)/2*h*h !+ g*vn%v(1)/2.d0*(hL**2-hL_reco**2) voir fluxnoncons
                     flux%v(3)=h*u%v(2)*un+gpes_phy*vn%v(2)/2*h*h !+ g*vn%v(2)/2.d0*(hL**2-hL_reco**2)voir fluxnoncons
                     flux%v(4)=un*(5.d-1*h*(u%v(1)**2+u%v(2)**2)+gpes_phy*h**2+gpes_phy*h*max(zL,zR)) ! flux d entropie
    !---------------------------------------------------
                 CASE DEFAULT
                     call print_err('fluxnum','Godunov is not implemented for this model')

              END SELECT

    !-------------------------------------------------------------------
          !~CASE (2) ! flux centrÃ© (instable sauf si integration en temps de rk4)
          !~    wstar%vprim%v=0.5d0*(wL%vprim%v+wR%vprim%v)
          !~    rL=wL%vprim%v(1)
          !~    rR=wR%vprim%v(1)
          !~    IF (rL.lt.rR) THEN
          !~       temp=rL
          !~       rL=rR
          !~       rR=temp
          !~    END IF
          !~    co=rR/rL
          !~    pe=6.
          !~    pm=pe-1.d0
          !~    IF (rR.ne.rL) wstar%vprim%v(1) = (co+1.D0)*rL/2.D0
          !~    CALL prim2bal(wstar)
          !~    CALL compflux(wstar,vn,flux)
          !~    visnu=0.d0
          !~    flux%v=flux%v-visnu*(wR%vbal%v-wL%vbal%v)

    !-------------------------------------------------------------------
          !~CASE (3) ! schÃ©ma vfroe
          !~    CALL riemann_vfroe(wL,wR,vn,0.d0,wstar)
          !~    CALL compflux(wstar,vn,flux)
      END SELECT
    !===========================   FIN DE LA ROUTINE    ====================
    END SUBROUTINE fluxnum