I’ve spent some of my spare time reading through Toro^{[1]} again, which provides a broad coverage of many different Riemann Solvers. I was keen to examine the solvers in depth so I could make a good choice for developing my upcoming 2D supersonic axisymmetric structured grid solver. In addition, it seems that the computational astrophysics community makes extensive use of the HLLC solver for relativistic hydrodynamics. Hence I wanted to try and implement a basic version in one-dimension, as an extension to the Roe solver presented in the previous article.

This morning I actually managed to implement the HLL solver, but not as yet the HLLC, as described in Toro^{[1]}, using Python. In addition I decided to modify the code base to include a new `utils.py`

file. This contains the `w2u`

and `u2w`

primitive/conservative variable conversion functions, as well as the specification of the Euler fluxes in conservative form. Since these are used across all Riemann solvers, it made sense to separate them out.

In addition to this I created a separate `hll_flux.py`

file, which contains the functions necessary to calculate the HLL fluxes. This keeps the main `1d_euler_hll.py`

file relatively clean. I will be continually adjusting the code-base such that code is refactored into its own modules or classes. I also modified the plotting function to create the classic “four-pane” of density, velocity, pressure and internal energy for a 1D problem, using markers as opposed to thick lines (as in the previous post). As yet I have not developed the exact Riemann solver solution, so I have not overlaid this onto the charts. That will come soon!

When developing the HLL solver I initially tried to make use of the entropy-based wave-speed estimates for calculating the largest wave speeds, for use in determining the HLL fluxes. The particle-speed estimate, $\tilde{u}$, and speed of sound estimate, $\tilde{a}$, needed for the wave speed estimates $S_L$ and $S_R$ are given below:

\begin{eqnarray}

\tilde{u} = \frac{\sqrt{\rho_L} u_L + \sqrt{\rho_R} u_R}{\sqrt{\rho_L} + \sqrt{\rho_R}}

\end{eqnarray}

And

\begin{eqnarray}

\tilde{a} = \left[ (\gamma - 1)(\tilde{H} - \frac{1}{2} \tilde{u}^2) \right]^{\frac{1}{2}}

\end{eqnarray}

Where the entropy estimate, $\tilde{H}$, is given by:

\begin{eqnarray}

\tilde{H} = \frac{\sqrt{\rho_L} H_L + \sqrt{\rho_R} H_R}{\sqrt{\rho_L} + \sqrt{\rho_R}}

\end{eqnarray}

With the enthalpy, $H = (E+p)/\rho$. This leads to estimates for the wave speeds as:

\begin{eqnarray}

S_L = \tilde{u} – \tilde{a},\;\; S_R = \tilde{u} + \tilde{a}

\end{eqnarray}

However, the output produced distinct gaps in the solutions of each of the quantities during the large rarefaction present in the modified Sod Test^{[1]}. This could be due to a bug on my part, but I checked the code twice. A third check is probably in order! You can see the four-pane output below:

Thus I decided that I would use the direct wave speeds instead, by directly sampling the left/right velocity and speeds of sound. Once I applied these wave speeds, the gaps disappeared. The output of this HLL solver on the Sod Test is given below:

Thus, the next set of tasks are as follows:

- Create the exact Riemann solver solution and apply it to all future four-panes
- Create the remaining shock-tube tests from Toro
^{[1]}and apply to HLL solver - Research the HLLC method (i.e. with added contact wave) and use approximate algorithm in Toro
^{[1]}to calculate wave speeds

In addition, I want to start looking at higher-order methods. In particular, I want to use a MUSCL-Hancock slope limiter approach, as my supervisor implemented for the GRP-based solver we had in our research group. This will provide second-order accuracy in time and space, which is a necessary requirement for making a robust flow simulator in the presence of subsonic and supersonic regions.

Once the higher order method is developed, tested and verified in 1D I can make use of Strang’s dimensional splitting technique to develop a 2D solver. In a previous post I mentioned CFD Books, which provides a 2D supersonic diffraction test case, with which I can compare to my own solver. At this stage I won’t be far off in modifying the solver to use the axisymmetric Euler equations and thus simulating a semi-cone intake for a supersonic air-breathing engine.

For reference, the code to produce these plots is given below. Firstly, `utils.py`

:

import numpy as np def w2u(w, gamma): """ Convert the primitive to conservative variables. """ u = np.zeros(3) u[0] = w[0] u[1] = w[0]*w[1] u[2] = w[2]/(gamma-1.0)+0.5*w[0]*w[1]*w[1] return u def u2w(u, gamma): """ Convert the conservative to primitive variables. """ w = np.zeros(3) w[0] = u[0] w[1] = u[1]/u[0] w[2] = (gamma-1.0)*( u[2] - 0.5*w[0]*w[1]*w[1] ) return w def internal_energy_e(w, gamma): """ Internal energy, e. See Toro, 2nd Ed., pg 3. """ return w[2]/(w[0]*(gamma - 1.0)) def total_energy_E(w, gamma): """ Total energy per unit volume, E. See Toro, 2nd Ed., pg 3. """ e = internal_energy_e(w, gamma) return w[0]*(0.5*w[1]*w[1] + e) def enthalpy_H(w, gamma): """ Calculate enthalpy, H = (E+p)/rho See Toro, 2nd Ed., pg 324. """ E = total_energy_E(w, gamma) return (E + w[2])/w[0] def euler_flux(w, gamma): """ Calculate the conservative Euler fluxes. """ rho = w[0] u = w[1] p = w[2] a2 = gamma*p/rho f_1 = rho*u f_2 = rho*u*u + p f_3 = rho*u*( a2/(gamma-1.0) + 0.5*u*u ) return np.array([f_1, f_2, f_3])

Secondly, the `hll_flux.py`

:

from math import sqrt import numpy as np from utils import total_energy_E, enthalpy_H, euler_flux, w2u def enthalpy_H_tilde(wL, wR, gamma): """ Calculate enthalpy estimate, H_tilde. See Toro, 2nd Ed., pg 324. """ H_L = enthalpy_H(wL, gamma) H_R = enthalpy_H(wR, gamma) num = sqrt(wL[0])*H_L + sqrt(wR[0])*H_R denom = sqrt(wL[0]) + sqrt(wR[0]) return num/denom def calc_u_tilde(wL, wR): """ Roe-average particle wave-speed estimate. See Toro, 2nd Ed., pg 324. """ num = sqrt(wL[0])*wL[1] + sqrt(wR[0])*wR[1] denom = sqrt(wL[0]) + sqrt(wR[0]) return num/denom def calc_a_tilde(wL, wR, gamma): """ Roe-average sound speed estimate. See Toro, 2nd Ed., pg 324. """ u_tilde = calc_u_tilde(wL, wR) H_tilde = enthalpy_H_tilde(wL, wR, gamma) return sqrt( (gamma - 1.0) * (H_tilde - 0.5 * u_tilde**2) ) def calc_wave_speeds_tilde(wL, wR, gamma): """ Calculate wave-speed estimates, S_L and S_R. See Toro, 2nd Ed., pg 324. """ u_tilde = calc_u_tilde(wL, wR) a_tilde = calc_a_tilde(wL, wR, gamma) S_L = u_tilde - a_tilde S_R = u_tilde + a_tilde return S_L, S_R def calc_wave_speeds_direct(wL, wR, gamma): """ Calculate wave-speeds directly. Returns S_L and S_R. See Toro, 2nd Ed., pg 324. """ S_L = wL[1] - sqrt(gamma*wL[2]/wL[0]) S_R = wR[1] + sqrt(gamma*wR[2]/wR[0]) return S_L, S_R def hll_flux(wL, wR, gamma): """ Use the HLL - Harten, Lax and van Leer solver to obtain the Euler fluxes See Toro, 2nd Ed., pg 320. """ 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) if 0.0 <= S_L: return F_L elif S_L <= 0.0 and 0.0 <= S_R: num = S_R*F_L - S_L*F_R + S_L*S_R*(U_R - U_L) denom = S_R - S_L return num/denom elif 0.0 >= S_R: return F_R

Thirdly, the `1d_euler.py`

:

import os.path import numpy as np from hll_flux import hll_flux from utils import w2u, u2w, internal_energy_e def test_case_sod(U, bcells, gamma): """ Populate the initial data with the Sod test case. """ mid_point = int(bcells/10.0*3.0) sod_L = np.array([1.0, 0.75, 1.0]) sod_R = np.array([0.125, 0.0, 0.1]) for i in range(0, mid_point): U[i,:] = w2u(sod_L, gamma) for i in range(mid_point, bcells): U[i,:] = w2u(sod_R, gamma) def calc_time_step(cfl, dx, gamma, bcells, U): """ Calculates the maximum wavespeeds and thus the timestep via an enforced CFL condition. """ max_speed = -1.0 for i in range(1,bcells-1): w = u2w(U[i], gamma) u = w[1] c = np.sqrt(gamma*w[2]/w[0]) max_speed = max(max_speed, abs(u)+c) dt = cfl*dx/max_speed # CFL condition return dt 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] = hll_flux(wL, wR, gamma) # Update solution for i in range(1, bcells-1): U[i] = U[i] + (dt/dx) * (fluxes[i-1]-fluxes[i]) # BCs U[0] = U[1] U[bcells-1] = U[bcells-2] if __name__ == "__main__": gamma = 1.4 cells = 100 bcells = cells + 2 dx = 1.0/cells cfl = 0.9 t = 0.0 tf = 0.2 nsteps = 0 U = np.zeros((bcells,3)) fluxes = np.zeros((bcells,3)) test_case_sod(U, bcells, gamma) for n in range(1, 50000): if (t==tf): break dt = calc_time_step(cfl, dx, gamma, bcells, U) if (t+dt > tf): dt = tf - t update_solution(U, fluxes, dt, dx, gamma, bcells) t += dt nsteps += 1 out_csv = open("output/1d_euler_hll.csv", "w") for elem in U: new_elem = u2w(elem, gamma) e = internal_energy_e(new_elem, gamma) out_csv.write("%s,%s,%s,%s\n" % ( new_elem[0], new_elem[1], new_elem[2], e) )

And finally, the updated plotting file `plot_1d_euler.py`

:

import sys import matplotlib.pyplot as plt import pandas as pd if __name__ == "__main__": filename = sys.argv[1] data = pd.io.parsers.read_csv( filename, delimiter=",", header=None, names=["rho", "u", "p", "e"], ) # Plot three charts fig = plt.figure() fig.patch.set_facecolor('white') # Set the outer colour to white # Plot the density ax1 = fig.add_subplot(221, ylabel='Density', ylim=[0.0, 1.2]) data['rho'].plot(ax=ax1, marker='o', markersize=4, color="black", linestyle='') # Plot the velocity ax2 = fig.add_subplot(222, ylabel='Velocity', ylim=[0.0, 1.5]) data['u'].plot(ax=ax2, marker='o', markersize=4, color="black", linestyle='') # Plot the pressure ax3 = fig.add_subplot(223, ylabel='Pressure', ylim=[0.0, 1.2]) data['p'].plot(ax=ax3, marker='o', markersize=4, color="black", linestyle='') # Plot the pressure ax4 = fig.add_subplot(224, ylabel='Internal Energy', ylim=[1.8, 3.8]) data['e'].plot(ax=ax4, marker='o', markersize=4, color="black", linestyle='') # Plot the figure plt.show()

References:

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