Skip to content

Commit

Permalink
Merge pull request #159 from rjleveque/rpt2_geoclaw
Browse files Browse the repository at this point in the history 
Modifications to rpt2_geoclaw
  • Loading branch information
@rjleveque
rjleveque authored Nov 8, 2020
2 parents 9ccef5b + e326696 commit 78c2236
Showing 1 changed file with 143 additions and 153 deletions.
296 changes: 143 additions & 153 deletions src/rpt2_geoclaw.f
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,191 +2,181 @@
subroutine rpt2(ixy,imp,maxm,meqn,mwaves,maux,mbc,mx,
& ql,qr,aux1,aux2,aux3,asdq,bmasdq,bpasdq)
! =====================================================
!
! Riemann solver in the transverse direction using
! Jacobian matrix from left cell (if imp==1) or right cell (if imp==2).
!
! Note this has been modified from the version used in v5.7.x and
! earlier, where Roe averages based on ql and qr were used, which
! is not correct. In addition:
! - a bug in the second component of the eigenvectors was fixed.
! - when s(2) is close to zero this component of flux difference
! is split equally between bmasdq and bpasdq to improve symmetry.
!
! Further modified to clean up and avoid a lot of work in dry cells.

!-----------------------last modified October 2020 ----------------------

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(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

integer, intent(in) :: ixy,maxm,meqn,maux,mwaves,mbc,mx,imp

real(kind=8), intent(in) :: ql(meqn,1-mbc:maxm+mbc)
real(kind=8), intent(in) :: qr(meqn,1-mbc:maxm+mbc)
real(kind=8), intent(in) :: asdq(meqn,1-mbc:maxm+mbc)
real(kind=8), intent(in) :: aux1(maux,1-mbc:maxm+mbc)
real(kind=8), intent(in) :: aux2(maux,1-mbc:maxm+mbc)
real(kind=8), intent(in) :: aux3(maux,1-mbc:maxm+mbc)

real(kind=8), intent(out) :: bmasdq(meqn,1-mbc:maxm+mbc)
real(kind=8), intent(out) :: bpasdq(meqn,1-mbc:maxm+mbc)

! local:
real(kind=8) :: s(mwaves), r(meqn,mwaves), beta(mwaves)
real(kind=8) :: h,u,v
real(kind=8) :: delf1,delf2,delf3
real(kind=8) :: dxdcm,dxdcp,topo1,topo3,eta

integer :: i,mw,mu,mv


if (ixy == 1) then
! normal solve was in x-direction
mu = 2
mv = 3
else
mu = 3
mv = 2
! normal solve was in y-direction
mu = 3
mv = 2
endif

! initialize all components of result to 0:
bmasdq(:,:) = 0.d0
bpasdq(:,:) = 0.d0


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
if (imp==1) then
h = qr(1,i-1)
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
h = ql(1,i)
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)
if (h <= tol) then
! fluctuation going into a dry cell, don't know how to split,
! so leave bmadsq(:,i)=bpasdq(:,i)=0 and go on to next i:
cycle
endif
s1l=vl-dsqrt(g*hl)
s3r=vr+dsqrt(g*hr)
! compute velocities in relevant cell, and other quantities:
s1=dmin1(roe1,s1l)
s3=dmax1(roe3,s3r)
if (imp==1) then
! fluctuation being split is left-going
u = qr(mu,i-1)/h
v = qr(mv,i-1)/h
eta = h + aux2(1,i-1)
topo1 = aux1(1,i-1)
topo3 = aux3(1,i-1)
else
! fluctuation being split is right-going
u = ql(mu,i)/h
v = ql(mv,i)/h
eta = h + aux2(1,i)
topo1 = aux1(1,i)
topo3 = aux3(1,i)
endif
s2=0.5d0*(s1+s3)
! check if cell that transverse waves go into are both too high:
! Note: prior to v5.8.0 this checked against max rather than min
if (eta < min(topo1,topo3)) cycle ! go to next i
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)
! if we get here, we want to do the splitting (no dry cells),
! so compute the necessary quantities:
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 =================================================
if (coordinate_system == 2) then
! on the sphere:
if (ixy == 2) then
dxdcp=(earth_radius*deg2rad)
dxdcm = dxdcp
else
if (imp == 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
else
! coordinate_system == 1 means Cartesian:
dxdcp = 1.d0
dxdcm = 1.d0
endif
c=====================Set-up eigenvectors===================================
c Determine some speeds necessary for the Jacobian
! In v5.7.x and prior versions,
! 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 possibly use q values in cells above/below,
! but these aren't available (only aux values).
! At any rate, there is no clear justification for using cells
! on the other side of the normal-solve interface.

! v5.8.0: modified to use left or right state alone in defining
! Jacobian, based on imp:

s(1) = v - dsqrt(g*h)
s(2) = v
s(3) = v + dsqrt(g*h)

c Determine asdq decomposition (beta)

delf1 = asdq(1,i)
delf2 = asdq(mu,i)
delf3 = asdq(mv, i)

! v5.8.0: fixed bug in beta(2): u in place of s(2)=v
beta(1) = (s(3)*delf1 - delf3) / (s(3) - s(1))
beta(2) = -u*delf1 + delf2
beta(3) = (delf3 - s(1)*delf1) / (s(3) - s(1))

c Set-up eigenvectors
r(1,1) = 1.d0
r(2,1) = s2
r(3,1) = s1
r(2,1) = u ! v5.8.0: fixed bug, u not s(2)=v
r(3,1) = s(1)

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
r(2,3) = u ! v5.8.0: fixed bug, u not s(2)=v
r(3,3) = s(3)


! compute fluctuations

do mw=1,3
if ((s(mw) < 0.d0) .and. (eta >= topo1)) 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
elseif ((s(mw) > 0.d0) .and. (eta >= topo3)) 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
endif
enddo ! loop on mw


c
enddo ! loop on i

return
end

0 comments on commit 78c2236

Please sign in to comment.