From e1911dce0e6da0b3586fc00f789680b38672c051 Mon Sep 17 00:00:00 2001 From: Randy LeVeque Date: Sat, 4 Jun 2016 11:10:31 -0700 Subject: [PATCH 1/3] change sonic test to use different criticaltol as proposed by Wenwen Li --- src/geoclaw_riemann_utils.f | 25 +++++++++++++++---------- 1 file changed, 15 insertions(+), 10 deletions(-) diff --git a/src/geoclaw_riemann_utils.f b/src/geoclaw_riemann_utils.f index eba41c1b..3cae9711 100644 --- a/src/geoclaw_riemann_utils.f +++ b/src/geoclaw_riemann_utils.f @@ -36,6 +36,7 @@ subroutine riemann_aug_JCP(maxiter,meqn,mwaves,hL,hR,huL,huR, 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 @@ -109,7 +110,10 @@ subroutine riemann_aug_JCP(maxiter,meqn,mwaves,hL,hR,huL,huR, c !determine the steady state wave ------------------- - criticaltol = 1.d-6 + !criticaltol = 1.d-6 + ! MODIFIED: + criticaltol = max(drytol*g, 1d-6) + criticaltol_2 = dsqrt(criticaltol) deldelh = -delb deldelphi = -g*0.5d0*(hR+hL)*delb @@ -145,15 +149,16 @@ subroutine riemann_aug_JCP(maxiter,meqn,mwaves,hL,hR,huL,huR, 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. + ! MODIFIED from 5.3.1 version + sonic = + & ((abs(s1s2bar).le.criticaltol) .or. + & (s1s2bar*s1s2tilde.le.criticaltol*criticaltol) .or. + & (s1s2bar*sE1*sE2.le.criticaltol*criticaltol) .or. + & (min(abs(sE1),abs(sE2)).lt.criticaltol_2) .or. + & (sE1.lt.criticaltol_2.and.s1m.gt. -criticaltol_2) .or. + & (sE2.gt.-criticaltol_2.and.s2m.lt. criticaltol_2) .or. + & ((uL+dsqrt(g*hL))*(uR+dsqrt(g*hR)).lt.0.d0) .or. + & ((uL- dsqrt(g*hL))*(uR- dsqrt(g*hR)).lt.0.d0)) c !find jump in h, deldelh if (sonic) then From 9c02ede506c8a4ea43578b1eea1c38f5f9fade60 Mon Sep 17 00:00:00 2001 From: Randy LeVeque Date: Sat, 4 Jun 2016 11:16:17 -0700 Subject: [PATCH 2/3] modify fw(3,:) in riemann_aug_JCP geoclaw solver to put transverse velocity jump into fw(3,1) or fw(3,3) depending on interface velocity rather than into fw(3,2), as suggested by Wenwen Li. This avoids some cases where transverse velocity does not propagate past jump in bathymetry, may improve some instability issues. --- src/geoclaw_riemann_utils.f | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/src/geoclaw_riemann_utils.f b/src/geoclaw_riemann_utils.f index 3cae9711..065d9502 100644 --- a/src/geoclaw_riemann_utils.f +++ b/src/geoclaw_riemann_utils.f @@ -259,9 +259,18 @@ subroutine riemann_aug_JCP(maxiter,meqn,mwaves,hL,hR,huL,huR, 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)= hR*uR*vR - hL*uL*vL - fw(3,1)- fw(3,3) + 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 From b3a3430f208a8338ba4ee59b1af4fc125d754874 Mon Sep 17 00:00:00 2001 From: Kyle Mandli Date: Fri, 23 Sep 2016 17:15:19 -0400 Subject: [PATCH 3/3] Make sonic check more efficient --- src/geoclaw_riemann_utils.f | 30 +++++++++++++++++++----------- 1 file changed, 19 insertions(+), 11 deletions(-) diff --git a/src/geoclaw_riemann_utils.f b/src/geoclaw_riemann_utils.f index 065d9502..c6b7ea3d 100644 --- a/src/geoclaw_riemann_utils.f +++ b/src/geoclaw_riemann_utils.f @@ -113,7 +113,7 @@ subroutine riemann_aug_JCP(maxiter,meqn,mwaves,hL,hR,huL,huR, !criticaltol = 1.d-6 ! MODIFIED: criticaltol = max(drytol*g, 1d-6) - criticaltol_2 = dsqrt(criticaltol) + criticaltol_2 = sqrt(criticaltol) deldelh = -delb deldelphi = -g*0.5d0*(hR+hL)*delb @@ -150,15 +150,24 @@ subroutine riemann_aug_JCP(maxiter,meqn,mwaves,hL,hR,huL,huR, c !find if sonic problem ! MODIFIED from 5.3.1 version - sonic = - & ((abs(s1s2bar).le.criticaltol) .or. - & (s1s2bar*s1s2tilde.le.criticaltol*criticaltol) .or. - & (s1s2bar*sE1*sE2.le.criticaltol*criticaltol) .or. - & (min(abs(sE1),abs(sE2)).lt.criticaltol_2) .or. - & (sE1.lt.criticaltol_2.and.s1m.gt. -criticaltol_2) .or. - & (sE2.gt.-criticaltol_2.and.s2m.lt. criticaltol_2) .or. - & ((uL+dsqrt(g*hL))*(uR+dsqrt(g*hR)).lt.0.d0) .or. - & ((uL- dsqrt(g*hL))*(uR- dsqrt(g*hR)).lt.0.d0)) + 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 @@ -637,4 +646,3 @@ subroutine riemanntype(hL,hR,uL,uR,hm,s1m,s2m,rare1,rare2, return end ! subroutine riemanntype---------------------------------------------------------------- -