Contents
| | | 6 Algorithm for numerical solution of flat-bottom Riemann problem by different finite volume methods of Godunov type | Contents | Index |
program godunov
c----67-------------------------------------------------------flowtec-|
c solves the flat-bottom st.venant equations numerically in
c one dimension (dam-break problem) by a finite-volume
c method using:
c ischeme fluxes
c 1 godunov (exact riemann solver)
c 3 roe (lin. riemann solver, flux-diff-split)
c 2 gallouet (lin. riemann solver, var-diff-split)
c----67-------------------------------------------------------flowtec-|
c 12.05.99
c BUGS: -?-
c----67---------------------------------------------------------------|
implicit double precision (a-h,o-z)
parameter(nmax=400,g=2.d0)
dimension q(nmax,2),rhs(nmax,2),flux(2)
parameter(cfl=.4d0)
c----67---------------------------------------------------------------|
parameter(hll=1.d0,ull=0.d0,hrr=5.d-1,urr=0.d0)
parameter(flx=1.d0,tfin=.1d0,nitmx=50000)
parameter(ischeme=3)
c----67---------------------------------------------------------------|
time=0.d0
nit=0
dx=flx/dfloat(nmax)
do i=1,nmax
rhs(i,1)=0.d0
rhs(i,2)=0.d0
enddo
c /* set the initial state (diaphragm at center)*/
do i=1,nmax/2
q(i,1)=hll
q(i,2)=ull*hll
enddo
do i=nmax/2+1,nmax
q(i,1)=hrr
q(i,2)=urr*hrr
enddo
c /* main time loop */
11 continue
nit=nit+1
dt=1.d35
c /* flux balance for each cell-interface */
do 20 icell=1,nmax-1
c /* solve interfacial riemann problem analytically */
if(ischeme.eq.1)then
call riemann(q(icell,1),q(icell,2)/q(icell,1),
$ q(icell+1,1),q(icell+1,2)/q(icell+1,1),
$ g,
$ fhc,uc,flamx)
elseif(ischeme.eq.2)then
call gallouet(q(icell,1),q(icell,2),
$ q(icell+1,1),q(icell+1,2),
$ g,
$ fhc,uc,flamx)
elseif(ischeme.eq.3)then
call roeflux(q(icell,1),q(icell,2),
$ q(icell+1,1),q(icell+1,2),
$ g,
$ flux,flamx)
else
write(*,*)'small problem here...',ischeme
endif
c /* calculate maximum GLOBAL time step */
dt=dmin1(dt,dx*cfl/flamx)
if(dt.eq.0.d0.or.flamx.eq.0.d0)then
write(*,*)'icell =',icell,dx,cfl,flamx
stop
endif
c /* accumulate rhs fluxes */
if (ischeme.eq.1.or.ischeme.eq.2)then
do ic=1,2
add=dflux(fhc,uc,g,ic)
rhs(icell,ic)=rhs(icell,ic)+add
rhs(icell+1,ic)=rhs(icell+1,ic)-add
enddo
else
do ic=1,2
rhs(icell,ic)=rhs(icell,ic)+flux(ic)
rhs(icell+1,ic)=rhs(icell+1,ic)-flux(ic)
enddo
endif
20 continue
time=time+dt
write(*,*)'time = ',time, ' dt= ',dt, ' tfin = ',tfin
if(time.gt.tfin.or.nit.gt.nitmx)then
call output(q,time,flx,nmax,g)
goto 10
endif
c /* update the variables */
do ival=2,nmax-1
do ic=1,2
q(ival,ic)=q(ival,ic)-dt/dx*rhs(ival,ic)
rhs(ival,ic)=0.d0
enddo
enddo
c /* check for non-physical values */
do i=1,nmax
if(q(i,1).lt.0d0)then
WRITE(*,*)'negative height: ',i,q(i,1),q(i,2)
stop
endif
enddo
goto 11
10 continue
stop
end
c----67-------------------------------------------------------flowtec-|
subroutine output(q,time,flx,ni,g)
implicit double precision (a-h,o-z)
dimension q(ni,2)
open(12,file='sol.dat')
write(12,100)time,ni
do i=1,ni
write(12,*)dfloat(i-1)/dfloat(ni-1)*flx,
$ q(i,1),q(i,2)/q(i,1),q(i,2),sqrt(g*q(i,1))
enddo
close(12)
100 format('# 1:x 2:h 3:u 4:h*u 5:c')
return
end
c----67-------------------------------------------------------flowtec-|
subroutine riemann(fhl,ul,fhr,ur,g,fhc,uc,flamx)
c----67-------------------------------------------------------flowtec-|
c solves the riemann problem for st.venant equations (flat bottom)
c analytically by considering the following zones:
c (e.g. for a rarefaction-shock configuration)
c
c (L) (1) (C) (2) (R)
c ------
c \\\\\\ ----------
c \\\\\ |
c \\\\ |
c \\\ | ^ state q
c \\ | |
c \ | |
c --------------------- -----> (x/t)
c----67---------------------------------------------------------------|
c 12.02.99
c BUGS: -?-
c----67---------------------------------------------------------------|
implicit double precision (a-h,o-z)
logical lvacuum,l1simple,l2simple,lldry,lrdry
parameter (f23=2.d0/3.d0,f13=1.d0/3.d0,f19=1.d0/9.d0)
c----67---------------------------------------------------------------|
c parameter (g=2.,time=1.d0,x0=2.5d0)
parameter (nnwtmx=20,small=1.d-6,width=1.2d0)
c----67---------------------------------------------------------------|
lrdry=.FALSE.
lldry=.FALSE.
cL=dsqrt(g*fhL)
cR=dsqrt(g*fhR)
if(fhl.ne.0.d0)then
if(fhR.eq.0.d0)then
c write(*,*)'dry bed on the right'
lrdry=.TRUE.
endif
B=fhR/fhL
C=(uR-uL)/cL
else
c write(*,*)'dry bed on the left'
lldry=.TRUE.
B=0.d0
C=0.d0
endif
c /* check for equal states on left and right */
if(fhr.eq.fhl.and.ur.eq.ul)then
fhc=fhr
uc=ur
flamx=dmax1(dabs(ur+sqrt(g*fhr)),dabs(ur-sqrt(g*fhr)),
$ dabs(ul+sqrt(g*fhl)),dabs(ul-sqrt(g*fhl)))
goto 112
endif
c /* check for vacuum (i.e. dry bed) */
lvacuum=.FALSE.
if(2.d0*(1.d0+dsqrt(B)).le.C)lvacuum=.TRUE.
c write(*,*)'appearance of vacuum:',lvacuum
c if(lvacuum)write(*,*)'will set velocity arbitrarily to zero'
c /* check of which type are the 1- and 2-families */
l1simple=.FALSE.
l2simple=.FALSE.
if(B.gt.0.d0)
$ check1=f1(dlog(B))*dsqrt(B)
if(check1.lt.C.or.lrdry)l1simple=.TRUE.
if(B.gt.0.d0)
$ check2=f1(-dlog(b))
if(check2.lt.C.or.lldry)l2simple=.TRUE.
c write(*,*)fhl,ul,fhr,ur
c write(*,*)'1-family:'
c if(l1simple)write(*,*)'simple wave'
c if(.not.l1simple.and..not.lldry)write(*,*)'shock'
c write(*,*)'2-family:'
c if(l2simple)write(*,*)'simple wave'
c if(.not.l2simple.and..not.lrdry)write(*,*)'shock'
c /* solve the 1-family */
if(l1simple)then
if(lvacuum.or.lrdry)then
fhC=0.d0
uc=0.d0
else
x1=-2.d0*dlog((1.d0+dsqrt(B))/2.-C/4.d0)
fhC=fhL*dexp(-x1)
CC=dsqrt(g*fhc)
uC=uL+cL*2.d0*(1.d0-dexp(-x1/2.d0))
c write(*,*)'1-rarefaction: hC=',fhC,' uC=',uC
endif
else
if(l2simple)then
if(lvacuum.or.lldry)then
fhC=0.d0
uc=0.d0
else
x2=-2.d0*dlog(-C/(dsqrt(B)*4.d0)+1.d0/(2.d0*dsqrt(B))+
$ .5d0)
fhC=fhR*dexp(-x2)
CC=dsqrt(g*fhc)
uC=uR-cC*2.d0*(dexp(x2/2.d0)-1.d0)
c write(*,*)'2-rarefaction: hC=',fhC,' uC=',uC
endif
else
c /* we have in fact two shocks: need to iterate */
z1=1.1d0
do n=1,nnwtmx
zold=z1
z1=zold-funcz1(zold,B,C)/dfuncz1(zold,B)
if(dabs(z1-zold).le.small)goto 111
enddo
write(*,*)'newton iteration not converged',z1,zold
stop
111 continue
fhC=z1*fhL
cc=dsqrt(fhc*g)
uC=uL-cL*(z1-1.d0)*dsqrt((z1+1.d0)/(2.d0*z1))
endif
endif
c /* we now know about the center state (c) */
c /* calculate limits of zone (1) and (2) */
c /* note: "positions" x** are actually given in "xi/t",i.e. velocities */
if(l1simple)then
xL1=uL-cL
if(.not.lvacuum.and..not.lrdry)then
x1C=uC-cC
else
c /* limit zone (1) such that vacuum is just reached through */
c /* the expansion */
x1C=uL+2.d0*cL
endif
else
c /* find shock position: calc shock speed from jump cond. */
c /* note: in this case, the zone has zero width */
if(.not.lldry)then
z=fhC/fhL
c write(*,*)'1-shock: hL/hC=',z
if(dabs(z-1.d0).gt.small)then
sigma=(uC*z-uL)/(z-1.d0)
else
sigma=ul-sqrt(g*fhl)
endif
xl1=sigma
x1C=sigma
else
c /* leave this bound undefined */
endif
endif
if(l2simple)then
x2R=uR+cR
if(.not.lvacuum.and..not.lldry)then
xC2=uC+cC
else
xC2=uR-2.d0*cR
endif
else
if(.not.lrdry)then
z=fhR/fhC
c write(*,*)'2-shock: hR/hC=',z
if(dabs(z-1.d0).gt.small)then
sigma=(uR*z-uC)/(z-1.d0)
else
sigma=ur+sqrt(g*fhr)
endif
x2R=sigma
xC2=sigma
else
c /* since nothing happens in the dry zone, extend it a bit from (1)(C) */
x2R=x1C+dabs(.1d0*x1c)
xC2=x2R
endif
endif
c /* get back to the case of a dry bed on the left */
if(lldry)then
x1C=xC2-dabs(.1d0*xC2)
xl1=x1c
endif
c /* compute maximum wave speed from limiting velocities xL1,x2R */
flamx=dmax1(dabs(xL1),dabs(x2R))
c /* finalize */
112 continue
return
end
c----67-------------------------------------------------------flowtec-|
subroutine gallouet(fhl,ql,fhr,qr,g,fhc,uc,flamx)
c----67-------------------------------------------------------flowtec-|
c solves the riemann problem for st.venant equations (flat bottom)
c in its linearized form, cf Buffard, Gallouet, Herard, CRAcadSci,
c t.326,Serie I,p.385-390, 1998
c----67---------------------------------------------------------------|
c 12.02.99
c BUGS: -?-
c----67---------------------------------------------------------------|
implicit double precision (a-h,o-z)
logical lvacuum
parameter (f23=2.d0/3.d0,f13=1.d0/3.d0,f19=1.d0/9.d0)
parameter (small=1.d-6)
c----67---------------------------------------------------------------|
if(fhl.ne.0.d0)then
ul=ql/fhl
else
ul=0.d0
endif
if(fhr.ne.0.d0)then
ur=qr/fhr
else
ur=0.d0
endif
cl=dsqrt(g*fhl)
cr=dsqrt(g*fhr)
um=(ul+ur)/2.d0
fhm=(fhr+fhl)/2.d0
cm=(cl+cr)/2.d0
c /* check for existence of dry bottom */
if((ul-ur).ge.2.d0*(cl+cr))then
lvacuum=.true.
uc=0.d0
cc=0.d0
flamx=small
goto 111
endif
c /* intermediate state in non-vacuum case */
uc=(um-(cr-cl))
cc=(cr+cl)/4.d0*(2.d0-(ur-ul)/(cr+cl))
flamx=dmax1(dabs(um+cm),dabs(um-cm))
111 continue
fhc=cc*cc/g
c
return
end
c----67---------------------------------------------------------------|
subroutine roeflux(fhl,ul,fhr,ur,g,rflux,flamx)
c----67-------------------------------------------------------flowtec-|
c solves the linearized riemann problem for st.venant equations
c (flat bottom) by using Roe's method (J.Comp.Physics, vol.43,1981)
c * fluxes are passed upon exit in vector rflux(1:2)
c----67---------------------------------------------------------------|
c 21.05.99
c BUGS: -?-
c----67---------------------------------------------------------------|
implicit double precision (a-h,o-z)
dimension a(2,2),cflux(2),uflux(2),diffv(2),rflux(2)
logical lvacuum
parameter (f23=2.d0/3.d0,f13=1.d0/3.d0,f19=1.d0/9.d0)
parameter (small=1.d-6)
c----67---------------------------------------------------------------|
cl=dsqrt(g*fhl)
cr=dsqrt(g*fhr)
c /* roe's average for shallow water */
fhc=(fhl+fhr)/2.d0
uc=(dsqrt(fhr)*ur+dsqrt(fhl)*ul)/(dsqrt(fhr)+dsqrt(fhl))
cc=dsqrt(g*fhc)
c /* diagonalization |A|=R*|lambda|*R1 */
call aroe(fhc,uc,cc,a)
c /* variable difference vector */
diffv(1)=fhr-fhl
diffv(2)=fhr*ur-fhl*ul
c /* centered flux */
do i=1,2
cflux(i)=dflux(fhr,ur,g,i)+dflux(fhl,ul,g,i)
enddo
c /* evaluate roe's flux: f=(fl+fr-|a|*(qr-ql))/2 */
do i=1,2
rflux(i)=cflux(i)
do j=1,2
rflux(i)=rflux(i)-a(i,j)*diffv(j)
enddo
rflux(i)=rflux(i)/2.d0
enddo
c /* estimate the maximum signal velocity */
flamx=dmax1(dabs(uc+cc),dabs(uc-cc))
c
return
end
c----67---------------------------------------------------------------|
subroutine aroe(hc,uc,cc,a)
c----67-------------------------------------------------------flowtec-|
c calculates roe's matrix
c----67---------------------------------------------------------------|
implicit double precision (a-h,o-z)
dimension a(2,2),r(2,2),r1(2,2),flambda(2)
c----67---------------------------------------------------------------|
c /* diagonalization matrix r */
r(1,1)=1.d0
r(1,2)=1.d0
r(2,1)=uc-cc
r(2,2)=uc+cc
c /* inverse r1 */
r1(1,1)=(uc+cc)/(2.d0*cc)
r1(1,2)=-1.d0/(2.d0*cc)
r1(2,1)=(cc-uc)/(2.d0*cc)
r1(2,2)=1.d0/(2.d0*cc)
c /* eigenvalues */
flambda(1)=uc-cc
flambda(2)=uc+cc
c /* matrix*diag*matrix multiply */
do i=1,2
do j=1,2
a(i,j)=0.d0
do k=1,2
a(i,j)=a(i,j)+r(i,k)*dabs(flambda(k))*r1(k,j)
enddo
enddo
enddo
c
return
end
c----67---------------------------------------------------------------|
double precision function f1(x)
implicit double precision (a-h,o-z)
c
if(x.lt.0.d0)then
f1=-(dexp(-x)-1.d0)*dsqrt((dexp(-x)+1)/(2.d0*dexp(-x)))
elseif(x.gt.0.d0)then
f1=2.d0*(1.d0-dexp(-x/2.d0))
else
f1=0.d0
endif
end
c----67---------------------------------------------------------------|
double precision function f2(x)
implicit double precision (a-h,o-z)
c
if(x.lt.0.d0)then
f2=(dexp(x)-1.d0)*dsqrt((dexp(x)+1.d0)/(2.d0*dexp(x)))
elseif(x.gt.0.d0)then
f2=2.d0*(dexp(x/2.d0)-1.d0)
else
f2=0.d0
endif
end
c----67---------------------------------------------------------------|
double precision function funcz1(z,B,C)
implicit double precision (a-h,o-z)
c
funcz1= -(z-1.D0)*dsqrt(2.D0)*dsqrt((z+1.D0)/z)/2.D0-
$ (z/B-1.D0)*dsqrt(2.D0)*dsqrt((z/B+1.D0)/z*B)*
$ dsqrt(B)/2.D0-C
c
end
c----67---------------------------------------------------------------|
double precision function dfuncz1(z,B)
implicit double precision (a-h,o-z)
c
dfuncz1= -dsqrt(2.D0)*dsqrt((z+1.D0)/z)/2.D0-
$ (z-1.D0)*dsqrt(2.D0)/dsqrt((z+1.D0)/z)*
$ (1.D0/z-(z+1.D0)/z**2)/4.D0-1.D0/dsqrt(B)*
$ dsqrt(2.D0)*dsqrt((z/B+1.D0)/z*B)/2.D0-
$ (z/B-1.D0)*dsqrt(2.D0)/dsqrt((z/B+1.D0)/z*B)*
$ dsqrt(B)*(1.D0/z-(z/B+1.D0)/z**2*B)/4.D0
c
end
c----67---------------------------------------------------------------|
double precision function dflux(h,u,g,i)
implicit double precision (a-h,o-z)
c
if (i.eq.1)then
dflux=h*u
elseif(i.eq.2)then
dflux=u*u*h+g*h*h/2.d0
else
write(*,*)'problem here:..',i
stop
endif
c
end
c----67---------------------------------------------------------------|
markus.uhlmann AT ciemat.es
| | | 6 Algorithm for numerical solution of flat-bottom Riemann problem by different finite volume methods of Godunov type | Contents | Index |