Applying Conjugate Gradient Method to Fluid Simulations
Author
Quote:Original post by grhodes_at_work
Have 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]
Quote:Original post by WinegumsQuote:Original post by grhodes_at_work
Have 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.
Author
Quote:Original post by grhodes_at_workQuote:Original post by WinegumsQuote:Original post by grhodes_at_work
Have 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.
Quote:Original post by Winegums
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.
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.
While it is true that A is indeed a sparse matrix, all you ever need to consider in terms of data structures is storing some extra bits per grid cell, so don't waste any energy trying to code some sparse matrix system.
In the case of stam, he solves Nabla^2 * q = Nabla w3. (A=Nabla^2)
The (Nabla^2) is the laplacian and depending on the boundary conditions it just boils down to multiplying the center cell with some sum and then subtracting cell neighbours. In the case of having no boundaries(and in 2D), we get a stencil like
..-1
-1 4 -1
..-1
Which means that taking the laplacian is the same thing as multing by A is the same thing as looping through your grid and multyplying the center of the stencil by 4 and then subtracting the neighbouring cells. But remember, the setup of the stencil relies on what boundary conditions you have. Try reading Mark Carlson's PhD thesis, chapter 3, where he discusses this. It's not really that hard, boils down to summing and removing bits depending on if you have air/solid boundary etc.
So it is quite enough to store something like(in 2D)
uint8 m_neighBours : 4; //one bit per face
uint8 m_numNeighbours : 3; // max 4
uint8 __pad : 1;
per cell and then just adjust your algorithms to make use of this instead.
Another thing you might consider skipping is the diffusion part of stams solver. Most of the diffusion is actually very undesirable since you always have too much numerical diffusion anyways.
[Edited by - mhennix on March 9, 2008 5:38:05 AM]
In the case of stam, he solves Nabla^2 * q = Nabla w3. (A=Nabla^2)
The (Nabla^2) is the laplacian and depending on the boundary conditions it just boils down to multiplying the center cell with some sum and then subtracting cell neighbours. In the case of having no boundaries(and in 2D), we get a stencil like
..-1
-1 4 -1
..-1
Which means that taking the laplacian is the same thing as multing by A is the same thing as looping through your grid and multyplying the center of the stencil by 4 and then subtracting the neighbouring cells. But remember, the setup of the stencil relies on what boundary conditions you have. Try reading Mark Carlson's PhD thesis, chapter 3, where he discusses this. It's not really that hard, boils down to summing and removing bits depending on if you have air/solid boundary etc.
So it is quite enough to store something like(in 2D)
uint8 m_neighBours : 4; //one bit per face
uint8 m_numNeighbours : 3; // max 4
uint8 __pad : 1;
per cell and then just adjust your algorithms to make use of this instead.
Another thing you might consider skipping is the diffusion part of stams solver. Most of the diffusion is actually very undesirable since you always have too much numerical diffusion anyways.
[Edited by - mhennix on March 9, 2008 5:38:05 AM]
Author
Quote:Original post by mhennix
Another thing you might consider skipping is the diffusion part of stams solver. Most of the diffusion is actually very undesirable since you always have too much numerical diffusion anyways.
I had noticed that since viscosity is 0 by default the diffusion step is quite pointless. Despite this, I have proposed implementing the CG method, and I'd really like to do this.
In saying that I'm thinking this is beyond me. I'm having issues understanding the problem let along the solution. grhodes_at_work you said x and b are transposed vectors, but they're matricies. The only way I could see this working is if you stretch the NxN matrix to be a N2x1 matrix (a vector). Does maths let you do this?
So I have the equation
Ax = b
where A is a matrix, and x and b are vectors.
I'm supposed to know A and b, and want to find x, which led me to:
So if I get A through the means suggested by Vorpy (which, interestingly, will just give me 1 along the diagonal when viscosity is 0), and we know x0 is b.
Actually, if viscosity (and, thus, a) is zero, won't that mean x = b?
So, do I just guess x, then multiply x and A and see what I get, and accept it if the value isn't too far off b? I'm guessing I'd start with x0 = b, since I'm advecting from b.
so error e would be (b - xn-1)
then, xn would be...xn-1 + e?
and if e < some threshold value, bail?
Conjugate gradient is a method for getting successively better guesses for x. So, yes, ultimately the iterative methods are methods of improving the guesses until they are good enough. Gauss-seidel relaxation is just a different way of accomplishing the same thing.
I think in Stam's fluid paper he uses the solver for performing the diffusion and projetion steps. It's for the projection step that he suggests the conjugate gradient method, since the projection is more important to the simulation quality.
In the paper, he doesn't even check to see if x is within an error tolerance, he just runs it for a certain number of iterations and it works well enough. Conjugate gradient can work the same way. Checking if the solution is within the required tolerance requires more calculation, but it can allow the iteration to cut off early if the solution is good enough, so it's a tradeoff.
As for viscosity affecting the matrix, I think you're misreading something. The solver is used for the diffusion step, where the diffusion constant is important, and in the projection step, which is independent of any physical constants. Maybe you mean that the diffusion is 0, and in that case, yes, there's no reason to perform the diffusion step since there would be no diffusion and the matrix would be the identity matrix. The relaxation method may be better for the diffusion step anyway since its accuracy isn't as important.
And about this N vs N^2 stuff: there are two different views here. There is the grid representing the fluid, in which you're taking the physical space and overlaying a grid, so each cell is a physical space. And then there is the mathematical matrix representing the system, where there is one row and one column in the matrix for every cell in the grid, and the matrix represents the interactions between the different cells, and most cells don't interact so the grid is very sparse. You do not explicitly build this matrix because it is so sparse and the simulation grid provides a much more convenient representation. The matrix is useful as a mathematical abstraction, not as the implementation.
I think in Stam's fluid paper he uses the solver for performing the diffusion and projetion steps. It's for the projection step that he suggests the conjugate gradient method, since the projection is more important to the simulation quality.
In the paper, he doesn't even check to see if x is within an error tolerance, he just runs it for a certain number of iterations and it works well enough. Conjugate gradient can work the same way. Checking if the solution is within the required tolerance requires more calculation, but it can allow the iteration to cut off early if the solution is good enough, so it's a tradeoff.
As for viscosity affecting the matrix, I think you're misreading something. The solver is used for the diffusion step, where the diffusion constant is important, and in the projection step, which is independent of any physical constants. Maybe you mean that the diffusion is 0, and in that case, yes, there's no reason to perform the diffusion step since there would be no diffusion and the matrix would be the identity matrix. The relaxation method may be better for the diffusion step anyway since its accuracy isn't as important.
And about this N vs N^2 stuff: there are two different views here. There is the grid representing the fluid, in which you're taking the physical space and overlaying a grid, so each cell is a physical space. And then there is the mathematical matrix representing the system, where there is one row and one column in the matrix for every cell in the grid, and the matrix represents the interactions between the different cells, and most cells don't interact so the grid is very sparse. You do not explicitly build this matrix because it is so sparse and the simulation grid provides a much more convenient representation. The matrix is useful as a mathematical abstraction, not as the implementation.
Maybe this is perfectly clear, but here goes:
You are correct in your assumption that x and b are N^2 x 1 vectors in theory. In practice however, you would still consider them to be grids(matrices).
A is a N^2 x N^2 matrix but in practice this boils down to accessing neighbours.
One row of A corresponds to the equation for one grid cell and in order for this to work out one needs to form N^2 vectors and N^2 x N^2 matrices. We want to form the laplacian for every single grid cell and in order to to this every single grid cell needs to be able to "see" all other cells.
in the end we want an equation like (per grid cell, depending on boundary conditions):
4*x[i,j] -x[i-1,j] -x[i,j-1] -x[i,j+1] -x[i+1,j] = b[i,j]
Which is exactly what you will get when multyplying the first row of A with x.
Here you can see that forming the big matrix is not necessary since we will always be working on grid cells that are x+-1,y+-1 away.
You are correct in your assumption that x and b are N^2 x 1 vectors in theory. In practice however, you would still consider them to be grids(matrices).
A is a N^2 x N^2 matrix but in practice this boils down to accessing neighbours.
One row of A corresponds to the equation for one grid cell and in order for this to work out one needs to form N^2 vectors and N^2 x N^2 matrices. We want to form the laplacian for every single grid cell and in order to to this every single grid cell needs to be able to "see" all other cells.
in the end we want an equation like (per grid cell, depending on boundary conditions):
4*x[i,j] -x[i-1,j] -x[i,j-1] -x[i,j+1] -x[i+1,j] = b[i,j]
Which is exactly what you will get when multyplying the first row of A with x.
Here you can see that forming the big matrix is not necessary since we will always be working on grid cells that are x+-1,y+-1 away.
Author
sorry for bumping this, and thanks again for the continued patience.
mhennix, given your example of
4*x[i,j] -x[i-1,j] -x[i,j-1] -x[i,j+1] -x[i+1,j] = b[i,j]
I'm a bit confused. I thought b corresponded to x0 in the diffusion function? Or are you saying this is Ax = b written out with the laplacian of A rather than the big matrix?
If so, would I have to just rearrange to get x[i,j], which is what I want to find?
I tried this and got
x[i,j] = (b[i,j] + [i-1,j] + x[i,j-1] + x[i,j+1] + x[i+1,j]) / 4
Which didn't work (everything blew up horribly). However changing that 4 to a 5 gave a stable, if kinda dodgy solver (I think density vanishes over time).
mhennix, given your example of
4*x[i,j] -x[i-1,j] -x[i,j-1] -x[i,j+1] -x[i+1,j] = b[i,j]
I'm a bit confused. I thought b corresponded to x0 in the diffusion function? Or are you saying this is Ax = b written out with the laplacian of A rather than the big matrix?
If so, would I have to just rearrange to get x[i,j], which is what I want to find?
I tried this and got
x[i,j] = (b[i,j] + [i-1,j] + x[i,j-1] + x[i,j+1] + x[i+1,j]) / 4
Which didn't work (everything blew up horribly). However changing that 4 to a 5 gave a stable, if kinda dodgy solver (I think density vanishes over time).
This topic is closed to new replies.
Advertisement
Popular Topics
Advertisement
