Exact Riemann Solver​ (no vacuum)

The Riemann problem is the solution of a system of hyperbolic conservation laws (in our case the Euler equations). This solver assumes two one-dimensional ideal gases separated at the origin x=0, with initial densities, velocities, and pressures given as input. It gives a solution to the set of equations

\rho + \rho u = 0$ $\rho u + \rho u^2 + p = 0$ $E + u(E + p) = 0$

with initial conditions

\rho(x,0) = \begin{cases} \rho_L, & x < 0 \\ \rho_R, & x > 0 \end{cases}$ $\rho(x,0)u(x,0) = \begin{cases} \rho_L u_L, & x < 0 \\ \rho_R u_R, & x > 0 \end{cases}$ $E(x,0) = \begin{cases} E_L, & x < 0 \\ E_R, & x > 0 \end{cases}$

The solutions to this IVP at later times are split into four regions - to the left, the initial left state, unperturbed by interaction with the right; between the left and center-left regions, there will either be a rarefaction fan or a shock (assumed ideal), across which the new equilibrium left state holds; the left and right fluids will be separated by a contact discontinuity, across which the new equilibrium right state holds (note that the right and left equilibrium states differ only in their densities); the center-right region will then be separated by a rarefaction fan or (again ideal) shock from the rightmost region, which is still in the initial right state.

The program reads input from a file called problem.data, and stores many of the relevant parameters in global variables.

NAMELIST/ProblemData/ rhoL,uL,pL,rhoR,uR,pR,gamma,tol,starttime,finaltime,frames,xL,xR,res

inquire(file='problem.data', exist=ex)
if (ex) then
 OPEN(UNIT=30, FILE='problem.data', STATUS="OLD")
 READ(30,NML=ProblemData)                               ! read into variables in namelist
 CLOSE(30)
else
 print *, 'Missing problem.data. Exiting.'
 stop           ! exit if missing
end if

The state vectors W_L and W_R are created with overloaded constructors that also calculate useful parameters A and B and the speed of sound c for use later in the program. (Note code between comment lines is excerpted from its actual location.)

L = region(rhoL,uL,pL)
R = region(rhoR,uR,pR)

!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
function ConstructRegion(rho,u,p)
        
ConstructRegion%W(irho) = rho
ConstructRegion%W(iu) = u
ConstructRegion%W(ip) = p
ConstructRegion%A = 2/((gamma + 1)*ConstructRegion%W(irho))
ConstructRegion%B = (gamma - 1)*ConstructRegion%W(ip)/(gamma + 1)
ConstructRegion%c = sqrt(gamma*ConstructRegion%W(ip)/ConstructRegion%W(irho))
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
        
call validate(L,R)

!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
validate(L,R)

critdiff = 2.*L%c/(gamma - 1.) + 2.*R%c/(gamma - 1.)      ! Difference between right and left velocities must be strictly less than this to prevent vacuum formation
        
if (critdiff <= R%W(iu) - L%W(iu)) then
 print *, 'Vacuum results from these initial states. Exiting.'
 stop
end if
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

These state vectors are then passed to the Newton-Raphson solver, which calls a subroutine to guess a good value for p₀ (currently not very effective - just uses tol) and then iterates until the tolerance for change in pressure between steps is met.

pstar = pguess(L,R)
        
! do while loop (exits once change is less than tolerance)
do
 pstar_old = pstar
 pstar = pstar_old - ((f(pstar_old,L) + f(pstar_old,R) + R%W(iu)-L%W(iu))/(df(pstar_old,L) + df(pstar_old,R)))          ! Newton-Raphson formula
 delta = 2.*abs((pstar - pstar_old)/(pstar + pstar_old))                                                                ! Relative change in pressure from previous iteration
 if (delta < tol) exit                                                                                                  ! Exit upon sufficient accuracy
end do

The functions f and df (f') for either the left or right regions are given by

f_{reg}(p_*) = \begin{cases} (p_* - p_{reg}) [\frac{A_{reg}}{p_* + B_{reg})}], & p_* > p_{reg} \\ f_{reg}(p_*) = \frac{2 c_{reg}}{\gamma - 1}[(\frac{p_*}{p_{reg}})^{\frac{\gamma - 1}{2 \gamma}} - 1], & p_* \leq p_{reg} \end{cases}$ $f'_{reg}(p_*) = \begin{cases} (\frac{A_{reg}}{B_{reg} + p_*})^{\frac{1}{2}}[1 - \frac{p_* - p_{reg}}{2(B_{reg} + p_*)}], & p_* > p_{reg} \\ \frac{1}{\rho_{reg} c_{reg}}(\frac{p_*}{p_{reg}})^{-\frac{\gamma + 1}{2 \gamma}}, & p_* \leq p_{reg} \end{cases}$

See below for derivation of function f.

The velocity and densities for the star regions can then be calculated:

! Calculates velocity in star region
exactVelocity(L,R,pstar)
        
exactVelocity = 0.5*(L%W(iu) + R%W(iu)) + 0.5*(f(pstar,R) - f(pstar,L)) ! Velocity in star region is the same on either side of CD

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

! Calculates densities in star regions
exactDensity(starReg,reg)
        
! Density is different for shocks and rarefactions - may differ on either side of CD
if (starReg%W(ip) > reg%W(ip)) then
 exactDensity = reg%W(irho)*((gamma - 1.)/(gamma + 1.) + (starReg%W(ip)/reg%W(ip)))/(((gamma - 1.)/(gamma + 1.))*(starReg%W(ip)/reg%W(ip)) + 1.)
else
 exactDensity = reg%W(irho)*(starReg%W(ip)/reg%W(ip))**(1./gamma)
end if

The final value of pstar is printed to console, then output is written to files, one for each time state requested. During output, the speeds of shocks and rarefaction boundaries are calculated:

SL = L%W(iu) - L%c*sqrt((gamma + 1)*Lstar%W(ip)/(2*gamma*L%W(ip)) + (gamma - 1)/(2*gamma))      ! Left shock speed
SHL = L%W(iu) - L%c                                                                             ! Head speed of left rarefaction
STL = Lstar%W(iu) - L%c*(Lstar%W(ip)/L%W(ip))**((gamma - 1)/(2*gamma))                          ! Tail speed of left rarefaction
SR = R%W(iu) + R%c*sqrt((gamma + 1)*Rstar%W(ip)/(2*gamma*R%W(ip)) + (gamma - 1)/(2*gamma))      ! right
SHR = R%W(iu) + R%c                                                                             ! right
STR = Rstar%W(iu) + R%c*(Lstar%W(ip)/L%W(ip))**((gamma - 1)/(2*gamma))                          ! right

Then a large conditional is entered to determine which region we are sampling at the given time and position:

! End once t > finaltime (last output will be t = finaltime)
do while (t <= finaltime)
 !print fileformat, "out/result",frame,".dat"
 write(filename,fileformat) "out/result",frame,".dat"   ! filename to output - goes to out/resultxxxxx.dat
 open(unit=40,file=filename,iostat=stat,status='replace')
 write (40,*) 'x rho u p'                               ! write column headers
 x = xL
 ! Calculate values for each position - similar to time above
 !!! could parallelize if I changed from while to for loop
 do while (x <= xR)
  if (t == 0) then      !initial state (so we don't divide by zero)
   if (x < 0) write (40,*) x, L%W
   if (x > 0) write (40,*) x, R%W
   if (x == 0) write (40,*) x, 0, 0, 0
  else
   if (x/t <= Lstar%W(iu)) then !left side of contact
   lright = .false.
    if (Lstar%W(ip) > L%W(ip)) then !shock on left side
     if (x/t <= SL) then !left side of shock
      write (40,*) x, L%W
     else
      write (40,*) x, Lstar%W
     end if
    else !rarefaction on left side
     if (x/t <= SHL) then !left of the rarefaction fan
      write (40,*) x, L%W
     else if (x/t >= STL) then !right of the rarefaction fan
      write (40,*) x, Lstar%W
     else !in rarefaction fan
      call rarefaction(L,x,t,Wfan,lright)

!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
rarefaction(reg,x,t,Wfan,lright)

! Rarefaction fan differs depending on whether it's a right or left rarefaction
if (.not. lright) then
 Wfan(irho) = reg%W(irho)*(2/(gamma + 1) + (gamma - 1)*(reg%W(iu) - x/t)/((gamma + 1)*reg%c))**(2/(gamma + 1))
 Wfan(iu) = (2/(gamma + 1))*(reg%c + (gamma - 1)*reg%W(iu)/2 + x/t)
 Wfan(ip) = reg%W(ip)*(2/(gamma + 1) + (gamma - 1)*(reg%W(iu) - x/t)/((gamma + 1)*reg%c))**(2*gamma/(gamma + 1))
else
 Wfan(irho) = reg%W(irho)*(2/(gamma + 1) - (gamma - 1)*(reg%W(iu) - x/t)/((gamma + 1)*reg%c))**(2/(gamma + 1))
 Wfan(iu) = (2/(gamma + 1))*(-reg%c + (gamma - 1)*reg%W(iu)/2 + x/t)
 Wfan(ip) = reg%W(ip)*(2/(gamma + 1) - (gamma - 1)*(reg%W(iu) - x/t)/((gamma + 1)*reg%c))**(2*gamma/(gamma + 1))
end if
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

      write (40,*) x, Wfan
     end if
    end if
   ! Repeat entire sequence above for right side (lright = .true.)
   end if
  end if
  x = x + dx            ! increment x to next position
 end do
 t = t + dt             ! increment t to next time
 frame = frame + 1      ! update frame counter
end do

These output files can then be processed by VisIt, Mathematica, etc. to make nice pictures.

View ​online, or download 1​, 2​, 3​, 4​, 5​

Deriving the function f and the Newton-Raphson procedure

f

From the shock jump conditions

\rho_L \hat{u_L} = \rho_{*L} \hat{u_*}$ $\rho_L \hat{u_L}^2 + p_L = \rho_{*L} \hat{u_*}^2 + p_*$ $\hat{u_L}(\hat{E_L} + p_L) = \hat{u_*}(\hat{E_{*L}} + p_*)$

and the energy and EOS of an ideal gas

e = \frac{p}{(\gamma - 1) \rho}$ $p = \rho R T$

where \hat{u} = u_L - S_L is the velocity of the fluid in the frame of the shock, the function f_L(p_*) = u_L - u_*,

f_L(p_*) = (p_* - p_L) [\frac{2}{(\gamma + 1) \rho_L (p_* + \frac{\gamma - 1}{\gamma + 1} p_L)}]$

can be derived to connect the left state with the pressure in the central, interaction region in the case of a shock. Similarly, in the case of a rarefaction, instead of the shock jump conditions, the isentropic equation

\frac{p_*}{\rho_*} = \frac{p_L}{\rho_L}$

and the Generalized Riemann Invariant

u_L + \frac{2 c_L}{\gamma - 1} = u_* + \frac{2 c_{*L}}{\gamma - 1}$

where c = \sqrt{\frac{\gamma p}{\rho}} is the speed of sound in the region, can be used to construct

f_L(p_*) = \frac{2 c_L}{\gamma - 1}[(\frac{p_*}{p_L})^{\frac{\gamma - 1}{2 \gamma}} - 1]$

The same analysis can be performed on the right side of the system to yield the same results.

Newton-Raphson

The Riemann problem can then be solved by finding the value p = p_* such that

f_L(p_*) + f_R(p_*) + u_R - u_L = 0$

The function and its first and second derivatives are well-behaved and monotonic, so we can use a simple Newton-Raphson iterative method, derived from a Taylor series expansion of f(p_*):

If our initial guess is p_k and f(p_k) \ne 0, we want to change p_k so that f(p_k + \delta) = f(p_{k+1}) = 0; the Taylor expansion is given by

f(p_k + \delta) = f(p_k) + \delta f'(p_k) + O(\delta^2)$

Dropping second order and higher terms, since we have f(p_k + \delta) = 0, we know f(p_k) + \delta f'(p_k) = 0, and solving for \delta,

\delta = -\frac{f(p_k)}{f'(p_k)}$

and

p_{k + 1} = p_k + \delta = p_k - \frac{f(p_k)}{f'(p_k)}$

The iteration procedure is terminated at the term such that

|\frac{2 (p_{k + 1} - p_k)}{p_{k + 1} + p_k}| < Tol$

where Tol is the (generally small) tolerance passed to the solver.

Source: See Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics