Convective Flux Jacobians

The Euler equations are an important system of equations in computational fluid dynamics (CFD). They are basically the full continuum equations of motion without the viscosity. Adding in the viscosity gives you the more realistic Navier-Stokes equations.

The conserved variable vector is

And the flux vector is

The flux Jacobian is simply the partial derivative of the flux with respect to the conserved variables (mass, momentum, energy).

The maxima to calculate the Jacobian is straightforward:

U : matrix([u1], [u2], [u3]);

f1 : u2;
f2 : u2^2 / u1 + p_exp;
f3 : (u2 / u1) * (u3 + p_exp);
F : matrix([f1], [f2], [f3]);

A : zeromatrix(nd, nd);

for j : 1 thru nd do
(
for i : 1 thru nd do
(
A[i,j] : diff(F[i,1], U[j,1])
)
);

Diagonalization

Using the chain rule lets you write the Euler equations in a quasi-linear form in terms of the flux Jacobian, now called A.

The flux Jacobian can be diagonalized into a matrix of right eigenvectors, diagonal matrix of eigenvalues, and a matrix of left eigenvectors.

The code to accomplish this in maxima is a little more involved, because of some substitutions needed to get everything in terms of u,gamma,a.

load("eigen");
VLV : eigenvectors(simp_A);

a_exp : (g * (g - 1) / u1) * (u3 - u2^2 / (2 * u1));
a_exp : sqrt(a_exp);
simp_a : ratsimp(ev(a_exp, u1 = rho, u2 = rho * u, u3 = rho * E));
a_eqn : a^2 = simp_a^2;
c_exp : g*(2*E-u^2);
a_eqn : facsum(a_eqn, c_exp);

/* put the eigenvectors into a matrix: */
V : columnvector( VLV[2] );
V : addcol( V,VLV[3],VLV[4] );
V : ratsimp( V );

/* put the eigenvalues into a diagonal matrix: */
load("diag");
L : diag( VLV[1][1] );

radsubstflag : true;
/* substitute in the speed of sound: */
V : ratsimp(ratsubst(a, simp_a, ratsimp(V)));
V : factor(ratsimp(ratsubst(a, simp_a, rootscontract(V))));
/* factor out gammas */
V : facsum(V,g);
c_a : part(solve(ratsubst(c,c_exp,a_eqn),c)[1],2);
V : ratsubst(c,c_exp,V);
V : ratsubst(ev(c_a),c,V);
V : facsum(V,u);

/* substitute in the speed of sound: */
L : ratsimp(ratsubst(a, simp_a, ratsimp(L)));

/* left eigenvectors: */
Vinv : factor( ratsimp( invert( V ) ) );
/* factor out gammas */
Vinv : facsum(Vinv,g);
d_exp : (g-1)/a^2;
Vinv : ratsubst(d,d_exp,Vinv);

Conclusion

The thing to note about this result is that there are three distinct wave speeds in the one dimensional Euler equations, u-a, u+a, and u, where a is the speed of sound. There are lots of possible ways to write the eigenvectors, in terms of the enthalpy, or the energy, or the pressure, but the form above with only the flow velocity, the speed of sound, and the ratio of specific heats is pretty elegant. It really illustrates what a fundamental quantity the speed of sound is for this set of equations.

Chapter 3 of Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction has a good work-up of this development and shows some alternative formulations for the flux Jacobian.

Kronecker Product of Sparse Matrices

Gilbert Strang discusses the use of a neat built-in function in Octave, kron(), that's pretty useful for doing multidimensional finite difference methods on PDEs:


kron() takes the Kroenecker product of two matrices. This is especially useful when constructing the large, sparse matrices needed for finite difference approximations.

K is the tridiagonal matrix created for a central difference approximation to the second derivative:

n = 5;
K = spdiag(2*ones(n,1),0) + spdiag(-ones(n-1,1),-1) + spdiag(-ones(n-1,1),1);
spy(K)

K is the differences for a 1D row of points, now we can create the 2-dimensional version for a plane of points with kron():

K2D = kron(speye(n),K) + kron(K,speye(n));
spy(K2D);

The extension to 3 dimensions is exactly as before, now we have a matrix describing the differences between points on the same line in the first sub and superdiagonal, differences across lines in the next bands out from those, and differences across planes in the furthest out bands:

K3D = kron(speye(n),K2D) + kron(K2D,speye(n));
spy(K3D);

Pretty simple! Now you're ready to solve multidimensional partial differential equations.

Fibonacci Numbers

Gilbert Strang has an interesting example in one of his linear algebra lectures on eigenvalues that is relevant to the articles I've written on speeding up calculations in Octave.

You might be familiar with the Fibonacci numbers, they are calculated by a two term recurrence:


We've already seen that recasting a calculation in terms of a matrix vector multiply or matrix vector solve is an important skill for writing good fast code in Octave. Strang demonstrates this at about halfway through the above lecture. The first order system that describes the generation of Fibonacci numbers is:


The nth term of the sequence is calculated by taking the nth power of A and applying it to the initial conditions:


When you diagonalize a matrix the powers of the matrix become very easy to calculate:

If you aren't familiar with the concept of eigenvalues and eigenvectors and diagonalizing a matrix I would highly recommend the linear algebra lectures on MIT's OCW website.

Here's the Octave code to accomplish the above. If you want, you can use tic() and toc() for timing. On my laptop the matrix method is about three times faster than the for loop method. What kind of results do you get? This advantage will grow as you go further into the sequence because the amount of work stays constant for the matrix method, but grows linearly for the recurrence method.

%% Our 'fibonacci' operator:
A = [1,1;1,0]
%% our initial conditions:
u_0 = [0;1];

%% get the matrix of eignevectors, V, and the matrix of eigenvalues, lambda
[V,lambda] = eig(A)

%% A = V lambda V^{-1}

Vinv = inverse(V)

%% This gives us an idea of the error introduced by calculating the
%% inverse with finite precision arithmetic:
Error = A - V*lambda*Vinv

n = 1400;
F = zeros(n,2);
%% Find the 100th Fibonaci number:
F(100,:) = V * (lambda^100) * Vinv * u_0;

%% check our answer with a simple loop:
fibo = zeros(n,1);
fibo(1) = 0;
fibo(2) = 1;
for i=3:n
fibo(i) = fibo(i-1) + fibo(i-2);
endfor

for i=2:n
F(i,:) = V * (lambda^i) * Vinv *u_0;
endfor

relative_error = abs( ( F(3:n,2) - fibo(3:n) ) ./ fibo(3:n) );

One thing to note is that one of the eigenvalues of our matrix is greater than one. This has implications for the numerical stability of the above method for calculating terms in the Fibonacci sequence. It means that the initial tiny round-off error in our eigenvector matrix inverse will be amplified as we continue to calculate higher and higher terms in the sequence. This behaviour is shown in the plot below, luckily the error starts out very small, and doesn't grow too large before we run out of range in the 32bit floating point numbers on the machine.

What's the answer? The 100th Fibonacci number is approximately 2.1892e20.

How to Speed Up Octave

Speed seems like a continuing theme on the Octave mailing lists, a recent thread on the list of "comparing the executive speed with Matlab" demonstrates that there really is very little speed overhead when you use Octave smartly. If you read down far enough, it also demonstrates that the Intel Fortran compiler is awesome!

There are four basic methods of speeding up a chunk of Octave code (in order of increasing speed):

1. Allocate memory smartly (using for loops)


2. Vectorize your code (avoiding for loops)
   


3. Use sparse matrix operators in place of for loops you can't vectorize


4. Write a compiled function in C/C++/Fortran to call from Octave

The first three options keep you in pure Octave. The fourth option, while requiring you to use a compiled language is really not too bad because you can use all of the types (vectors, matrices, etc) defined in the Octave libraries to write your function. That's almost as easy as doing it in Octave.

Speeding Up For Loops

This one really amounts to not being foolish. The only way to speed up for loops is to make sure you allocate all of your memory (using zeros() for instance) before you enter the for loop. The Octave Manual covers this topic here. The better thing to do is vectorize your code and avoid for loops all together.

Vectorizing

What you really want to do is cast your calculation as a series of vector and matrix operations, or use the built-in (usually compiled) functions that operate directly on vectors. The post Is Octave Slow? demonstrates the sorts of speed up you can get when you get rid of for loops.

Sparse Matrices for Data Dependancy

What about for loops you can't vectorize? The ones where the result at iteration n depends on what you calculated in iteration n-1, n-2, etc. You can define a sparse matrix that represents your for loop and then do a sparse matrix solve. This is much faster because you are now using compiled functions. The drawback is that you are trading some increase (sometimes substantial) in memory usage to buy this speed up.

• Is Octave Slow? covers this approach as well.


• More on Speeding up Octave covers the details on sparse matrix allocation.

Compile It

The manual has a good intro on using mkoctfile to compile your function and make it callable from Octave. Probably the most useful thing is that you can use the Octave types, so you don't have to code up your own matrix or vector type. At this point the speed really has little to do with Octave and everything to do with the quality of the compiler you use and how efficiently you code up your calculation.

Conclusion

The real power of Octave isn't that you can get your calculations to be nearly as fast as if you wrote them all in a compiled language. It's that you can prototype a calculation quickly and make sure that it is correct, and then there is a fairly well defined path to follow to speed up your correct calculation.

Thresholding for Edge Detection

The real world is noisy, and the differencing operations of edge detection tend to magnify that noise. We can apply some statistical methods to our edge detection to hopefully overcome that. A very basic approach is to threshold based on a histogram of the gradient magnitude. Taking a look at the distribution of gradient magnitude for our trusty old electricity meter (shown below), it looks like a pixel level below about 235 puts us safely out into the tail of the distribution. Most pixels (in this image) are not edge pixels, so they are clustered on the high end of the scale.

This is pretty easy to do with numpy and Image (PIL):

import numpy
import Image
pix = numpy.array( img.getdata() )
for i in xrange(0,pix.size):
if pix[i] < 235: pix[i] = 0 else: pix[i] = 255 img_out.putdata( pix )

Of course the histogram shape will depend on how much noise and how many edges are in the particular picture. This one had black dials on a white background, and not much noise so thresholding is a pretty fruitful strategy for really pulling out the edges.

You can make an intuitive argument that gradient magnitude distributions will have roughly similar shapes to that shown above. Especially as image resolutions steadily increase, the vast majority of pixels will not be on an edge. Sort of the curse of dimensionality in reverse, I guess that makes it a blessing.

Incompatibility between Matplotlib and Numpy

If you are running Fedora 9 you may have run in to this problem (at least one other guy did). Taking a look at the bugzilla over at RedHat shows that this has been a recurring problem over the past couple years (first cropping up in FC5). It seems that the packagers for numpy and matplotlib are having trouble synchronizing the updates for numpy and matplotlib. They are fast moving projects so one occasionally changes the API quietly and breaks the other. Of course there are other dependencies that probably play in to the decision of when to update a particular package, but this is the one that breaks my workflow.

The easiest fix for me was to download the latest matplotlib sources and compile/install them manually:

$ tar -zxvf matplotlib-0.98.5.2.tar.gz
$ cd matplotlib-0.98.5.2
$ python setup.py build
$ su
$ python setup.py install

This is less than ideal, because now I'm not running a "stock" Fedora where all of my software updates are taken care of automagically, but the increase in work-load is negligible.

A work-around for the incompatibility between numpy and matplotlib that requires editing some of the python sources:
http://www.mail-archive.com/numpy-discussion@scipy.org/msg13753.html

Or you can just upgrade to Fedora 10. Unfortunately that's not an option for me because of a closed proprietary device driver that relies on < 2.6.27 kernels. I haven't compiled my own kernels since FC1 days, and don't really feel like getting back into the habit.

A Treatise on the Calculus of Finite Differences

Boole, George. "A Treatise on the Calculus of Finite Differences", Macmillan and Co., 1872.

The public domain and Google books is a great combination. You can download the classic by George Boole, A Treatise on the Calculus of Finite Differences (TCFD) and read it at your leisure, for free!

Even though it was written in 1860, this book probably even more relevant today than it was then. Why? Because it provides the mathematical basis for solving higher and higher resolution partial differential equations with faster and faster computers for the complex domains of modern engineering problems. The faster, more accurate and more efficient solution to these types of problems provides the foundation for improving quality of life for everyone.

I like his definition of the calculus of finite differences, it says so much, succinctly:

The Calculus of Finite Differences may be strictly defined as the science which is occupied about the ratios of the simultaneous increments of quantities mutually dependant.

This book has to be understood as a second volume to accompany Boole's first book on Differential Equations. He stresses the importance of the relationship between the two in his preface to TCFD. The discussion of the distinction between a calculus of limits and the calculus of differences highlights this connection.

Like most old math or science books, it is the discussion of important applications that feels dated. The ideas are timeless, but the important applications change with the problems and the technology of the time.

Another timeless idea that bears repeating around all practitioners of the numerical art (from Chapter 8):

... we shall very often need to use the method of Finite Differences for the purpose of shortening numerical calculation, and here the mere knowledge that the series obtained are convergent will not suffice; we must also know the degree of approximation.

To render our results trustworthy and useful we must find the limits of the error produced by taking a given number of terms of the expansion instead of calculating the exact value of the function that gave rise thereto.

Treating Image Edges

There are some good examples floating around of creating customized image processing kernels in the python imaging library (PIL). It's very convenient to simply define your own mask, and the library takes care of the convolution.

One drawback is that the edges don't get filtered in this process. There are various ad-hoc things you can do. You can just drop a line of pixels all the way around the picture, so the gradient image of a 800x600 picture would be 798x598 pixels (if you're using a compact stencil), or you could just copy the edge pixels from the old picture to the new picture. So in effect, the edges don't get filtered. For most applications this is probably no big deal, but it might be, and we can do better.

All we have to do is change our difference operator for the edge pixels, use a forward or backward difference on the edges rather than the central differences used on the interior. Then you get gradient info for your edge pixels too!

Customizing the PIL kernel class has been done, so I'm going to use Fortran90 and F2Py for my difference operators. I'll still use the PIL for all of the convenient image reading/writing utilities though. No need to implement that in Fortran!

We can use our little Maxima script to generate the differences we need on the edges. If you've taken a basic numerical methods course you'll probably recognize them as the three point forward and backward difference approximations for the first derivative.



The Python to do just that is here:

import sys, math, Image, numpy
from filters import *


if __name__ == '__main__':
in_img = Image.open(sys.argv[1])
in_img = in_img.convert("L")
out_img = Image.new("L", in_img.size)
pix = in_img.getdata()
# put the pixels into a fortran contiguous array:
pix_array = numpy.array(pix)
pix_array = pix_array.reshape( in_img.size, order='F' )
# calculate the gradient
g = first_order_gradient(pix_array,in_img.size[0],in_img.size[1])
gmag = 255 - numpy.sqrt( g[0,:,:]*g[0,:,:] + g[1,:,:]*g[1,:,:] )
gmag.astype("int")
# pack the gradient magnitude into a new image
out_img.putdata( gmag.reshape( gmag.size, order='F' ) )
out_img.save(sys.argv[2])

Where I've used F2Py to generate the filters module from a Fortran90 file. The module defines the subroutine first_order_gradient(), which implements the ideas discussed above. The gradients on the interior points are calculated using a central difference:

! apply the central operator to the interior pixels
do j = 2, nj-1
do i = 2, ni-1
g(1,i,j) = (in_image(i+1,j) - in_image(i-1,j)) / 2.0
g(2,i,j) = (in_image(i,j+1) - in_image(i,j-1)) / 2.0
end do
end do

The gradients along the edges are calculated using the forward and backward differences (only x edges shown for brevity):

! loop over the y dimension and apply the asymmetric operators to the
! x min and x max edges
do j = 2, nj-1 ! leave out the corners
i = 1 ! minimum edge, go forward
g(1,i,j) = (-3.0*in_image(i,j) + 4.0*in_image(i+1,j) - in_image(i+2,j))/2.0
i = ni ! maximum edge, go backward
g(1,i,j) = (3.0*in_image(i,j) - 4.0*in_image(i-1,j) + in_image(i-2,j))/2.0
! no change for the y derivative
i = 1
g(2,i,j) = (in_image(i,j+1) - in_image(i,j-1)) / 2.0
i = ni
g(2,i,j) = (in_image(i,j+1) - in_image(i,j-1)) / 2.0
end do

Finally the four corners are calculated separately, because each one has a unique combination of stencils required. Putting the loops in Fortran speeds things up considerably compared to the pure Python implementations I did before. Here is the gradient magnitude of the dials from Poor Man's ADC calculated with the above technique:

Sobel, Scharr and Truncation Error

In image processing the discrete difference operators to calculate gradients and second derivatives are called "masks". They are often referred to as "stencils" in other applications of finite differences for solving PDEs, estimating sensitivities, etc.

Some of the kernels (a mask plus particular values) for calculating gradients are Sobel:

[ 1 2 1]
Gy = [ 0 0 0]
[-1 -2 -1]

[ 1 0 -1]
Gx = [ 2 0 -2]
[ 1 0 -1]

or Sharr:

[ 3 10 3]
Gy = [ 0 0 0]
[-3 -10 -3]

[ 3 0 -3]
Gx = [ 10 0 -10]
[ 3 0 -3]

or the simpler stencil we used in the previous posts, which only used neighbouring pixels in one direction to estimate the gradient. What if we want to design our own masks? A useful tool is the multidimensional Taylor series:

In two dimensions we have:

which includes up to third order terms. We can't actually include the two pure third order terms because we only have three points in either the x or y direction, and we need a minimum of four. We can include the mixed third order terms though, so that gives us 8 constraints on our coefficients. Instead of including all 9 points in the stencil we have to drop one, so in the Maxima code below you'll notice that the central point is commented out.

Lets examine the truncation error for the two masks above. First we define a vector of delta x and delta y corresponding to each point in the stencil:

/* define the points in the stencil or mask: */
steps : matrix( /* dx, dy */
[-1, 1],
[0, 1],
[1, 1],
[-1, 0],
/* [0,0], */
[1, 0],
[-1, -1],
[0,-1],
[1, -1] );

Then vectors describing the coefficients in the two masks (we'll just use the y-component for this example):

/* coefficients in sobel mask */
sobel_y : matrix(
[1],
[2],
[1],
[0],
/*[0],*/
[0],
[-1],
[-2],
[-1]
);
/* coefficients in sharr mask */
sharr_y : matrix(
[3],
[10],
[3],
[0],
/*[0],*/
[0],
[-3],
[-10],
[-3]
);

To evaluate the truncation error we also need to create an operator that is based on the multidimensional Taylor series shown above. This will let us simply multiply those vectors by the operator to see what the masks are actually estimating. Here's the operator:

n : matrix_size( steps );
A : zeromatrix(n[1],n[1]);

for j : 1 thru n[1] do ( A[1,j] : 1 );/* f */
for j : 1 thru n[1] do ( A[2,j] : steps[j,1] );/* df/dx */
for j : 1 thru n[1] do ( A[3,j] : steps[j,2] );/* df/dy */
for j : 1 thru n[1] do ( A[4,j] : steps[j,1]^2/2! );/* d2f/dx2 */
for j : 1 thru n[1] do ( A[5,j] : steps[j,2]^2/2! );/* d2f/dy2 */
for j : 1 thru n[1] do ( A[6,j] : 2*steps[j,1]*steps[j,2]/2! );/* d2f/dxdy */
for j : 1 thru n[1] do ( A[7,j] : 3*steps[j,1]*steps[j,2]^2/3! );/* d3f/dxdy2 */
for j : 1 thru n[1] do ( A[8,j] : 3*steps[j,1]^2*steps[j,2]/3! );/* d3f/dx2dy */

Now the truncation error terms are shown by a couple of matrix vector multiplies:

sobel_te : A.sobel_y;
sharr_te : A.sharr_y;

Rather than just getting a single non-zero entry in the third place of the resulting vector we get a non-zero entry in the last place as well. Looking back at our operator, this is the row for the third order mixed derivative term, (d^3 f)/(dx^2 dy). In other words both the Sobel and Sharr operators introduce an error in their estimation of the gradient that is proportional to this mixed derivative.

How to avoid this error? Well, we can choose a right hand side vector that only has a 1 in the third position and zeros everywhere else, and then invert our operator to solve for what the values in our mask should be.

/* define which derivative we want to estimate, constrain the unwanted
ones to be zero: */
derivs : matrix(
[0], /* f */
[0], /* df/dx */
[1], /* df/dy */
[0], /* d2f/dx2 */
[0], /* d2f/dy2 */
[0], /* d2f/dxdy */
[0], /* d3f/dx2dy */
[0] /* d3f/dxdy2 */
);
inv_A : invert(A);
coeffs : inv_A.derivs;

This gives a result with only non-zero coefficients for the points at [0,1] and [0,-1], in other words, the simple estimate we've already been using! There's really no need to include all those other points for the gradient calculation.