Hi I am working with a code that uses Roe solver with entropy fix.I checked the Roe paper,most of this code matches with the paper but I can not figure out 6 line of this code. I guess it is related to entropy fix issue.

I read the following papers: 1-P.L. Roe, Journal of Computational Physics,43,(1981),"Approximate Riemann solvers,Parameter vectors and difference schemes" 2-P.L. Roe, J. Pike, Computing methods in Applied Science and Engineering,(1984),(North-Holland),"Efficient Construction and Utilization of Approximate Riemann Solution" 3-P.L. Roe, SIAM J. Sci. Stat. Comput.,Vol 13,No 2,1992,"Sonic Flux Formulae" 4-the book by Toro,"Riemann solvers"

I wonder if any expert help me or refer me to a reference regarding this issue.

And this is the code(I specified the line that I have problem with ******):

SUBROUTINE riemann(nhat,ds,Qleft,Qright,flux) !................................................. ......... ! ! compute the Roe flux with entropy fix. !................................................. .........

USE physics !

IMPLICIT DOUBLE PRECISION (a-h,o-z)

DOUBLE PRECISION, DIMENSION(5) :: Qleft,Qright,flux

DOUBLE PRECISION, DIMENSION(3) :: nhat !

rho = Qleft(1)

rhou = Qleft(2)

rhov = Qleft(3)

rhow = Qleft(4)

rhoe = Qleft(5)

rhon = Qright(1)

rhoun = Qright(2)

rhovn = Qright(3)

rhown = Qright(4)

rhoen = Qright(5)

ul = rhou/rho

vl = rhov/rho

wl = rhow/rho

pleft = (gamma-1.d0)*(rhoe - 0.5d0/rho*(rhou**2 + rhov**2 + rhow**2 )) !

ur = rhoun/rhon

vr = rhovn/rhon

wr = rhown/rhon

pright = (gamma-1.d0)*(rhoen - 0.5d0/rhon*(rhoun**2 + rhovn**2+ rhown**2)) !

ql = nhat(1)*ul + nhat(2)*vl + nhat(3)*wl

qr = nhat(1)*ur + nhat(2)*vr + nhat(3)*wr

hl = 0.5d0*(ul*ul + vl*vl + wl*wl) + gamma/(gamma-1.d0)*pleft/rho

hr = 0.5d0*(ur*ur + vr*vr + wr*wr) + gamma/(gamma-1.d0)*pright/rhon

!--------------------- !square root averaging !---------------------

rtd = sqrt(rho*rhon)

betal = rho/(rho + rtd)

betar = 1.d0 - betal

utd = betal*ul + betar*ur

vtd = betal*vl + betar*vr

wtd = betal*wl + betar*wr

htd = betal*hl + betar*hr

atd2 = (gamma-1.d0)*(htd - 0.5d0*(utd*utd + vtd*vtd + wtd*wtd))

atd = sqrt(atd2)

qtd = utd*nhat(1) + vtd*nhat(2) + wtd*nhat(3) !

IF(qtd >= 0.0d0) THEN

dw1 = 0.5d0*((pright - pleft)/atd2 - (qr - ql)*rtd/atd) ************************************************* I could not figure out the next 5 line *************************************************

sp1 = qtd - atd

sp1m = min(sp1,0.0d0)

hd1m = ((gamma+1.d0)/4.d0*atd/rtd)*dw1

eta1 = max(-abs(sp1) - hd1m,0.0d0)

udw1 = dw1*(sp1m - 0.5d0*eta1)

rql = rho*ql

flux(1) = ds*(rql + udw1)

flux(2) = ds*(rql*ul + pleft*nhat(1) + udw1*(utd - atd*nhat(1)))

flux(3) = ds*(rql*vl + pleft*nhat(2) + udw1*(vtd - atd*nhat(2)))

flux(4) = ds*(rql*wl + pleft*nhat(3) + udw1*(wtd - atd*nhat(3)))

flux(5) = ds*(rql*hl + udw1*(htd - qtd*atd))

ELSE

dw4 = 0.5d0*((pright - pleft)/atd2 + (qr - ql)*rtd/atd)

sp4 = qtd + atd

sp4p = max(sp4,0.0d0)

hd4 = ((gamma+1.d0)/4.d0*atd/rtd)*dw4

eta4 = max(-abs(sp4) + hd4,0.0d0)

udw4 = dw4*(sp4p + 0.5d0*eta4)

rqr = rhon*qr

flux(1) = ds*(rqr - udw4)

flux(2) = ds*(rqr*ur + pright*nhat(1) - udw4*(utd + atd*nhat(1)))

flux(3) = ds*(rqr*vr + pright*nhat(2) - udw4*(vtd + atd*nhat(2)))

flux(4) = ds*(rqr*wr + pright*nhat(3) - udw4*(wtd + atd*nhat(3)))

flux(5) = ds*(rqr*hr - udw4*(htd + qtd*atd))

ENDIF

RETURN

END SUBROUTINE riemann

Thanks

Mehdi

Yes, those lines of code are indeed the entropy fix.

I'm at home so don't have access to papers etc, but the paper you want is that by Harten and Hyman (probably in J. Comp. Physics) which explains an entropy fix and its application.