# Applying Conjugate Gradient Method to Fluid Simulations

This topic is 3543 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[i] = x[i] - a * (x[i-1] + x[i+1] + x[i-width] + x[i+width] - 4*x[i])

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[i] = (1+4*a)*x[i] - 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

##### Share on other sites
Quote:
 Original post by VorpyMultiplying 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.

But isn't x just a 2 or 3 dimensional point? and A is a huge matrix, so the equation would be irrelevant after the first 2 or three points?...

1 0 0 0 0
0 1 0 0 0
0 0 1 0 0 . [3,1]
0 0 0 1 0
0 0 0 0 1

Even if I changed [3,1] to be [3,1,0,0,0] it would still be pointless calculating anything beyond the first two rows.

EDIT: Or I could be an idiot. x is of length N, so is x a 2xN matrix? With N rows of 2 points?

Also if the non zero elements are always close to the diagonal, what happens at early/late points when there isn't space to the left or right? do they just get pushed to the other edge of the row?

uttumuttu:

Thanks for your post, though I'm having some trouble understanding the laplacian stuff (relying on Wikipedia here!). I looked at discrete laplacian stuff, but wouldn't that just mean I take the data in grid form then average each cell from its peers?

##### Share on other sites
x is of length N. There is one entry in x for each cell in the grid. Ax=b is just a concise mathematical notation for a system of linear equations.

For example this system:
x + y + z = 1x + y - z = 9x - y + z = -5
could be written as Ax=b, with A as:
1  1  11  1 -11 -1  1
x as the column vector [x y z]
and b as [1 9 -5]

The methods being discussed, like gauss seidel relaxation and conjugate gradient, are just efficient methods for numerically solving systems of equations. The systems can also be solved by inverting the matrix, or by using gaussian elimination. The iterative methods tend to be recommended for this purpose because they are robust and an exact solution is not required for fluid simulation in games(more iterations get you closer to the solution, so you can trade off quality and speed).

##### Share on other sites
Quote:
 Original post by Vorpyx is of length N. There is one entry in x for each cell in the grid. Ax=b is just a concise mathematical notation for a system of linear equations.

For some fluids simulations (e.g., full/unreduced Euler, Navier-Stokes formulations) a cell's state value (x) can actually be a vector, e.g., density/mass, scalar momentum for each coordinate direction, some representation of energy. Of course you can "unroll" a cell to end up with multiple rows/unknowns per grid cell. For a 2D Euler solver, for example, you'd have 4 entries per cell (mass, 2 momentums, 1 energy), so a 100 x 100 grid would have not 10000, but 40000 unknowns.

Quote:
 Original post by VorpyThe methods being discussed, like gauss seidel relaxation and conjugate gradient, are just efficient methods for numerically solving systems of equations. The systems can also be solved by inverting the matrix, or by using gaussian elimination. The iterative methods tend to be recommended for this purpose because they are robust and an exact solution is not required for fluid simulation in games(more iterations get you closer to the solution, so you can trade off quality and speed).

To expand on Vorpy's comments...I can't really recommend ever using direct inversion for many problems. First, direct matrix inversion + back substitution can simply be too expensive as the number of unknown goes up, e.g., iterative methods can often produce a solution faster, as Vorpy said. Not only faster, but MUCH faster. Even if you do converge fully to a machine-epsilon accurate solution that is the best you can expect with floating point math. In particular, techniques such as multigrid and some conjugate gradient implementations can converge in approximately O(n) time, giving very accurate results in far shorter time then Gaussian Elimination or other direct inversion techniques. Second, direct inversion methods tend to accumulate floating point roundoff error as the number of unknowns goes up, and thus may produce a very inaccurate result, while iterative methods do not have quite the same roundoff problem, though they suffer from their own stability/convergence problems that are very much a function of the choice of governing equation as well as the choice of discretization/computational molecule. (Our guideline, when I was studying CFD as an undergraduate, was that with more than 1000 unknowns, direct inversion would not produce an accurate enough result for engineering purposes. 1000 unknowns isn't much at all for fluids...its a 32x32 grid in 2D (which can be okay enough for some effects), or 10x10x10 in 3D (maybe bare minimum for any effect at all). You definitely would want to avoid O(n2) or greater for semi-large or very larger problems...and Gaussian Elimination is O(n3)....not inexpensive at all. Only use for very small/tiny problems in real-time, as iterative techniques are almost always going to be faster.

##### Share on other sites
Hi Winegums,

Out of curiosity, what is the aim of your solver? Real-time or offline?
Are you planning on doing free-surface simulations(i.e using level-sets to track the surface) or just smoke simulation?

If you are doing realtime stuff (meaning surfaces are off the table), I wouldn't bother with CG at all. The crappy Gauss-Seidel solver given in Stam's paper actually works pretty damn good. Instead of worrying about the solver part, try focusing on the advection part of your algorithm as it will significantly improve the visual quality way more than the introduction of a CG solver ever will.
I would recommend the MacCormack method presentein this paper:
http://physbam.stanford.edu/~fedkiw/papers/stanford2006-09.pdf

If you are, however, doing off-line simulations you better be prepared to bang your head in the wall at full tilt over and over again over a few months time. I know I did :). First off, you need to switch to a staggered("MAC") grid since this improves simulation quality and stability quite a lot.

The best thing is to start with Robert Bridson's excellent Siggraph fluid course as it does a very good job of breaking down academic papers into graspable parts: http://www.cs.ubc.ca/~rbridson/fluidsimulation/

Furthermore, I would recommend reading Mark Carlson's PhD thesis as it explains boundary conditions pretty well in a practical manner: http://www.cc.gatech.edu/~carlson/papers/carlson-thesis.pdf
(Chapter 3 in particular).

Let me know what specific questions you have and I'll try to answer them..

Cheers and Good Luck!

/M

##### Share on other sites
I did some fluid sim work a while back and I have a Sparse Matrix solver if you need it.

http://ksk1.net/dl/SMatrix.h
http://ksk1.net/dl/SMatrix.cpp

Also this is a very cool paper from BYU explaining a lot of the math. (esp the MAC grid and the general Navier Stokes solver for the interface)

FluidFlowForTheRestOfUs.pdf

---
Note that I did not write this solver. I forgot who I got it from. But I know it to work pretty well.
Also please do not use it in any commercial setting, I'm just providing it as a stop-gap till you write your own or find a GPL-ish library.

[modified by grhodes_at_work to put in a working link to the BYU paper]

[Edited by - Purple-Haze on February 25, 2008 2:15:17 PM]

##### Share on other sites
mhennix makes a very good point. Although fancier methods such as CG (and the multigrid method I mentioned) can give blazingly fast convergence times...they are often very difficult to tune, and might not provide fast convergence, or even any convergence, for general problems. The more basic iterative methods are more well-rounded in this regard.

##### Share on other sites
thanks for all the replies, I've been off trying to do other bits of my dissertation and totally forgot about this thread :/

Quote:
 Original post by mhennixHi Winegums,Out of curiosity, what is the aim of your solver? Real-time or offline?Are you planning on doing free-surface simulations(i.e using level-sets to track the surface) or just smoke simulation?

The aim was to get a real time fluid simulation (in 2d and 3d), while looking at how realistic looking it can be made without the frame rate being impacted too heavily.

I don't want to flake, as I think I need to see this through to prove that the last few weeks of research haven't been a waste, but are there easier methods to implement? I think learning numerical approximations as well as fluid simulations AND trying to tie those into code all at once has proven too much for me. I looked at runge-kutta, which we've done in class before, and I think that might be easier to impliment.

Thanks for all the links too, I think I have a busy evening of reading ahead of me..

##### Share on other sites
Quote:
 Original post by WinegumsI looked at runge-kutta, which we've done in class before, and I think that might be easier to impliment.

Hmm...well, there are two things to consider in fluids simulation. There is the spatial problem...distributing/diffusing properties in space in such a way that the governing equations and boundary/initial conditions are satisfied at an instant in time. And then there is the time problem, updating properties over time---time marching. In reality, they are coupled of course. But in numerical simulation they often are separated. Part of Stam's method is a spatial only step, solving for velocity on a fixed grid. This part of his technique is formulated in the Eulerian point-of-view (note that this has nothing to do with the Euler equations of fluid flow and it has nothing to do with the two-step Euler integrators). The other part is the time marching part, in Stam's case using a Lagrangian point-of-view to move the particles that carry properties via an advection step. (The time marching part could be done in a Eulerian view also, and would simply apply advection at each grid cell instead of at particles.) The RK methods (much like Verlet integration popular in rigid body simulation) are formulated for use in solving the time marching part. The methods you brought up in your original post, CG, the thing I mentioned---multigrid, the other iterative methods such as Gauss Siedel, are applied to solving the spatial problem. And for fluids that have a basis in the Eulerian view, this is the way things are done, e.g., solve the spatial problem at an instant in time using methods suitable for solving systems of simultaneous equations, then use time integrators to advect properties in time, rinse and repeat. So, you can't really just say RK vs. CG. They are applied to different aspects of the problem.

So, there is another method that you might find more intuitive, and it is a purely particle-based/Lagrangian view approach. No grids directly involved (though they are involved to help speed things up). That method, one variation of this approach, is called Smoothed Particle Hydrodynamics (SPH). It uses kernel functions to represent the spatial problem, and the issue of simulation the spatial problem is embedded in the kernel functions...no need to do CG or Gauss-Siedel or whatever to solve the spatial problem. The main loop just does the time marching and the kernel functions are re-evaluated prior to each time step. Although this can be more intuitive, it probably isn't much easier in practice...just different. The kernel functions require that you know some N nearest neighbor particles for every particle and so you have to track that, which is where you might introduce some grid or spatial partition scheme to avoid searching through all particles for every step to find the neighbors. Further, for real-time you will have to have a limited number of particles and if you want to render something that looks like a fluid, you may have to generate a surface that emcompasses the particles, using some implicit surface technique such as marching cubes. The surface changes every time step, or every time you choose to reevaluate it. Messy. There are some advantages to purely Lagrangian-based fluids simulation (e.g., handling mixed types of fluids/fluid-solid interaction/splashing more readily), but it too can be really messy to do in real-time.

[Edited by - grhodes_at_work on March 3, 2008 12:45:18 AM]

##### Share on other sites
I did look at SPH when I started research for my dissertation proposal, but I couldn't really get into it and wrap my head around it, whereas the explanations of grid based stuff seemed to have an easier on ramp.

And now time is becoming a factor, and I don't think I have time to switch methods and really understand them. The particle method seemed more intuitive to me, but everything I saw was harder to disassemble than the grid stuff.

I found some sort of runge-kutta implimentation. However, I didn't realise until a wee while ago that I needed a function to use with the method, then I read your reply which concreted this realisation.

##### Share on other sites
after talking to my technical supervisor I decided to give this another crack.

So my first problem will be building the matrix. Lets ignore that for now.

The next problem I see is multiplying an N^2 matrix by a 2 * N matrix, how can this resolve?

Sorry for all the questions.

##### Share on other sites
Have you taken a linear algebra class?

##### Share on other sites
Quote:
 Original post by grhodes_at_workHave you taken a linear algebra class?

Yup. Believe it or not (and I can see why you wouldn't) I've been doing maths modules in university for the last four years. Last thing we did was some work on eigenvectors and inertia tensors.

EDIT: I suck, in my mind I was trying to do N^2 x N^2 by an N x N matrix.

DOUBLE EDIT: No, wait, I'm right. I still don't see how you can multiply uneven matricies. I think I've filled up a 3x3 with 0's to multiply it by a 4x4, but multiplying an N^2 x N^2 by a N x 2 seems a lot harder.

[Edited by - Winegums on March 7, 2008 9:11:20 AM]

##### Share on other sites
Quote:
Original post by Winegums
Quote:
 Original post by grhodes_at_workHave you taken a linear algebra class?

Yup. Believe it or not (and I can see why you wouldn't) I've been doing maths modules in university for the last four years. Last thing we did was some work on eigenvectors and inertia tensors.

EDIT: I suck, in my mind I was trying to do N^2 x N^2 by an N x N matrix.

DOUBLE EDIT: No, wait, I'm right. I still don't see how you can multiply uneven matricies. I think I've filled up a 3x3 with 0's to multiply it by a 4x4, but multiplying an N^2 x N^2 by a N x 2 seems a lot harder.

I didn't mean to offend. Your initial question about multiplying matrices was how to multiply an N^2 by a 2*N, written just like that. I assumed, since the notation wasn't standard, that by N2 you meant N x N, and by 2*N you meant 2 x N. Which isn't possible if you intended the standard form of rows x columns. So my analysis is what led me to ask the question.

And the same is true of your latest post, how to multiply an N2 x N2 times an N x 2. Not possible. The number of columns in the left-most matrix must match the number of rows in the right-most matrix.

I think you are not formulating the problem correctly. Which leads to the question, where are you getting your N2? A normal finite difference style fluids solver would have either an N x N, with each cell being a vector and with matching vector boundary/initial conditions, or a kN x kN, with k being the number of scalar properties associated with a cell (so for 2D, Euler governing equations, k would be 4, for mass, x momentum, y momentum, and energy). And the boundary conditions/initial conditions would match that size, e.g., kN per cell.

##### Share on other sites
Quote:
Original post by grhodes_at_work
Quote:
Original post by Winegums
Quote:
 Original post by grhodes_at_workHave you taken a linear algebra class?

Yup. Believe it or not (and I can see why you wouldn't) I've been doing maths modules in university for the last four years. Last thing we did was some work on eigenvectors and inertia tensors.

EDIT: I suck, in my mind I was trying to do N^2 x N^2 by an N x N matrix.

DOUBLE EDIT: No, wait, I'm right. I still don't see how you can multiply uneven matricies. I think I've filled up a 3x3 with 0's to multiply it by a 4x4, but multiplying an N^2 x N^2 by a N x 2 seems a lot harder.

I didn't mean to offend. Your initial question about multiplying matrices was how to multiply an N^2 by a 2*N, written just like that. I assumed, since the notation wasn't standard, that by N2 you meant N x N, and by 2*N you meant 2 x N. Which isn't possible if you intended the standard form of rows x columns. So my analysis is what led me to ask the question.

And the same is true of your latest post, how to multiply an N2 x N2 times an N x 2. Not possible. The number of columns in the left-most matrix must match the number of rows in the right-most matrix.

I think you are not formulating the problem correctly. Which leads to the question, where are you getting your N2? A normal finite difference style fluids solver would have either an N x N, with each cell being a vector and with matching vector boundary/initial conditions, or a kN x kN, with k being the number of scalar properties associated with a cell (so for 2D, Euler governing equations, k would be 4, for mass, x momentum, y momentum, and energy). And the boundary conditions/initial conditions would match that size, e.g., kN per cell.

Sorry I think I've misexplained myself.

I thought, based on some things said earlier in the thread, that to get matrix A you take each value of your initial matrix as a row (with five nonzero elements per row). Since the initial matrix is NxN (N2) I thought A would therefore be an N2 by N2 matrix.

##### Share on other sites
Quote:
 Original post by WinegumsSorry I think I've misexplained myself.I thought, based on some things said earlier in the thread, that to get matrix A you take each value of your initial matrix as a row (with five nonzero elements per row). Since the initial matrix is NxN (N2) I thought A would therefore be an N2 by N2 matrix.

Ah, :). I'm re-reading posts myself..so... Actually, looking back at Vorpy's post, if the idea is that N is the # of cells in each direction, and you have a 2D grid (e.g., 2 directions) making the grid of cells an NxN grid and N2 total cells, then in fact it is correct that the matrix is a spare N2 x N2. And using two point differences then as Vorpy said, and I think I earlier confirmed, 5 nonzero entries per row. I didn't pay enough attention to the earlier posts, and got into the idea that N was the total number of cells. My bad, I apologize.

So, then, in this case, x and b will be transposed vectors of length N2, which multiply perfectly well with the A matrix of size N2 x N2.

I will again state that depending on the particular form of the governing equations you're solving, each entry in A might not be a scalar, e.g., a cell could represent 4 sub-equations for a 2D problem based on the Euler's formulation of unsteady, rotational, inviscid, compressible flow with some energy equation to close the system, and in this case each cell of A would be a sub 4x4 matrix, with each entry of x and b being vectors of length 4. If you're using Stam's approach, he has reformulated the problem so A represents velocity, e.g., it is velocity he is diffusing, and for 2D that is a 2 component vector. I believe he solves the two components separately, so you end up with two solves instead of one big 2N2 x 2N2 problem.
I hope this little anectdote doesn't confuse too much.

##### Share on other sites

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

## Create an account

Register a new account

• ### Forum Statistics

• Total Topics
628714
• Total Posts
2984357

• 23
• 11
• 10
• 13
• 14