# Applying Conjugate Gradient Method to Fluid Simulations

This topic is 3850 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

## Recommended Posts

Hi, I'm trying to write a solver for a grid based fluid system using the conjugate gradient method. I've read about it on wikpedia (http://en.wikipedia.org/wiki/Conjugate_gradient) and then read a better article here (http://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf) Now I'm at a point where I want to apply the maths to code, but I'm stuck at the point of mapping my code and variables to the maths. Given Ax = b, What's A? My thinking is that it's the density array that holds all the forces acting on the fluid grid, but arrays and matrices are different aren't they? I'm guessing x is the current position of the element in question as a vector. The problem then is that I don't see how you can then multiply a two element vector (asuming a 2d fluid) against an NxN grid (where N is the number of grid elements). I've found general mathematical implementations of the CG method, but these use things like a general Matrix class, and I'm having problems seeing how to apply these to the problem at hand.

##### Share on other sites
Quote:
 Original post by WinegumsHi, I'm trying to write a solver for a grid based fluid system using the conjugate gradient method.

How are you trying to use the CG method? It's not normal to use it for either fluids (i.e. complete 3D internal motion), or fluid surfaces, as far as I know.

Quote:
 Given Ax = b,

A is a square (NxN) matrix that you know in advance. It's typically sparse - i.e. is mostly empty, and has some other properties

x is an unknown vector of length N.

b is a known vector of length N

Maybe if you explain a bit more about what you're trying to do then someone can give more advice, since it's not clear why/how you're trying to use the CG method, and it might not be appropriate anyway. What kind of fluid do you want to simulate?

##### Share on other sites

I'm using CG as per a recommended extension for the fluid system proposed in:

http://www.dgp.toronto.edu/people/stam/reality/Research/pdf/GDC03.pdf

The paper says this should be reasonably straightforward to implement, but without knowing how to match up the data in the code to the variables in the equation, I'm sort of lost.

##### Share on other sites
It says to use it in place of the gauss-siedel solver. There is one row and one column in the matrix for every cell in the grid. Most of the entries in the matrix are 0. There are only 5 non-zero entries in each row.

x is actually called x in the paper. This equation from the paper:
x0[IX(i,j)] = x[IX(i,j)] - a*(x[IX(i-1,j)]+x[IX(i+1,j)]+x[IX(i,j-1)]+x[IX(i,j+1)]-4*x[IX(i,j)]);
describes the matrix. x0 is b.

##### Share on other sites
Quote:
 Original post by VorpyIt says to use it in place of the gauss-siedel solver. There is one row and one column in the matrix for every cell in the grid. Most of the entries in the matrix are 0. There are only 5 non-zero entries in each row.

How is it you've concluded that the matrix only has 5 non zero entities in each row?

Quote:
 x is actually called x in the paper. This equation from the paper:x0[IX(i,j)] = x[IX(i,j)] - a*(x[IX(i-1,j)]+x[IX(i+1,j)]+x[IX(i,j-1)]+x[IX(i,j+1)]-4*x[IX(i,j)]);describes the matrix. x0 is b.

Is that equation explicitly from the paper? It doesn't seem to be in it, or have you rearranged the one in the diffusion function?

either way, I'm a bit confused. Are you saying that
Ax = b is the same as
x0[IX(i,j)] = x[IX(i,j)] - a*(x[IX(i-1,j)]+x[IX(i+1,j)]+x[IX(i,j-1)]+x[IX(i,j+1)]-4*x[IX(i,j)]);

with x0 being b?

and the x referred to in the functions is the x used in the equation? this part really confuses me, as I thought x was a vector?

Assuming x is x, should I then use that with b to find A?

##### Share on other sites
The equation is from the diffusion equation given in the paper; I just copied and pasted it. It is given in code form, but it is also given as part of the description of the algorithm outside of the code segments. Each row of the matrix is described by the equation. Each row in the matrix has a different corresponding (i,j) coordinate. The elements in the equation correspond to non-zero elements of the matrix.

x0 is b.
x is x.
The matrix is a representation of the equations relating the elements of x0 and x. (one equation for every cell in the fluid simulation)

x0[IX(i,j)] = x[IX(i,j)] - a*(x[IX(i-1,j)]+x[IX(i+1,j)]+x[IX(i,j-1)]+x[IX(i,j+1)]-4*x[IX(i,j)]);

translated to look more like a row of a matrix equation:

b = x - a * (x[i-1] + x[i+1] + x[i-width] + x[i+width] - 4*x)

and then distributing out the constants to make it even more clear that the right hand side is the result of multiplying a row of constants in A by the vector of unknowns x:

b = (1+4*a)*x - a*x[i-1] - a*x[i+1] - a*x[i-width] - a*x[i+width]

For example, in a 3 by 3 matrix, with elements 0 to 8, the equation for the center element would be:

b[4] = (1+4*a)*x[4] - a*x[3] - a*x[5] - a*x[1] - a*x[7]

As a dot product, this would look like:

b[4] = (0, -a, 0, -a, 1+4*a, -a, 0, -a, 0) * x

The vector (0, -a, 0, -a, 1+4*a, -a, 0, -a, 0) is one row of A, for this 3 by 3 cell simulation.

The equations for boundary cells may need to be slightly different to address the boundary conditions.

every x[n] not in the equation corresponds to a 0 in that row of the matrix

[Edited by - Vorpy on February 15, 2008 12:52:34 AM]

##### Share on other sites
Quote:
Original post by Winegums
Quote:
 Original post by VorpyIt says to use it in place of the gauss-siedel solver. There is one row and one column in the matrix for every cell in the grid. Most of the entries in the matrix are 0. There are only 5 non-zero entries in each row.

How is it you've concluded that the matrix only has 5 non zero entities in each row?

That arises directly from the choice of discretizations. For a 2D flow, you'd get 5 nonzero entries based on a finite difference/finite volume discretization of the governing equations in which the derivatives at a given cell are computed via 2 point differences in either direction....2 nonzero values in x for x derivatives, 2 nonzero values in y for y derivatives, plus 1 nonzero value for the cell itself = 5 total. This collection of nodes, the ones that have nonzero coefficients in each row of the matrix, is sometimes called the "computational molecule," or "stencil."

How to Understand CFD Jargon

Practical CFD Simulations on Programmable Graphics
Hardware using SMAC
(This paper has some excellent graphics illustrating some more complex computational molecules.)

##### Share on other sites
Ok thanks for your replies guys. I read what you wrote a few times over, and I may be getting it, albeit slowly.

I still need to find A individually don't I?
Does this mean A will be a matrix of dimensions N^2? (with N being the width/height of the system)
Also I'm a bit unsure about which elements in the matrix will be non-zero, will it be ones centred around the diagonal?

And then, assuming I have A, how do you multiply a vector (which only has an X and Y component) by a matrix? That is asuming that b is meant to just be a vector?

##### Share on other sites
If the fluid simulation is N cells by N cells, then the matrix will have N^2 rows and N^2 columns. Since there are only 5 nonzero entries in each row, the matrix is very sparse, and you'll want to use a sparse matrix algorithm (an algorithm optimized for matrices which contain mostly zeros). A sparse matrix representation usually uses linked lists or something similiar to avoid having to store all the zero entires, and to quickly locate the nonzero entries.

The matrix in this problem has a very regular structure, and the non-zero entries are centered around the main diagonal, although they are not necessarily close to the main diagonal. I'm not sure what sparse conjugate gradient algorithms look like, but it may be possible to optimize such a method for this particular problem because of the regular structure of the matrix.

You probably don't want to use a general purpose implementation that is designed for a dense matrix.

Multiplying a column vector by a matrix, as in Ax, is equivalent to creating a new column vector which contains the dot products of the rows of the matrix and the column vector.

##### Share on other sites
Actually, you don't need to construct a matrix
of any kind. If you're using CG to enforce
mass conservation, then "multiplication of A and x"
simply means computing the laplacian of x.

- Mikko

1. 1
2. 2
Rutin
25
3. 3
4. 4
5. 5

• 11
• 10
• 13
• 20
• 14
• ### Forum Statistics

• Total Topics
632950
• Total Posts
3009377
• ### Who's Online (See full list)

There are no registered users currently online

×