diff --git a/src/__init__.py b/src/__init__.py index 10705bf4..3a9bf4d3 100755 --- a/src/__init__.py +++ b/src/__init__.py @@ -63,6 +63,7 @@ from . import euler_hlle_1D from . import nonlinear_elasticity_fwave_1D from . import reactive_euler_with_efix_1D + from . import shallow_hlle_1D from . import shallow_roe_with_efix_1D from . import shallow_bathymetry_fwave_1D from . import shallow_roe_tracer_1D @@ -82,6 +83,7 @@ from . import euler_hlle_with_walls_2D from . import kpp_2D from . import psystem_2D + from . import shallow_hlle_2D from . import shallow_roe_with_efix_2D from . import shallow_bathymetry_fwave_2D from . import sw_aug_2D diff --git a/src/rp1_shallow_hlle.f90 b/src/rp1_shallow_hlle.f90 new file mode 100644 index 00000000..afc00424 --- /dev/null +++ b/src/rp1_shallow_hlle.f90 @@ -0,0 +1,78 @@ +! ========================================================= +subroutine rp1(maxmx,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,wave,s,amdq,apdq) +! ========================================================= + +! solve Riemann problem for the 1D shallow water equations using the HLLE +! approximate Riemann solver. + +! waves: 2 +! equations: 2 + +! Conserved quantities: +! 1 depth (h) +! 2 momentum (hu) + + implicit none + + integer, intent(in) :: maxmx, meqn, mwaves, mbc, mx, maux + double precision, dimension(meqn,1-mbc:maxmx+mbc), intent(in) :: ql, qr + double precision, dimension(maux,1-mbc:maxmx+mbc), intent(in) :: auxl, auxr + double precision, dimension(meqn, mwaves, 1-mbc:maxmx+mbc), intent(out) :: wave + double precision, dimension(meqn, 1-mbc:maxmx+mbc), intent(out) :: amdq, apdq + double precision, dimension(mwaves, 1-mbc:maxmx+mbc), intent(out) :: s + + double precision :: u_l, u_r, h_l, h_r, c_l, c_r + double precision :: hsqrt_l, hsqrt_r, u_hat, h_hat, c_hat, grav + double precision :: h_m, hu_m + double precision :: rho_m, rhou_m, E_m, s1, s2 + integer :: m, i, mw + + common /cparam/ grav + + do i=2-mbc,mx+mbc + h_l = qr(1,i-1) + h_r = ql(1,i) + u_l = qr(2,i-1) / qr(1,i-1) + u_r = ql(2,i ) / ql(1,i ) + c_l = dsqrt(grav*h_l) + c_r = dsqrt(grav*h_r) + + ! Roe averages + hsqrt_l = dsqrt(qr(1,i-1)) + hsqrt_r = dsqrt(ql(1,i)) + h_hat = 0.5*(h_l + h_r) + u_hat = (hsqrt_l*u_l + hsqrt_r*u_r) / (hsqrt_l + hsqrt_r) + c_hat = dsqrt(grav*h_hat) + + s1 = min(u_l - c_l, u_hat - c_hat) + s2 = max(u_r + c_r, u_hat + c_hat) + + ! Middle state + h_m = (ql(2,i) - qr(2,i-1) - s2*h_r + s1*h_l) / (s1-s2) + hu_m = (ql(2,i)*(u_r-s2) - qr(2,i-1)*(u_l-s1) + 0.5*grav*(h_r**2-h_l**2)) / (s1-s2) + + wave(1,1,i) = h_m - h_l + wave(2,1,i) = hu_m - qr(2,i-1) + s(1,i) = s1 + + wave(1,2,i) = h_r - h_m + wave(2,2,i) = ql(2,i) - hu_m + s(2,i) = s2 + + end do + + do m=1, meqn + do i=2-mbc, mx+mbc + amdq(m,i) = 0.d0 + apdq(m,i) = 0.d0 + do mw=1, mwaves + if (s(mw,i) < 0.d0) then + amdq(m,i) = amdq(m,i) + s(mw,i)*wave(m,mw,i) + else + apdq(m,i) = apdq(m,i) + s(mw,i)*wave(m,mw,i) + endif + end do + end do + end do + +end subroutine rp1 diff --git a/src/rpn2_shallow_hlle.f90 b/src/rpn2_shallow_hlle.f90 new file mode 100644 index 00000000..88bbed12 --- /dev/null +++ b/src/rpn2_shallow_hlle.f90 @@ -0,0 +1,102 @@ +! ===================================================== +subroutine rpn2(ixy,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,wave,s,amdq,apdq) +! ===================================================== + +! HLLE solver for the 2D shallow water equations. + +! waves: 2 +! equations: 3 + +! Conserved quantities: +! 1 depth +! 2 x-momentum +! 3 y-momentum + + implicit none + + integer, intent(in) :: ixy, maxm, meqn, mwaves, maux, mbc, mx + double precision, dimension(meqn,1-mbc:maxm+mbc), intent(in) :: ql, qr + double precision, dimension(maux,1-mbc:maxm+mbc), intent(in) :: auxl, auxr + double precision, dimension(meqn, mwaves, 1-mbc:maxm+mbc), intent(out) :: wave + double precision, dimension(mwaves, 1-mbc:maxm+mbc), intent(out) :: s + double precision, dimension(meqn, 1-mbc:maxm+mbc), intent(out) :: amdq, apdq + + double precision :: u_l, u_r, c_l, c_r, u_hat, c_hat, v_l, v_r, v_hat, h_hat + double precision :: h_l, h_r, hsqrt_l, hsqrt_r, hsq2 + double precision :: grav + double precision :: h_m, hu_m, hv_m, s1, s2 + integer :: depth, mu, mv + integer :: i, m, mw + + common /cparam/ grav + +! # set mu to point to the component of the system that corresponds +! # to momentum in the direction of this slice, mv to the orthogonal +! # momentum: +! + + depth = 1 + if (ixy.eq.1) then + mu = 2 + mv = 3 + else + mu = 3 + mv = 2 + endif + + do i=2-mbc,mx+mbc + h_l = qr(depth,i-1) + h_r = ql(depth,i) + ! Velocity + u_l = qr(mu,i-1) / qr(depth,i-1) + u_r = ql(mu,i ) / ql(depth,i ) + v_l = qr(mv,i-1) / qr(depth,i-1) + v_r = ql(mv,i ) / ql(depth,i ) + ! Sound speed + c_l = dsqrt(grav*h_l) + c_r = dsqrt(grav*h_r) + + hsqrt_l = dsqrt(qr(depth,i-1)) + hsqrt_r = dsqrt(ql(depth,i)) + hsq2 = hsqrt_l + hsqrt_r + h_hat = 0.5*(h_l + h_r) + u_hat = (u_l*hsqrt_l + u_r*hsqrt_r) / hsq2 + v_hat = (v_l*hsqrt_l + v_r*hsqrt_r) / hsq2 + c_hat = dsqrt(grav*h_hat) + + ! Speeds of non-shear waves + s1 = min(u_l - c_l, u_hat - c_hat) + s2 = max(u_r + c_r, u_hat + c_hat) + + ! middle" state + h_m = (ql(mu,i) - qr(mu,i-1) - s2*ql(depth,i) + s1*qr(depth,i-1))/(s1-s2) + hu_m = (ql(mu,i)*(u_r-s2) - qr(mu,i-1)*(u_l-s1) + 0.5*grav*(h_r**2 - h_l**2) ) / (s1-s2) + hv_m = (ql(mv,i)*u_r - qr(mv,i-1)*u_l - s2*ql(mv,i) + s1*qr(mv,i-1))/(s1-s2) + + wave(depth,1,i) = h_m - h_l + wave(mu,1,i) = hu_m - qr(mu,i-1) + wave(mv,1,i) = hv_m - qr(mv,i-1) + s(1,i) = s1 + + wave(depth,2,i) = h_r - h_m + wave(mu,2,i) = ql(mu,i) - hu_m + wave(mv,2,i) = ql(mv,i) - hv_m + s(2,i) = s2 + end do + + + do m=1, meqn + do i=2-mbc, mx+mbc + amdq(m,i) = 0.d0 + apdq(m,i) = 0.d0 + do mw=1, mwaves + if (s(mw,i) .lt. 0.d0) then + amdq(m,i) = amdq(m,i) + s(mw,i)*wave(m,mw,i) + else + apdq(m,i) = apdq(m,i) + s(mw,i)*wave(m,mw,i) + endif + end do + end do + end do + +end subroutine rpn2 diff --git a/src/setup.py b/src/setup.py index 9ca3fdfe..4f4aef0d 100644 --- a/src/setup.py +++ b/src/setup.py @@ -19,6 +19,7 @@ 'euler_hlle', 'nonlinear_elasticity_fwave', 'reactive_euler_with_efix', + 'shallow_hlle', 'shallow_roe_with_efix', 'shallow_bathymetry_fwave', 'shallow_roe_tracer'] @@ -91,6 +92,8 @@ def configuration(parent_package='',top_path=None): special_target_list = \ [{'ext' :'kpp_2D', 'srcs':['rpn2_kpp.f90','rpt2_dummy.f90']}, + {'ext' :'shallow_hlle_2D', + 'srcs':['rpn2_shallow_hlle.f90','rpt2_dummy.f90']}, {'ext' :'euler_hlle_2D', 'srcs':['rpn2_euler_hlle.f90','rpt2_dummy.f90']}, {'ext' :'euler_hlle_with_walls_2D', diff --git a/src/static.py b/src/static.py index f59cb469..89e5e7b3 100644 --- a/src/static.py +++ b/src/static.py @@ -21,6 +21,7 @@ 'shallow_roe_tracer_1D' : 3, 'shallow_roe_1D' : 2, 'shallow_hll_1D' : 2, + 'shallow_hlle_1D' : 2, 'shallow_bathymetry_fwave_1D' : 2, 'shallow_1D_py' : 2, 'vc_advection_1D' : 1, @@ -35,6 +36,7 @@ 'euler_hlle_with_walls_2D' : 4, 'kpp_2D' : 1, 'psystem_2D' : 3, + 'shallow_hlle_2D' : 3, 'shallow_roe_with_efix_2D' : 3, 'shallow_bathymetry_fwave_2D' : 3, 'sw_aug_2D' : 3, @@ -69,6 +71,7 @@ 'shallow_roe_tracer_1D' : 3, 'shallow_roe_1D' : 2, 'shallow_hll_1D' : 2, + 'shallow_hlle_1D' : 2, 'shallow_bathymetry_fwave_1D' : 2, 'shallow_1D_py' : 2, 'vc_advection_1D' : 1, @@ -83,6 +86,7 @@ 'euler_hlle_with_walls_2D' : 2, 'kpp_2D' : 1, 'psystem_2D' : 2, + 'shallow_hlle_2D' : 2, 'shallow_roe_with_efix_2D' : 3, 'shallow_bathymetry_fwave_2D' : 3, 'sw_aug_2D' : 3,