Can someone point me to resources on implementing NSE solvers? Tutorials on the web would be prefferable, but books and and papers are good too. I spent a lot of time looking, and couldn''t really find anything on google.

Here are three nice papers:

1.) Stable Fluids by Jos Stam (1999) - presents an elegant scheme for solving the equations robustly and with speed. A great starting point.

2.) Visual simulation of Smoke (2001) by Ron Fedkiw, Jos Stam and Henrik Jensen. Pretty much equal to the first paper, but it explains some tehniques for improving the accuracy.

3.) The Motion of Hot, Turbulent Gas (1997) by Nick Foster and Dimitri Metaxas. Shows a typical explicit scheme for solving the equations. IMHO not as good as Stam's method (especially for starters) but it should made clear that this is pretty much the "traditional" view of solving ordinary differential equations.

How to find these papers? Just use google.
Btw, you'll probably make a great game with these techniques ;P

- Mikko Kauppila

[edited by - uutee on October 31, 2002 12:20:38 PM]

Thank you very much for your responce. The papers are *very* useful. They are, however, a little above my head May be I''ll take differential equations class next semester. Anyway, I''m sure if I stare at them enough I will eventually get there.

For now, if anyone has some source code of an NSE solver I would really like to take a look at it. Again, thank you for your help.

// here is the prototype for the modulus function, it should be// quite straightforward to codeinline int IMod( int a, int b );   // integer modulus, [0,b-1]// here are the prototypes for the interpolators. the source arrays// are of the size N*NVector2D InterpolateV( Vector2D * V, int N, real xx, real yy );real InterpolateS( real * s, int N, real xx, real yy );// V0 and S0 are the actual simulation grids of size N*N. V0 is the// velocity field and S0 is the "substance" field which you want to// advect through the liquid (for example smoke densities).// the old fields are initially stored in these grids and the// new fields are returned in these grids.// V1 and S1 are temporary grids of size N*N, they can be initialized// with any value. they return nothing in them.// F is the force grid of size N*N, user supplied.// N is the size of the simulation grids (each grid is N*N)// dt is the time step// FluidMatrixMultiplyvoid FluidSimulate( Vector2D * V0, Vector2D * V1, real * S0, real * S1,                    Vector2D * F, int N, real dt,                    void (*FluidMatrixMultiply)(real * src, real * dst) ){    int i, x, y;/************************ STEP 1 - ADD FORCES ******************************/    for( i=0; i        V0 += F*dt;<br><br>/*********************** STEP 2 - ADVECT FIELDS ****************************/<br><br>    for( y=0; y<N; y++ )<br>    for( x=0; x<N; x++ )<br>    {<br>        real xx = (real)x - V0[y*N+x].x*dt;  // trace back in time<br>        real yy = (real)y - V0[y*N+x].y*dt;<br><br>        V1[y*N+x] = InterpolateV( V0, N, xx, yy ); // interpolate to the<br>        S1[y*N+x] = InterpolateS( S0, N, xx, yy ); // temp. grids<br>    }<br><br>    for( i=0; i<N*N; i++ )<br>    {<br>        V0 = V1;   // move temp. grids to the actual <br>        S0 = S1;   // simulation grids<br>    }<br><br>/***************** STEP 3 - ENFORCE INCOMPRESSIBILITY **********************/<br><br>    real * mb = new real[N*N];   // gradient dot velocity<br>    real * mx = new real[N*N];   // pressure which we want to solve for<br><br>    for( y=0; y<N; y++ )<br>    for( x=0; x<N; x++ )<br>    {<br>        // zero the pressure…<br>        mx.n[y*X+x] = 0;<br>        // compute gradient dot velocity for this point<br>        mb.n[y*N+x] = (V0[y*N+IMod(x+1,N)].x - V0[y*N+IMod(x-1,N)].x)/2.0+<br>                      (V0[IMod(y+1,N)*N+x].y - V0[IMod(y-1,N)*N+x].y)/2.0;<br>    }<br><br>    // solve the linear system<br>    CG<real>( FluidMatrixMultiply, mx, mb, 1e-4, N*N );<br><br>    // subtract the gradient of pressure from the <br>    for( y=0; y<N; y++ )<br>    for( x=0; x<N; x++ )<br>    {<br>        Vector2D g;    // compute the gradient of pressure<br>        g.x = (mx.n[ y*N+IMod(x+1,N) ] - mx.n[ y*N+IMod(x-1,N) ])/2.0;<br>        g.y = (mx.n[ IMod(y+1,N)*N+x ] - mx.n[ IMod(y-1,N)*N+x ])/2.0;<br>        V0[ y*N+x ] -= g;    // subtract it from the velocity<br>    }<br>}<br> </pre> <br><br><SPAN CLASS=editedby>[edited by - uutee &#111;n November 5, 2002 12:25:20 PM]</SPAN>  </i>   

Once again the lovely forum screwed up my indents! But yes, I tailored a specially simplified version of the simulator for you to study. There are two "black-box" thing not explained in the code: the other one is the CG (Conjugate Gradient) routine, which has the prototype:

int CG( void (*MatrixMultiply)(real * src, real * dst), real * X, real * B, real eps, int N );

This solves the linear system A*X=B, where the multiplication by A is handled through the user supplied function MatrixMultiply, X is initially filled with zeroes and B is also user supplied. N is the size of the vectors, it''s NOT equal to the N used in the other parts of the code (Ncg = Nfluid*Nfluid, it''s the square).

The other black-box is the MatrixMultiply routine, or more specifically the matrix used in it. The matrix has a following form: for each point (x,y) in the simulation grid we assign one row of the matrix (row = y*N+x, so there are N*N rows).

Each row has a following structure: the diagonal is set to -4, and all the adjacent grid points [those are (x+1,y) (x-1,y) (x,y+1) (x,y-1)] are assigned one column from that row (column = y*N+x, where x,y are the coordinates of the adjacent grid point). The values at those matrix nodes are set to 1. All other nodes on the row are set to zero (if x-1 < 0 or x=N or whatever, we warp the coordinate around the grid).

At this point it should be obvious that we must use a special data structure for our matrix, which is of the size (N*N)x(N*N) (where N is the fluid grid size)! It''s easy to notice that each row consists of only five nonzero entries (the diagonal and the adjacent points). Therefore, matrix storage can be handled in just O(n) amount of memory.

It''s not really the easiest thing at the first glance, but once you get used to it, you''ll find it a very elegant and powerful technique.

- Mikko Kauppila

uutree: Thank you very much for taking the time to post source code and explain everything in detail! It makes things a lot clearer. Again, thank you very much for your time.

uutree,

If you wrap your source code with "" and "" (removing the spaces) then your indents will not be lost, .

Graham Rhodes
Senior Scientist
Applied Research Associates, Inc.