Update to JCP GeoClaw Riemann solver

[rjleveque]

Some modifications proposed by Wenwen Li to deal with sonic point better and to properly handle transverse velocities. I've been using this version on several problems recently and it seems fine, and works better in some cases.

Note: These changes will cause some geoclaw regression tests to fail, so the regression data should be updated in geoclaw.

This may help with certain instability issues as in clawpack/geoclaw#209.

Some notes from her email to me:

---

GeoClaw uses the same threshold value ‘criticaltol’ for all criteria; however, ‘s1s2bar’ and ‘s1s2tilde’ have the unit of (m/s)^2, and sE1, sE2, s1m, s2m, uL and uR have the unit of m/s. So I think it might be better to use different threshold values for different lines. I propose to set the threshold values based on the value of sqrt(drytol*g)

---

I had been noticing that there were horizontal/vertical lines in velocity and momentum plots for some tsunami scenarios --

When waves propagate from a deep cell towards a shallow cell (assuming huL > 0, huR > 0, hvL ≠ 0, hvR = 0), at the topo-jump interface, it is possible that s(1) < 0, s(2) < 0 and s(3) > 0. In GeoClaw Riemann solver, the hv component of all three fluxes Fw(1), Fw(2) and Fw(3) were calculated as below (a):

(a)
fw(3,1)=fw(3,1)*vL                            ! hv component of Fw(1)
fw(3,3)=fw(3,3)*vR                           ! hv component of Fw(3)
fw(3,2)= hR*uR*vR - hL*uL*vL - fw(3,1)- fw(3,3)      ! hv component of Fw(2)

If vR==0, then the hv component of Fw(3) will always be zero, which means the hv component of the left cell cannot affect the right cell even if there should be mass transfer from left to right (since uL, vL > 0). Furthermore, as long as s(2) < 0, the value of hvL won’t affect the hv component of Fw(3), so it won’t affect the hv component of the right cell at all, even if both uL and uR could be rightward. Thus the hv component will be trapped and accumulated at the left cell since there is incoming hv flux from the left but no outgoing hv flux.

The above issue is due to s(2) being not always in the same direction of real mass flux at the interface. Since s(2) is closely related to the calculations of flux for h and hu components, instead of changing s(2) itself, we could give a better estimation of real mass flux at the interface then distribute the hv flux component between fw(1) and fw(3) only.

In GeoClaw Riemann solver, the h flux component is zero for mid-wave (fw(2).h = 0) and Fw(1).h + F2(3).h = huR – huL. (Fw(1).h is the mass growth rate at the left cell.) According to the mass conservation, the real mass excessive flux at the interface should be
Fw(real_mass) = huL + Fw(1).h = huR - Fw(3).h
So we may change the code (a) to (b):

(b)
fw(3,1)=fw(3,1)*vL
fw(3,3)=fw(3,3)*vR
fw(3,2)=0d0

hustar_interface = huL + fw(1,1)   ! or 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

[mandli]

I tested this out and as expected we get slightly different results and the tests no longer pass with
[3.092074, 2.04987] vs. [ 3.090622, 2.053757] and [1.322509158600000e-02, 1.054953019800000e-02] vs. [1.323738840200000e-02, 1.053218228800000e-02] (full exceptions below). On inspection I do not see anything different so we may just need to regenerate the regression data.

======================================================================
FAIL: Test bowl-slosh example
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/Users/mandli/src/clawpack/geoclaw/tests/bowl_slosh/regression_tests.py", line 58, in runTest
    self.check_gauges(save=save, gauge_id=1, indices=(2, 3))
  File "/Users/mandli/src/clawpack/clawpack/clawutil/test.py", line 355, in check_gauges
    raise e
AssertionError: ('\nNot equal to tolerance rtol=1e-14, atol=1e-08\n\n(mismatch 100.0%)\n x: array([ 3.092074,  2.04987 ])\n y: array([ 3.090622,  2.053757])', 'Sum Check Failed for gauge = 1', ['3.092073650000000e+00', '2.049870060000000e+00'], ['3.090621760000000e+00', '2.053756900000000e+00'])

======================================================================
FAIL: Test multi-layer basic plane-waves.
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/Users/mandli/src/clawpack/geoclaw/tests/multilayer/regression_tests.py", line 63, in runTest
    self.check_gauges(save=save, gauge_id=0, indices=(3, 6))
  File "/Users/mandli/src/clawpack/clawpack/clawutil/test.py", line 355, in check_gauges
    raise e
AssertionError: ('\nNot equal to tolerance rtol=1e-14, atol=1e-08\n\n(mismatch 100.0%)\n x: array([ 0.013225,  0.01055 ])\n y: array([ 0.013237,  0.010532])', 'Sum Check Failed for gauge = 0', ['1.322509158600000e-02', '1.054953019800000e-02'], ['1.323738840200000e-02', '1.053218228800000e-02'])

[rjleveque]

@mandli: let's wait until clawpack/geoclaw#225 is passing Travis tests and merged in and then update the geoclaw regression tests to merge this.

[mandli]

+1

[mandli]

The regression data in GeoClaw has been updated in clawpack/geoclaw#228.

[rjleveque]

I will merge this and then update regression data in geoclaw.