II.3 Matrix inversion methods

The solution step in CFD basically involves evaluating M-1 in order to advance the solution by one timestep. In theory, any valid technique could be used to find M-1. In practice there are constraints imposed by the practicalities of the computation. For instance, you may know about Gaussian elimination or Cramers rule for exactly inverting small (eg 3 × 3) matrices by hand. If we were to use these methods to invert matrices in a CFD code, we would quickly find problems. First, these methods require us to store all N × N components of the matrix, which is completely impractical - N is the number of cells in the mesh, which can be as large as 10 million. Second, these methods scale as N3. This means that if we double the number of cells (double N), the time taken to execute the algorithm increases by a factor of 23 = 8. For these reasons we use approximate, iterative algorithms for inverting matrices. The matrices we are trying to invert are sparse (ie. mostly contain zeros), and so only non-zero entries need to be stored. The methods used are often iterative, so although they are approximate, the iteration can be repeated until any desired accuracy has been achieved.

A good matrix algorithm needs to be

For example, a 1-d FV method links each cell to its 2 neighbours. This produces a matrix similar to that shown in figure 6.


PIC
Figure 6: Tri-diagonal matrix produced from a simple 1-d FV matrix.

Only the elements D, A and B need to be stored in memory (a considerable saving). Each row of the matrix generates a relation

Afi- 1 + Dfi + Bfi+1 =  Qi
This can be combined with the equation for the preceeding row to provide a forward elimination equation for fi given fi+1. The next 2 rows provide an equation for fi+1 given fi+2 etc. Having produced these relationships, we start at the Nth gridpoint (the forward elimination equation for this gridpoint involves the boundary condition, which we know) and back substitute to find the values at all the points.

This procedure forms the Tri-Diagonal Matrix or TDMA algorithm - probably the simplest matrix inversion algorithm available. It can be extended to 2 and 3d case. Other more sophisticated methods are also used - the exact algorithm used often depends on the exact type of equations being solved, the structure of the mesh, and other application specific details.