In the previous article I discussed the basic road-map for creating a TVD HLLC Riemann solver that would work in two-dimensions. I’ve now managed to implement both the exact Riemann solver, as given by Toro^{[1]} as well as the implementation of the first-order HLLC. Hence, at least for the modified Sod Test^{[1]} I can display the four-pane of density, velocity, pressure and internal energy.

The fluxes for the HLLC method are given by:

\begin{eqnarray}

{\bf F}^{hllc}_{i+\frac{1}{2}} = \left\{

\begin{array}{lr}

{\bf F}_L & 0 \le S_L \\

{\bf F}_{*L} = {\bf F}_L + S_L({\bf U}_{*L} – U_L) & S_L \le 0 \le S_{*} \\

{\bf F}_{*R} = {\bf F}_R + S_R({\bf U}_{*R} – U_R) & S_{*} \le 0 \le S_{R} \\

{\bf F}_{R} & 0 \ge S_{R}

\end{array}

\right.\\

\end{eqnarray}

In order to calculate them it is necessary to compute $S_{*}$, ${\bf U}_{*L}$ and ${\bf U}_{*R}$. Firstly, we calculate $S_{*}$ via:

\begin{eqnarray}

S_{*} = \frac{p_R – p_L + \rho_L u_L (S_L – u_L) – \rho_R u_R (S_R – u_R)}{\rho_L (S_L – u_L) – \rho_R (S_R – u_R)}

\end{eqnarray}

Secondly, we calculate each of the velocities in the starred region, on either side of the contact discontinuity, $K \in \{L, R \}$:

\begin{eqnarray}

{\bf U}_{*K} = \rho_K \left( \frac{S_K – u_K}{S_K – S_{*}} \right) \left[

\begin{array}{c}

1 \\

S_{*} \\

\frac{E_K}{\rho_K} + (S_{*} -u_K) \left[ S_{*} + \frac{p_K}{\rho_K (S_K - u_K)} \right] \\

\end{array}

\right]\\

\end{eqnarray}

Thus the fluxes are now fully determined.

The code to carry out these calculations has been added to `hll_fluxes.py`

, as discussed in the previous article:

def calc_S_star(wL, wR, S_L, S_R): """ Calculate the wave speed S_{*} in the starred region. See Toro, 2nd Ed., pg 327, eqn 10.58. """ dL = wL[0] uL = wL[1] pL = wL[2] dR = wR[0] uR = wR[1] pR = wR[2] num = pR - pL + dL*uL*(S_L - uL) - dR*uR*(S_R - uR) denom = dL*(S_L - uL) - dR*(S_R - uR) return num/denom def calc_U_k(wK, S_K, S_star, gamma): """ Calculate the velocities in the starred region on either side of the contact discontinuity, U_L_star and U_R_star. See Toro, 2nd Ed., pg 323, eqn 10.33. """ dK = wK[0] uK = wK[1] pK = wK[2] EK = total_energy_E(wK, gamma) coef = dK*((S_K - uK)/(S_K - S_star)) vec = np.array([ 1.0, S_star, EK/dK + (S_star - uK)*(S_star + (pK/(dK*(S_K-uK)))) ]) return coef*vec def hllc_flux(wL, wR, gamma): """ Use the HLLC - HLL extended with contact, via Toro, to obtain the Euler fluxes See Toro, 2nd Ed., pg 323. """ S_L, S_R = calc_wave_speeds_direct(wL, wR, gamma) F_L = euler_flux(wL, gamma) F_R = euler_flux(wR, gamma) U_L = w2u(wL, gamma) U_R = w2u(wR, gamma) S_star = calc_S_star(wL, wR, S_L, S_R) U_L_star = calc_U_k(wL, S_L, S_star, gamma) U_R_star = calc_U_k(wR, S_R, S_star, gamma) F_L_star = F_L + S_L*(U_L_star - U_L) F_R_star = F_R + S_R*(U_R_star - U_R) if 0.0 <= S_L: return F_L elif S_L <= 0.0 and 0.0 <= S_star: return F_L_star elif S_star <= 0.0 and 0.0 <= S_R: return F_R_star elif 0.0 >= S_R: return F_R

It only remained to change `1d_euler_hll.py`

to import the necessary flux calculation and adjust the flux to use the HLLC variant:

.. from hll_flux import hllc_flux .. .. def update_solution(U, fluxes, dt, dx, gamma, bcells): """ Updates the solution of the equation via the Godunov procedure. """ # Create fluxes for i in range(0, bcells-1): wL = u2w(U[i], gamma) wR = u2w(U[i+1], gamma) fluxes[i] = hllc_flux(wL, wR, gamma) .. ..

In the previous post I mentioned wanting to code up an exact Riemann solver so that I could compare its solution with the approximate solvers created here. I utilised the listing given in Toro^{[1]}, but coded it up in Python. Here is the full listing for `riemann_exact.py`

:

import sys from math import sqrt # Gammas gamma = 1.4 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 def test_pressure_positivity(cL, cR, uL, uR): if (g4*(cL+cR) <= (uR-uL)): print "Riemann Exact: Vacuum generated by data!" sys.exit() def guessp(dL, dR, uL, uR, pL, pR, cL, cR): quser = 2.0 cup = 0.25*(dL + dR)*(cL + cR) ppv = 0.5*(pL + pR) + 0.5*(uL - uR)*cup ppv = max(0.0, ppv) pmin = min(pL, pR) pmax = max(pL, pR) qmax = pmax/pmin if ( qmax <= quser and ( pmin <= ppv and ppv <= pmax ) ): # Select PVRS Riemann Solver return ppv else: # Select Two-Rarefaction Riemann Sover if ppv < pmin: pq = (pL/pR)**g1 uM = (pq*uL/cL + uR/cR + g4*(pq - 1.0))/(pq/cL + 1.0/cR) ptL = 1.0 + g7*(uL - uM)/cL ptR = 1.0 + g7*(uM - uR)/cR return 0.5*(pL*ptL**g3 + pR*ptR**g3) else: # Select Two-Shock Riemann Solver with # PVRS as estimate geL = sqrt((g5/dL)/(g6*pL + ppv)) geR = sqrt((g5/dR)/(g6*pR + ppv)) return (geL*pL + geR*pR - (uR - uL))/(geL + geR) def star_pu(dL, dR, uL, uR, pL, pR, cL, cR): tol_pre = 1.0e-6 nr_iter = 20 p_start = guessp(dL, dR, uL, uR, pL, pR, cL, cR) print "p_start: %s" % p_start p_old = p_start u_diff = uR - uL for i in range(1, nr_iter+1): fL, fLd = prefun(p_old, dL, pL, cL) fR, fRd = prefun(p_old, dR, pR, cR) p = p_old - (fL + fR + u_diff)/(fLd + fRd) change = 2.0*abs((p-p_old)/(p+p_old)) if change <= tol_pre: break if p < 0.0: p = tol_pre p_old = p u = 0.5*(uL + uR + fR - fL) return p, u def prefun(pT, dk, pk, ck): if pT <= pk: # Rarefaction wave prat = pT/pk f = g4*ck*(prat**g1 - 1.0) fd = (1.0/(dk*ck))*prat**(-g2) return f, fd else: # Shock wave ak = g5/dk bk = g6*pk qrt = sqrt(ak/(bk+pT)) f = (pT-pk)*qrt fd = (1.0 - 0.5*(pT - pk)/(bk + pT))*qrt return f, fd def sample(pM, uM, s): if s <= uM: if pM <= pL: shL = uL - cL if s <= shL: d = dL u = uL p = pL else: cmL = cL*(pM/pL)**g1 stL = uM - cmL if s > stL: d = dL*(pM/pL)**(1.0/gamma) u = uM p = pM else: u = g5*(cL + g7*uL + s) c = g5*(cL + g7*(uL - s)) d = dL*(c/cL)**g4 p = pL*(c/cL)**g3 else: pmL = pM/pL sL = uL - cL*sqrt(g2*pmL + g1) if s <= sL: d = dL u = uL p = pL else: d = dL*(pmL + g6)/(pmL*g6 + 1.0) u = uM p = pM else: if pM > pR: pmR = pM/pR sR = uR + cR*sqrt(g2*pmR + g1) if s >= sR: d = dR u = uR p = pR else: d = dR*(pmR + g6)/(pmR*g6 + 1.0) u = uM p = pM else: shR = uR + cR if s >= shR: d = dR u = uR p = pR else: cmR = cR*(pM/pR)**g1 stR = uM + cmR if s <= stR: d = dR*(pM/pR)**(1.0/gamma) u = uM p = pM else: u = g5*(-cR + g7*uR + s) c = g5*(cR - g7*(uR - s)) d = dR*(c/cR)**g4 p = pR*(c/cR)**g3 return d, u, p if __name__ == "__main__": dom_len = 1.0 diaph = 0.3 cells = 100 time_out = 0.2 # Sod Test dL = 1.0 uL = 0.75 pL = 1.0 dR = 0.125 uR = 0.0 pR = 0.1 # Sound speeds cL = sqrt(gamma*pL/dL) cR = sqrt(gamma*pR/dR) dS = 0.0 uS = 0.0 pS = 0.0 test_pressure_positivity(cL, cR, uL, uR) pM, uM = star_pu(dL, dR, uL, uR, pL, pR, cL, cR) print "pM, uM: %s, %s" % (pM, uM) dx = dom_len/cells riem_file = open("output/riemann_exact.csv", "w") for i in range(1, cells+1): x_pos = (float(i) - 0.5)*dx s = (x_pos - diaph)/time_out d, u, p = sample(pM, uM, s) riem_file.write("%s,%s,%s,%s\n" % (d, u, p, p/d/g8) ) riem_file.close()

After running both `riemann_exact.py`

and `1d_euler_hll.py`

I was able to produce the following four-pane:

It appears to be in complete agreement with the visualisation provided in Toro^{[1]} for the modified Sod Test.

I will be the first to admit that these implementations in Python are far from production ready in an execution speed sense. While they produce the correct code, they are not making use of any of the real optimisations provided by NumPy or SciPy. So what is the rationale for using Python? My reasons are given below:

- Python is quick to develop in. There is little "boilerplate" to worry about, as there is in C++ or Fortran.
- Python has mature and extensive plotting facilities, via Matplotlib, while C++ and Fortran do not. For the latter GNUPlot or TecPlot are often employed.
- Python is
*fast enough*to give me a baseline to work against when I finally come to make a production-ready, high-performance solution in C++ or Fortran.

In order to be fully happy with the HLLC implementation I am going to implement the remaining numerical tests and apply them to this solver, comparing with those in Toro^{[1]}. Subsequent to this (assuming no errors!) I will be extending the HLLC solver to a MUSCL-Hancock TVD approach, making the solver second order in both space and time. This will provide a robust solver that I can then use for a two-dimensional axisymmetric simulator.

**References**

- [1] - Toro, E.F.,
*Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction, 3rd Ed.*, Springer