# A Compressible Aerodynamics Library for Python

**Update 2021-02-18:** I have got around to writing proper documentation for
the `compflow`

library, which is available elsewhere on this
site. The below account is a more readable description of my
motivation for writing it and the development process.

## Problem statement

In modelling high-speed fluid flow, there are some frequently-used analytical relations which give a non-dimensional group as a function of Mach number, \(\mathit{Ma}\), and specific heat ratio, \(\gamma\). The relations are generally valid for any one-dimensional flow of a perfect gas. For example, the mass flow rate \(\dot{m}\), specific heat capacity at constant pressure \(c_p\), stagnation temperature \(T_0\), stagnation pressure \(p_0\), flow area \(A\) form the following non-dimensional group,

$$ \frac{\dot{m}\sqrt{c_p T_0}}{A p_0} = \frac{\gamma}{\sqrt{\gamma-1}}\mathit{Ma} \left(1+\frac{\gamma-1}{2}\mathit{Ma}^2\right) ^{-\frac{1}{2}\frac{\gamma+1}{\gamma-1}}\ . \tag{I} $$

Equation (I) is widely used for aerodynamic calculations involving compressible flows, particularly for internal flows such as in gas turbines. Evaluating the relation numerically given values for \(\mathit{Ma}\) and \(\gamma\) is straightforward, if a little fiddly to code and a good candidate for abstraction into a function. Equation (I) does not admit an analytical solution for \(\mathit{Ma}\), so the reverse problem, where values of \({\dot{m}\sqrt{c_p T_0}}/{A p_0}\) and \(\gamma\) are given, requires expensive iteration. This turned out to be a bottleneck in one of my programs.

## Implementation

I created the Python library `compflow`

to address these needs, after I was surprised to find that an implementation
did not already exist.

The first implementation was in pure `numpy`

, and only addressed the need for
abstraction — optimization of speed came afterwards. I took my undergraduate
engineering databook and coded up 9 of the compressible flow relations
tabulated therein. In Python, Equation (I) is simply:

```
def mcpTo_APo_from_Ma(Ma, ga):
return (
ga / np.sqrt(ga - 1.0) * Ma
* (1. + 0.5 * (ga - 1.0) * Ma ** 2.)
** (-0.5 * (ga + 1.0) / (ga - 1.0))
)
```

We can invert Equation (I) using the Newton solver from `scipy`

. The method
defaults to estimating the derivative by a finite difference approximation, but
it is faster to use an analytical form for the derivative if we have it. We
differentiate Equation (I), put it in another function `der_mcpTo_APo_from_Ma`

,
and feed that into the Newton solver:

```
from scipy.optimize import newton
def Ma_from_mcpTo_APo(mcpTo_APo, ga, supersonic=False):
# Function to find root of, and its derivative
def f(x):
return mcpTo_APo_from_Ma(x, ga) - mcpTo_APo
def fp(x):
return der_mcpTo_APo_from_Ma(x, ga)
# Choose inital guess for a sub or supersonic solution
if supersonic:
Ma_guess = 1.5* np.ones_like(mcpTo_APo)
else:
Ma_guess = 0.5* np.ones_like(mcpTo_APo)
# Do the iteration
return newton(f, Ma_guess , fprime=fp)
```

I also wrote some automated tests to check values returned by my functions against those written in the book, which picked up a few bugs! The tests also verify that the functions behave correctly in limiting cases, \(\mathit{Ma} \rightarrow 0\) and \(\mathit{Ma} \rightarrow \infty\), which is not guaranteed and depends on the how the algebraic form is evaluated numerically.

At this point, inverting Equation (I) using the Newton solver from `scipy`

was
still my bottleneck. Since the specific heat ratio is a constant for a
particular gas, it is feasible to cache a lookup table for a non-dimensional
quantity for each value of \(\gamma\). This would cost more up front, but
eventually win out if many calculations at the same \(\gamma\) are to be done, as
we don’t have to repeat the work involved in iterating. I implemented a simple
linear interpolation with adaptive sampling algorithm, where the spacing
between tabulated \(\mathit{Ma}\) values is bisected until the local error is
within a fixed tolerance.

Once a lookup table is created, querying it is about x50 as fast as iteratively
solving. Not bad, but I thought I could do better by dropping to a compiled
language for these calculations. `f2py`

, the Fortran to Python interface
generator seemed to be the simplest
option, compared to the reams of boilerplate code I would need to write for a C
extension to `numpy`

.

I like Fortran because fast code is easy to write. One can just write a do
loop and the compiler takes care of the rest. Having started programming with
`BASIC`

and proceeding on to `MATLAB`

, the old-school syntax comes naturally
for me. In Fortran, Equation (I) becomes:

```
SUBROUTINE MCPTO_APO(Y,M,G,N)
C G/SQRT(G-1)*M*(1+(G-1)/2*M^2)^(-(G+1)/(G-1)/2)
C ARGUMENTS
INTEGER N
REAL*8 G
REAL*8 M(N)
REAL*8 Y(N)
Cf2py intent(in) m
Cf2py intent(in) g
Cf2py intent(out) y
Cf2py intent(hide) n
C INTERMEDIATE VARS
REAL*8 M_GP1_GM1_2
REAL*8 GM1_2
REAL*8 G_SQ_GM1
GM1_2 = (G-1.0D0)/2.0D0
M_GP1_GM1_2 = (G+1.0D0)/(G-1.0D0)/(-2.0D0)
G_SQ_GM1 = G / SQRT(G-1.0D0)
C MAIN LOOP
DO I=1,N
Y(I)=G_SQ_GM1*M(I)*(1.0D0+GM1_2*M(I)*M(I))**M_GP1_GM1_2
ENDDO
END
```

The comment lines starting `Cf2py`

are special directives to instruct `f2py`

to
produce the most efficient wrapper. The only optimisation I made is to move the
computation of some constants such as \((\gamma - 1)/2\) outside the loop,
because their values are the same for each \(\mathit{Ma}\) in the input
vector. The Fortran routines for inversion are lengthier but conceptually
simple and straightforward to code, just basic Newton-Rhapson iterations.

## Benchmarks

To quantify the speed improvements, I coded up a simple benchmark that compares
the execution time of native `numpy`

, lookup table, and Fortran
implementations. The speedup depends on the size of the input \(\mathit{Ma}\)
array, \(n\).

For forward calculations, i.e. evaluation of Equation (I) for a given
\(\mathit{Ma}\), the Fortran implementation yields a x18 speedup for scalar \(n=1\)
inputs. Although each of the builtin `numpy`

functions are well-optimised C,
which should be hard to beat individually, doing the entire calculation in one
compiled go has clear benefits. For arrays with \(n\ge1000\) elements, the
implementations are equal in speed. Where the overhead of going in and out of
compiled functions is relatively smaller, and the cost is dominated by the
calculations themselves, our advantage reduces.

For inverse calculations, solving Equation (I) for \(\mathit{Ma}\) given
\({\dot{m}\sqrt{c_p T_0}}/{A p_0}\), the lookup table offers a speedup of x45 for
scalar inputs, while the compiled version has an enormous speedup of x400!
This is a great result. While the `scipy`

Newton solver is very general and
powerful, a simple, bespoke compiled version can be much faster. However, for
input arrays with \(n\ge20\), querying the lookup table wins, probably because
`scipy`

is clever enough to reuse some work when evaluating multiple table
entries at once.

## Conclusion

I am fairly happy with the library as it stands, and any future work will be directed at improving the documentation and encouraging other people to use it. The library is available on sourcehut or PyPI for the reader to install. Let me know if you have found it useful!