HLLC-Based Riemann Solver and Exact Solution

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:

1D Euler HLLC-Based Riemann Solver Sod Test, First Order

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

Riemann Solver Update and Relativistic Hydrodynamics

I thought I’d provide a brief update of where I’m currently at with the Riemann solver learning and implementation.

HLLC Riemann Solver

I’ve been studying Toro[1] extensively over the last week or so. In particular, I’ve been looking at how to implement the HLLC solver in one-dimension, in a first order setting. This is essentially a small step-up from the HLL solver, which introduces a contact wave between the left/right shock/expansions. Doing so makes the physics more realistic. At this stage the HLLC will be my preferred choice for approximate Riemann solver going forward.

Once the HLLC is implemented, I want to create a MUSCL-Hancock Total Variation Diminishing (TVD) second order (time and space) solver using HLLC, still in one dimension. This will make use of slope limiters to create the TVD aspect of the code. In order to verify the accuracy I can compare it to the first order HLLC solver using the list of numerical tests provided in Toro. If successful, we should see minimal spurious oscillations (Gibbs phenomena) and less smearing of contact and shock discontinuities.

At this stage we will have a robust second-order 1D Riemann solver using the Godunov method.

I’ve also been looking into multi-dimensional methods using such a solver. Toro has a nice discussion on using Strang’s dimensional splitting technique as well as a Finite Volume Method (FVM). Reading this chapter actually clarified a few points for me about such solvers are constructed.

Since my PhD supervisor used a FVM approach, I will try and replicate that for the TVD HLLC in two-dimensions using the standard Euler equations. Eventually I will need to solve the axisymmetric Euler equations, which believe slightly differently. Using a FVM approach allows me to relatively straightforwardly extend the solver to use non-Cartesian geometries, which will lead naturally to a block-structured approach.

I’m concerned that a pure Python implementation, as I’m running now, will be too slow for a useful simulator. Hence I will likely reprogram the method in Fortran or C++. Having the Python implementation will give me a baseline from which to work from and an automatic test case to compare it with.

Relativistic Hydrodynamics

I’ve also purchased Relativistic Hydrodynamics[2] by Rezzolla and Zanotti. It is an absolutely fantastic book, which covers a significant amount of ground. It starts with some basic review of general relativity concepts and then discusses statistical physics and kinetic theory in a non-relativistic setting, leading to the derivation of the Euler equations. This is extremely useful because once relativistic considerations are taken into account the derivation becomes somewhat familiar.

I have not studied statistical physics in any formal setting before so I’m finding the chapter quite hard going. Hence I’m very much looking forward to taking the Theoretical Minimum course on Statistical Mechanics, which I’m hoping will clear up a lot of the issues! Also, Landau & Lifshitz have a volume on Statistical Physics (V)[3] so I will be purchasing that at some stage.

I’ve also had a brief review of the chapter in Relativistic Hydrodynamics on high-resolution shock capturing schemes (i.e. everything discussed above!), but it has also made me aware of “high resolution methods”, which it defines as methods with order greater or equal to 3. I have not come across these methods before (ENO, WENO and Discontinuous Galerkin) so I’m very excited to start browsing the literature on such techniques.

References

  • [1] – Toro, E.F., Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction, 3rd Ed., Springer
  • [2] – Rezzolla, L., Zanotti, O., Relativistic Hydrodynamics, CUP Oxford
  • [3] – Landau, L. D., Lifshitz, E. M., Statistical Physics: Volume 5 (Course on Theoretical Physics), Butterworth-Heinemann

HLL Riemann Solver Implementation in Python

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:

1D Euler HLL-Based Riemann Solver Sod Test With Estimated Wave Speeds

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:

1D Euler HLL-Based Riemann Solver Sod Test With Estimated Wave Speeds

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

Theoretical Physics Update

Unfortunately, I’ve been traveling a fair bit in the last couple of months, as well as finishing off the writing of a book for one of my e-commerce sites. Hence I’ve not had as much time as I would like to work towards theoretical physics learning.

However, I’m going to be restarting the Theoretical Minimum courses within the next week or so, beginning with the Special Relativity and Electrodynamics course. I’m very much looking forward to this because it will be my first direct study of electromagnetism in a formal (i.e. lecture-based) setting. Not only that, but I’m keen to see how the Lagrangian approach is utilised for classical fields.

In addition, the special relativity aspect will be extremely useful when I come to expand the Riemann solvers discussed previously to a relativistic setting.

Certainly a lot to look forward to!

A One-Dimensional Python-Based Riemann Solver

I’ve not updated the site in a month or so. Unfortunately, despite my prior article about carving out a study day, life sometimes gets in the way! I’ve been spending a great deal of time finishing off the writing of a book for QuantStart, an entrepreneurial venture of mine in the quantitative finance realm. Despite a busy writing schedule I have spent some time on physics, if only reading in a more sedate manner. I’ve been trying to follow Quantum Field Theory For The Gifted Amateur, which is an absolutely fantastic higher level, albeit still quite mathematically technical, overview of Quantum Field Theory. More on that in later posts.

My Compressible Computational Fluid Dynamics Research

This post is concerned with Computational Fluid Dynamics (CFD), which I believe to be relevant here because it is fundamentally computational physics, and thus deserves space on the blog! In particular it is a continuation of my postgraduate aeronautics research at Imperial College. I’ve not looked at my PhD work for some time but I thought I would take some of the content I researched and wrap it up into a paper, which I will attempt to submit to a journal. This particular work concerned idealised two-dimensional axisymmetric supersonic air-breathing engine (“ramjet”) intake simulations under conditions that are “off-design”, namely where the compression process is not optimised. This leads to some interesting physical phenomena in the inviscid Euler regime, including strong-shock Mach stems and shock-expansion fan interactions.

Of particular interest to me is a hysteresis effect that occurs in a certain domain of intake cone-angle and free-stream Mach number. My supervisor and I were able to show that such a hysteresis effect can indeed be simulated numerically. Further, we were able to map out the regime in the cone-angle/Mach number parameter space where such hysteresis occurred. The domain was bounded by the so-called Detachment Criterion and Von-Neumann Criterion.

The goal of this article is to describe how I plan to replicate this result with more modern tools, e.g. Python and C++, and then extend the result to multiple engine geometries, potentially including curved surfaces within the intakes. This would allow an isentropic intake cone as well as a curved cowl to replicate the shock-on-lip and shock-on-shoulder conditions.

There are two components to the simulation. The first will be the design and implementation of a finite-volume method (FVM) compressible inviscid solver, using a Godunov-type method and a Riemann solver. My supervisor previously had a sophisticated Fortran RANS-based turbulence code that utilised the Generalised Riemann Problem of Ben Artzi and Falcovitz[1]. I plan to build my own in Python, due to reasons of development speed and testing. Once I am happy with the results, I will replicate the code in C++ and optimise, making sure to compare results with the Python code at every step in order to ensure validity. This code will solve the two-dimensional axisymmetric Euler equations in an inhomogeneous “stretched” block-structured fashion, which makes the assumption that the tangent free-stream flow vector is parallel to the cone axis of symmetry, i.e. zero angle of attack. The output will consist of a velocity, density and pressure fields. A suitable visualisation tool will also be needed in the absence of the commercial TecPlot code!

The second code will utilise the Method of Characteristics (MoC) to produce a highly performant analytic solution to flow-field values at suitable domains within an engine. I had previously written such a code in C++, known as the Axisymmetric Method of Characteristics Solver, or ARMOCS. I aim to re-implement this in Python, since execution speed will be more than sufficient for this task. The goal is then to produce the domain of duality that sits between the Detachment Criterion and the Von Neumann Criterion in free-stream Mach number/cone-angle space. Once these boundaries have been established, the FVM code can perform the actual simulation to check whether the hysteresis effect manifests itself. I have already established that this is the case for the particular geometry I studied for my PhD, but I plan to consider more sophisticated regimes and make the code more general.

Much has changed since I carried out the research in 2008 (six years ago!). In the interim more Riemann Solvers have been developed. Hence I have more scope to “swap out” the underlying solver being used in the FVM simulation. In addition, computational speed and random-access memory have moved on substantially. Thus I will able to try more highly resolved meshes, albeit with problems identified by Quirk[3]. I have not yet carried out an up-to-date literature review on the latest shockwave/expansion-fan literature, so the first step will be to bring myself up to speed.

Python FVM Godunov-Type Roe-Based Euler Solver

I’ve already begun the process by constructing a one-dimensional Roe-based Godunov solver for the inviscid Euler equations in one-dimension. An interval $[0,1] \in \mathbb{R}$ is divided into $N$ cells, each of which contains a discretised pressure, velocity and density field. Such a grouping of the three fields is known as a primitive formulation. The Godunov procedure of updating the relations is recast into the conservative formulation, which essentially provides a (discretised) density field, a momentum field and an internal energy field. The conservative formulation of the Euler equations is given by:

\begin{eqnarray}
\frac{\partial U}{\partial t} = \frac{\partial F(U)}{\partial x}
\end{eqnarray}

Where $F(U) = (\rho, \rho u, )$ are the fluxes. The numerical evolution of the system is governed by the Godunov updating procedure:

\begin{eqnarray}
U_i^{n+1} = \frac{\Delta t}{\Delta x} (F(U_{i+\frac{1}{2}}^n) – F(U_{i-\frac{1}{2}}^n))
\end{eqnarray}

For this simple situation we are considering a modified version of the classic Sod Test, which is found in Chapter 6 of Toro[4]. We choose a CFL number of 0.9, a diaphragm point of the “Riemann shock tube” of $x = 0.3$, for the interval $[0,1]$, a primitive left state $W_L$ of $(\rho_L, u_L, p_L) = (1.0, 0.75, 1.0)$ and a primitive right state $W_R$ of $(\rho_R, u_R, p_R) = (0.1, 0.125, 0)$.

The first step in the code is to import the necessary Python libraries. We simply require NumPy and os.path:

import os.path

import numpy as np

It is useful to have functions that can convert between primitive states $W$ and conservative variable states $U$:

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

The flux of the conservative variables is calculated using the Euler equations:

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])

The majority of the code is concerned with an implementation of a Roe-based approximate Riemann solver:

def roe_flux(wL, wR, gamma):
    """
    Use the Roe approximate Riemann solver to calculate fluxes.
    """
    uL = w2u(wL, gamma)
    uR = w2u(wR, gamma)

    # Primitive and other variables.
    # Left state
    rhoL = wL[0]
    vL = wL[1]
    pL = wL[2]
    aL = np.sqrt(gamma*pL/rhoL)
    HL = ( uL[2] + pL ) / rhoL

    # Right state
    rhoR = wR[0]
    vR = wR[1]
    pR = wR[2]
    aR = np.sqrt(gamma*pR/rhoR)
    HR = ( uR[2] + pR ) / rhoR

    # First compute the Roe Averages
    RT = np.sqrt(rhoR/rhoL);
    rho = RT*rhoL
    v = (vL+RT*vR)/(1.0+RT)
    H = (HL+RT*HR)/(1.0+RT)
    a = np.sqrt( (gamma-1.0)*(H-0.5*v*v) )

    # Differences in primitive variables.
    drho = rhoR - rhoL
    du = vR - vL
    dP = pR - pL

    # Wave strength (Characteristic Variables).
    dV = np.array([0.0,0.0,0.0])
    dV[0] = 0.5*(dP-rho*a*du)/(a*a)
    dV[1] = -( dP/(a*a) - drho )
    dV[2] = 0.5*(dP+rho*a*du)/(a*a)

    # Absolute values of the wave speeds (Eigenvalues)
    ws = np.array([0.0,0.0,0.0])
    ws[0] = abs(v-a)
    ws[1] = abs(v)
    ws[2] = abs(v+a)

    # Modified wave speeds for nonlinear fields (the so-called entropy fix, which
    # is often implemented to remove non-physical expansion shocks).
    # There are various ways to implement the entropy fix. This is just one
    # example. Try turn this off. The solution may be more accurate.
    Da = max(0.0, 4.0*((vR-aR)-(vL-aL)) )
    if (ws[0] < 0.5*Da):
        ws[0] = ws[0]*ws[0]/Da + 0.25*Da
    Da = max(0.0, 4.0*((vR+aR)-(vL+aL)) )
    if (ws[2] < 0.5*Da):
        ws[2] = ws[2]*ws[2]/Da + 0.25*Da

    # Right eigenvectors
    R = np.zeros((3,3))

    R[0][0] = 1.0
    R[1][0] = v - a
    R[2][0] = H - v*a

    R[0][1] = 1.0
    R[1][1] = v
    R[2][1] = 0.5*v*v

    R[0][2] = 1.0
    R[1][2] = v + a
    R[2][2] = H + v*a

    # Compute the average flux.
    flux = 0.5*( euler_flux(wL, gamma) + euler_flux(wR, gamma) )

    # Add the matrix dissipation term to complete the Roe flux.
    for i in range(0,3):
        for j in range(0,3):
            flux[i] = flux[i] - 0.5*ws[j]*dV[j]*R[i][j]
    return flux

We need to enforce a CFL condition to maintain numerical stability:

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

We then carry out the Godunov temporal updating procedure:

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] = roe_flux(wL, wR, gamma)

    # Update solution
    for i in range(1, bcells-1):
        U[i] = U[i] + (dt/dx) * (fluxes[i-1]-fluxes[i])

    # Boundary Conditions
    U[0] = U[1]
    U[bcells-1] = U[bcells-2]

Finally we wrap the previous functions in a __main__ block and output the primitive variables to disk:

if __name__ == "__main__":
    gamma = 1.4
    cells = 1000
    bcells = cells + 2
    dx = 1.0/cells

    cfl = 0.8
    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_roe_godunov_euler_sod.csv", "w")
    for elem in U:
        new_elem = u2w(elem, gamma)
        out_csv.write("%s,%s,%s\n" % tuple(new_elem))

We now plot a three-pane of density, velocity and pressure. Comparison with the Sod Test in Toro[4] gives full agreement for $N=1000$. The matplotlib code for the output is as follows:

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"],
    )

    # Plot three charts
    fig = plt.figure()
    fig.patch.set_facecolor('white')     # Set the outer colour to white
    
    # Plot the density
    ax1 = fig.add_subplot(311, ylabel='Density', ylim=[0.0, 1.2])
    data['rho'].plot(ax=ax1, color="black", lw=2.)

    # Plot the velocity
    ax2 = fig.add_subplot(312, ylabel='Velocity', ylim=[0.0, 1.5])
    data['u'].plot(ax=ax2, color="black", lw=2.)

    # Plot the pressure
    ax3 = fig.add_subplot(313, ylabel='Pressure', ylim=[0.0, 1.2])
    data['p'].plot(ax=ax3, color="black", lw=2.)

    # Plot the figure
    plt.show()

The Sod Test visualisation is given below:

1D Euler Roe-Based Riemann Solver Sod Test Scalar Quantities

While this is a "good start" there are many improvements needed before being able to simulate a full isentropic supersonic intake:

  • I need to consider the two-dimensional axisymmetric Euler equations on a homogeneous grid. For this I will employ Strang's operator splitting technique, that allows a one-dimensional "sweep" in each spatial direction.
  • I need to use a more robust second-order Riemann solver, such as the GRP. This will take some time to implement as it is a non-trivial algorithm. Thankfully the Ben-Artzi paper[] and monograph provides a detailed overview of how to carry it out. I would also like to try the HLLE/HLLC solvers of Harten, Lax, van Leer, Einfeldt and Toro.
  • While I will maintain a structured grid, I would like to add grid inhomogeneity, such that certain areas can receive a higher resolution than others. This is appropriate at the cowl lip and cowl shoulder. Eventually I would like to use a fully unstructured adaptive mesh refinement code, but this is some time away!
  • The grid will of course need to have a non-Cartesian geometry to simulate the isometric intake and curved cowl.

In the next post I will discuss the development of the 2D axisymmetric solver using operator splitting.