Intro

The Riemann problem is an IVP for the Euler Equations, which consists of 2 constant initial data states separated by a discontinuity between them, say at x=0. For x<x=0, we say the data state is XL, for which there are initial variables WL=(rho_L, p_L, u_L). Similarly for the initial right data state, XR=x>x=0, the initial variables are WR=(rho_R, p_R, u_R). The solution of this IVP consists of 3 nonlinear waves, a left wave that is either a shock or a rarefaction, a center contact discontinuity, and a right wave that is either a shock or a rarefaction wave. Depending on which type of L- or R-wave is present, different expressions exist that describe the change in variables across them. The contact discontinuity is special in that the pressure (p) and velocity (u) are constant across it.

Program Outline for Exact Riemann Solver

Given the above properties of the Riemann problem, an algebraic expression can be derived which gives p in the central "star-region" (denoted by '*'). The overall structure of the Riemann problem then is to solve this algebraic equation for p*. Once p* is known, u* follows immediately. The remaining rho*_L and rho*_R follow from expressions valid for the specific L- or R- wave present.

Specifically, this means that the code will determine at every point (x,t), that point's relative position to the different waves present. Once the position is determined, the solution is given by analytic expressions for the following five possibilities: pre- or post- shock, ahead of rarefaction head, behind rarefaction tail, or within the rarefaction fan. The position of each sampling point (x,t) is determined by a characteristic speed in the grid given by s = dx/t, where dx is the distance from the initial discontinuity to the sampling point (x,t), and t is the simulation end time. This s for every point on the grid can be compared to the known present waves, as their speeds are known exactly. In this way, the relative position of the sampling point with respect to the waves on the grid is determined.

Exact Riemann Solver

PROGRAM ExactRiemannSolver

--———————————————————————————————————————————————- ! This program finds the solution of p*, u*, and rho*_L and rho*_R given initial data states wL ! and wR of the Riemann problem ! Author: Erica Kaminski, 12.17.12 --———————————————————————————————————————————————-

IMPLICIT NONE

REAL

DomLen, TimeOut, DisPos

INTEGER

cells, i

REAL

rho_L, p_L, u_L !Left data vars

REAL

rho_R, p_R, u_R !Right data vars

REAL, dimension(3)

wL !Left data state

REAL, dimension(3)

wR !Right data state

REAL

Cs_L, Cs_R !Sound speed Left and Right states

REAL, PARAMETER

gamma = 1.4 ! Ratio specific heats

REAL

p_Star, u_Star ! Output of program

REAL

mpa ! Normalization constant for P

REAL

u_diff, dx, xpos, s

REAL

rhos, us, ps

NAMELIST /InitialData/ DomLen, TimeOut, DisPos, cells, rho_L, p_L, u_L, rho_R, p_R, u_R, mpa

--———————————————————————————————————————————————— ! Main body of program here --————————————————————————————————————————————————

OPEN(UNIT = 3, File = "initial.data", status = 'old') !Initialize data READ(3, NML = InitialData) CLOSE(3)

wL = (/rho_L, p_L, u_L/) !Populate w arrays wR = (/rho_R, p_R, u_R/)

Cs_L = Sqrt(Gamma*p_L/rho_L) Cs_R = Sqrt(Gamma*p_R/rho_R)

U_DIFF = u_R - u_L

!PRESSURE POSITIVITY CONDITION —-

IF (2.0*Cs_L/(gamma-1.0) + 2.0*Cs_R/(gamma-1.0).LE.U_DIFF) THEN

PRINT *, 'PRESSURE POSITIVITY CONDITION VIOLATED!' PRINT *, 'VACUUM STATE IS CREATED BY INITIAL CONDITIONS' PRINT *, 'THIS RIEMANN SOLVER IS NOT CORRECT METHOD FOR SOLUTION' PRINT *, 'PROGRAM STOPPED' STOP

END IF

CALL Find_Pstar(p_Star, u_Star) !check that p_Star = p from subprpgram with a print statement - I did, it does

OPEN(UNIT=2, FILE='exact.out', STATUS='OLD', POSITION='APPEND')

dx = DomLen/REAL(cells) WRITE(2,*)'xpos rho u P iE' DO 10 i = 1, cells

xpos = dx*(REAL(i)-0.5) !cell-centered s = (xpos - DisPos)/TimeOut !characteristic speed CALL SampleExact(p_Star, u_Star, s, rhos, us, ps) WRITE(2, 20)xpos, rhos, us, ps/mpa, ps/rhos/(gamma-1.0)/mpa !the last may be internal energy…?

10 Continue

CLOSE(2)

20 FORMAT(5(F14.6,2X))

CONTAINS

SUBROUTINE SampleExact(p_Star, u_Star, s, rhos, us, ps)

REAL,INTENT(IN)

p_Star, u_Star, s REAL,INTENT(OUT) ::rhos, us, ps !sampled density, velocity, pressure at point (x,t)

REAL

SHL, STL !Speed of left rarefaction head and tail

REAL

SHR, STR !Speed of right rarefaction head and tail

REAL

CsL_Star, CsR_Star !Sound speed behind left and right rarefactions

REAL

SL, SR !Left and right shock speeds

REAL

G1, G2, G3, G4, G5, G6, ppl, ppr

CsL_Star = Cs_L*(p_Star/p_L)(G3/G2) CsR_Star = Cs_R*(p_Star/p_R)(G3/G2)

G1 = gamma + 1.0 G2 = 2.0*gamma G3 = gamma - 1.0 ppl = p_Star/p_L ppr = p_Star/p_R G4 = G3/G1 G5 = G1/G2 G6 = G3/G2

IF(s.LE.u_Star)THEN

! (x,t) is to the left of contact IF(p_L.GE.p_Star)THEN

!Left wave is a rarefaction, based on initial data !3 possible left data states then - WL, WL_Fan, W_Star - next we determine which

! PRINT*,'IT"S A RAREFACTION'

!Are these negative quantities:? SHL = u_L - Cs_L STL = u_Star - CsL_Star

IF(S.LE.SHL)THEN

! PRINT *, 'LEFT DATA STATE'

!Left data state rhos = rho_L us = u_L ps = p_L

ELSE

! IF((SHL.GE.S).AND.(S.LE.STL))THEN IF((SHL.LE.S).AND.(S.LE.STL)) THEN

! PRINT*, 'iNSIDE RARE WAVE'

!Inside rarefaction fan rhos = rho_L*(2/G1 + G3/(G1*Cs_L)*(u_L - s))(2/G3) us = 2/G1*(Cs_L + (G3/2)*u_L + s) ps = p_L*( 2/G1 + G3/(G1*Cs_L)*(u_L - s) )(G2/G3)

ELSE

!In left star region

! PRINT *, 'IN STAR REGION'

rhos = rho_L*ppl(1.0/gamma) !follows from isentropic law us = u_Star ps = p_Star

END IF

END IF

ELSE

!Left wave is a shock !Two states possible - pre-shock, post-shock

! PRINT*,'IT"S A SHOCK ! ! !'

SL = (u_L - Cs_L)*SQRT(G1/G2*ppl + G3/G2) IF(s.LE.SL)THEN

!Data state is left state/pre-shock rhos = rho_L us = u_L ps = p_L

ELSE

!Data state is post-shock rhos = rho_L*((ppl + G3/G1)/((G3/G1)*ppl + 1.0)) us = u_Star ps = p_Star

END IF

END IF

ELSE

!(x,t) is on right of contact IF(p_Star.LE.p_R)THEN

!Rarefaction, and 3 wave states possible: wR, WRFan, W*R SHR = Cs_R + u_R STR = CsR_Star + u_Star IF(S.LE.STR)THEN

!In r star region rhos = rho_R*ppr(1.0/gamma) !follows from isentropic law — CHECK? us = u_Star ps = p_Star

ELSE

IF((STR.LE.S).AND.(S.LE.SHR))THEN

!Inside right rarefaction fan rhos = rho_R*(2/G1 - G3/(G1*Cs_R)*(u_R - s))(2/G3) us = 2/G1*(-Cs_R + (G3/2)*u_R + s) ps = p_R*( 2/G1 - G3/(G1*Cs_R)*(u_R - s) )(G2/G3) !Check these expressions

ELSE

!In right data state rhos = rho_R us = u_R ps = p_R

END IF

END IF

ELSE

!Right wave is a shock SR = (u_R + Cs_R)*SQRT(G1/G2*ppr + G3/G2) IF(s.LE.SR)THEN

!In post-shock region rhos = rho_R*((ppr + G3/G1)/((G3/G1)*ppr + 1.0)) us = u_Star ps = p_Star

ELSE

!In pre-shock region rhos = rho_R us = u_R ps = p_R

END IF

END IF

END IF

END SUBROUTINE SampleExact

SUBROUTINE Find_Pstar(p, u)

!Does mpa NEED to be read into this subroutine, if it is already initialized in main body of program? --——————————————————————————————————————————————- ! Implements Newton's iterative method to solve for p_star --——————————————————————————————————————————————-

REAL, INTENT(OUT)

p, u

REAL

p0 ! Initial estimate of P_star

REAL, PARAMETER

tol = 0.0000001

REAL

change ! Relative change in pressure due to iteration

REAL

pOld, pNew! Cs_L, Cs_R, rho_L, rho_R, p_L, p_R, u_L, u_R — these should be all accessible as they are global vars

INTEGER

i, NmrItr

REAL

f_p, f_pD, FL, FLD, FR, FRD ! pressure function in Toro

REAL

Delta NAMELIST /Itr/ NmrItr

!statement labels are not global: OPEN(UNIT = 3, FILE = 'initial.data') READ(3, NML = Itr) CLOSE(3)

!Get estimate for P_star CALL P_estimate(p0)

pOld = p0 ! Initialize pOld

!Check that the vars are global:

! PRINT *, 'u_L, p_L, rho_L, cs_L =', u_L, p_L, rho_L, cs_L

OPEN(UNIT = 2, FILE = 'exact.out', STATUS = 'UNKNOWN') WRITE(2,*) '—————————————————————————————————————' WRITE(2,*) ' Iteration number Change ' WRITE(2,*) '—————————————————————————————————————'

!Need define Cs_L here too or main program enough? DO 10 i = 1,NmrItr

CALL PressFunct(FL, FLD, pOld, rho_L, p_L, Cs_L) CALL PressFunct(FR, FRD, pOld, rho_R, p_R, Cs_R) f_p = FL + FR + U_diff f_pD = FLD + FRD Delta = f_p/f_pD pNew = pOld - Delta ! This is corrected p. change = 2.0*ABS((pNew - pOld)/(pNew + pOld)) WRITE(2,*) i, change IF (change.LE.tol) GOTO 20 IF (pNew.LT.0.0) pNew = tol pOld = pNew

10 Continue

WRITE(2,*) 'Newton Rhapson Method Diverged!'

20 WRITE(2,*) 'Pressure in Star Region and U in star region'

p = pNew/MPA u = (1.0/2.0)*(u_L + u_R) + (1.0/2.0)*(fR - fL) WRITE(2,*) 'p_Star = ', p WRITE(2,*) 'u_Star = ', u

CLOSE(2)

END SUBROUTINE Find_Pstar

SUBROUTINE P_estimate(pEst)

!Employs the subroutine from Toro for estimating initial P_star

REAL,INTENT(OUT)

pEst

REAL

G1, G2, G3, G4, G5, G6, G7, G8

REAL

CUP, GEL, GER, pMax, pMin, PPV, PQ, PTL, PTR, QMAX, QUSER, UM

G1 = (gamma - 1.0)/(2.0*gamma) G2 = (gamma + 1.0)/(2.0*gamma) G3 = 2.0*gamma/(gamma - 1.0) G4 = 2.0/(gamma - 1.0) G5 = 2.0/(gamma + 1.0) G6 = (gamma - 1.0)/(gamma + 1.0) G7 = (gamma - 1.0)/2.0 G8 = (gamma - 1.0)

QUSER = 2.0

! Guess P using the PVRS solver

CUP = 0.25*(rho_L + rho_R)*(Cs_L + Cs_R) PPV = 0.5*(p_L + p_R) + 0.5*(u_L - u_R)*CUP PPV = MAX(0.0, PPV) PMIN = MIN(p_L, p_R) PMAX = MAX(p_L, p_R) QMAX = PMAX/PMIN

IF (QMAX.LE.QUSER.AND.(PMIN.LE.PPV.AND.PPV.LE.PMAX)) THEN

!SELECT PVRS SOLVER pEst = PPV

ELSE

IF(PPV.LT.PMIN) THEN

!SELECT 2-RAREFACTION RS PQ = (p_L/p_R)G1 UM = (PQ*u_L/Cs_L + u_R/Cs_R + G4*(PQ - 1.0))/(PQ/Cs_L + 1.0/Cs_R) PTL = 1.0 + G7*(u_L - UM)/Cs_L PTR = 1.0 + G7*(UM - u_R)/Cs_R pEst = 0.5*(p_L*PTLG3 + p_R*PTRG3)

ELSE

!Select the 2-shock RS with PVRS as estimate

GEL = SQRT((G5/rho_L)/(G6*p_L + PPV)) GER = SQRT((G5/rho_R)/(G6*p_R + PPV)) pEst = (GEL*p_L + GER*p_R - (u_R - u_L))/(GEL + GER)

END IF

END IF

END SUBROUTINE P_estimate

SUBROUTINE PressFunct(f, fD, pItr, rho, p, Cs)

REAL, INTENT(OUT)

f, fd ! Pressure function and its derivative in Toro

REAL, INTENT(IN)

pItr ! Pressure correction as computed by Find_Pstar

REAL, INTENT(IN)

rho, p, Cs !Initial data

REAL

A, B ! Constants in pressure function

A = 2.0/((gamma+1.0)*rho) B = ((gamma-1.0)/(gamma+1.0))*p

IF(pItr.LT.p)THEN !Rarefaction wave

f = (2.0*Cs/(gamma - 1.0))*((pItr/p)((gamma-1.0)/(2.0*gamma)) - 1.0) fd = (1.0/(rho*Cs))*(pItr/p)(-(gamma+1)/(2*gamma))

ELSE

!Shock f = (pItr - p)*SQRT(A/(pItr+B)) fd = Sqrt(A /(B + pItr))*(1 - (pItr - p)/(2*(B+pItr)))

END IF

END SUBROUTINE PressFunct

END PROGRAM ExactRiemannSolver

Results

1. 

2. 

3. 

4.