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:

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

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

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

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

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.

Quantum Mechanics and The Next Steps

This week’s study day consisted of Lectures 8 and 9 from the Theoretical Minimum Quantum Mechanics course. Susskind has heavily emphasised Entanglement, since it is clearly a concept that causes significant confusion! I myself have been working through both the lecture courses and the accompanying book. The book is significantly more detailed than the lecture courses and includes many topics not covered in the videos. In particular, there is a brief treatment of the Quantum Harmonic Oscillator as well as Creation/Annhilation Operators.

Lecture 9 began the discussion of the (quantum) continuous particle in one dimension. Thus we are finally beginning to discuss position and momentum operators, as well as the time-dependent Schrödinger equation. Susskind spent a good deal of time discussing the Dirac Delta Function (and briefly Distribution Theory) and Fourier Series. Thankfully, these are concepts that I’m quite familiar with having studied Linear Analysis and Fourier Analysis the past in a mathematical setting. It was particularly pleasing to see the Fourier treatment applied to quantum mechanics.

I met up with a practising String Theorist a few nights ago, who happened to be working on the Ads/CFT Correspondance. It was fascinating to discuss Quantum Field Theory and how it relates to String Theory, at a level which was well beyond my current understanding! It was recommended that I take a look at the seminal two-volume work by Cohen-Tannoudji on Quantum Mechanics. I’ve also just received the first two volumes of Landau & Lifshitz in the post today, namely Mechanics and The Classical Theory of Fields.

The key question now is how to proceed upon completion of the Quantum Mechanics lectures on the Theoretical Minimum. The Special Relativity and Electrodynamics is an obvious next step, as is Landau & Lifshitz Vol II. However, I would like some slightly less daunting texts on electromagnetism and electrodynamics, prior to L&L Vol II. Griffiths is supposed to be a good undergraduate introduction, but from what I gather shies away from the Lagrangian/Hamiltonian approach, which is exactly what I need in order to tackle more advanced material.

Irrespective of the subject matter being studied, I definitely need to start tackling some actual questions. There are plenty of Classical Mechanics and Quantum Mechanics questions to consider. Also, I really want to start putting together my own notes into Classical Mechanics, that will ultimately become a course that will live on this site. I think I’m going to devote some of the next study day to putting together a brief syllabus.

Quantum Gravity or Numerical Relativity for Research?

For my study day this week I spent a lot of time working through Lecture 7 of the Theoretical Minimum Quantum Mechanics course. The theme of the lecture is Quantum Entanglement. At this stage we’re limiting ourselves to considering a singlet state, as I described in the previous post.

The thrust of the lecture involved Density Matrices and ultimately Bell’s Theorem. These are deep results and are not something which come easily, hence the need to repeat the material and carry it out over two lectures. It is interesting that having a bit of a background in stochastic processes and financial engineering certainly helps with the probabilistic concepts.

I’m glad that I’m finally beginning to get a solid grasp of elementary QM concepts. Once the course is over I’m going to work my way through the Special Relativity and Electrodynamics course. At this stage I really want to make sure I’m up to scratch on Classical Mechanics, QM, SR and EM by virtue of working through some textbooks and questions on the respective topics. Essentially, if I’m happy with the Landau & Lifshitz topics for Vols I (Classical Mechanics), II (Classical Theory of Fields) and III (Non-Relativistic Quantum Mechanics), then I feel I am in a good position to have a go at more complicated topics.

I am currently deciding on whether to pursue a deeper study of Quantum Field Theory or General Relativity, as a means of getting closer towards independent research. At this stage I am drawn towards both Quantum Gravity and Numerical Relativity as potential research areas. Perhaps I need not restrict myself to either. Nevertheless I am going to require a solid education in QFT, String Theory and GR in order to continue.

The next logical step is to “parallelise” learning of QFT and GR via the Theoretical Minimum lecture series, which has a specific GR module, and study QFT via some of the introductory textbooks. This is because there isn’t a specific QFT course at the Theoretical Minimum, although a lot of the necessary material is covered in the Advanced Quantum Mechanics course.

However, once I cover General Relativity via a number of textbooks (and even a workbook), I will have to decide whether to push hard on Numerical Relativity. This is a logical continuation of my PhD research into compressible fluid dynamics, so I should be able to begin researching more rapidly. But I am still very interested in the prospect of theoretical research into Quantum Gravity. Watch this space!