Merge pull request #145 from ketch/sw_hlle

Shallow water HLLE solvers

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

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

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

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',

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,