ketch
2/21/2017 - 6:27 PM
Comparison of some Riemann solvers for the variable speed-limit LWR traffic model
Comparison of some Riemann solvers for the variable speed-limit LWR traffic model
! 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
! 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
! 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
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