c ========================================================= subroutine rp1(maxmx,meqn,mwaves,mbc,mx,ql,qr,auxl,auxr, & wave,s,amdq,apdq) c ========================================================= c c # solve Riemann problems for the 1D shallow water equations c # (h)_t + (u h)_x = 0 c # (h u)_t + ( u^2 h + 0.5 grav*h^2 )_x = 0 c c # using an exact Riemann solver with the possibility to handle c # the dry-bed case h=0 in the data or in the solution. c c # On input, ql contains the state vector at the left edge of each cell c # qr contains the state vector at the right edge of each cell c # On output, wave contains the waves, c # s the speeds, c # amdq the left-going flux difference A^- \Delta q c # apdq the right-going flux difference A^+ \Delta q c c # Note that the i'th Riemann problem has left state qr(i-1,:) c # and right state ql(i,:) c # From the basic clawpack routine step1, rp is called with ql = qr = q. c c # Here meqn=mwaves=2 should be passed from the calling routine c # Written by James Rossmanith and Randy LeVeque c # Based on the procedure in Section 13.10 of c # "Finite Volume Methods for Hyperbolic Problems" by R. J. LeVeque c c # NOTATION: For p = 1 or 2: c # iwavep = 1 if p-wave is a shock c # = 2 if p-wave is a rarefaction c # if p-wave is a shock: c # sl(p) = sr(p) = shock speed c # if p-wave is a rarefaction: c # sl(p) < sr(p) are the speeds of left and right edge c implicit double precision (a-h,o-z) dimension ql(1-mbc:maxmx+mbc, meqn) dimension qr(1-mbc:maxmx+mbc, meqn) dimension s(1-mbc:maxmx+mbc, mwaves) dimension wave(1-mbc:maxmx+mbc, meqn, mwaves) dimension amdq(1-mbc:maxmx+mbc, meqn) dimension apdq(1-mbc:maxmx+mbc, meqn) dimension auxl(1-mbc:maxmx+mbc, *) dimension auxr(1-mbc:maxmx+mbc, *) dimension sl(1:2) dimension sr(1:2) dimension fl(1:2), fr(1:2), fstar(1:2) c c # local storage c --------------- dimension delta(2) common /comrp/ grav logical verbose data verbose /.false./ c # set verbose to true only for debugging -- gives lots of output! c MaxIters = 20 !# for Newton's method TOL = 1.0d-14 !# convergence criterion for Newton critH = 1.d-20 !# Critical value of h: set to zero below this. do 30 i=2-mbc,mx+mbc c # Define left and right states for each i hl = qr(i-1,1) ! left state height hr = ql(i,1) ! right sate height hul = qr(i-1,2) ! left state momentum hur = ql(i,2) ! right state momentum if (verbose) then write(6,*) ' ' write(6,*) 'i = ', i write(6,*) 'hl, hul:', hl,hul write(6,*) 'hr, hur:', hr,hur endif c # check for negative depth and reset to zero: if (hl .lt. 0.d0) then write(6,*) ' *** resetting hl to zero from ',hl hl = 0.d0 endif if (hr .lt. 0.d0) then write(6,*) ' *** resetting hr to zero from ',hr hr = 0.d0 endif c # Compute ul and ur, and check for dry states in data: c ------------------------------------------------------ if (hl.gt.critH) then ul = hul/hl else c # dry state to left hl = 0.d0 ul = 0.d0 iwave1 = 1 !# 1-wave is a zero-strength shock iwave2 = 2 !# 2-wave is a rarefaction endif c if (hr.gt.critH) then ur = hur/hr else c # dry state to right hr = 0.d0 ur = 0.d0 iwave1 = 2 !# 1-wave is a rarefaction iwave2 = 1 !# 2-wave is a zero-strength shock if (hl.eq.0.d0) then c # if both states are dry, both waves are zero-strength shocks iwave1 = 1 endif endif c c # wave speeds: cl = dsqrt(grav*hl) !# wave speed in left state cr = dsqrt(grav*hr) !# wave speed in right state c c --------------------------------------- c STEP #1: Compute intermediate height hm c --------------------------------------- if (hl.eq.0.d0 .or .hr.eq.0.d0) then c # if h=0 in one state or both, then hm = 0: hm = 0.d0 c # um is determined by connecting to the wet state with a c # rarefaction, or um=0 if both states are dry: um = ur + ul - 2.d0*(dsqrt(grav*hr) - dsqrt(grav*hl)) cm = 0.d0 c # jump to Step 2: go to 200 endif c c # check for dry state appearing in solution: c -------------------------------------------- c c # uml is velocity at right edge of a 1-rarefaction connected to c # ql and extended to hm = 0: uml = ul + 2.d0*cl c # umr is velocity at left edge of a 2-rarefaction connected to c # qr and extended to hm = 0: umr = ur - 2.d0*cr if (uml .le. umr) then c # dry state appears in solution if (verbose) write(6,*) 'dry state appears' hm = 0.d0 cm = 0.d0 c c # both waves are rarefactions: iwave1 = 2 iwave2 = 2 sl(1) = ul - cl sr(1) = uml - cm sl(2) = umr + cm sr(2) = ur + cr c # um is not unique in this case, but only hm*um is used later c # and this will be zero since hm=0, so might as well set um=0: um = 0.d0 c # jump to Step 3: go to 300 endif c # Solve Riemann problem for hm, unless hm=0 has already been found c ------------------------------------------------------------------ c # initial guess for hm based on two rarefaction waves, from (13.55) c # in "Finite Volume Methods for Hyperbolic Problems": h = (1.d0/(16.d0*grav))*(ul - ur + 2.d0*(cl + cr))**2 c c if (h.le.hl .and. h.le.hr) then c # two rarefaction waves is correct if (verbose) write(6,*) 'two rarefactions' iwave1 = 2 iwave2 = 2 hm = h um = ul + 2.d0*(dsqrt(grav*hl) - dsqrt(grav*hm)) cm = dsqrt(grav*hm) go to 200 endif c c c # Two rarefaction waves is not correct. Use hm computed c # above as initial guess for Newton itertion to find correct hm, c # by solving phi_l(h) - phi_r(h) = 0 as described in Section 13.10. c # dphil and dphir denote derivatives phi'_l(h) and phi'_r(h) c c c c # Newton iteration: c if (verbose) write(6,*) 'beginning Newton, h = ',h do 50 iter=1,MaxIters c # left wave if (h.gt.hl) then iwave1 = 1 sqterm = dsqrt(0.5d0*grav*(1.d0/h + 1.d0/hl)) phil = ul - (h-hl)*sqterm dphil = -sqterm + grav*(h-hl) / (4.d0*h*h*sqterm) else iwave1 = 2 phil = ul + 2.d0*(dsqrt(grav*hl) - dsqrt(grav*h)) dphil = -dsqrt(grav/h) endif c # right wave if (h.gt.hr) then iwave2 = 1 sqterm = dsqrt(0.5d0*grav*(1.d0/h + 1.d0/hr)) phir = ur + (h-hr)*sqterm dphir = sqterm - grav*(h-hr) / (4.d0*h*h*sqterm) else iwave2 = 2 phir = ur - 2.d0*(dsqrt(grav*hr) - dsqrt(grav*h)) dphir = dsqrt(grav/h) endif c c # Newton update: deltah = (phil-phir) / (dphil-dphir) h = h - deltah if (verbose) write(6,*) 'iter,h:',iter,h if (dabs(deltah) .lt. TOL) go to 100 50 continue c write(6,599) iter 599 format('nonconvergence of Newton iteration, iterations = ', & i4) stop c 100 continue c if (verbose) then write(6,601) iter 601 format('Newton iteration converged in ',i4,' iterations') write(6,*) ' ' endif c c # After convergence of Newton, we should have phil = phir = um, c # the intermediate velocity in the Riemann solution, c # and h is the value of hm, the intermediate depth. c hm = h um = phil cm = dsqrt(grav*hm) c if (dabs(phil-phir) .gt. 1d-4) then write(6,*) 'phil and phir do not agree:',phil,phir write(6,*) 'dphil, dphir:',dphil,dphir write(6,*) 'delta h:',deltah stop endif 200 continue c c --------------------------------------------------------- c STEP #2: Compute wave speeds c --------------------------------------------------------- c if (iwave1.eq.1) then c # 1-wave is a shock: if (dabs(hl-hm).le.1d-10) then s1 = ul - cl else s1 = (hl*ul - hm*um)/(hl-hm) endif sl(1) = s1 sr(1) = s1 else c # 1-wave is a rarefaction: sl(1) = ul - cl sr(1) = um - cm endif c if (iwave2.eq.1) then c # 2-wave is a shock: if (dabs(hr-hm).le.1d-10) then s2 = ur + cr else s2 = (hr*ur - hm*um)/(hr-hm) endif sl(2) = s2 sr(2) = s2 else c # 2-wave is a rarefaction: sl(2) = um + cm sr(2) = ur + cr endif c 300 continue c c ---------------------------------------------------------------- c STEP #3: Set hstar, ustar to Godunov values along x/t = 0 c ---------------------------------------------------------------- c if (sl(1).lt.0.d0 .and. sr(1).gt.0.d0) then c # 1-wave: transonic rarefaction if (verbose) write(6,*) '1-wave is transonic' hstar = (1.d0/(9.d0*grav))*(ul + 2.d0*dsqrt(grav*hl))**2 ustar = ul - 2.d0*( dsqrt(grav*hstar) - dsqrt(grav*hl) ) elseif (sl(2).lt.0.d0 .and. sr(2).gt.0.d0) then c # 2-wave: transonic rarefaction if (verbose) write(6,*) '2-wave is transonic' hstar = (1.d0/(9.d0*grav))*(ur - 2.d0*dsqrt(grav*hr))**2 ustar = ur + 2.d0*( dsqrt(grav*hstar) - dsqrt(grav*hr) ) elseif (sl(1).gt.0.d0) then c # fully supersonic and rightgoing if (verbose) write(6,*) 'fully supersonic rightgoing' hstar = hl ustar = ul elseif (sr(2).lt.0.d0) then c # fully supersonic and leftgoing if (verbose) write(6,*) 'fully supersonic leftgoing' hstar = hr ustar = ur else c # 1-wave is leftgoing and c # 2-wave is rightgoing if (verbose) write(6,*) 'qstar = qm' hstar = hm ustar = um endif cstar = dsqrt(grav*hstar) c 400 continue c c ---------------------------------------------- c STEP #4: Compute flux fstar and fluctuations c ---------------------------------------------- c c # amdq = SUM s*wave over left-going waves c # apdq = SUM s*wave over right-going waves c c # Flux f at ql, qr, and qstar fl(1) = ul*hl fr(1) = ur*hr fstar(1) = ustar*hstar fl(2) = hl*ul**2 + 0.5d0*grav*hl**2 fr(2) = hr*ur**2 + 0.5d0*grav*hr**2 fstar(2) = hstar*ustar**2 + 0.5d0*grav*hstar**2 c c # fluctuations: apdq(i,1) = fr(1) - fstar(1) apdq(i,2) = fr(2) - fstar(2) amdq(i,1) = fstar(1) - fl(1) amdq(i,2) = fstar(2) - fl(2) c c c --------------------------------------- c STEP #5: Compute waves and their speeds c --------------------------------------- if (verbose) then write(6,*) 'sl(1), sr(1): ', sl(1),sr(1) write(6,*) 'sl(2), sr(2): ', sl(2),sr(2) endif c # 1-wave wave(i,1,1) = hm - hl wave(i,2,1) = hm*um - hl*ul s(i,1) = 0.5d0*( sl(1) + sr(1) ) c # 2-wave wave(i,1,2) = hr - hm wave(i,2,2) = ur*hr - um*hm s(i,2) = 0.5d0*( sl(2) + sr(2) ) if (verbose) then write(6,*) 'hstar, ustar:',hstar,ustar write(6,*) 'amdq:',amdq(i,1),amdq(i,2) write(6,*) 'apdq:',apdq(i,1),apdq(i,2) write(6,*) 's1 and 1-wave:',s(i,1),wave(i,1,1),wave(i,2,1) write(6,*) 's2 and 2-wave:',s(i,2),wave(i,1,2),wave(i,2,2) endif c 30 continue c return end