TABLE OF CONTENTS
- 1. PHY/fluxnum
PHY/fluxnum [ 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