TABLE OF CONTENTS
PHY/svriemann [ Modules ]
NOM
svriemann(hg,ug,hd,ud,xi,h,u)
DESCRIPTION
Resolution du probleme de Riemann pour Saint Venant
Solving the Riemann problem for Saint Venant
ENTREES / INPUT
hg : hauteur d'eau gauche / water level left
ug : vitesse gauche / velocity left
hd : hauteur d'eau droit / water level right
ud : vitesse droit / velocity right
xi : x/t
SORTIES / OUTPUT
h : hauteur d'eau interface / water level at the interface
u : vitesse interface / interface velocity
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) :: hg !! hauteur d eau gauche / water level left real(kind=kind(0.d0)),intent(in) :: hd !! vitesse gauche / velocity left real(kind=kind(0.d0)),intent(in) :: ug !! hauteur d eau droit / water level right real(kind=kind(0.d0)),intent(in) :: ud !! vitesse droit / velocity right real(kind=kind(0.d0)),intent(in) :: xi !! x/t real(kind=kind(0.d0)),intent(out) :: h !! hauteur d eau a l interface / water level at the interface real(kind=kind(0.d0)),intent(out) :: u !! vitesse a l interface / interface velocity real(kind=kind(0.d0)) :: g,eps integer,parameter :: maxiter=100 real(kind=kind(0.d0)) :: cg,cd,cstar,ustar,hstar,dh,res real(kind=kind(0.d0)) :: xlambda1m,xlambda1p,xlambda2m,xlambda2p integer :: config,iter !=========================== DEBUT DU CODE EXECUTABLE ================== eps=hzero_phy g=gpes_phy ! test type configuration IF ((ud-ug) .GE. (2.d0*dsqrt(g)*(dsqrt(hg)+dsqrt(hd)))) THEN ! dry zone centrale config = 1 ELSEIF (hg .LE. eps .and. hd .GT. eps) THEN ! dry zone gauche config = 2 ELSEIF (hd .LE. eps .and. hg .GT. eps) THEN ! dry zone droite config = 3 ELSEIF (hg .LE. eps .and. hd .LE. eps) THEN !dry zone a droite et a gauche config = 4 ELSE ! pas de dry zone config = 0 ENDIF select case(config) ! pas de dry zone case(0) ! !------- Initialisations ! hstar=.5*(hg+hd) ! newton-raphson pour calculer h* iter=0 res=1.d0 dh=1.d0 do while((dabs(dh).gt.eps).or.(dabs(res).gt.eps)) iter=iter+1 ! calcul du residu (on veut que res=0) res= ud-ug+(hstar-hd)*zsv(hstar,hd)+(hstar-hg)*zsv(hstar,hg) ! calcul de l accroissement dh=zsv(hstar,hd) +(hstar-hd)*dzsv(hstar,hd) & + zsv(hstar,hg) +(hstar-hg)*dzsv(hstar,hg) dh=-res/dh ! mise a jour hstar=hstar+dh ! test d arret if (iter.ge.maxiter) then stop 'Newton iteration max atteinte!' endif enddo ! calcul de u* ustar=ud+(hstar-hd)*zsv(hstar,hd) ! celerite du "son" a gauche droite et milieu cg=dsqrt(g*hg) cd=dsqrt(g*hd) cstar=dsqrt(g*hstar) ! calcul de lambda(i,plus) et lambda(i,moins) ! onde 1 if (hstar.gt.hg) then ! choc xlambda1m=(hg*ug-hstar*ustar)/(hg-hstar) xlambda1p=xlambda1m else ! detente xlambda1m=ug-cg xlambda1p=ustar-cstar endif ! onde 2 if (hstar.gt.hd) then ! choc xlambda2m=(hd*ud-hstar*ustar)/(hd-hstar) xlambda2p=xlambda2m else ! detente xlambda2p=ud+cd xlambda2m=ustar+cstar endif ! h_interface et u_interface if (xi .lt. xlambda1m) then ! etat gauche h=hg u=ug elseif ((xi .ge. xlambda1m).and.(xi .lt. xlambda1p)) then ! 1-detente h=(2.d0*cg+ug-xi)**2/(9.d0*g) u=(2.d0*(cg+xi)+ug)/3.d0 elseif ((xi .ge. xlambda1p).and.(xi .lt. xlambda2m)) then ! etat milieu h=hstar u=ustar elseif ((xi .ge. xlambda2m).and.(xi .lt. xlambda2p)) then ! 2-detente h=(2.d0*cd-ud+xi)**2/(9.d0*g) u=(-2.d0*cd+2.d0*xi+ud)/3.d0 elseif (xi .ge. xlambda2p) then ! etat droit h=hd u=ud endif ! dry zone centrale case(1) ! etat milieu hstar=0.d0 ustar=0.d0 ! celerite du "son" a gauche droite et milieu cg=dsqrt(g*hg) cd=dsqrt(g*hd) cstar=0.d0 ! calcul de lambda(i,plus) et lambda(i,moins) ! onde 1 xlambda1m=ug-cg xlambda1p=ug+2.d0*cg ! onde 2 xlambda2m=ud-2.d0*cd xlambda2p=ud+cd ! h_interface et u_interface if (xi .lt. xlambda1m) then ! etat gauche h=hg u=ug elseif ((xi .ge. xlambda1m).and.(xi .lt. xlambda1p)) then ! 1-detente h=(2.d0*cg+ug-xi)**2/(9.d0*g); u=(2.d0*cg+ug+2*xi)/3.d0 elseif ((xi .ge. xlambda1p).and.(xi .lt. xlambda2m)) then ! etat milieu h=hstar u=ustar elseif ((xi .ge. xlambda2m).and.(xi .lt. xlambda2p)) then ! 2-detente h=(2.d0*cd-ud+xi)**2/(9.d0*g); u=(-2.d0*cd+ud+2*xi)/3.d0 elseif (0.d0 .ge. xlambda2p) then ! etat droit h=hd u=ud endif ! dry zone gauche case(2) cd=dsqrt(g*hd) ! calcul de lambda(i,plus) et lambda(i,moins) ! onde 2 xlambda2m=ud-2.d0*cd xlambda2p=ud+cd ! h_interface et u_interface if (xi .lt. xlambda2m) then ! etat milieu h=0.d0 u=0.d0 elseif ((xi .ge. xlambda2m).and.(xi .lt. xlambda2p)) then ! 2-detente h=(2.d0*cd-ud+xi)**2/(9.d0*g); u=(-2.d0*cd+ud+2*xi)/3.d0 elseif (xi .ge. xlambda2p) then ! etat droit h=hd u=ud endif ! dry zone droite case(3) cg=dsqrt(g*hg) ! calcul de lambda(i,plus) et lambda(i,moins) ! onde 1 xlambda1m=ug-cg xlambda1p=ug+2*cg ! h_interface et u_interface if (xi .lt. xlambda1m) then ! etat gauche h=hg u=ug elseif ((xi .ge. xlambda1m).and.(xi .lt. xlambda1p)) then ! 1-detente h=(2*cg+ug-xi)**2/(9.d0*g); u=(2*cg+ug+2*xi)/3.d0 elseif (xi .ge. xlambda1p) then ! etat milieu h=0.d0 u=0.d0 endif ! dry zone a droite et a gauche case(4) h=0.d0 u=0.d0 end select return ! !=========================== FIN DE LA ROUTINE ==================== contains function zsv(h1,h2) implicit none real(kind=kind(0.d0)) :: gg,h1,h2,zsv,zsol gg=9.81d0 zsol=0.d0 if (h1.gt.h2) then zsol=dsqrt(gg/2.d0)*dsqrt((h1+h2)/(h1*h2)) else zsol=2.d0*dsqrt(gg)/(dsqrt(h1)+dsqrt(h2)) endif zsv=zsol end function zsv function dzsv(h1,h2) implicit none real(kind=kind(0.d0)) :: gg,h1,h2,dzsv,dzsol gg=9.81d0 dzsol=0.d0 if (h1.gt.h2) then dzsol= dsqrt(2.d0)*dsqrt(gg)/4.d0*sqrt(h1*h2/(h1+h2))* & (1.d0/(h1*h2)-(h1+h2)/(h1**2*h2)) else dzsol=-dsqrt(gg)/(dsqrt(h1)+dsqrt(h2))**2/dsqrt(h1) endif dzsv=dzsol end function dzsv end subroutine svriemann