Normal and transverse Riemann solvers for shallow water equations on a quadrilateral grid.

rpn2.f90

! =====================================================
subroutine rpn2(ixy,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,wave,s,amdq,apdq)
! =====================================================
!
!   # Roe-solver for the 2D shallow water equations
!   #  on a quadrilateral grid
!
!   # solve Riemann problems along one slice of data.
!
!   # On input, ql contains the state vector at the left edge of each cell
!   #           qr contains the state vector at the right edge of each cell
!
!   # This data is along a slice in the x-direction if ixy=1
!   #                            or the y-direction if ixy=2.
!   # On output, wave contains the waves, s the speeds,
!   # and amdq, apdq the decomposition of the flux difference
!   #   f(qr(i-1)) - f(ql(i))
!   # into leftgoing and rightgoing parts respectively.
!   # With the Roe solver we have
!   #    amdq  =  A^- \Delta q    and    apdq  =  A^+ \Delta q
!   # where A is the Roe matrix.  An entropy fix can also be incorporated
!   # into the flux differences.
!
!   # Note that the i'th Riemann problem has left state qr(:,i-1)
!   #                                    and right state ql(:,i)
!   # From the basic clawpack routines, this routine is called with ql = qr

    implicit double precision (a-h,o-z)

    dimension wave(meqn, mwaves, 1-mbc:maxm+mbc)
    dimension    s(mwaves, 1-mbc:maxm+mbc)
    dimension   ql(meqn, 1-mbc:maxm+mbc)
    dimension   qr(meqn, 1-mbc:maxm+mbc)
    dimension apdq(meqn, 1-mbc:maxm+mbc)
    dimension amdq(meqn, 1-mbc:maxm+mbc)
    dimension auxl(7, 1-mbc:maxm+mbc)
    dimension auxr(7, 1-mbc:maxm+mbc)

!   # Gravity constant set in the shallow1D.py file
    common /cparam/  grav

!   # Roe averages quantities of each interface
    parameter (maxm2 = 1800)
    double precision u(-6:maxm2+7),v(-6:maxm2+7),a(-6:maxm2+7),h(-6:maxm2+7)

!   local arrays
!   ------------
    dimension delta(3)
    logical :: efix

    dimension unorl(-6:maxm2+7), unorr(-6:maxm2+7)
    dimension utanl(-6:maxm2+7), utanr(-6:maxm2+7)
    dimension alf(-6:maxm2+7)
    dimension beta(-6:maxm2+7)

    data efix /.true./    !# Use entropy fix for transonic rarefactions

!   # rotate the velocities q(2) and q(3) so that it is aligned with grid
!   # normal.  The normal vector for the face at the i'th Riemann problem
!   # is stored in the aux array
!   # in locations (1,2) if ixy=1 or (4,5) if ixy=2.  The ratio of the
!   # length of the cell side to the length of the computational cell
!   # is stored in aux(3) or aux(6) respectively.

    if (ixy == 1) then
        inx = 1
        iny = 2
        ilenrat = 3
    else
        inx = 4
        iny = 5
        ilenrat = 6
    endif

!     # determine rotation matrix
!             [ alf  beta ]
!             [-beta  alf ]
!
!     # note that this reduces to identity on standard cartesian grid

    do i=2-mbc,mx+mbc
        alf(i) = auxl(inx,i)
        beta(i) = auxl(iny,i)
        unorl(i) = alf(i)*ql(2,i) + beta(i)*ql(3,i)
        unorr(i-1) = alf(i)*qr(2,i-1) + beta(i)*qr(3,i-1)
        utanl(i) = -beta(i)*ql(2,i) + alf(i)*ql(3,i)
        utanr(i-1) = -beta(i)*qr(2,i-1) + alf(i)*qr(3,i-1)
    enddo

!   # compute the Roe-averaged variables needed in the Roe solver.
!   # These are stored in the common block comroe since they are
!   # later used in routine rpt2 to do the transverse wave splitting.

    do 10 i = 2-mbc, mx+mbc
        h(i) = (qr(1,i-1)+ql(1,i))*0.50d0
        hsqrtl = dsqrt(qr(1,i-1))
        hsqrtr = dsqrt(ql(1,i))
        hsq2 = hsqrtl + hsqrtr
        u(i) = (unorr(i-1)/hsqrtl + unorl(i)/hsqrtr) / hsq2
        v(i) = (utanr(i-1)/hsqrtl + utanl(i)/hsqrtr) / hsq2
        a(i) = dsqrt(grav*h(i))
    10 continue

!   # now split the jump in q at each interface into waves
!
!   # find a1 thru a3, the coefficients of the 3 eigenvectors:
    do 20 i = 2-mbc, mx+mbc
        delta(1) = ql(1,i) - qr(1,i-1)
        delta(2) = unorl(i) - unorr(i-1)
        delta(3) = utanl(i) - utanr(i-1)
        a1 = ((u(i)+a(i))*delta(1) - delta(2))*(0.50d0/a(i))
        a2 = -v(i)*delta(1) + delta(3)
        a3 = (-(u(i)-a(i))*delta(1) + delta(2))*(0.50d0/a(i))

!       # Compute the waves.

        wave(1,1,i) = a1
        wave(2,1,i) = alf(i)*a1*(u(i)-a(i)) - beta(i)*a1*v(i)
        wave(3,1,i) = beta(i)*a1*(u(i)-a(i)) + alf(i)*a1*v(i)
        s(1,i) = (u(i)-a(i)) * auxl(ilenrat,i)

        wave(1,2,i) = 0.0d0
        wave(2,2,i) = -beta(i)*a2
        wave(3,2,i) = alf(i)*a2
        s(2,i) = u(i) * auxl(ilenrat,i)

        wave(1,3,i) = a3
        wave(2,3,i) = alf(i)*a3*(u(i)+a(i)) - beta(i)*a3*v(i)
        wave(3,3,i) = beta(i)*a3*(u(i)+a(i)) + alf(i)*a3*v(i)
        s(3,i) = (u(i)+a(i)) * auxl(ilenrat,i)
    20 continue

!    # compute flux differences amdq and apdq.
!    ---------------------------------------

    if (efix) go to 110

!   # no entropy fix
!   ----------------
!
!   # amdq = SUM s*wave   over left-going waves
!   # apdq = SUM s*wave   over right-going waves

    do 100 m=1,3
        do 100 i=2-mbc, mx+mbc
            amdq(m,i) = 0.d0
            apdq(m,i) = 0.d0
        do 90 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
        90 continue
    100 continue
    go to 900

!-----------------------------------------------------

  110 continue

!   # With entropy fix
!   ------------------
!
!   # compute flux differences amdq and apdq.
!   # First compute amdq as sum of s*wave for left going waves.
!   # Incorporate entropy fix by adding a modified fraction of wave
!   # if s should change sign.

    do 200 i=2-mbc,mx+mbc

!       check 1-wave

        him1 = qr(1,i-1)
        s0 =  (unorr(i-1)/him1 - dsqrt(grav*him1)) * auxl(ilenrat,i)

!       check for fully supersonic case :

        if (s0 > 0.0d0 .AND. s(1,i) > 0.0d0) then
            do 60 m=1,3
                amdq(m,i)=0.0d0
            60 continue
            goto 200
        endif

        h1 = qr(1,i-1)+wave(1,1,i)
        hu1= unorr(i-1)+ alf(i)*wave(2,1,i) + beta(i)*wave(3,1,i)
        s1 = (hu1/h1 - dsqrt(grav*h1))* auxl(ilenrat,i)

!       speed just to right of 1-wave

        if (s0 < 0.0d0 .AND. s1 > 0.0d0) then
!       transonic rarefaction in 1-wave
            sfract = s0*((s1-s(1,i))/(s1-s0))
        else if (s(1,i) < 0.0d0) then
!       1-wave is leftgoing
            sfract = s(1,i)
        else
!       1-wave is rightgoing
            sfract = 0.0d0
        endif
        do 120 m=1,3
            amdq(m,i) = sfract*wave(m,1,i)
        120 continue

!       check 2-wave

        if (s(2,i) > 0.0d0) then
!       #2 and 3 waves are right-going
            go to 200 
        endif

        do 140 m=1,3
            amdq(m,i) = amdq(m,i) + s(2,i)*wave(m,2,i)
        140 continue

!       check 3-wave

        hi = ql(1,i)
        s03 = (unorl(i)/hi + dsqrt(grav*hi)) * auxl(ilenrat,i)
        h3=ql(1,i)-wave(1,3,i)
        hu3=unorl(i)- (alf(i)*wave(2,3,i) + beta(i)*wave(3,3,i))
        s3=(hu3/h3 + dsqrt(grav*h3)) * auxl(ilenrat,i)
        if (s3 < 0.0d0 .AND. s03 > 0.0d0) then
!       transonic rarefaction in 3-wave
            sfract = s3*((s03-s(3,i))/(s03-s3))
        else if (s(3,i) < 0.0d0) then
!       3-wave is leftgoing
            sfract = s(3,i)
        else
!       3-wave is rightgoing
            goto 200
        endif

        do 160 m=1,3
            amdq(m,i) = amdq(m,i) + sfract*wave(m,3,i)
        160 continue
    200 continue

!   compute rightgoing flux differences :

    do 220 m=1,3
        do 220 i = 2-mbc,mx+mbc
            df = 0.0d0
            do 210 mw=1,mwaves
                df = df + s(mw,i)*wave(m,mw,i)
            210 continue
            apdq(m,i)=df-amdq(m,i)
    220 continue

    900 continue

    return
    end subroutine rpn2

rpt2.f90

! =====================================================
subroutine rpt2(ixy,imp,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,aux1,aux2,aux3,asdq,bmasdq,bpasdq)
! =====================================================
    implicit double precision (a-h,o-z)

!   # Riemann solver in the transverse direction for the shallow water
!   # equations  on a quadrilateral grid.
!
!   # Split asdq (= A^* \Delta q, where * = + or -)
!   # into down-going flux difference bmasdq (= B^- A^* \Delta q)
!   #    and up-going flux difference bpasdq (= B^+ A^* \Delta q)

    dimension     ql(meqn, 1-mbc:maxm+mbc)
    dimension     qr(meqn, 1-mbc:maxm+mbc)
    dimension   asdq(meqn, 1-mbc:maxm+mbc)
    dimension bmasdq(meqn, 1-mbc:maxm+mbc)
    dimension bpasdq(meqn, 1-mbc:maxm+mbc)
    dimension   aux1(7, 1-mbc:maxm+mbc)
    dimension   aux2(7, 1-mbc:maxm+mbc)
    dimension   aux3(7, 1-mbc:maxm+mbc)

!   # Gravity constant set in the shallow1D.py file
    common /cparam/ grav

!   # Roe averages quantities of each interface
    parameter (maxm2 = 1800)
    double precision u(-6:maxm2+7),v(-6:maxm2+7),a(-6:maxm2+7),h(-6:maxm2+7)

!   local arrays
!   ------------
    dimension   alf(-6:maxm2+7)
    dimension  beta(-6:maxm2+7)
    dimension  wave(meqn, mwaves, -6:maxm2+7)
    dimension     s(mwaves, -6:maxm2+7)
    dimension delta(3)

    if (ixy == 1) then
        inx = 4
        iny = 5
        in = 3
    else
        inx = 1
        iny = 2
        in = 6
    endif

!   # imp is used to flag whether wave is going to left or right,
!   # since states and grid are different on each side

    if (imp == 1) then
!   # asdq = amdq, moving to left
        ix1 = 2-mbc
        ixm1 = mx+mbc
    else
!   # asdq = apdq, moving to right
        ix1 = 1-mbc
        ixm1 = mx+mbc
    endif

!   --------------
!   # up-going:
!   --------------

!   # determine rotation matrix for interface above cell, using aux3
!          [ alf  beta ]
!          [-beta  alf ]

    do i=ix1,ixm1
        if (imp == 1) then
            i1 = i-1
        else
            i1 = i
        endif
        alf(i) = aux3(inx,i1)
        beta(i) = aux3(iny,i1)
        h(i) = ql(1,i1)
        u(i) = (alf(i)*ql(2,i1) + beta(i)*ql(3,i1)) / h(i)
        v(i) = (-beta(i)*ql(2,i1) + alf(i)*ql(3,i1)) / h(i)
        a(i) = dsqrt(grav*h(i))
    enddo

!   # now split asdq into waves:
!
!   # find a1 thru a3, the coefficients of the 3 eigenvectors:

    do 20 i = ix1,ixm1
        delta(1) = asdq(1,i) 
        delta(2) = alf(i)*asdq(2,i) + beta(i)*asdq(3,i)
        delta(3) = -beta(i)*asdq(2,i) + alf(i)*asdq(3,i)
        a1 = ((u(i)+a(i))*delta(1) - delta(2))*(0.50d0/a(i))
        a2 = -v(i)*delta(1) + delta(3)
        a3 = (-(u(i)-a(i))*delta(1) + delta(2))*(0.50d0/a(i))

!       # Compute the waves.

        wave(1,1,i) = a1
        wave(2,1,i) = alf(i)*a1*(u(i)-a(i)) - beta(i)*a1*v(i)
        wave(3,1,i) = beta(i)*a1*(u(i)-a(i)) + alf(i)*a1*v(i)
        s(1,i) = (u(i)-a(i)) / aux2(in,i)

        wave(1,2,i) = 0.0d0
        wave(2,2,i) = -beta(i)*a2
        wave(3,2,i) = alf(i)*a2
        s(2,i) = u(i) / aux2(in,i)

        wave(1,3,i) = a3
        wave(2,3,i) = alf(i)*a3*(u(i)+a(i)) - beta(i)*a3*v(i)
        wave(3,3,i) = beta(i)*a3*(u(i)+a(i)) + alf(i)*a3*v(i)
        s(3,i) = (u(i)+a(i)) / aux2(in,i)
    20 continue

!    # compute flux difference bpasdq
!    --------------------------------

    do 40 m=1,3
        do 40 i=ix1,ixm1
            bpasdq(m,i) = 0.d0
            do 30 mw=1,mwaves
                bpasdq(m,i) = bpasdq(m,i) + dmax1(s(mw,i),0.d0)*wave(m,mw,i)
            30 continue
    40 continue

!   --------------
!   # down-going:
!   --------------

!   # determine rotation matrix for interface below cell, using aux2
!           [ alf  beta ]
!           [-beta  alf ]

    do i=ix1,ixm1

        if (imp == 1) then
            i1 = i-1
        else
            i1 = i
        endif

        alf(i) = aux2(inx,i1)
        beta(i) = aux2(iny,i1)
        u(i) = (alf(i)*ql(2,i1) + beta(i)*ql(3,i1)) / h(i)
        v(i) = (-beta(i)*ql(2,i1) + alf(i)*ql(3,i1)) / h(i)
    enddo

!   # now split asdq into waves:
!
!   # find a1 thru a3, the coefficients of the 3 eigenvectors:
    do 80 i = ix1,ixm1
       delta(1) = asdq(1,i) 
       delta(2) = alf(i)*asdq(2,i) + beta(i)*asdq(3,i)
       delta(3) = -beta(i)*asdq(2,i) + alf(i)*asdq(3,i)
       a1 = ((u(i)+a(i))*delta(1) - delta(2))*(0.50d0/a(i))
       a2 = -v(i)*delta(1) + delta(3)
       a3 = (-(u(i)-a(i))*delta(1) + delta(2))*(0.50d0/a(i))

!      # Compute the waves.

       wave(1,1,i) = a1
       wave(2,1,i) = alf(i)*a1*(u(i)-a(i)) - beta(i)*a1*v(i)
       wave(3,1,i) = beta(i)*a1*(u(i)-a(i)) + alf(i)*a1*v(i)
       s(1,i) = (u(i)-a(i)) / aux2(in,i)

       wave(1,2,i) = 0.0d0
       wave(2,2,i) = -beta(i)*a2
       wave(3,2,i) = alf(i)*a2
       s(2,i) = u(i) / aux2(in,i)

       wave(1,3,i) = a3
       wave(2,3,i) = alf(i)*a3*(u(i)+a(i)) - beta(i)*a3*v(i)
       wave(3,3,i) = beta(i)*a3*(u(i)+a(i)) + alf(i)*a3*v(i)
       s(3,i) = (u(i)+a(i)) / aux2(in,i)
    80 continue

!   # compute flux difference bmasdq
!   --------------------------------

    do 100 m=1,3
        do 100 i=ix1,ixm1
            bmasdq(m,i) = 0.d0
            do 90 mw=1,mwaves
                bmasdq(m,i) = bmasdq(m,i) + dmin1(s(mw,i), 0.d0)*wave(m,mw,i)
            90 continue
    100 continue

    return
    end subroutine rpt2