As part of my “hobbyist” foray into Numerical Relativity (NR) I wanted to refamiliarise myself with the research I carried out at the Hypersonics Group of Imperial College as part of my doctorate. I was investigating ramjet intakes in a high supersonic/low hypersonic inviscid regime using a Generalised Riemann Problem (GRP) based Godunov Scheme Riemann Solver numerical code written in Fortran.

It didn’t take me long to realise that there was a significant hydrodynamic component to Numerical Relativity, which involved some very sophisticated computational fluid dynamics (CFD). In particular the field is known as relativistic hydrodynamics and encompasses modelling of not only compressible fluids but also those subject to magnetic field interactions at relativistic velocities. Hence we are in a relativistic compressible magnetohydrodynamic (MHD) regime with dynamically varying spacetime geometry. As you can imagine this makes the simulations rather complex.

If I am to be usefully contributing to this field in the future it is absolutely essential that I “continue where I left off” in regards to aerodynamic modelling and bridge the gap with NR. Thus I have set myself an additional goal of creating a complex CFD solver, for the particular field I was looking at, as well as gaining additional insight into the MHD regime.

In particular I have been working towards constructing a two-dimensional axisymmetric simulation of an isentropic ramjet intake with external cowl using Fortran as the implementation language and Python as the visualisation tool. The finite-volume method simulation will ultimately have the following features:

- Utilises the latest approximate Riemann solvers in the literature, which include the aforementioned GRP and the Rotated-RHLL
^{[1]} due to Nishikawa and Kitamura.
- Flexible triangular-element unstructured grid, useful for more general geometries.
- Second-order accuracy via slope limiters.
- Will be an explicit Euler simulator, rather than Navier-Stokes (i.e. viscosity) at this stage. Later implementations will introduce laminar, and ultimately turbulent, viscosity.

To this end I have been researching around the internet to see what open source code and tutorials are already available in order to save time. I actually managed to find an exceptional CFD teaching resource called CFD Books by Katate Masatsuka (which is actually a pen name!). The site has many extremely well commented Fortran 90 codes for various compressible CFD implementations using Riemann solvers.

In particular there is an example of a 2D unstructured explicit Euler solver using Roe/Rotated-RHLL approximate Riemann solvers, a Van Albada slope limiter and a two-step Runge-Kutta timestepper. The example uses a classic supersonic compressible flow test known as a shock-diffraction problem. It is a great way to learn how such CFD codes work as all of the Fortran code is clearly commented and very well written.

There is also a one-dimensional example of Sod’s test (a classic Riemann solver test case, which has an exact solution) using a Roe approximate solver, a Minmod slope limiter and a two-step Runge-Kutta timestepper. I actually spent yesterday coding up this example myself, trying to follow along the whole way in order to see if I can replicate the results. It was a very beneficial experience, since I had previously not looked at some Fortran code in earnest in approximately 5 years!

## Installing a Python Scientific Stack

I thought it would also be instructive to show you how I managed to install a Python and Fortran scientific stack on my Ubuntu 14.04 Desktop machine so that you could also replicate the results.

Firstly we update and upgrade Ubuntu:

sudo apt-get update && sudo apt-get upgrade

Then we install the necessary packages for working with Python and Fortran:

sudo apt-get install build-essential git-core python2.7-dev python-dev python-pip liblapack-dev libblas-dev libatlas-base-dev gfortran libpng-dev libjpeg8-dev libfreetype6-dev libqt4-core libqt4-gui libqt4-dev libzmq-dev

This is my approach to creating a set of directories with a git repository:

mkdir -p ~/science/aerocfd/projroot
mkdir ~/science/aerocfd/data
mkdir ~/science/aerocfd/tmp
mkdir ~/science/aerocfd/test
cd ~/science/aerocfd/projroot
git init .

Then a space to store the 2D Unstructured Euler Shock Diffraction problem from CFDBooks:

cd ~/science/aerocfd/test
mkdir -p idolikecfd/2d_unstructured_euler_shock_diffraction
cd idolikecfd/2d_unstructured_euler_shock_diffraction

Then we download all of the Fortran files from CFDBooks:

wget http://www.cfdbooks.com/cfdcodes/OSSAN_Euler2D/data_package_v0.f90
wget http://www.cfdbooks.com/cfdcodes/OSSAN_Euler2D/euler_solver_v0.f90
wget http://www.cfdbooks.com/cfdcodes/OSSAN_Euler2D/main_v0.f90
wget http://www.cfdbooks.com/cfdcodes/OSSAN_Euler2D/readme_v0.txt
wget http://www.cfdbooks.com/cfdcodes/OSSAN_Euler2D/twod_rectangular_grid_v0.f90

This compiles the software and executes it to produce the TecPlot output files:

source readme

The following is only necessary if you want to visualise the output using a Python (open source) stack, rather than using the commercial TecPlot product. This creates a Python virtual environment (“virtualenv”):

mkdir -p ~/science/venv/aerocfd
source ~/science/venv/aerocfd/bin/activate

This installs all of the necessary scientific packages for Python:

pip install numpy
pip install scipy
pip install matplotlib
pip install ipython
pip install pandas

I then wrote the following script to plot the output of the 2D Euler shock diffraction problem:

import matplotlib
matplotlib.use('TkAgg')
from matplotlib.mlab import griddata
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import numpy.ma as ma
from numpy.random import uniform
if __name__ == "__main__":
filename = "idolikecfd/2d_unstructured_euler_shock_diffraction/test_results.dat"
data = pd.io.parsers.read_csv(filename, delimiter="\t",
header=None, names=["x","y","rho","u","v","p","M"],
delim_whitespace=True)
x = np.array(data["x"])
y = np.array(data["y"])
z = np.array(data["rho"])
xi = np.linspace(0.0,1.0,500)
yi = np.linspace(0.0,1.0,500)
zi = griddata(x,y,z,xi,yi)
CS = plt.contour(xi,yi,zi,50,linewidths=0.5,colors='k')
plt.xlim(0,1)
plt.ylim(0,1)
plt.title('Density Plot')
plt.show()

This is the output from the Python plotter, which I think compares relatively favourably to TecPlot:

You can clearly see the expansion fan forming at the “corner” of the shock diffraction input tube (middle left on the image) and the associated circular shock wave forming into the region.

In the next set of articles I hope to try and construct a 1D relativistic Riemann solver problem, although I’ll have to brush up on my references first!

## References

- [1] – Nishikawa, H.; Kitamura, K. (2008), “Very simple, carbuncle-free, boundary-layer-resolving, rotated-hybrid Riemann solvers”, J. Comput. Phys. 227 (4): 2560–2581