Comparison of some Riemann solvers for the variable speed-limit LWR traffic model

Makefile

all: traffic_fwave_old.so traffic_fwave.so traffic_qwave.so

traffic_fwave_old.so: rp1_traffic_vc_tracer_fwave_old.f90
	f2py -c rp1_traffic_vc_tracer_fwave_old.f90 -m traffic_fwave_old

traffic_fwave.so: rp1_traffic_vc_tracer_fwave.f90
	f2py -c rp1_traffic_vc_tracer_fwave.f90 -m traffic_fwave

traffic_qwave.so: rp1_traffic_vc_tracer_qwave.f90
	f2py -c rp1_traffic_vc_tracer_qwave.f90 -m traffic_qwave

clean:
	rm *.so

rp1_traffic_vc_tracer_fwave.f90

! Riemann solver for the variable speed limit LWR traffic model
! with a tracer:

! q_t + (u_max q(1-q))_x = 0
! p_t + u_max(1-q) p_x = 0

! Here the speed limit u_max (stored in aux(1,i) may vary with x and t.

! waves: 2
! equations: 2
! aux fields: 1

! Conserved quantities:
!       1 q
!       2 p

! Auxiliary variables:
!       1 u_max

! See http://www.clawpack.org/riemann.html for a detailed explanation
! of the Riemann solver API.

subroutine rp1(maxm,num_eqn,num_waves,num_aux,num_ghost,num_cells, &
               ql,qr,auxl,auxr,fwave,s,amdq,apdq)
    implicit none
    ! Inputs
    integer, intent(in) :: maxm, num_eqn, num_waves, num_aux, num_ghost, num_cells
    double precision, intent(in), dimension(num_eqn, 1-num_ghost:maxm+num_ghost) :: ql,qr
    double precision, intent(in), dimension(num_aux, 1-num_ghost:maxm+num_ghost) :: auxl,auxr
    ! Outputs
    double precision, intent(out) :: s(num_waves,1-num_ghost:num_cells+num_ghost)
    double precision, intent(out) :: fwave(num_eqn,num_waves,1-num_ghost:num_cells+num_ghost)
    double precision, intent(out), dimension(num_eqn, 1-num_ghost:maxm+num_ghost) :: amdq,apdq
    ! Locals
    double precision :: fim1, fi, sim1, si, f0
    integer :: i

    do i = 2-num_ghost, num_cells+num_ghost
        ! Compute the wave, speeds, and flux difference

        ! compute characteristic speed in each cell:
        sim1 = auxr(1,i-1)*(1.d0 - 2.d0*qr(1,i-1))
        si   = auxl(1,i  )*(1.d0 - 2.d0*ql(1,i  ))
        s(2,i) = auxl(1,i) * (1.d0 - ql(1,i))

        ! compute flux in each cell and flux difference:
        fim1 = auxr(1,i-1)*qr(1,i-1)*(1.d0 - qr(1,i-1))
        fi   = auxl(1,i  )*ql(1,i  )*(1.d0 - ql(1,i  ))

        fwave(1,1,i) = fi - fim1
        fwave(2,1,i) = 0.d0
        fwave(1,2,i) = 0.d0
        fwave(2,2,i) = s(2,i)*(ql(2,i) - qr(2,i-1))

        apdq(2,i) = fwave(2,2,i)
        amdq(2,i) = 0.d0 ! Traffic always moves right

        if (sim1 .lt. 0.d0 .and. si .le. 0.d0) then
            ! left-going
            s(1,i) = sim1
            amdq(1,i) = fwave(1,1,i)
            apdq(1,i) = 0.d0
        else if (sim1 .ge. 0.d0 .and. si .gt. 0.d0) then
            ! right-going
            s(1,i) = si
            amdq(1,i) = 0.d0
            apdq(1,i) = fwave(1,1,i)
        else if (sim1 .lt. 0.d0 .and. si .gt. 0.d0) then
            ! transonic rarefaction
            s(1,i) = 0.5d0*(sim1 + si)

            ! entropy fix:  (perhaps doesn't work for all cases!!!)
            ! This assumes the flux in the transonic case should
            ! correspond to q=0.5 on the side with the smaller umax value.
            f0 = dmin1(auxr(1,i-1),auxl(1,i))*0.25d0
            ! split fwave between amdq and apdq:
            amdq(1,i) = f0 - fim1
            apdq(1,i) = fi - f0

        else
            ! transonic shock
            s(1,i) = 0.5d0*(sim1 + si)
            if (fi-fim1 .lt. 0.d0) then
                amdq(1,i) = fwave(1,1,i)
                apdq(1,i) = 0.d0
            else if (fi-fim1 .gt. 0.d0) then
                amdq(1,i) = 0.d0
                apdq(1,i) = fwave(1,1,i)
            else
                amdq(1,i) = 0.5d0 * fwave(1,1,i)
                apdq(1,i) = 0.5d0 * fwave(1,1,i)
            endif
        endif

    enddo
    return
    end

rp1_traffic_vc_tracer_fwave_old.f90

! Riemann solver for the variable speed limit LWR traffic model
! with a tracer:

! q_t + (u_max q(1-q))_x = 0
! p_t + u_max(1-q) p_x = 0

! Here the speed limit u_max (stored in aux(1,i) may vary with x and t.

! waves: 2
! equations: 2
! aux fields: 1

! Conserved quantities:
!       1 q
!       2 p

! Auxiliary variables:
!       1 u_max

! See http://www.clawpack.org/riemann.html for a detailed explanation
! of the Riemann solver API.

subroutine rp1(maxm,num_eqn,num_waves,num_aux,num_ghost,num_cells, &
               ql,qr,auxl,auxr,fwave,s,amdq,apdq)
    implicit none
    ! Inputs
    integer, intent(in) :: maxm, num_eqn, num_waves, num_aux, num_ghost, num_cells
    double precision, intent(in), dimension(num_eqn, 1-num_ghost:maxm+num_ghost) :: ql,qr
    double precision, intent(in), dimension(num_aux, 1-num_ghost:maxm+num_ghost) :: auxl,auxr
    ! Outputs
    double precision, intent(out) :: s(num_waves,1-num_ghost:num_cells+num_ghost)
    double precision, intent(out) :: fwave(num_eqn,num_waves,1-num_ghost:num_cells+num_ghost)
    double precision, intent(out), dimension(num_eqn, 1-num_ghost:maxm+num_ghost) :: amdq,apdq
    ! Locals
    double precision :: fim1, fi, sim1, si, f0
    integer :: i

    do i = 2-num_ghost, num_cells+num_ghost
        ! Compute the wave, speeds, and flux difference

        ! compute characteristic speed in each cell:
        sim1 = auxr(1,i-1)*(1.d0 - 2.d0*qr(1,i-1))
        si   = auxl(1,i  )*(1.d0 - 2.d0*ql(1,i  ))
        s(2,i) = auxl(1,i) * (1.d0 - ql(1,i))

        ! compute flux in each cell and flux difference:
        fim1 = auxr(1,i-1)*qr(1,i-1)*(1.d0 - qr(1,i-1))
        fi   = auxl(1,i  )*ql(1,i  )*(1.d0 - ql(1,i  ))

        fwave(1,1,i) = fi - fim1
        fwave(2,1,i) = 0.d0
        fwave(1,2,i) = 0.d0
        fwave(2,2,i) = s(2,i)*(ql(2,i) - qr(2,i-1))

        apdq(2,i) = fwave(2,2,i)
        amdq(2,i) = 0.d0 ! Traffic always moves right

        if (sim1 .lt. 0.d0 .and. si .le. 0.d0) then
            ! left-going
            s(1,i) = sim1
            amdq(1,i) = fwave(1,1,i)
            apdq(1,i) = 0.d0
        else if (sim1 .ge. 0.d0 .and. si .gt. 0.d0) then
            ! right-going
            s(1,i) = si
            amdq(1,i) = 0.d0
            apdq(1,i) = fwave(1,1,i)
        else if (sim1 .lt. 0.d0 .and. si .gt. 0.d0) then
            ! transonic rarefaction
            s(1,i) = 0.5d0*(sim1 + si)

            ! entropy fix:  (perhaps doesn't work for all cases!!!)
            ! This assumes the flux in the transonic case should
            ! correspond to q=0.5 on the side with the smaller umax value.
            f0 = dmin1(auxr(1,i-1),auxl(1,i))*0.25d0
            ! split fwave between amdq and apdq:
            amdq(1,i) = f0 - fim1
            apdq(1,i) = fi - f0

        else
            ! transonic shock
            s(1,i) = 0.5d0*(sim1 + si)
            if (s(1,i) .lt. 0.d0) then
                amdq(1,i) = fwave(1,1,i)
                apdq(1,i) = 0.d0
            else if (s(1,i) .gt. 0.d0) then
                amdq(1,i) = 0.d0
                apdq(1,i) = fwave(1,1,i)
            else
                amdq(1,i) = 0.5d0 * fwave(1,1,i)
                apdq(1,i) = 0.5d0 * fwave(1,1,i)
            endif
        endif

    enddo
    return
    end

rp1_traffic_vc_tracer_qwave.f90

! Riemann solver for the variable speed limit LWR traffic model
! with a tracer:

! q_t + (u_max q(1-q))_x = 0
! p_t + u_max(1-q) p_x = 0

! Here the speed limit u_max (stored in aux(1,i) may vary with x and t.

! waves: 2
! equations: 2
! aux fields: 1

! Conserved quantities:
!       1 q
!       2 p

! Auxiliary variables:
!       1 u_max

! See http://www.clawpack.org/riemann.html for a detailed explanation
! of the Riemann solver API.

subroutine rp1(maxm,num_eqn,num_waves,num_aux,num_ghost,num_cells, &
               ql,qr,auxl,auxr,wave,s,amdq,apdq)
    implicit none
    ! Inputs
    integer, intent(in) :: maxm, num_eqn, num_waves, num_aux, num_ghost, num_cells
    double precision, intent(in), dimension(num_eqn, 1-num_ghost:maxm+num_ghost) :: ql,qr
    double precision, intent(in), dimension(num_aux, 1-num_ghost:maxm+num_ghost) :: auxl,auxr
    ! Outputs
    double precision, intent(out) :: s(num_waves,1-num_ghost:num_cells+num_ghost)
    double precision, intent(out) :: wave(num_eqn,num_waves,1-num_ghost:num_cells+num_ghost)
    double precision, intent(out), dimension(num_eqn, 1-num_ghost:maxm+num_ghost) :: amdq,apdq
    ! Locals
    double precision :: f_l, f_r, s_l, s_r, f0, q_l, q_r, v_l, v_r
    integer :: i

    do i = 2-num_ghost, num_cells+num_ghost
        q_r = ql(1,i)
        q_l = qr(1,i-1)
        v_r = auxl(1,i)
        v_l = auxr(1,i-1)

        ! compute characteristic speed in each cell:
        s_l = v_l*(1.d0 - 2.d0*q_l)
        s_r = v_r*(1.d0 - 2.d0*q_r)
        s(2,i) = v_r * (1.d0 - q_r)

        ! compute flux in each cell and flux difference:
        f_l = v_l*q_l*(1.d0 - q_l)
        f_r = v_r*q_r*(1.d0 - q_r)

        wave(1,1,i) = q_r - q_l ! Trying this even though it's not right
        wave(2,1,i) = 0.d0
        wave(1,2,i) = 0.d0
        wave(2,2,i) = s(2,i)*(ql(2,i) - qr(2,i-1))

        if (q_r .ne. q_l) then
            s(1,i) = (f_r-f_l)/(q_r - q_l)
        else
            s(1,i) = 0.d0
        endif

        apdq(2,i) = wave(2,2,i)
        amdq(2,i) = 0.d0 ! Traffic always moves right


        if ((f_l .ge. 0.25d0*v_r) .and. (s_r .gt. 0.d0)) then
            ! left-going shock, right-going rarefaction
            f0 = 0.25d0*v_r
        elseif ((f_r .ge. 0.25d0*v_l) .and. (s_l .lt. 0.d0)) then
            ! right-going shock, left-going rarefaction
            f0 = 0.25d0*v_l
        elseif ((s_r .le. 0.d0) .and. (f_l .gt. f_r)) then
                ! left-going shock
                f0 = f_r
        elseif ((s_l .ge. 0.d0) .and. (f_r .gt. f_l)) then
                ! right-going shock
                f0 = f_l
        elseif ((s_l .le. 0.d0) .and. (s_r .ge. 0.d0)) then
            ! Transonic rarefaction
            if (v_r .le. v_l) then
                f0 = 0.25d0*v_r
            else
                f0 = 0.25d0*v_l
            endif
        elseif (f_l .le. f_r) then
            ! left-going rarefaction
            f0 = f_r
        else
            ! right-going rarefaction
            f0 = f_l
        endif
        amdq(1,i) = f0 - f_l
        apdq(1,i) = f_r - f0


    enddo
    return
    end

Variable_speed_riemann_solver_comparison.ipynb