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]

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.

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?

>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

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.

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]

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;}

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.