Archived

This topic is now archived and is closed to further replies.

pickups

fluid dynamics - projection step.

Recommended Posts

hello. im working on a fluid solver, and i implemented it in 2D successfuly. but, when i try it in 3D, every thing works except for the projection step, what practically means the flui''s velocity explodes. i tried to fix it, but i dont know what to do next. did anybody of you ever try this stuff, or even succeeded? i would really appriciate your help, thanx.

Share this post


Link to post
Share on other sites
First of all, you didn''t really tell how you''re solving the projection step. Do you use a relaxation scheme (as in Foster&Metaxas and Stam''s GDC paper) or a (much better) sparse linear system solver (Stam and Fedkiw papers)?

Furthermore something needs to be known about the boundary conditions you''re using. Periodic or solid walls, perhaps?

Something needs to be known about your discretization, too. Are you using a MAC staggered grid or a cell-centered grid?

- Mikko Kauppila

Share this post


Link to post
Share on other sites
On the other hand, the fact that your solver actually did work in 2D suggests of some bug or a misinterpreted formula in the 3D version.

If you''re forming a sparse linear system, what kind of linear equation are you using in 3d?

- Mikko Kauppila

Share this post


Link to post
Share on other sites
all right, ill fill the details.

im using a relaxation solver, gauss-seidel solver. i tried to find some info about the conjugate gradient method but couldnt find anything, and i dont have access to university books.

i use a cell centerd grid, and my boundary conditions are fixed walls, i think it is called free-slip, becuase on the horizonal walls i keep the horizonal velocity and on the vertical walls i keep the vertical velocity.
i used the code from stams GDC 2003 paper, and turned it to 3D.
i compute the velocity divergent with the following code, from stams paper:

Div[i,j,k] = 0.5*h*(U[i+1,j,k] - U[i-1,j,k] +
V[i,j+1,k] - V[i,j-1,k] + W[i,j,k+1] - W[i,j,k-1])

altough it should be: 0.5*(...)/h. stams doesnt do it like this, so i followed him.

then, i compute the pressure filed like this:

P[i,j,k] = (P[i+1,j,k] + P[i-1,j,k] +
P[i,j+1,k] + P[i,j-1,k] + P[i,j,k+1] + P[i,j,k-1])-Div[i,j,k])/6

this is what stam does, although in the math it should be :
....-h*h*Div[])/6.

what do you think?

Share this post


Link to post
Share on other sites
Get ghostscript and GSView here, for reading .PS files:
http://www.cs.wisc.edu/~ghost/doc/AFPL/get800.htm

The "canned algorithms" section in Shewchuk''s paper gives straight pseudocodes for implementing the methods, so it is only a few minutes'' job.

I really suggest implementing a version with periodic boundary conditions at first if you''re doing the CG method for solving Poissons''s equation on a cell-centered grid.

- Mikko Kauppila

Share this post


Link to post
Share on other sites
allright thank you.
i will read this paper, im sure ill understand it, it looks good.
but what do you mean periodic boundary conditions? why is it better?
and one more thing, what should i expect after implementing CG, visually speaking? i mean, i actually have a working 3D solver, with a gauss-seidel solver for the projection. it works nice, but i can bearly see small scale swirls. will CG improve it, or i should implement vorticity confinement?

thanks for your help.

Share this post


Link to post
Share on other sites
>but what do you mean periodic boundary conditions?

Periodic boundary conditions is when the fluid wraps around...
In semi-Lagrangian integration this means that if your particle "escapes" the fluid domain, it enters from the opposite side.

Also when constructing the Poisson linear equation it means that you should wrap your coordinates... it should be quite trivial, really.

>why is it better?

It''s not better. It''s just easy to implement (on a cell-centered grid).

You see, if your matrix in the Poisson equation (Ax=b when discretized) is singular, as it usually is, the sum of the nodes in the vector b should be zero. That is most easy to achieve with periodic boundary conditions. In a MAC staggered grid it''s also easy to achieve even with solid walls.

>and one more thing, what should i expect after implementing
>CG, visually speaking?

You should see more turbulence, that''s for sure. Computationally speaking you''ll also see a huge speedup, as only 20 iterations produces a near-perfect solution.

- Mikko Kauppila

Share this post


Link to post
Share on other sites
My previous post was done in a hurry, so maybe I should clarify a bit:

When a fluid wraps around the coordinate x is supposed to be enforced to be in the bounds [0,N[ (for integers that's of course [0,N-1]). Normal modulo operator (%) isn't too good in this because it can't handle negative values in the way we'd like to (if I remember right). A better solution is:

inline int IMod( int a, int N )
{
while( a < 0 ) a += N;
while( a >= N ) a -= N;
return a;
}

// real is either float or double
inline real RMod( real a, int N )
{
while( a < 0 ) a += N;
while( a >= N ) a -= N;
return a;
}

This is of course not readily usable for the grid that Stam's GDC2k3 paper uses, because it has extra "ghost" values at points 0 and N+1, but I'm sure handling this won't be a big problem for wrapping points around.

Btw, if you find a nice (elegant) way to encapture boundary conditions in the cell-centered grid of the GDC2k3 paper, please tell I'm currently constructing a MAC staggered velocity field out of the cell-centered velocities (by averaging), because that allows me to reuse the old (convenient) code to do boundary conditions in a MAC staggered grid.

And more about what you should except visually: you should see more things than turbulence. That is global interactions - the Poisson equation is naturally of the kind that if you apply a force to one point (P1) on the grid, this force is felt at *all* points over the grid - not just on the grid points near P1. If you use a relaxation method, which relies on local interactions only, you'll need an awfully large amount of iterations to encapture global interactions well. In that huge number of iterations accuracy might be lost too, so CG will surely pay you back the time spent implementing it.

- Mikko Kauppila

[edited by - uutee on July 31, 2003 3:49:23 PM]

Share this post


Link to post
Share on other sites
Ok. ill try to implement it. i only have high-school math knolegde, becuase i just finished high-school, but i can handle it.

i dont think i understand what you mean about encapturing boundary-conditions. and also, to wrap around the coordinates, even for negative values, you can do this:
x'' = (n + (x%n))%n.

and maybe this also, although im not sure about negative values:
x'' = x & (n - 1).
this one at least works when n is a power of 2, which is OK.

Share this post


Link to post
Share on other sites
Hi.

i read the paper, and i understand it more or less. but one thing remains unclear.
my goal is to solve the poisson equation:

laplace(P) = divergent(U)

the CG method is used to solve Ax = b.

i just dont see where the two meat
i guess its pretty simple, but i dont get it. maybe divergent(U) is b and laplace(P) is A? so what is X? what does the poisson equation has to do with minimization?

Share this post


Link to post
Share on other sites
>i only have high-school math knolegdee

So have I. Stam''s solver is very simple, some high schools even teach linear algebra.

>i dont think i understand what you mean about encapturing >boundary-conditions

I meant encapturing some non-trivial boundary conditions, such as having solid walls for the boundaries. Of course you could talk about encapturing boundary conditions in a fluid that wraps around, but computationally it''s trivial.

>laplace(P) = divergent(U)
>the CG method is used to solve Ax = b.
>i just dont see where the two meat
>i guess its pretty simple, but i dont get it.
>maybe divergent (U) is b and laplace(P) is A? so what is X?

You''re right. divergence U is b, laplacian is A, and P is X (that is, we''re solving for a pressure field). You need to discretize divergence and laplacian operators to form a linear system of equations (Ax=b).

@ means differentiation operator, and I really have to assume you understand it (with high-school math one should). You should also know the basic central-difference discretizations for first- and second-order partial derivatives on a uniform grid (I''ll give them below).

laplacian P = @^2P/@x^2 + @^2P/@y^2 + @^2P/@z^2

@^2P/@x^2 = [P(x+1,y,z) - 2*p(x,y,z) + P(x-1,y,z)]/(h^2)
and the same for @^2P/@y^2 and @^2P/@z^2. Sum them up to get the laplacian.

divergence U(x,y,z) = @u/@x + @v/@y + @w/@z
(u,v,w are the x,y,z components of the velocity, respectively)
this is discretized as: @u/@x = [u(x+1,y) - u(x-1,y)]/(2h) (and same for @v/@y and @w/@z).

The final linear system will be of the form (both sides are multiplied by h^2 here):

P(x+1,y,z) + P(x,y+1,z) + P(x,y,z+1) - 6*P(x,y,z) + P(x-1,y,z) + P(x,y-1,z) + P(x,y,z-1) = (h/2)*[ u(x+1,y,z) - u(x-1,y,z) + v(x,y+1,z) - v(x,y-1,z) + w(x,y,z+1) - w(x,y,z-1)]

This is then very trivial to turn into a sparse matrix. Each diagonal entry is -6, and each non-diagonal column entry Y on a row X is 1 if and only if the point X is adjacent to Y on the grid (otherwise it''s a zero, which we don''t store).

- Mikko Kauppila

Share this post


Link to post
Share on other sites
Hi.
actually, where i live, they dont teach these things at all, this unfortunate, because i thinks it interesting.
anyway, i know all this differenciation stuff, so i know the equations you wrote below, i wrote them myself a few posts ago.
but i dont understand the part of how to converts the final formula t a spares matrix.

first of all, how can the divergent be a vector? i can think of it as a matrix, becuase for every point on the grid there is a divergent scalar.
and how can the laplacian be a matrix? it is just an operator.

Share this post


Link to post
Share on other sites
I guess, the key is to understand how many entries your sparse matrix contains. Suppose we have a simulation grid with N cells in each direction, that is: N*N*N cells for the whole grid. Our matrix for the linear system has (N*N*N)*(N*N*N) entries. So, in the case where you have a uniform cartesian grid, where dx=dy=dz you get a band matrix (3 bands). Since all entries, except those in the bands, are zero, you only use those nonzero values. This results in one equation per cell. So in the Gauss-Seidel Solver you rearrange this equation (laplace(p)=div(u)), in order to implicitly compute the inverse of the matrix, (when all grid cells are traversed).
So in your Gauss-Seidel-Solver you compute one solution vector of dimension (N*N*N).

Hope that helps. ;-)

Cheers,
Oliver

[edited by - pilarsor on August 1, 2003 7:21:11 AM]

Share this post


Link to post
Share on other sites
quote:
Original post by uutee
My previous post was done in a hurry, so maybe I should clarify a bit:

When a fluid wraps around the coordinate x is supposed to be enforced to be in the bounds [0,N[ (for integers that''s of course [0,N-1]). Normal modulo operator (%) isn''t too good in this because it can''t handle negative values in the way we''d like to (if I remember right). A better solution is:

inline int IMod( int a, int N )
{
while( a < 0 ) a += N;
while( a >= N ) a -= N;
return a;
}

// real is either float or double
inline real RMod( real a, int N )
{
while( a < 0 ) a += N;
while( a >= N ) a -= N;
return a;
}

This is of course not readily usable for the grid that Stam''s GDC2k3 paper uses, because it has extra "ghost" values at points 0 and N+1, but I''m sure handling this won''t be a big problem for wrapping points around.

Btw, if you find a nice (elegant) way to encapture boundary conditions in the cell-centered grid of the GDC2k3 paper, please tell I''m currently constructing a MAC staggered velocity field out of the cell-centered velocities (by averaging), because that allows me to reuse the old (convenient) code to do boundary conditions in a MAC staggered grid.

And more about what you should except visually: you should see more things than turbulence. That is global interactions - the Poisson equation is naturally of the kind that if you apply a force to one point (P1) on the grid, this force is felt at *all* points over the grid - not just on the grid points near P1. If you use a relaxation method, which relies on local interactions only, you''ll need an awfully large amount of iterations to encapture global interactions well. In that huge number of iterations accuracy might be lost too, so CG will surely pay you back the time spent implementing it.

- Mikko Kauppila

[edited by - uutee on July 31, 2003 3:49:23 PM]


Instead of using "real".. why not just...


template <class T>
inline T RMod(T a, int N )
{
while( a < 0 ) a += N;
while( a >= N ) a -= N;
return a;
}

Share this post


Link to post
Share on other sites
Ok, so the vector b is just a "streched" grid, storing the veolcity divergent at each point, which i calculate once beforce the projection.
about the matrix, i really dont get it . what is a matrix band? i still dont understand what to put where in this big matrix.

btw. thank you very much for your help.

Share this post


Link to post
Share on other sites