diff --git a/src/2d/bouss/Makefile.bouss b/src/2d/bouss/Makefile.bouss index 225e7a539..07a4e6c01 100644 --- a/src/2d/bouss/Makefile.bouss +++ b/src/2d/bouss/Makefile.bouss @@ -167,6 +167,6 @@ COMMON_SOURCES += \ $(BOUSSLIB)/resetBoussStuff.f \ $(BOUSSLIB)/simpleBound.f90 \ $(BOUSSLIB)/umfpack_support.f \ - $(BOUSSLIB)/rpn2_geoclaw.f \ - $(BOUSSLIB)/rpt2_geoclaw_sym.f \ - $(BOUSSLIB)/geoclaw_riemann_utils.f \ + $(CLAW)/riemann/src/rpn2_geoclaw.f \ + $(CLAW)/riemann/src/rpt2_geoclaw.f \ + $(CLAW)/riemann/src/geoclaw_riemann_utils.f \ diff --git a/src/2d/bouss/geoclaw_riemann_utils.f b/src/2d/bouss/geoclaw_riemann_utils.f deleted file mode 100644 index 4745cf43d..000000000 --- a/src/2d/bouss/geoclaw_riemann_utils.f +++ /dev/null @@ -1,658 +0,0 @@ -c----------------------------------------------------------------------- - subroutine riemann_aug_JCP(maxiter,meqn,mwaves,hL,hR,huL,huR, - & hvL,hvR,bL,bR,uL,uR,vL,vR,phiL,phiR,pL,pR,sE1,sE2,drytol,g,rho, - & 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,pL,pR - double precision drytol,g,rho - - - !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 criticaltol_2, hustar_interface - double precision s1s2bar,s1s2tilde,hbar,hLstar,hRstar,hustar - double precision huRstar,huLstar,uRstar,uLstar,hstarHLL - double precision deldelh,deldelphi,delP - 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 - delP = pR - pL - 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 - ! MODIFIED: - criticaltol = max(drytol*g, 1d-6) - criticaltol_2 = sqrt(criticaltol) - deldelh = -delb - deldelphi = -0.5d0 * (hR + hL) * (g * delb + delp / rho) - -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 - ! MODIFIED from 5.3.1 version - sonic = .false. - if (abs(s1s2bar) <= criticaltol) then - sonic = .true. - else if (s1s2bar*s1s2tilde <= criticaltol**2) then - sonic = .true. - else if (s1s2bar*sE1*sE2 <= criticaltol**2) then - sonic = .true. - else if (min(abs(sE1),abs(sE2)) < criticaltol_2) then - sonic = .true. - else if (sE1 < criticaltol_2 .and. s1m > -criticaltol_2) then - sonic = .true. - else if (sE2 > -criticaltol_2 .and. s2m < criticaltol_2) then - sonic = .true. - else if ((uL+dsqrt(g*hL))*(uR+dsqrt(g*hR)) < 0.d0) then - sonic = .true. - else if ((uL- dsqrt(g*hL))*(uR- dsqrt(g*hR)) < 0.d0) then - sonic = .true. - end if - -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 ! adjust deldelh for well-balancing of atmospheric pressure difference - deldelh = deldelh - delP/(rho*g) - -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)) - deldelphi = deldelphi - hbar * delp / rho - - 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 - - fw(:,:) = 0.d0 ! initialize before setting some waves - - 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). - ! MODIFIED from 5.3.1 version - fw(3,1)=fw(3,1)*vL - fw(3,3)=fw(3,3)*vR - fw(3,2)= 0.d0 - - hustar_interface = huL + fw(1,1) ! = huR - fw(1,3) - if (hustar_interface <= 0.0d0) then - fw(3,1) = fw(3,1) + (hR*uR*vR - hL*uL*vL - fw(3,1)- fw(3,3)) - else - fw(3,3) = fw(3,3) + (hR*uR*vR - hL*uL*vL - fw(3,1)- fw(3,3)) - end if - - - 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,pL,pR,sE1,sE2,drytol,g, - & rho,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,pL,pR - double precision drytol,g,rho - - !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,delP - 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 - delP = pR - pL - - convergencetol= 1.d-16 - criticaltol = 1.d-99 - - deldelh = -delb - deldelphi = -0.5d0 * (hR + hL) * (g * delb + delP / rho) - -! !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,pL,pR,s1,s2,drytol,g,rho, - & 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,pL,pR - double precision drytol,g,rho - - double precision sw(mwaves) - double precision fw(meqn,mwaves) - - !local - double precision delh,delhu,delphi,delb,delhdecomp,delphidecomp - double precision deldelh,deldelphi,delP - double precision beta1,beta2 - - - !determine del vectors - delh = hR-hL - delhu = huR-huL - delphi = phiR-phiL - delb = bR-bL - delP = pR - pL - - deldelphi = -0.5d0 * (hR + hL) * (g * delb + delP / rho) - 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/2d/bouss/rpn2_geoclaw.f b/src/2d/bouss/rpn2_geoclaw.f deleted file mode 100644 index 076fcd79f..000000000 --- a/src/2d/bouss/rpn2_geoclaw.f +++ /dev/null @@ -1,305 +0,0 @@ -c====================================================================== - subroutine rpn2(ixy,maxm,meqn,mwaves,maux,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, rho - use geoclaw_module, only: earth_radius, deg2rad - use amr_module, only: mcapa - - use storm_module, only: pressure_forcing, pressure_index - - 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(meqn,3) - double precision sw(3) - - double precision hR,hL,huR,huL,uR,uL,hvR,hvL,vR,vL,phiR,phiL,pL,pR - 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 :: debug - - logical rare1,rare2 - - ! In case there is no pressure forcing - pL = 0.d0 - pR = 0.d0 - debug = .false. - - ! initialize all components to 0 - fw(:,:) = 0.d0 - fwave(:,:,:) = 0.d0 - s(:,:) = 0.d0 - amdq(:,:) = 0.d0 - apdq(:,:) = 0.d0 - - !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 - - -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) - if (pressure_forcing) then - pL = auxr(pressure_index, i-1) - pR = auxl(pressure_index, i) - end if - - 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,meqn,mwaves,hL,hR,huL, - & huR,hvL,hvR,bL,bR,uL,uR,vL,vR,phiL,phiR,pL,pR,sE1,sE2, - & drytol,g,rho,sw,fw) - -C call riemann_ssqfwave(maxiter,meqn,mwaves,hL,hR,huL,huR, -C & hvL,hvR,bL,bR,uL,uR,vL,vR,phiL,phiR,pL,pR,sE1,sE2,drytol,g, -C & rho,sw,fw) - -C call riemann_fwave(meqn,mwaves,hL,hR,huL,huR,hvL,hvR, -C & bL,bR,uL,uR,vL,vR,phiL,phiR,pL,pR,sE1,sE2,drytol,g,rho,sw, -C & 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============================================= - - 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 - - if (debug) then - 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("+++ ampdq ",2i4,2e25.15) - 152 format("+++ fwave ",8x,3e25.15) - enddo - enddo - endif - - return - end subroutine diff --git a/src/2d/bouss/rpt2_geoclaw_sym.f b/src/2d/bouss/rpt2_geoclaw_sym.f deleted file mode 100644 index ce51c0d3c..000000000 --- a/src/2d/bouss/rpt2_geoclaw_sym.f +++ /dev/null @@ -1,260 +0,0 @@ -! ===================================================== - subroutine rpt2(ixy,imp,maxm,meqn,mwaves,maux,mbc,mx, - & ql,qr,aux1,aux2,aux3,asdq,bmasdq,bpasdq) -! ===================================================== - use geoclaw_module, only: g => grav, tol => dry_tolerance - use geoclaw_module, only: coordinate_system,earth_radius,deg2rad - - implicit none -! -! # Riemann solver in the transverse direction using an einfeldt -! Jacobian. - -!-----------------------last modified 1/10/05---------------------- - - integer ixy,maxm,meqn,maux,mwaves,mbc,mx,imp - - 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(mwaves) - double precision r(meqn,mwaves) - double precision beta(mwaves) - 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,s1down,s3up - double precision delf1,delf2,delf3,dxdcd,dxdcu - double precision dxdcm,dxdcp,topo1,topo3,eta - - integer i,m,mw,mu,mv - - !write(83,*) 'rpt2, imp = ',imp - - ! initialize all components to 0: - bmasdq(:,:) = 0.d0 - bpasdq(:,:) = 0.d0 - - 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) - -!===========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 (imp.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 (imp.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================= -c vhat=(vr*dsqrt(hr))/(dsqrt(hr)+dsqrt(hl)) + -c & (vl*dsqrt(hl))/(dsqrt(hr)+dsqrt(hl)) -c -c uhat=(ur*dsqrt(hr))/(dsqrt(hr)+dsqrt(hl)) + -c & (ul*dsqrt(hl))/(dsqrt(hr)+dsqrt(hl)) -c hhat=(hr+hl)/2.d0 - - ! Note that we used left right states to define Roe averages, - ! which is consistent with those used in rpn2. - ! But now we are computing upgoing, downgoing waves either in - ! cell on left (if imp==1) or on right (if imp==2) so we - ! should perhaps use Roe averages based on values above or below, - ! but these aren't available. - - !roe1=vhat-dsqrt(g*hhat) - !roe3=vhat+dsqrt(g*hhat) - - ! modified to at least use hl,vl or hr,vr properly based on imp: - ! (since s1 and s3 are now vertical velocities, - ! it made no sense to use h,v in cell i-1 for downgoing - ! and cell i for upgoing) - - if (imp == 1) then - ! asdq is leftgoing, use q from cell i-1: - if (hl <= tol) go to 90 - s1 = vl-dsqrt(g*hl) - s3 = vl+dsqrt(g*hl) - s2 = vl - uhat = ul - hhat = hl - else - ! asdq is rightgoing, use q from cell i: - if (hr <= tol) go to 90 - s1 = vr-dsqrt(g*hr) - s3 = vr+dsqrt(g*hr) - s2 = vr - uhat = ur - hhat = hr - endif - - ! don't use Roe averages: - !s1=dmin1(roe1,s1down) - !s3=dmax1(roe3,s3up) - - !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) - - ! fixed bug in beta(2): uhat in place of s(2) - beta(1) = (s3*delf1/(s3-s1))-(delf3/(s3-s1)) - beta(2) = -uhat*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) = uhat ! fixed bug, uhat not 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) = uhat ! fixed bug, uhat not 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 ((abs(s(mw)) > 0.d0) .and. - & (abs(s(mw)) < 0.001d0*dsqrt(g*hhat))) then - ! split correction symmetrically if nearly zero - ! Note wave drops out if s(mw)==0 exactly, so don't split - bmasdq(1,i) =bmasdq(1,i) + - & 0.5d0*dxdcm*s(mw)*beta(mw)*r(1,mw) - bmasdq(mu,i)=bmasdq(mu,i)+ - & 0.5d0*dxdcm*s(mw)*beta(mw)*r(2,mw) - bmasdq(mv,i)=bmasdq(mv,i)+ - & 0.5d0*dxdcm*s(mw)*beta(mw)*r(3,mw) - bpasdq(1,i) =bpasdq(1,i) + - & 0.5d0*dxdcp*s(mw)*beta(mw)*r(1,mw) - bpasdq(mu,i)=bpasdq(mu,i)+ - & 0.5d0*dxdcp*s(mw)*beta(mw)*r(2,mw) - bpasdq(mv,i)=bpasdq(mv,i)+ - & 0.5d0*dxdcp*s(mw)*beta(mw)*r(3,mw) - elseif (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 - - !if ((i>3) .and. (i<6)) then - if (.false.) then - ! DEBUG - write(83,*) 'i = ',i - write(83,831) s(1),s(2),s(3) - write(83,831) beta(1),beta(2),beta(3) - do m=1,3 - write(83,831) r(m,1),r(m,2),r(m,3) - enddo - do m=1,3 - write(83,831) asdq(m,i), bmasdq(m,i), bpasdq(m,i) - 831 format(3d20.12) - enddo - endif -c======================================================================== - enddo ! do i loop -c - - -c - - return - end