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.

Finishing Quantum Mechanics at the Theoretical Minimum

Yesterday I finished the final lecture in the Theoretical Minimum Quantum Mechanics course. I very much enjoyed the final two lectures as we finally began discussing the one-dimensional particle case. It was fascinating to see the interplay between the position space and momentum space, which was simply a manifestation of the Fourier Transform. I’m also gaining a good ability to calculate using the logic of QM, namely that of infinite-dimensional vector spaces. This isn’t surprising, since this was a big part of my undergraduate mathematics degree, but it was rather exciting to see it applied to less abstract situations.

While I’ve now finished the lectures, I still have to finish the accompanying text. I’m actually now onto the material covered in the final two lectures so am not too far behind the video content. The book does include some discussion on the quantum harmonic oscillator, as I alluded to in previous posts, so this will hopefully prepare me somewhat for tackling -very basic- introductory Quantum Field Theory texts.

Again, as I’ve alluded to in previous posts I’m now going to move onto the Special Relativity and Electrodynamics course, which goes into significant detail on Classical Field Theory via the Lagrangian/Hamiltonian approach. Prior to tackling this, however, I might brush up on Lagrangian Mechanics for continuous systems (at least for one dimension). The usual example given here is the derivation of the wave equation via the one-dimensional constrained vibrating string. There’s a good discussion of it in A Student’s Guide to Lagrangians and Hamiltonians.

Concurrently to the Theoretical Minimum lectures I am going to try working my way through Landau & Lifshitz Vol II (The Classical Theory of Fields). I’ve had some suggestions as to which sections are worthwhile to pursue. In addition I’ve been recommended the Feynmann Lectures on Physics as a more intuitive guide when L&L becomes a bit terse! Although the latter is somewhat pricey, so it may have to wait.

How to Create a Dedicated Study Day

My current employment situation means that I have a reasonable degree of flexibility when it comes to allocating my time. I’m going to assume that the majority of you have one or two free days every weekend, and/or a number of evenings over the week, that are not directly taken up with employment.

One of the ways in which I have rapidly accelerated my physics learning recently was to carve out a dedicated study day per week. In particular, I have chosen to work on physics every Monday, without fail (up to emergencies, of course!). Monday was chosen simply because it gave me something highly enjoyable to look forward to after the weekend. I cannot emphasise enough how useful this approach has been in being able to actually make significant progress.

The most important part of attempting to create dedicated study time is to make a commitment to yourself that you will actually carry out the study every week, without fail, at that particular time. Consistency is the key to making good progress. Obviously, this is not straightforward when one has a busy lifestyle and family obligations! I have found that because I look forward to this particular day every week, that I never postpone it. I guess this is simply a question of motivation!

So what should you do on a study day? This will clearly depend upon your personal situation, but I have found that watching one or two lectures (between 1-2 hours in length) and making extensive notes is extremely helpful. I don’t overwhelm myself with trying to learn too much in a day and I make consistent progress over the longer term. For instance, at the time of writing, I’ve almost finished my second course on the Theoretical Minimum after only a couple of months.

The key is really to be consistent. Even a single day taken every week, over the course of a few years, can add up to a lot of study.