void Grid::DrawDensity(float width, float height)
{
float scaleX = width/(float)n;
float scaleY = height/(float)n;
glBegin(GL_QUADS);
for( int i = 0; i <= n; i++ )
{
// Subtract 0.5f because when forming a QUAD we need the lower left hand corner, not the middle
float x = (i)*scaleX;//(1/h);
for( int j = 0; j <= n; j++ )
{
float y = (j)*scaleY;//(1/h);
// We can sample the density field because it is practically just a height field of values between 0-1
float lowerLeft = dens[Index(i,j)];
float upperLeft = dens[Index(i,j+1)];
float lowerRight = dens[Index(i+1,j)];
float upperRight = dens[Index(i+1,j+1)];
//std::cout << lowerLeft << std::endl;
glColor3f(lowerLeft, lowerLeft, lowerLeft);
glVertex2f(x, y);
glColor3f(upperLeft, upperLeft, upperLeft);
glVertex2f(x, y+scaleY);
glColor3f(lowerRight, lowerRight, lowerRight);
glVertex2f(x+scaleX, y);
glColor3f(upperRight, upperRight, upperRight);
glVertex2f(x+scaleX, y+scaleY);
}
}
glEnd();
}
void Grid::DrawVelocity(float width, float height)
{
int i, j;
float x, y, h;
h = 1.0f/n;
int scaleX = width/n;
int scaleY = height/n;
glColor3f ( 1.0f, 1.0f, 1.0f );
glLineWidth ( 1.0f );
glBegin ( GL_LINES );
for ( i=1 ; i<=n ; i++ ) {
x = (i)*scaleX;
for ( j=1 ; j<=n ; j++ ) {
y = (j)*scaleY;
glVertex2f ( x, y );
glVertex2f ( x+vx[Index(i,j)]*scaleX, y+vy[Index(i,j)]*scaleY );
}
}
glEnd ();
}
Navier-Stokes based on Jos Stam's work
Author
Hello,
I'm trying to build a rather small fluid simulator based on Jos Stam's paper "Real-Time Fluid Dynamics for Games" from GDC 03. My problem is that my simulator seems to work irregularly (obviously :-P). When I add a source of density the density seems to diffuse outward but it disappear from where it was created and it just keep getting bigger. Then, it'll start to disappear from the inside out, creating "circular" ring shape (not that it's really all that circular). The velocity also doesn't seem to be acting normal either (but it's hard to tell).
As an extra note, the drawing routines included with the code of the paper didn't work at all for me, so I had to change them to scale the values by the ScreenWidth/NumGridCells and ScreenHeight/NumGridCells. Here's the code for the drawing, maybe it might be a clue for what I'm doing wrong?
Can you guys tell what might be wrong from the behavior I described above? I'll post code if no one can think of what might be wrong from above.
Thanks a million!
I didn't notice anything during a look over the code, although I didn't examine it in depth.
What I'm inclined to suggest is first testing your drawing functions. Disable the water simulation, and populate the cells with known values - some simple function or other is probably a good idea (perhaps val = c*dist_from_origin, or c/dist_for_origin, where 'c' is some constant, and, in the latter case, with an appropriate special case for the cell, if any, that lies at the origin). This should help you to check that the drawing functions are working properly, and, if not, hopefully lead you to correcting them.
That said, based on your description, could it not simply be that you're adding too much density, and ending up with a water crest that is higher than your system is managing to render, leading to the appearance of a flat-topped ring?
Perhaps a screenshot would help, in order to (hopefully) better get a sense of the bug.
What I'm inclined to suggest is first testing your drawing functions. Disable the water simulation, and populate the cells with known values - some simple function or other is probably a good idea (perhaps val = c*dist_from_origin, or c/dist_for_origin, where 'c' is some constant, and, in the latter case, with an appropriate special case for the cell, if any, that lies at the origin). This should help you to check that the drawing functions are working properly, and, if not, hopefully lead you to correcting them.
That said, based on your description, could it not simply be that you're adding too much density, and ending up with a water crest that is higher than your system is managing to render, leading to the appearance of a flat-topped ring?
Perhaps a screenshot would help, in order to (hopefully) better get a sense of the bug.
Author
Here's a screenshots:
http://img149.imageshack.us/img149/8274/picture1gi3.png
and if anyone has a mac, an app to actually try it out (press v t switch between drawing the velocity and density): http://dl.getdropbox.com/u/20237/Navier-Stokes%20Sim.zip
Obviously the drawing is wrong, but so is the simulation.
http://img149.imageshack.us/img149/8274/picture1gi3.png
and if anyone has a mac, an app to actually try it out (press v t switch between drawing the velocity and density): http://dl.getdropbox.com/u/20237/Navier-Stokes%20Sim.zip
Obviously the drawing is wrong, but so is the simulation.
Author
Well, I got my drawing routines to work properly, and the velocity it seems to work just like Stam's, so I'm pretty confident that the velocity portion of the code is correct, but since the diffusion uses the same diffusion and advection functions as the velocity I can't understand why it isn't working correctly... Here's my code, I hope someone might be able to find what's going wrong :)
.cpp
.h
.cpp
/* * Grid.cpp * Navier-Stokes * * Created by Jedd Haberstro on 5/22/08. * Copyright 2008 Jedd Haberstro. All rights reserved. * */#include <OpenGL/OpenGL.h>#include <OpenGL/gl.h>#include <iostream>#include <cmath>#include "Grid.h"#define SWAP(x0,x) {float * tmp=x0;x0=x;x=tmp;}void Grid::SimStep(){ // Add forces /*AddSource(vx, vxOld); AddSource(vy, vyOld); SWAP(vxOld, vx); Diffuse(1, vx, vxOld, viscosity); SWAP(vyOld, vy); Diffuse(2, vy, vyOld, viscosity); Project(vx, vy, vxOld, vyOld); SWAP(vxOld, vx); SWAP(vyOld, vy); Advect(1, vx, vxOld, vxOld, vyOld); Advect(2, vy, vyOld, vxOld, vyOld); Project(vx, vy, vxOld, vyOld); AddSource(dens, densOld); SWAP(densOld, dens); Diffuse(0, dens, densOld, diffusionRate); SWAP(densOld, dens); Advect(0, dens, densOld, vx, vy);*/ //Velocity step AddSource(vx, vxOld); AddSource(vy, vyOld); Diffuse(1, vxOld, vx, viscosity); Diffuse(2, vyOld, vy, viscosity); Project(vxOld, vyOld, vx, vy); Advect(1, vx, vxOld, vxOld, vyOld); Advect(2, vy, vyOld, vxOld, vyOld); Project(vx, vy, vxOld, vyOld); //Density step AddSource(dens, densOld); Diffuse(0, densOld, dens, diffusionRate); Advect(0, dens, densOld, vx, vy);}void Grid::DrawDensity(float width, float height){ int i, j; float px, py; float scaleX = width/(float)n; float scaleY = height/(float)n; /*for( j = 0; i < n-1; j++ ) { glBegin(GL_TRIANGLE_STRIP); i = 0; px = scaleX + i * scaleX; py = scaleY + j * scaleY; glColor3f(dens[Index(i,j)], dens[Index(i,j)], dens[Index(i,j)]); glVertex2f(px, py); for( i = 0; i < n-1; i++ ) { px = scaleX + i * scaleX; py = scaleY + (j + 1) * scaleY; glColor3f(dens[Index(i,j)], dens[Index(i,j)], dens[Index(i,j)]); glVertex2f(px, py); px = scaleX + (i + 1) * scaleX; py = scaleY + j * scaleY; glColor3f(dens[Index(i,j)], dens[Index(i,j)], dens[Index(i,j)]); glVertex2f(px, py); } px = scaleX + (n - 1) * scaleX; py = scaleY + (j + 1) * scaleY; glColor3f(dens[Index(i,j)], dens[Index(i,j)], dens[Index(i,j)]); glVertex2f(px, py); glEnd(); }*/ glBegin(GL_TRIANGLE_STRIP); for( int i = 0; i <= n; i++ ) { // Subtract 0.5f because when forming a QUAD we need the lower left hand corner, not the middle for( int j = 0; j <= n; j++ ) { // We can sample the density field because it is practically just a height field of values between 0-1 float lowerLeft = dens[Index(i,j)]; float upperLeft = dens[Index(i,j+1)]; float lowerRight = dens[Index(i+1,j)]; float upperRight = dens[Index(i+1,j+1)]; //std::cout << lowerLeft << std::endl; //glColor3f(1,0,0); glColor3f(lowerLeft, lowerLeft, lowerLeft); glVertex2f(i*scaleX, j*scaleY); glColor3f(upperLeft, upperLeft, upperLeft); glVertex2f(i*scaleX, (j+1)*scaleY); glColor3f(lowerRight, lowerRight, lowerRight); glVertex2f((i+1)*scaleX, j*scaleY); glColor3f(upperRight, upperRight, upperRight); glVertex2f((i+1)*scaleX, (j+1)*scaleY); } } glEnd();}void Grid::DrawVelocity(float width, float height){ int i, j; int scaleX = width/n; int scaleY = height/n; glColor3f ( 1.0f, 1.0f, 1.0f ); glLineWidth ( 1.0f ); glBegin(GL_LINES); for (i = 0; i < n; i++) { for (j = 0; j < n; j++) { glColor3f(1, 0, 0); glVertex2f(scaleX + i * scaleX, scaleY + j * scaleY); glVertex2f((scaleX + i * scaleX) + 1000 * vx[Index(i,j)], (scaleY + j * scaleY) + 1000 * vy[Index(i,j)]); } } glEnd();} // Add a source of density for a single cellvoid Grid::AddDensity(float s, int i, int j){ if( i > n || i < 0 || j > n || j < 0 ) return; densOld[Index(i,j)] += s*dt;} void Grid::AddForce(float fx, float fy, int i, int j){ if( i > n || i < 0 || j > n || j < 0 ) return; vxOld[Index(i,j)] += fx*dt; vyOld[Index(i,j)] += fy*dt;}void Grid::AddGlobalForce(float fx, float fy){ for( int i = 0; i < size; i++ ) { vxOld += fx*dt; vyOld += fy*dt; }}int Grid::ConvertPosToGrid(float mousePos, float size){ return ((mousePos/size)*(n+1));} //Gauss-Seidel relaxationvoid Grid::LinearSolve( int b, float* x, float* x1, float a, float c ){ for( int k = 0; k < 20; k++ ) { for( int i = 1; i < n; i++ ) { for( int j = 1; j < n; j++ ) { x[Index(i,j)] = ( x1[Index(i,j)] + a * ( x[Index(i-1,j)] + x[Index(i+1,j)] + x[Index(i,j-1)] + x[Index(i,j+1)] ) )/c; } } SetBoundary( b, x ); }} void Grid::Diffuse( int b, float* x, float* x1, float diff ){ float a = dt*diff*n*n; LinearSolve(b, x, x1, a, 1+4*a);} void Grid::Advect( int b, float* d, float* d0, float* u, float* v ){ int i0, j0, i1, j1; float x, y, s0, t0, s1, t1, dt0; dt0 = dt*n; for( int i = 1; i <= n; i++ ) { for( int j = 1; j <= n; j++ ) { x = i-dt0*u[Index(i,j)]; y = j-dt0*v[Index(i,j)]; if( x < 0.5f ) x = 0.5f; if( x > n+0.5f ) x = n + 0.5f; i0 = (int)x; i1 = i0 + 1; if( y < 0.5f ) y = 0.5f; if( y > n + 0.5f ) y = n + 0.5f; j0 = (int)y; j1 = j0+1; s1 = x-i0; s0 = 1-s1; t1 = y-j0; t0 = 1-t1; d[Index(i,j)] = s0 * (t0 * d0[Index(i0,j0)] + t1*d0[Index(i0,j1)]) + s1 * (t0 * d0[Index(i1,j0)] + t1*d0[Index(i1,j1)]); } } SetBoundary( b, d );} void Grid::Project( float* x, float* y, float* p, float* div ){ for( int i = 1; i <= n; i++ ) { for( int j = 1; j <= n; j++ ) { div[Index(i,j)] = -0.5f * (x[Index(i+1,j)] - x[Index(i-1,j)] + y[Index(i,j+1)] - y[Index(i,j-1)])/n; p[Index(i,j)] = 0; } } SetBoundary( 0, div ); SetBoundary( 0, p ); LinearSolve( 0, p, div, 1, 4 ); for( int i = 1; i <= n; i++ ) { for( int j = 1; j <= n; j++ ) { x[Index(i,j)] -= 0.5f * n * (p[Index(i+1,j)] - p[Index(i-1,j)]); y[Index(i,j)] -= 0.5f * n * (p[Index(i,j+1)] - p[Index(i,j-1)]); } } SetBoundary( 1, x ); SetBoundary( 2, y );} void Grid::SetBoundary( int b, float* x ){ for( int i = 1; i <= n; i++ ) { x[Index(0 ,i)] = b==1 ? -x[Index(1,i)] : x[Index(1,i)]; x[Index(n+1,i)] = b==1 ? -x[Index(n,i)] : x[Index(n,i)]; x[Index(i,0 )] = b==2 ? -x[Index(i,1)] : x[Index(i,1)]; x[Index(i,n+1)] = b==2 ? -x[Index(i,n)] : x[Index(i,n)]; } x[Index(0 ,0 )] = 0.5f*(x[Index(1,0 )] + x[Index(0 ,1)]); x[Index(0 ,n+1)] = 0.5f*(x[Index(1,n+1)] + x[Index(0 ,n)]); x[Index(n+1,0 )] = 0.5f*(x[Index(n,0 )] + x[Index(n+1,1)]); x[Index(n+1,n+1)] = 0.5f*(x[Index(n,n+1)] + x[Index(n+1,n)]);}void Grid::AddSource(float* x, float* source){ for( int i = 0; i < size; i++ ) { x += source*dt; }}void Grid::Reset(){ for( int i = 0; i < size; i++ ) vx = vy = vxOld = vyOld = dens = densOld = 0.0f;}.h
/* * Grid.h * Navier-Stokes * * Created by Jedd Haberstro on 5/22/08. * Copyright 2008 Jedd Haberstro. All rights reserved. * */#ifndef _GRID_H_#define _GRID_H_struct Grid{ Grid() : n(100), size((n+2)*(n+2)), h(1.0f/(float)n), dt(1.0f/60.0f), diffusionRate(0.0f), viscosity(0.0f) { vx = new float[size]; vy = new float[size]; vxOld = new float[size]; vyOld = new float[size]; dens = new float[size]; densOld = new float[size]; } Grid(int n, float dt, int diff, int visc) { this->n = n; size = (n+2)*(n+2); h = 1.0f/(float)n; this->dt = dt; diffusionRate = diff; viscosity = visc; vx = new float[size]; vy = new float[size]; vxOld = new float[size]; vyOld = new float[size]; dens = new float[size]; densOld = new float[size]; } void SimStep(); void DrawDensity(float width, float height); void DrawVelocity(float width, float height); // Add a source of density for a single cell void AddDensity(float s, int i, int j); void AddForce(float fx, float fy, int i, int j); void AddGlobalForce(float fx, float fy); int ConvertPosToGrid(float mousePos, float size); void Reset(); private: inline int Index(int i, int j) { return (i)+(n+2)*(j); } //Gauss-Seidel relaxation void LinearSolve( int b, float* x, float* x1, float a, float c ); void Diffuse( int b, float* x, float* x1, float diff ); void Advect( int b, float* d, float* d0, float* u, float* v ); void Project( float* x, float* y, float* p, float* div ); void SetBoundary( int b, float* x ); void AddSource(float* x, float* source); private: float diffusionRate; // Rate of diffusion, when > 0 it density will diffuse to other cells float viscosity; // Viscosity of the fluid float dt; // The (constant) timestep int n; // 1D size of grid float h; // Physical length of each side of the grid int size; // Total number of gridcells in the grid //float* densSources; // Sources of density float* vy; // Array of the current x velocity at each cell float* vx; // Array of the current y velocity at each cell float* vxOld; // Array of the old x velocity at each cell float* vyOld; // Array of the old y velocity at each cell float* dens; // Array of the current density at each cell float* densOld; // Array of the old density at each cell};#endif /* _GRID_H_ */
Well, what is it doing now? Is it much the same thing? (I don't know, you see, how much of your previous picture was a result of the drawing function and how much was a result of the simulation function. ^^; )
I'm afraid that I don't think that I've worked with this particular system before (it looks like a different system to that which I used for my own water simulation recently), so I fear that I'll be of little use at the moment, but I do have a few more questions:
Are you sure that the forces and inputs that you're applying are not too high? Have you tried lower (perhaps even much lower) values? Are you sure that you're supposed to be using the same functions for both velocity and density?
Having looked at the code, are you sure that, in LinearSolve, the function in the innermost loop shouldn't use only x1 on the right? As it is, it looks as though already-transformed values (that is, values from the density map that you're writing into, densOld, if I'm not much mistaken) are being used in calculating the current value.
I'm afraid that I don't think that I've worked with this particular system before (it looks like a different system to that which I used for my own water simulation recently), so I fear that I'll be of little use at the moment, but I do have a few more questions:
Are you sure that the forces and inputs that you're applying are not too high? Have you tried lower (perhaps even much lower) values? Are you sure that you're supposed to be using the same functions for both velocity and density?
Having looked at the code, are you sure that, in LinearSolve, the function in the innermost loop shouldn't use only x1 on the right? As it is, it looks as though already-transformed values (that is, values from the density map that you're writing into, densOld, if I'm not much mistaken) are being used in calculating the current value.
Author
The density is much of the same thing. I add more density to it, it acts like liquid, but then it quickly grows to fill most of the screen, then eats its self from the inside out.
I'm sure that the force and density values I'm adding aren't to high. When rendering the velocity by itself, the velocity part is acting exactly how it's suppose to be.
Yes I am sure I should be using the same functions.
I'm sure LinearSolve is correct. The method is semi-implicit. Having look at the paper, it is described here:
"The basic idea behind our method is to find the densities which when diffused backward
in time yield the densities we started with. In code:
This is a linear system for the unknowns x[IX(i,j)]. We could build the matrix for this linear
system and then call a standard matrix inversion routine. However, this is overkill for this
problem because the matrix is very sparse: only very few of its elements are non-zero.
Consequently we can use a simpler iterative technique to invert the matrix. The simplest
iterative solver which works well in practice is Gauss-Seidel relaxation."
What would be causing the density to blow up and not dissipate??
[Edited by - bronxbomber92 on May 26, 2008 12:12:20 AM]
I'm sure that the force and density values I'm adding aren't to high. When rendering the velocity by itself, the velocity part is acting exactly how it's suppose to be.
Yes I am sure I should be using the same functions.
I'm sure LinearSolve is correct. The method is semi-implicit. Having look at the paper, it is described here:
"The basic idea behind our method is to find the densities which when diffused backward
in time yield the densities we started with. In code:
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)]); This is a linear system for the unknowns x[IX(i,j)]. We could build the matrix for this linear
system and then call a standard matrix inversion routine. However, this is overkill for this
problem because the matrix is very sparse: only very few of its elements are non-zero.
Consequently we can use a simpler iterative technique to invert the matrix. The simplest
iterative solver which works well in practice is Gauss-Seidel relaxation."
What would be causing the density to blow up and not dissipate??
[Edited by - bronxbomber92 on May 26, 2008 12:12:20 AM]
All right, I believe that I see how you came to your equation. However, I seem to get 5a instead of 1-4a for my value of 'c':
[edit] Argh, sorry about the broken h-scroll. >_< I've changed the box to use source instead of code tags, and have manually broken the long, ironically-begun "for the sake of convenience" line. [/edit]
Have I made a mistake somewhere?
Another thing occurs to me, spurred, I think, by your mention of the density blowing up and not dissipating - is your timestep perhaps too large? That, as I recall, can have some unfortunate consequences...
(If you're using a variable timestep, then perhaps consider building in caps to the amount of density that may be moved about, and/or caps in other calculations.)
[Edited by - Thaumaturge on May 26, 2008 11:24:04 AM]
[edit] Argh, sorry about the broken h-scroll. >_< I've changed the box to use source instead of code tags, and have manually broken the long, ironically-begun "for the sake of convenience" line. [/edit]
The initial function: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)]);I'm guessing that what's given as x0 above is what you have as x1 in your code.For the sake of convenience, I'll call x[IX(i-1,j)], x[IX(i+1,j)], x[IX(i,j-1)] and x[IX(i,j-1)] 'e', 'f', 'g' and 'h', x[IX(i,j)] will be 'x' and x0[IX(i,j)] will be 'x0', giving us:x0 = x - a*(e + f + g + h - 4x)x0 = x - ae - af - ag - ah + 4axx0 = 5ax - ae - af - ag - ah5ax = x0 + ae + af + ag + ah = x0 + a(e + f + g + h)x = (x0 + a(e + f + g + h))/(5a)Have I made a mistake somewhere?
Another thing occurs to me, spurred, I think, by your mention of the density blowing up and not dissipating - is your timestep perhaps too large? That, as I recall, can have some unfortunate consequences...
(If you're using a variable timestep, then perhaps consider building in caps to the amount of density that may be moved about, and/or caps in other calculations.)
[Edited by - Thaumaturge on May 26, 2008 11:24:04 AM]
Author
Hmm, changing it to 5a makes it utterly uaeless :/. Same thing goes with lowering the timestep from 0.1 to something like 1/30.
What type of fluid simulation have you made?
What type of fluid simulation have you made?
I used the method described in this lecture by YannL.
For the most part it worked well for me, I think. My implementation doesn't deal terribly nicely with sudden drop-offs in the underlying surface level, for example, and, since I was using a non-fixed timestep I ended up placing some limitations on the amounts that the system was allowed to transfer (which were partially successful) but it did ripple nicely.
If you're interested, the thread that I created on the topic when I was having trouble should be here. The screenshots are old, however - I think that the final result looked better (especially since it included colour and transparency ;)). I also managed to deal with one of the problems that I was having in that thread, of the simulation looking horrible when confined, again by limiting the amounts that the system was allowed to pass around, as I recall.
For the most part it worked well for me, I think. My implementation doesn't deal terribly nicely with sudden drop-offs in the underlying surface level, for example, and, since I was using a non-fixed timestep I ended up placing some limitations on the amounts that the system was allowed to transfer (which were partially successful) but it did ripple nicely.
If you're interested, the thread that I created on the topic when I was having trouble should be here. The screenshots are old, however - I think that the final result looked better (especially since it included colour and transparency ;)). I also managed to deal with one of the problems that I was having in that thread, of the simulation looking horrible when confined, again by limiting the amounts that the system was allowed to pass around, as I recall.
This topic is closed to new replies.
Advertisement
Popular Topics
Advertisement
