diff --git a/src/geoclaw_riemann_utils.f b/src/geoclaw_riemann_utils.f new file mode 100644 index 00000000..eba41c1b --- /dev/null +++ b/src/geoclaw_riemann_utils.f @@ -0,0 +1,626 @@ +c----------------------------------------------------------------------- + subroutine riemann_aug_JCP(maxiter,meqn,mwaves,hL,hR,huL,huR, + & hvL,hvR,bL,bR,uL,uR,vL,vR,phiL,phiR,sE1,sE2,drytol,g,sw,fw) + + ! solve shallow water equations given single left and right states + ! This solver is described in J. Comput. Phys. (6): 3089-3113, March 2008 + ! Augmented Riemann Solvers for the Shallow Equations with Steady States and Inundation + + ! To use the original solver call with maxiter=1. + + ! This solver allows iteration when maxiter > 1. The iteration seems to help with + ! instabilities that arise (with any solver) as flow becomes transcritical over variable topo + ! due to loss of hyperbolicity. + + + + implicit none + + !input + integer meqn,mwaves,maxiter + double precision fw(meqn,mwaves) + double precision sw(mwaves) + double precision hL,hR,huL,huR,bL,bR,uL,uR,phiL,phiR,sE1,sE2 + double precision hvL,hvR,vL,vR + double precision drytol,g + + + !local + integer m,mw,k,iter + double precision A(3,3) + double precision r(3,3) + double precision lambda(3) + double precision del(3) + double precision beta(3) + + double precision delh,delhu,delphi,delb,delnorm + double precision rare1st,rare2st,sdelta,raremin,raremax + double precision criticaltol,convergencetol,raretol + double precision s1s2bar,s1s2tilde,hbar,hLstar,hRstar,hustar + double precision huRstar,huLstar,uRstar,uLstar,hstarHLL + double precision deldelh,deldelphi + double precision s1m,s2m,hm + double precision det1,det2,det3,determinant + + logical rare1,rare2,rarecorrector,rarecorrectortest,sonic + + !determine del vectors + delh = hR-hL + delhu = huR-huL + delphi = phiR-phiL + delb = bR-bL + delnorm = delh**2 + delphi**2 + + call riemanntype(hL,hR,uL,uR,hm,s1m,s2m,rare1,rare2, + & 1,drytol,g) + + + lambda(1)= min(sE1,s2m) !Modified Einfeldt speed + lambda(3)= max(sE2,s1m) !Modified Eindfeldt speed + sE1=lambda(1) + sE2=lambda(3) + lambda(2) = 0.d0 ! ### Fix to avoid uninitialized value in loop on mw -- Correct?? ### + + + hstarHLL = max((huL-huR+sE2*hR-sE1*hL)/(sE2-sE1),0.d0) ! middle state in an HLL solve + +c !determine the middle entropy corrector wave------------------------ + rarecorrectortest=.false. + rarecorrector=.false. + if (rarecorrectortest) then + sdelta=lambda(3)-lambda(1) + raremin = 0.5d0 + raremax = 0.9d0 + if (rare1.and.sE1*s1m.lt.0.d0) raremin=0.2d0 + if (rare2.and.sE2*s2m.lt.0.d0) raremin=0.2d0 + if (rare1.or.rare2) then + !see which rarefaction is larger + rare1st=3.d0*(sqrt(g*hL)-sqrt(g*hm)) + rare2st=3.d0*(sqrt(g*hR)-sqrt(g*hm)) + if (max(rare1st,rare2st).gt.raremin*sdelta.and. + & max(rare1st,rare2st).lt.raremax*sdelta) then + rarecorrector=.true. + if (rare1st.gt.rare2st) then + lambda(2)=s1m + elseif (rare2st.gt.rare1st) then + lambda(2)=s2m + else + lambda(2)=0.5d0*(s1m+s2m) + endif + endif + endif + if (hstarHLL.lt.min(hL,hR)/5.d0) rarecorrector=.false. + endif + +c ## Is this correct 2-wave when rarecorrector == .true. ?? + do mw=1,mwaves + r(1,mw)=1.d0 + r(2,mw)=lambda(mw) + r(3,mw)=(lambda(mw))**2 + enddo + if (.not.rarecorrector) then + lambda(2) = 0.5d0*(lambda(1)+lambda(3)) +c lambda(2) = max(min(0.5d0*(s1m+s2m),sE2),sE1) + r(1,2)=0.d0 + r(2,2)=0.d0 + r(3,2)=1.d0 + endif +c !--------------------------------------------------- + + +c !determine the steady state wave ------------------- + criticaltol = 1.d-6 + deldelh = -delb + deldelphi = -g*0.5d0*(hR+hL)*delb + +c !determine a few quanitites needed for steady state wave if iterated + hLstar=hL + hRstar=hR + uLstar=uL + uRstar=uR + huLstar=uLstar*hLstar + huRstar=uRstar*hRstar + + !iterate to better determine the steady state wave + convergencetol=1.d-6 + do iter=1,maxiter + !determine steady state wave (this will be subtracted from the delta vectors) + if (min(hLstar,hRstar).lt.drytol.and.rarecorrector) then + rarecorrector=.false. + hLstar=hL + hRstar=hR + uLstar=uL + uRstar=uR + huLstar=uLstar*hLstar + huRstar=uRstar*hRstar + lambda(2) = 0.5d0*(lambda(1)+lambda(3)) +c lambda(2) = max(min(0.5d0*(s1m+s2m),sE2),sE1) + r(1,2)=0.d0 + r(2,2)=0.d0 + r(3,2)=1.d0 + endif + + hbar = max(0.5d0*(hLstar+hRstar),0.d0) + s1s2bar = 0.25d0*(uLstar+uRstar)**2 - g*hbar + s1s2tilde= max(0.d0,uLstar*uRstar) - g*hbar + +c !find if sonic problem + sonic=.false. + if (abs(s1s2bar).le.criticaltol) sonic=.true. + if (s1s2bar*s1s2tilde.le.criticaltol) sonic=.true. + if (s1s2bar*sE1*sE2.le.criticaltol) sonic = .true. + if (min(abs(sE1),abs(sE2)).lt.criticaltol) sonic=.true. + if (sE1.lt.0.d0.and.s1m.gt.0.d0) sonic = .true. + if (sE2.gt.0.d0.and.s2m.lt.0.d0) sonic = .true. + if ((uL+sqrt(g*hL))*(uR+sqrt(g*hR)).lt.0.d0) sonic=.true. + if ((uL-sqrt(g*hL))*(uR-sqrt(g*hR)).lt.0.d0) sonic=.true. + +c !find jump in h, deldelh + if (sonic) then + deldelh = -delb + else + deldelh = delb*g*hbar/s1s2bar + endif +c !find bounds in case of critical state resonance, or negative states + if (sE1.lt.-criticaltol.and.sE2.gt.criticaltol) then + deldelh = min(deldelh,hstarHLL*(sE2-sE1)/sE2) + deldelh = max(deldelh,hstarHLL*(sE2-sE1)/sE1) + elseif (sE1.ge.criticaltol) then + deldelh = min(deldelh,hstarHLL*(sE2-sE1)/sE1) + deldelh = max(deldelh,-hL) + elseif (sE2.le.-criticaltol) then + deldelh = min(deldelh,hR) + deldelh = max(deldelh,hstarHLL*(sE2-sE1)/sE2) + endif + +c !find jump in phi, deldelphi + if (sonic) then + deldelphi = -g*hbar*delb + else + deldelphi = -delb*g*hbar*s1s2tilde/s1s2bar + endif +c !find bounds in case of critical state resonance, or negative states + deldelphi=min(deldelphi,g*max(-hLstar*delb,-hRstar*delb)) + deldelphi=max(deldelphi,g*min(-hLstar*delb,-hRstar*delb)) + + del(1)=delh-deldelh + del(2)=delhu + del(3)=delphi-deldelphi + +c !Determine determinant of eigenvector matrix======== + det1=r(1,1)*(r(2,2)*r(3,3)-r(2,3)*r(3,2)) + det2=r(1,2)*(r(2,1)*r(3,3)-r(2,3)*r(3,1)) + det3=r(1,3)*(r(2,1)*r(3,2)-r(2,2)*r(3,1)) + determinant=det1-det2+det3 + +c !solve for beta(k) using Cramers Rule================= + do k=1,3 + do mw=1,3 + A(1,mw)=r(1,mw) + A(2,mw)=r(2,mw) + A(3,mw)=r(3,mw) + enddo + A(1,k)=del(1) + A(2,k)=del(2) + A(3,k)=del(3) + det1=A(1,1)*(A(2,2)*A(3,3)-A(2,3)*A(3,2)) + det2=A(1,2)*(A(2,1)*A(3,3)-A(2,3)*A(3,1)) + det3=A(1,3)*(A(2,1)*A(3,2)-A(2,2)*A(3,1)) + beta(k)=(det1-det2+det3)/determinant + enddo + + !exit if things aren't changing + if (abs(del(1)**2+del(3)**2-delnorm).lt.convergencetol) exit + delnorm = del(1)**2+del(3)**2 + !find new states qLstar and qRstar on either side of interface + hLstar=hL + hRstar=hR + uLstar=uL + uRstar=uR + huLstar=uLstar*hLstar + huRstar=uRstar*hRstar + do mw=1,mwaves + if (lambda(mw).lt.0.d0) then + hLstar= hLstar + beta(mw)*r(1,mw) + huLstar= huLstar + beta(mw)*r(2,mw) + endif + enddo + do mw=mwaves,1,-1 + if (lambda(mw).gt.0.d0) then + hRstar= hRstar - beta(mw)*r(1,mw) + huRstar= huRstar - beta(mw)*r(2,mw) + endif + enddo + + if (hLstar.gt.drytol) then + uLstar=huLstar/hLstar + else + hLstar=max(hLstar,0.d0) + uLstar=0.d0 + endif + if (hRstar.gt.drytol) then + uRstar=huRstar/hRstar + else + hRstar=max(hRstar,0.d0) + uRstar=0.d0 + endif + + enddo ! end iteration on Riemann problem + + do mw=1,mwaves + sw(mw)=lambda(mw) + fw(1,mw)=beta(mw)*r(2,mw) + fw(2,mw)=beta(mw)*r(3,mw) + fw(3,mw)=beta(mw)*r(2,mw) + enddo + !find transverse components (ie huv jumps). + fw(3,1)=fw(3,1)*vL + fw(3,3)=fw(3,3)*vR + fw(3,2)= hR*uR*vR - hL*uL*vL - fw(3,1)- fw(3,3) + + return + + end !subroutine riemann_aug_JCP------------------------------------------------- + + +c----------------------------------------------------------------------- + subroutine riemann_ssqfwave(maxiter,meqn,mwaves,hL,hR,huL,huR, + & hvL,hvR,bL,bR,uL,uR,vL,vR,phiL,phiR,sE1,sE2,drytol,g,sw,fw) + + ! solve shallow water equations given single left and right states + ! steady state wave is subtracted from delta [q,f]^T before decomposition + + implicit none + + !input + integer meqn,mwaves,maxiter + + double precision hL,hR,huL,huR,bL,bR,uL,uR,phiL,phiR,sE1,sE2 + double precision vL,vR,hvL,hvR + double precision drytol,g + + !local + integer iter + + logical sonic + + double precision delh,delhu,delphi,delb,delhdecomp,delphidecomp + double precision s1s2bar,s1s2tilde,hbar,hLstar,hRstar,hustar + double precision uRstar,uLstar,hstarHLL + double precision deldelh,deldelphi + double precision alpha1,alpha2,beta1,beta2,delalpha1,delalpha2 + double precision criticaltol,convergencetol + double precision sL,sR + double precision uhat,chat,sRoe1,sRoe2 + + double precision sw(mwaves) + double precision fw(meqn,mwaves) + + !determine del vectors + delh = hR-hL + delhu = huR-huL + delphi = phiR-phiL + delb = bR-bL + + convergencetol= 1.d-16 + criticaltol = 1.d-99 + + deldelh = -delb + deldelphi = -g*0.5d0*(hR+hL)*delb + +! !if no source term, skip determining steady state wave + if (abs(delb).gt.0.d0) then +! + !determine a few quanitites needed for steady state wave if iterated + hLstar=hL + hRstar=hR + uLstar=uL + uRstar=uR + hstarHLL = max((huL-huR+sE2*hR-sE1*hL)/(sE2-sE1),0.d0) ! middle state in an HLL solve + + alpha1=0.d0 + alpha2=0.d0 + +! !iterate to better determine Riemann problem + do iter=1,maxiter + + !determine steady state wave (this will be subtracted from the delta vectors) + hbar = max(0.5d0*(hLstar+hRstar),0.d0) + s1s2bar = 0.25d0*(uLstar+uRstar)**2 - g*hbar + s1s2tilde= max(0.d0,uLstar*uRstar) - g*hbar + + +c !find if sonic problem + sonic=.false. + if (abs(s1s2bar).le.criticaltol) sonic=.true. + if (s1s2bar*s1s2tilde.le.criticaltol) sonic=.true. + if (s1s2bar*sE1*sE2.le.criticaltol) sonic = .true. + if (min(abs(sE1),abs(sE2)).lt.criticaltol) sonic=.true. + +c !find jump in h, deldelh + if (sonic) then + deldelh = -delb + else + deldelh = delb*g*hbar/s1s2bar + endif +! !bounds in case of critical state resonance, or negative states + if (sE1.lt.-criticaltol.and.sE2.gt.criticaltol) then + deldelh = min(deldelh,hstarHLL*(sE2-sE1)/sE2) + deldelh = max(deldelh,hstarHLL*(sE2-sE1)/sE1) + elseif (sE1.ge.criticaltol) then + deldelh = min(deldelh,hstarHLL*(sE2-sE1)/sE1) + deldelh = max(deldelh,-hL) + elseif (sE2.le.-criticaltol) then + deldelh = min(deldelh,hR) + deldelh = max(deldelh,hstarHLL*(sE2-sE1)/sE2) + endif + +c !find jump in phi, deldelphi + if (sonic) then + deldelphi = -g*hbar*delb + else + deldelphi = -delb*g*hbar*s1s2tilde/s1s2bar + endif +! !bounds in case of critical state resonance, or negative states + deldelphi=min(deldelphi,g*max(-hLstar*delb,-hRstar*delb)) + deldelphi=max(deldelphi,g*min(-hLstar*delb,-hRstar*delb)) + +!---------determine fwaves ------------------------------------------ + +! !first decomposition + delhdecomp = delh-deldelh + delalpha1 = (sE2*delhdecomp - delhu)/(sE2-sE1)-alpha1 + alpha1 = alpha1 + delalpha1 + delalpha2 = (delhu - sE1*delhdecomp)/(sE2-sE1)-alpha2 + alpha2 = alpha2 + delalpha2 + + !second decomposition + delphidecomp = delphi - deldelphi + beta1 = (sE2*delhu - delphidecomp)/(sE2-sE1) + beta2 = (delphidecomp - sE1*delhu)/(sE2-sE1) + + if ((delalpha2**2+delalpha1**2).lt.convergencetol**2) then + exit + endif +! + if (sE2.gt.0.d0.and.sE1.lt.0.d0) then + hLstar=hL+alpha1 + hRstar=hR-alpha2 +c hustar=huL+alpha1*sE1 + hustar = huL + beta1 + elseif (sE1.ge.0.d0) then + hLstar=hL + hustar=huL + hRstar=hR - alpha1 - alpha2 + elseif (sE2.le.0.d0) then + hRstar=hR + hustar=huR + hLstar=hL + alpha1 + alpha2 + endif +! + if (hLstar.gt.drytol) then + uLstar=hustar/hLstar + else + hLstar=max(hLstar,0.d0) + uLstar=0.d0 + endif +! + if (hRstar.gt.drytol) then + uRstar=hustar/hRstar + else + hRstar=max(hRstar,0.d0) + uRstar=0.d0 + endif + + enddo + endif + + delhdecomp = delh - deldelh + delphidecomp = delphi - deldelphi + + !first decomposition + alpha1 = (sE2*delhdecomp - delhu)/(sE2-sE1) + alpha2 = (delhu - sE1*delhdecomp)/(sE2-sE1) + + !second decomposition + beta1 = (sE2*delhu - delphidecomp)/(sE2-sE1) + beta2 = (delphidecomp - sE1*delhu)/(sE2-sE1) + + ! 1st nonlinear wave + fw(1,1) = alpha1*sE1 + fw(2,1) = beta1*sE1 + fw(3,1) = fw(1,1)*vL + ! 2nd nonlinear wave + fw(1,3) = alpha2*sE2 + fw(2,3) = beta2*sE2 + fw(3,3) = fw(1,3)*vR + ! advection of transverse wave + fw(1,2) = 0.d0 + fw(2,2) = 0.d0 + fw(3,2) = hR*uR*vR - hL*uL*vL -fw(3,1)-fw(3,3) + !speeds + sw(1)=sE1 + sw(2)=0.5d0*(sE1+sE2) + sw(3)=sE2 + + return + + end subroutine !------------------------------------------------- + + +c----------------------------------------------------------------------- + subroutine riemann_fwave(meqn,mwaves,hL,hR,huL,huR,hvL,hvR, + & bL,bR,uL,uR,vL,vR,phiL,phiR,s1,s2,drytol,g,sw,fw) + + ! solve shallow water equations given single left and right states + ! solution has two waves. + ! flux - source is decomposed. + + implicit none + + !input + integer meqn,mwaves + + double precision hL,hR,huL,huR,bL,bR,uL,uR,phiL,phiR,s1,s2 + double precision hvL,hvR,vL,vR + double precision drytol,g + + double precision sw(mwaves) + double precision fw(meqn,mwaves) + + !local + double precision delh,delhu,delphi,delb,delhdecomp,delphidecomp + double precision deldelh,deldelphi + double precision beta1,beta2 + + + !determine del vectors + delh = hR-hL + delhu = huR-huL + delphi = phiR-phiL + delb = bR-bL + + deldelphi = -g*0.5d0*(hR+hL)*delb + delphidecomp = delphi - deldelphi + + !flux decomposition + beta1 = (s2*delhu - delphidecomp)/(s2-s1) + beta2 = (delphidecomp - s1*delhu)/(s2-s1) + + sw(1)=s1 + sw(2)=0.5d0*(s1+s2) + sw(3)=s2 + ! 1st nonlinear wave + fw(1,1) = beta1 + fw(2,1) = beta1*s1 + fw(3,1) = beta1*vL + ! 2nd nonlinear wave + fw(1,3) = beta2 + fw(2,3) = beta2*s2 + fw(3,3) = beta2*vR + ! advection of transverse wave + fw(1,2) = 0.d0 + fw(2,2) = 0.d0 + fw(3,2) = hR*uR*vR - hL*uL*vL -fw(3,1)-fw(3,3) + return + + end !subroutine ------------------------------------------------- + + + + + +c============================================================================= + subroutine riemanntype(hL,hR,uL,uR,hm,s1m,s2m,rare1,rare2, + & maxiter,drytol,g) + + !determine the Riemann structure (wave-type in each family) + + + implicit none + + !input + double precision hL,hR,uL,uR,drytol,g + integer maxiter + + !output + double precision s1m,s2m + logical rare1,rare2 + + !local + double precision hm,u1m,u2m,um,delu + double precision h_max,h_min,h0,F_max,F_min,dfdh,F0,slope,gL,gR + integer iter + + + +c !Test for Riemann structure + + h_min=min(hR,hL) + h_max=max(hR,hL) + delu=uR-uL + + if (h_min.le.drytol) then + hm=0.d0 + um=0.d0 + s1m=uR+uL-2.d0*sqrt(g*hR)+2.d0*sqrt(g*hL) + s2m=uR+uL-2.d0*sqrt(g*hR)+2.d0*sqrt(g*hL) + if (hL.le.0.d0) then + rare2=.true. + rare1=.false. + else + rare1=.true. + rare2=.false. + endif + + else + F_min= delu+2.d0*(sqrt(g*h_min)-sqrt(g*h_max)) + F_max= delu + + & (h_max-h_min)*(sqrt(.5d0*g*(h_max+h_min)/(h_max*h_min))) + + if (F_min.gt.0.d0) then !2-rarefactions + + hm=(1.d0/(16.d0*g))* + & max(0.d0,-delu+2.d0*(sqrt(g*hL)+sqrt(g*hR)))**2 + um=sign(1.d0,hm)*(uL+2.d0*(sqrt(g*hL)-sqrt(g*hm))) + + s1m=uL+2.d0*sqrt(g*hL)-3.d0*sqrt(g*hm) + s2m=uR-2.d0*sqrt(g*hR)+3.d0*sqrt(g*hm) + + rare1=.true. + rare2=.true. + + elseif (F_max.le.0.d0) then !2 shocks + +c !root finding using a Newton iteration on sqrt(h)=== + h0=h_max + do iter=1,maxiter + gL=sqrt(.5d0*g*(1/h0 + 1/hL)) + gR=sqrt(.5d0*g*(1/h0 + 1/hR)) + F0=delu+(h0-hL)*gL + (h0-hR)*gR + dfdh=gL-g*(h0-hL)/(4.d0*(h0**2)*gL)+ + & gR-g*(h0-hR)/(4.d0*(h0**2)*gR) + slope=2.d0*sqrt(h0)*dfdh + h0=(sqrt(h0)-F0/slope)**2 + enddo + hm=h0 + u1m=uL-(hm-hL)*sqrt((.5d0*g)*(1/hm + 1/hL)) + u2m=uR+(hm-hR)*sqrt((.5d0*g)*(1/hm + 1/hR)) + um=.5d0*(u1m+u2m) + + s1m=u1m-sqrt(g*hm) + s2m=u2m+sqrt(g*hm) + rare1=.false. + rare2=.false. + + else !one shock one rarefaction + h0=h_min + + do iter=1,maxiter + F0=delu + 2.d0*(sqrt(g*h0)-sqrt(g*h_max)) + & + (h0-h_min)*sqrt(.5d0*g*(1/h0+1/h_min)) + slope=(F_max-F0)/(h_max-h_min) + h0=h0-F0/slope + enddo + + hm=h0 + if (hL.gt.hR) then + um=uL+2.d0*sqrt(g*hL)-2.d0*sqrt(g*hm) + s1m=uL+2.d0*sqrt(g*hL)-3.d0*sqrt(g*hm) + s2m=uL+2.d0*sqrt(g*hL)-sqrt(g*hm) + rare1=.true. + rare2=.false. + else + s2m=uR-2.d0*sqrt(g*hR)+3.d0*sqrt(g*hm) + s1m=uR-2.d0*sqrt(g*hR)+sqrt(g*hm) + um=uR-2.d0*sqrt(g*hR)+2.d0*sqrt(g*hm) + rare2=.true. + rare1=.false. + endif + endif + endif + + return + + end ! subroutine riemanntype---------------------------------------------------------------- + diff --git a/src/rpn2_geoclaw.f b/src/rpn2_geoclaw.f new file mode 100644 index 00000000..31a1cd11 --- /dev/null +++ b/src/rpn2_geoclaw.f @@ -0,0 +1,288 @@ +c====================================================================== + subroutine rpn2(ixy,maxm,meqn,maux,mwaves,mbc,mx,ql,qr,auxl,auxr, + & fwave,s,amdq,apdq) +c====================================================================== +c +c Solves normal Riemann problems for the 2D SHALLOW WATER equations +c with topography: +c # h_t + (hu)_x + (hv)_y = 0 # +c # (hu)_t + (hu^2 + 0.5gh^2)_x + (huv)_y = -ghb_x # +c # (hv)_t + (huv)_x + (hv^2 + 0.5gh^2)_y = -ghb_y # + +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 +c This data is along a slice in the x-direction if ixy=1 +c or the y-direction if ixy=2. + +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 routines, this routine is called with +c ql = qr +c +c +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! ! +! # This Riemann solver is for the shallow water equations. ! +! ! +! It allows the user to easily select a Riemann solver in ! +! riemannsolvers_geo.f. this routine initializes all the variables ! +! for the shallow water equations, accounting for wet dry boundary ! +! dry cells, wave speeds etc. ! +! ! +! David George, Vancouver WA, Feb. 2009 ! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + use geoclaw_module, only: g => grav, drytol => dry_tolerance + use geoclaw_module, only: earth_radius, deg2rad + use amr_module, only: mcapa + + implicit none + + !input + integer maxm,meqn,maux,mwaves,mbc,mx,ixy + + double precision fwave(meqn, mwaves, 1-mbc:maxm+mbc) + double precision s(mwaves, 1-mbc:maxm+mbc) + double precision ql(meqn, 1-mbc:maxm+mbc) + double precision qr(meqn, 1-mbc:maxm+mbc) + double precision apdq(meqn,1-mbc:maxm+mbc) + double precision amdq(meqn,1-mbc:maxm+mbc) + double precision auxl(maux,1-mbc:maxm+mbc) + double precision auxr(maux,1-mbc:maxm+mbc) + + !local only + integer m,i,mw,maxiter,mu,nv + double precision wall(3) + double precision fw(3,3) + double precision sw(3) + + double precision hR,hL,huR,huL,uR,uL,hvR,hvL,vR,vL,phiR,phiL + double precision bR,bL,sL,sR,sRoe1,sRoe2,sE1,sE2,uhat,chat + double precision s1m,s2m + double precision hstar,hstartest,hstarHLL,sLtest,sRtest + double precision tw,dxdc + + logical rare1,rare2 + + !loop through Riemann problems at each grid cell + do i=2-mbc,mx+mbc + +!-----------------------Initializing----------------------------------- + !inform of a bad riemann problem from the start + if((qr(1,i-1).lt.0.d0).or.(ql(1,i) .lt. 0.d0)) then + write(*,*) 'Negative input: hl,hr,i=',qr(1,i-1),ql(1,i),i + endif + + !Initialize Riemann problem for grid interface + do mw=1,mwaves + s(mw,i)=0.d0 + fwave(1,mw,i)=0.d0 + fwave(2,mw,i)=0.d0 + fwave(3,mw,i)=0.d0 + enddo + +c !set normal direction + if (ixy.eq.1) then + mu=2 + nv=3 + else + mu=3 + nv=2 + endif + + !zero (small) negative values if they exist + if (qr(1,i-1).lt.0.d0) then + qr(1,i-1)=0.d0 + qr(2,i-1)=0.d0 + qr(3,i-1)=0.d0 + endif + + if (ql(1,i).lt.0.d0) then + ql(1,i)=0.d0 + ql(2,i)=0.d0 + ql(3,i)=0.d0 + endif + + !skip problem if in a completely dry area + if (qr(1,i-1) <= drytol .and. ql(1,i) <= drytol) then + go to 30 + endif + + !Riemann problem variables + hL = qr(1,i-1) + hR = ql(1,i) + huL = qr(mu,i-1) + huR = ql(mu,i) + bL = auxr(1,i-1) + bR = auxl(1,i) + + hvL=qr(nv,i-1) + hvR=ql(nv,i) + + !check for wet/dry boundary + if (hR.gt.drytol) then + uR=huR/hR + vR=hvR/hR + phiR = 0.5d0*g*hR**2 + huR**2/hR + else + hR = 0.d0 + huR = 0.d0 + hvR = 0.d0 + uR = 0.d0 + vR = 0.d0 + phiR = 0.d0 + endif + + if (hL.gt.drytol) then + uL=huL/hL + vL=hvL/hL + phiL = 0.5d0*g*hL**2 + huL**2/hL + else + hL=0.d0 + huL=0.d0 + hvL=0.d0 + uL=0.d0 + vL=0.d0 + phiL = 0.d0 + endif + + wall(1) = 1.d0 + wall(2) = 1.d0 + wall(3) = 1.d0 + if (hR.le.drytol) then + call riemanntype(hL,hL,uL,-uL,hstar,s1m,s2m, + & rare1,rare2,1,drytol,g) + hstartest=max(hL,hstar) + if (hstartest+bL.lt.bR) then !right state should become ghost values that mirror left for wall problem +c bR=hstartest+bL + wall(2)=0.d0 + wall(3)=0.d0 + hR=hL + huR=-huL + bR=bL + phiR=phiL + uR=-uL + vR=vL + elseif (hL+bL.lt.bR) then + bR=hL+bL + endif + elseif (hL.le.drytol) then ! right surface is lower than left topo + call riemanntype(hR,hR,-uR,uR,hstar,s1m,s2m, + & rare1,rare2,1,drytol,g) + hstartest=max(hR,hstar) + if (hstartest+bR.lt.bL) then !left state should become ghost values that mirror right +c bL=hstartest+bR + wall(1)=0.d0 + wall(2)=0.d0 + hL=hR + huL=-huR + bL=bR + phiL=phiR + uL=-uR + vL=vR + elseif (hR+bR.lt.bL) then + bL=hR+bR + endif + endif + + !determine wave speeds + sL=uL-sqrt(g*hL) ! 1 wave speed of left state + sR=uR+sqrt(g*hR) ! 2 wave speed of right state + + uhat=(sqrt(g*hL)*uL + sqrt(g*hR)*uR)/(sqrt(g*hR)+sqrt(g*hL)) ! Roe average + chat=sqrt(g*0.5d0*(hR+hL)) ! Roe average + sRoe1=uhat-chat ! Roe wave speed 1 wave + sRoe2=uhat+chat ! Roe wave speed 2 wave + + sE1 = min(sL,sRoe1) ! Eindfeldt speed 1 wave + sE2 = max(sR,sRoe2) ! Eindfeldt speed 2 wave + + !--------------------end initializing...finally---------- + !solve Riemann problem. + + maxiter = 1 + + call riemann_aug_JCP(maxiter,3,3,hL,hR,huL, + & huR,hvL,hvR,bL,bR,uL,uR,vL,vR,phiL,phiR,sE1,sE2, + & drytol,g,sw,fw) + +c call riemann_ssqfwave(maxiter,meqn,mwaves,hL,hR,huL,huR, +c & hvL,hvR,bL,bR,uL,uR,vL,vR,phiL,phiR,sE1,sE2,drytol,g,sw,fw) + +c call riemann_fwave(meqn,mwaves,hL,hR,huL,huR,hvL,hvR, +c & bL,bR,uL,uR,vL,vR,phiL,phiR,sE1,sE2,drytol,g,sw,fw) + +c !eliminate ghost fluxes for wall + do mw=1,3 + sw(mw)=sw(mw)*wall(mw) + + fw(1,mw)=fw(1,mw)*wall(mw) + fw(2,mw)=fw(2,mw)*wall(mw) + fw(3,mw)=fw(3,mw)*wall(mw) + enddo + + do mw=1,mwaves + s(mw,i)=sw(mw) + fwave(1,mw,i)=fw(1,mw) + fwave(mu,mw,i)=fw(2,mw) + fwave(nv,mw,i)=fw(3,mw) +! write(51,515) sw(mw),fw(1,mw),fw(2,mw),fw(3,mw) +!515 format("++sw",4e25.16) + enddo + + 30 continue + enddo + + +c==========Capacity for mapping from latitude longitude to physical space==== + if (mcapa.gt.0) then + do i=2-mbc,mx+mbc + if (ixy.eq.1) then + dxdc=(earth_radius*deg2rad) + else + dxdc=earth_radius*cos(auxl(3,i))*deg2rad + endif + + do mw=1,mwaves +c if (s(mw,i) .gt. 316.d0) then +c # shouldn't happen unless h > 10 km! +c write(6,*) 'speed > 316: i,mw,s(mw,i): ',i,mw,s(mw,i) +c endif + s(mw,i)=dxdc*s(mw,i) + fwave(1,mw,i)=dxdc*fwave(1,mw,i) + fwave(2,mw,i)=dxdc*fwave(2,mw,i) + fwave(3,mw,i)=dxdc*fwave(3,mw,i) + enddo + enddo + endif + +c=============================================================================== + + +c============= compute fluctuations============================================= + amdq(1:3,:) = 0.d0 + apdq(1:3,:) = 0.d0 + do i=2-mbc,mx+mbc + do mw=1,mwaves + if (s(mw,i) < 0.d0) then + amdq(1:3,i) = amdq(1:3,i) + fwave(1:3,mw,i) + else if (s(mw,i) > 0.d0) then + apdq(1:3,i) = apdq(1:3,i) + fwave(1:3,mw,i) + else + amdq(1:3,i) = amdq(1:3,i) + 0.5d0 * fwave(1:3,mw,i) + apdq(1:3,i) = apdq(1:3,i) + 0.5d0 * fwave(1:3,mw,i) + endif + enddo + enddo +!-- do i=2-mbc,mx+mbc +!-- do m=1,meqn +!-- write(51,151) m,i,amdq(m,i),apdq(m,i) +!-- write(51,152) fwave(m,1,i),fwave(m,2,i),fwave(m,3,i) +!--151 format("++3 ampdq ",2i4,2e25.15) +!--152 format("++3 fwave ",8x,3e25.15) +!-- enddo +!-- enddo + + return + end subroutine diff --git a/src/rpt2_geoclaw.f b/src/rpt2_geoclaw.f new file mode 100644 index 00000000..ffebbeaf --- /dev/null +++ b/src/rpt2_geoclaw.f @@ -0,0 +1,195 @@ +c +c +c ===================================================== + subroutine rpt2(ixy,maxm,meqn,maux,mwaves,mbc,mx, + & ql,qr,aux1,aux2,aux3, + & ilr,asdq,bmasdq,bpasdq) +c ===================================================== + use geoclaw_module, only: g => grav, tol => dry_tolerance + use geoclaw_module, only: coordinate_system,earth_radius,deg2rad + + implicit none +c +c # Riemann solver in the transverse direction using an einfeldt +c Jacobian. + +c-----------------------last modified 1/10/05---------------------- + + integer ixy,maxm,meqn,maux,mwaves,mbc,mx,ilr + + double precision ql(meqn,1-mbc:maxm+mbc) + double precision qr(meqn,1-mbc:maxm+mbc) + double precision asdq(meqn,1-mbc:maxm+mbc) + double precision bmasdq(meqn,1-mbc:maxm+mbc) + double precision bpasdq(meqn,1-mbc:maxm+mbc) + double precision aux1(maux,1-mbc:maxm+mbc) + double precision aux2(maux,1-mbc:maxm+mbc) + double precision aux3(maux,1-mbc:maxm+mbc) + + double precision s(3) + double precision r(3,3) + double precision beta(3) + double precision abs_tol + double precision hl,hr,hul,hur,hvl,hvr,vl,vr,ul,ur,bl,br + double precision uhat,vhat,hhat,roe1,roe3,s1,s2,s3,s1l,s3r + double precision delf1,delf2,delf3,dxdcd,dxdcu + double precision dxdcm,dxdcp,topo1,topo3,eta + + integer i,m,mw,mu,mv + + abs_tol=tol + + if (ixy.eq.1) then + mu = 2 + mv = 3 + else + mu = 3 + mv = 2 + endif + + do i=2-mbc,mx+mbc + + hl=qr(1,i-1) + hr=ql(1,i) + hul=qr(mu,i-1) + hur=ql(mu,i) + hvl=qr(mv,i-1) + hvr=ql(mv,i) + +c===========determine velocity from momentum=========================== + if (hl.lt.abs_tol) then + hl=0.d0 + ul=0.d0 + vl=0.d0 + else + ul=hul/hl + vl=hvl/hl + endif + + if (hr.lt.abs_tol) then + hr=0.d0 + ur=0.d0 + vr=0.d0 + else + ur=hur/hr + vr=hvr/hr + endif + + do mw=1,mwaves + s(mw)=0.d0 + beta(mw)=0.d0 + do m=1,meqn + r(m,mw)=0.d0 + enddo + enddo + dxdcp = 1.d0 + dxdcm = 1.d0 + + if (hl <= tol .and. hr <= tol) go to 90 + +* !check and see if cell that transverse waves are going in is high and dry + if (ilr.eq.1) then + eta = qr(1,i-1) + aux2(1,i-1) + topo1 = aux1(1,i-1) + topo3 = aux3(1,i-1) +c s1 = vl-sqrt(g*hl) +c s3 = vl+sqrt(g*hl) +c s2 = 0.5d0*(s1+s3) + else + eta = ql(1,i) + aux2(1,i) + topo1 = aux1(1,i) + topo3 = aux3(1,i) +c s1 = vr-sqrt(g*hr) +c s3 = vr+sqrt(g*hr) +c s2 = 0.5d0*(s1+s3) + endif + if (eta.lt.max(topo1,topo3)) go to 90 + + if (coordinate_system.eq.2) then + if (ixy.eq.2) then + dxdcp=(earth_radius*deg2rad) + dxdcm = dxdcp + else + if (ilr.eq.1) then + dxdcp = earth_radius*cos(aux3(3,i-1))*deg2rad + dxdcm = earth_radius*cos(aux1(3,i-1))*deg2rad + else + dxdcp = earth_radius*cos(aux3(3,i))*deg2rad + dxdcm = earth_radius*cos(aux1(3,i))*deg2rad + endif + endif + endif + +c=====Determine some speeds necessary for the Jacobian================= + vhat=(vr*dsqrt(hr))/(dsqrt(hr)+dsqrt(hl)) + + & (vl*dsqrt(hl))/(dsqrt(hr)+dsqrt(hl)) + + uhat=(ur*dsqrt(hr))/(dsqrt(hr)+dsqrt(hl)) + + & (ul*dsqrt(hl))/(dsqrt(hr)+dsqrt(hl)) + hhat=(hr+hl)/2.d0 + + roe1=vhat-dsqrt(g*hhat) + roe3=vhat+dsqrt(g*hhat) + + s1l=vl-dsqrt(g*hl) + s3r=vr+dsqrt(g*hr) + + s1=dmin1(roe1,s1l) + s3=dmax1(roe3,s3r) + + s2=0.5d0*(s1+s3) + + s(1)=s1 + s(2)=s2 + s(3)=s3 +c=======================Determine asdq decomposition (beta)============ + delf1=asdq(1,i) + delf2=asdq(mu,i) + delf3=asdq(mv, i) + + beta(1) = (s3*delf1/(s3-s1))-(delf3/(s3-s1)) + beta(2) = -s2*delf1 + delf2 + beta(3) = (delf3/(s3-s1))-(s1*delf1/(s3-s1)) +c======================End ================================================= + +c=====================Set-up eigenvectors=================================== + r(1,1) = 1.d0 + r(2,1) = s2 + r(3,1) = s1 + + r(1,2) = 0.d0 + r(2,2) = 1.d0 + r(3,2) = 0.d0 + + r(1,3) = 1.d0 + r(2,3) = s2 + r(3,3) = s3 +c============================================================================ +90 continue +c============= compute fluctuations========================================== + + bmasdq(1,i)=0.0d0 + bpasdq(1,i)=0.0d0 + bmasdq(2,i)=0.0d0 + bpasdq(2,i)=0.0d0 + bmasdq(3,i)=0.0d0 + bpasdq(3,i)=0.0d0 + do mw=1,3 + if (s(mw).lt.0.d0) then + bmasdq(1,i) =bmasdq(1,i) + dxdcm*s(mw)*beta(mw)*r(1,mw) + bmasdq(mu,i)=bmasdq(mu,i)+ dxdcm*s(mw)*beta(mw)*r(2,mw) + bmasdq(mv,i)=bmasdq(mv,i)+ dxdcm*s(mw)*beta(mw)*r(3,mw) + elseif (s(mw).gt.0.d0) then + bpasdq(1,i) =bpasdq(1,i) + dxdcp*s(mw)*beta(mw)*r(1,mw) + bpasdq(mu,i)=bpasdq(mu,i)+ dxdcp*s(mw)*beta(mw)*r(2,mw) + bpasdq(mv,i)=bpasdq(mv,i)+ dxdcp*s(mw)*beta(mw)*r(3,mw) + endif + enddo +c======================================================================== + enddo +c + +c + + return + end