<previous: MUSCL Differencing	up to Table of Contents	next: References >
open User's Guide (in this window) open Applet Page (in new window)

Matrix Solution Techniques

    The one remaining item to consider is that the implicit time integration algorithms result in a set of linear matrix equations of the form:

(75)

The implementation of the time and spatial construction algorithms results in a system of equations in which both the left and right hand side matrices are completely known and the vector of conservative variables, q, remains as the only unknown.  In general 2D and 3D CFD calculations, the matrix L would be a sparse matrix of high bandwidth which would almost certainly require iterative methods such as the Gauss-Seidel method or the Conjugate Gradient method in order to solve.  Fortunately, in 1-D the matrix L is completely composed of zero elements except for the terms along the matrix diagonal and elements adjacent to the diagonal.  This form is a special form which allows the use of some simpler methods of solution.  A variation of the well-known Thomas algorithm may be employed to generate a solution with is both more efficient and exact when compared to iterative methods.
    For the situation here, the 1-D finite folume formulation generates a system of equations which consists of N cells with 3 equations and unknowns per cell.  This results in 3N equations and unknowns which is best treated as a system of block matrices.  Thus, there are N block equations (with each "block" consisting of 3 unknowns).  The dependency of the implicit algorithms generate a Jacobian matrix for the diagonal block and the two blocks immediately on either side.  This results in a pentadiagonal matrix system which has the form shown in eqn. (76).  Due to the large size of the matrix system, the scaled down image presented can be selected for a link to a larger, more readable version.

(click image to enlarge scale - opens in a new window)

(76)

where each of the implicit block matrices are given in the form (a time step coefficient actually resides in front of the matrix, and the diagonal entry has an extra term of +1 in it -- but these are small details):

(77)

    The traditional Thomas algorithm works for a scalar set of equations in tri-diagonal form.  Here, the equation set is in block form and is pentadiagonal, but the method may be extended to include this case as well without significant effort.  Since the L matrix is comprised almost entirely of zero blocks, it is prudent to only store the necessary information. Therefore, in Gryphon, the L matrix is a matrix of N rows by 5 columns, with each element itself being a 3 x 3 matrix.  This requires storage space of 45N double precision numbers, which is an incredible savings over storing the entire matrix.  It is understood that the columns represent the second left, left, diagonal, right, and second right blocks at each row.  It is hence more efficient to simplify the notation a bit to make things more clear.  In each row, the C matrix is always the diagonal entry, with A and B representing the two left-most components and D and E representing the two rightmost components.  The right hand side term M remains as well.

(78)

    The modified Thomas algorithm proceeds according to the following algorithm:

• FOR rows 3 to N, eliminate the second left block A in each row by subtracting the previous row times the target A matrix times the inverse of the previous row's B matrix.  This means carrying out the following operation on both the left and right hand sides of the system:
  

  (79)

  
  
• FOR rows 2 to N, eliminate the left block B (which was modified in the previous step) in each row by subtracting the previous row times the target B matrix times the inverse of the previous row's C matrix.  This means carrying out the following operation on both the left and right hand sides of the system:
  

  (80)

  
  
• The last row is completely diagonalized a this point.  Therefore, the final block may be solved by taking the inverse of the modified C matrix time the right hand side as shown:
  

  (81)

  
  
• The second to last row, can next be solved for since the solution for cell N is known.  It is calculated as shown:
  

  (82)

  
  
• FOR rows N-2 to 1, working backwards, each conservative variable vector j can be found because the j+1 and j+2 vectors are already known, leading to the following general formula for backward substitution: 
  

  (83)

  
  

    It is important to realize that the above procedure is a set of matrix operations, not scalar operations.  The inverse operations are full matrix inversions of a 3 by 3 matrix.  In general numerical mathematics, iterative methods are typically used to determine the inverse of a large matrix due to the high cost of calculating an inverse exactly.  Fortunately, the exact procedure is useful here since the matrix operations involve only inverting a 3 x 3 matrix (and no larger).  The inverse of any invertable 3 x 3 matrix is given in eqn. (84).  This formula must be used every time a -1 exponent appears in the above algorithm.

(click image to enlarge scale - opens in a new window)

(84)

    The Thomas algorithm amounts to a very specialized form of Gaussian elimination, which only works for diagonally structured matrices of this type.

<previous: MUSCL Differencing

up to Table of Contents

next: References >