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!