# Fluid Simulation

This topic is 2787 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

## The full package can be found here:http://servut.us/snaileri/fluid/

##### Share on other sites
Hi there, I've had the same trouble getting started on fluid simulations, but after reading a number of articles on the topic, I managed to get some very good simulations up and running. Although It has some bugs and problems, I found this article on SPH very helpful:

http://cs.clemson.edu/~bpelfre/sph_tutorial.pdf

And here's some code:

http://sacredsoftware.net/svn/misc/watertower2/prototypes/keith/brandonpelfrey_sph/brandonpelfrey.cc

Cheers,
Mike

##### Share on other sites
Thanks, I will check it out.

##### Share on other sites
Okay, the code you gave uses SPH, and while I though it would be similar to MPM, it's not.
I find it extremely difficult to read the MPM code as I don't know where the abbreviations come from.
I would be very thankful if someone with more skills on fluid simulations could look at the code and tell what 'particle.x', 'particle.u', 'particle.cx', 'particle.px[3]','particle.gx[3]' might be. Same with the nodes' variables: 'd', 'u', 'ax'...
It seems like x is a pair for y and u is a pair for v, so once I know what properties these variables store, it would be wise to convert the pairs into 2D vectors.

I saw that in the SPH code particles check their neighbors by looping through every other particle, which is quite a slow way. Line 178:
[source]// Now look at every other particle
int particles_size = particles.size();
for(int j = i + 1; j < particles_size; ++j) {
[/source]
I'm not sure how exactly the neighbor-check is done in MPM, but I can't find a similar line in that.
Can anyone find any similarities between the two codes?

##### Share on other sites
If you don't understand how fluids work in general, trying to make sense of source code is going to be a nightmare. To that end, I'd suggest a little light reading that I found useful when I was starting down the fluid-sim path:
http://www.cs.ubc.ca/~rbridson/fluidsimulation/ - These are a nice introduction to fluids and how to simulate them. Very clearly written and easy to read!
http://www.dgp.toronto.edu/people/stam/reality/Research/pdf/GDC03.pdf - Jos Stam's paper on video game smoke. Lots of code in the paper, but you'll have to move it over to some execution environment before it'll run.
http://www.dgp.toronto.edu/people/stam/reality/Research/pdf/smoke.pdf - Jos Stam, Ronald Fedkiw, and Henrik Wann Jensen's classic paper.

After reading these three a time or few, you should be able to teach somebody waiting in the checkout line with you more than they ever wanted to know about fluids. At that point, you should be much more able to look at the code and make sense of it.

If you've got a handle on all that already, then I apologize! However, I still suggest that you read these: they do a very good job of putting code with theory.

Once you've done that, you can get into fancier techniques. That's where the fun is, so stick with it!

- jouley

##### Share on other sites
Thank you!

EDIT:
So, there's a book and then there's a PDF called 'fluids_notes.pdf'.
Is the PDF sufficient enough or do you recommend that I get the book? Did you learn the stuff from the book or the PDF?

##### Share on other sites

So, there's a book and then there's a PDF called 'fluids_notes.pdf'.
Is the PDF sufficient enough or do you recommend that I get the book? Did you learn the stuff from the book or the PDF?

I found the PDF sufficient for my purposes. I haven't looked at the book, so I can't say how much useful material it adds.

Best of luck!

##### Share on other sites
Ok... here goes:

To solve the Navier Stokes equations; which represent a continuum, one must find a way to discretize it with a mesh or a grid, in which case, it's no longer a continuum. A couple of Astrophysicists back in the 80's figured out a way to treat the fluid as a bunch of particles that behave in a fluid-like way. This is really nice because you don't have to mesh the thing. (Which is a pain.)

The key things to know when dealing with SPH is that each property of the fluid is held in the particles:
every particle has mass, viscosity, density (weird, I know), velocity, position and any other property which is represented in the equations. Each property is calculated based upon a weighting function.

The weighting function for every particle is a function of the difference in the distance from particle p, and every other particle, over a smoothing length.

g_p(x,y) ~= f(x',y', h)

// Steps in SPH:
// (0) Initialize all variables. Apply Boundary Conditions.
// (1) Solve for density for all particles
// (2) Compute pressure from an Equation of State
// (3) Calculate Viscosity from Energy Equation.
// (4) Take all calculated quantities and plug them into the
// momentum equation, get the Accelerations of the particles.
// (Get Force)
// (5) Use Acceleration in time integration to get positions and velocities.
// (6) Send to rendering pipeline then repeat steps 1-5.

I've annotated the java code where I recognize structures. This fellow didn't explain what they did very well at all. We learn a lesson from this: Comment your code, so that the next poor fool who reads it doesn't feel the need to wet themselves.

// forked from zlaper's Liquid simulation //Is this demo using MPM algorithm? //What are the node's m, d, gx, gy, u, v, ax and ay variables? //What are particle's x, y, u, v, dudx, dudy, dvdx, dvdy, cx, cy, px[3], py[3], gx[3] and gy[3] variables? //What are the phases in the liquid simulation? What happens and in what order? What are some of the fluid dynamics functions that are used in each chunk? //Where are the code chunks that control the viscosity, smoothing, elasticity and gravity? If they don't exist, how to add them? // This looks like SPH. // m is mass, d is density, gx and gy... weighting function. // u - velocity in x, v - velocity in y // ax - acceleration in x, ay - acceleration in y // dudx --> du/dx, dudy --> du/dy // cx --> speed of sound in x dim, cy --> speed of sound in y dim // px --> pressure in x, py --> pressure in y // Steps in SPH: // (0) Initialize all variables. Apply Boundary Conditions. // (1) Solve for density for all particles // (2) Compute pressure from an Equation of State // (3) Calculate Viscosity from Energy Equation. // (4) Take all calculated quantities and plug them into the // momentum equation, get the Accelerations of the particles. // (Get Force) // (5) Use Acceleration in time integration to get positions and velocities. // (6) Send to rendering pipeline then repeat steps 1-5. package com.jasonpeinko.liquid; import java.util.ArrayList; import com.badlogic.gdx.Gdx; import com.badlogic.gdx.Input; import com.badlogic.gdx.Input.Buttons; import com.badlogic.gdx.graphics.GL10; import com.badlogic.gdx.graphics.Texture; import com.badlogic.gdx.graphics.g2d.SpriteBatch; import com.badlogic.gdx.graphics.glutils.ImmediateModeRenderer; import com.badlogic.gdx.math.Vector2; public class Liquid { boolean mousePressed = false; int mx,my = 0; int mxOld,myOld = 0; Vector2 mVec; private static final float TICK = 1/60.0f; private ArrayList<particle> particles = new ArrayList<particle>(); private int gsizeX = 160; //200 private int gsizeY = 150; //150 private float accum = 0.0f; private Node[][] grid = new Node[gsizeX][gsizeY]; private ArrayList<node> active = new ArrayList<node>(); private Material water = new Material(1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f); //Only the first parameter is used - Material.m ImmediateModeRenderer renderer; Texture texture; Texture texture2; SpriteBatch batch; public Liquid() { mVec = new Vector2(0,0); init(); } public void init() { int i, j; for (i = 0; i < gsizeX; i++) { //The grid is created for (j = 0; j < gsizeY; j++) { this.grid[j] = new Node(); } } Particle p; for (i = 0; i < 100; i++) { //10,000 Particles are spawned for (j = 0; j < 100; j++) { p = new Particle(water, (i) + 4f, (j) + 4f, 0.0f, 0.0f); particles.add(p); } } //System.out.println("Particles: " + particles.size()); //Print the quantity of particles renderer = new ImmediateModeRenderer(); texture = new Texture(Gdx.files.internal("data/water.png")); texture2 = new Texture(Gdx.files.internal("data/node.png")); batch = new SpriteBatch(); } public void render() { //Input mxOld = mx; //The MousePos on the previous frame myOld = my; mx = Gdx.input.getX()/4; //The current MousePos my = Gdx.input.getY()/4; mVec.set(mx-mxOld, my-myOld); //System.out.println(mVec.x + "," + mVec.y); //Print Force vector if (Gdx.input.isButtonPressed(Buttons.LEFT)) { mousePressed = true; } else mousePressed = false; accum += Gdx.graphics.getDeltaTime(); while(accum >= TICK) { accum -= TICK; simulate(TICK); } batch.begin(); //We render stuff for (int i = 0 ; i<grid.length; i++)="" {="" we="" don't="" render="" just="" particles,="" but="" part="" of="" the="" grid="" as="" well="" for="" (int="" j="0" ;="" j<grid.length;="" j++)="" if="" (grid[j].active)="" only="" active="" nodes="" are="" drawn="" batch.draw(texture2,="" i*4,="" j*4);="" }="" (particle="" p="" :="" particles)="" batch.draw(texture,="" p.x*4,="" p.y*4);="" particles="" and="" grid-nodes="" will="" be="" with="" a="" multiplier="" 4.0x,="" since="" simulation="" is="" so="" small="" batch.end();="" mumbo="" jumbo="" begins="" here="" public="" void="" simulate(float="" dt)="" delta="" not="" being="" used="" some="" reason="" i="0;" <="" gsizex;="" this="" to="" reset="" mouse="" affection="" status="" gsizey;="" grid[j].affectedbymouse="false;" (node="" n="" active)="" at="" beginning="" step,="" every="" grid-cell="" zeroed="" n.m="(n.d" =="" n.gx="n.gy" n.u="n.v" n.ax="n.ay" 0.0f);="" n.active="false;" active.clear();="" g="mx-10;" g<mx+10;g++)="" check="" what="" affected="" by="" h="my-10;" h<my+10;h++)="" (g<0="" ||="">=gsizeX){ continue; } if (h<0 || h>=gsizeY){ continue; } grid[g][149-h].affectedByMouse = true; } } // H-L Comment // This is comparing each particle with all of the other particles... actually this is inefficient. // It scales as O(N^2). He/She should have used a linked list. //H-L --> // Pressure is calculated. He/She got it from a model. for (Particle p : particles) { //What happens in this loop? p.cx = (int) (p.x - 0.5F); //We check on which grid-cell (node) the particle is standing on p.cy = (int) (p.y - 0.5F); float x = p.cx - p.x; p.px[0] = (0.5F*x*x + 1.5F*x + 1.125F); p.gx[0] = (x + 1.5F); x += 1.0F; p.px[1] = (-x*x + 0.75F); p.gx[1] = (-2.0F*x); x += 1.0F; p.px[2] = (0.5F*x*x - 1.5F*x + 1.125F); p.gx[2] = (x - 1.5F); float y = p.cy - p.y; p.py[0] = (0.5F*y*y + 1.5F*y + 1.125F); p.gy[0] = (y + 1.5F); y += 1.0F; p.py[1] = (-y*y + 0.75F); p.gy[1] = (-2.0F*y); y += 1.0F; p.py[2] = (0.5F*y*y - 1.5F*y + 1.125F); p.gy[2] = (y - 1.5F); // Wow... ok... speed of sound is dynamically calculated. for (int i = 0; i < 3; i++) { for (int j = 0; j < 3; j++) { int cxi = p.cx + i; int cyj = p.cy + j; Node n = grid[cxi][cyj]; if (!n.active) { active.add(n); n.active = true; } float phi = p.px*p.py[j]; n.m += phi*water.m; n.d += phi; float dx = p.gx*p.py[j]; float dy = p.px*p.gy[j]; n.gx += dx; n.gy += dy; } } } for (Particle p : particles) { //What happens in this loop? int cx = (int) p.x; int cy = (int) p.y; int cxi = cx + 1; //Area on a cell? int cyi = cy + 1; float p00 = grid[cx][cy].d; //The node's, on which the particle is standing, data is gathered float x00 = grid[cx][cy].gx; float y00 = grid[cx][cy].gy; float p01 = grid[cx][cyi].d; //Corner of a node? float x01 = grid[cx][cyi].gx; float y01 = grid[cx][cyi].gy; float p10 = grid[cxi][cy].d; //Corner of a node? float x10 = grid[cxi][cy].gx; float y10 = grid[cxi][cy].gy; float p11 = grid[cxi][cyi].d; //Corner of a node? float x11 = grid[cxi][cyi].gx; float y11 = grid[cxi][cyi].gy; float pdx = p10 - p00; float pdy = p01 - p00; float C20 = 3.0F*pdx - x10 - 2.0F*x00; float C02 = 3.0F*pdy - y01 - 2.0F*y00; float C30 = -2.0F*pdx + x10 + x00; float C03 = -2.0F*pdy + y01 + y00; float csum1 = p00 + y00 + C02 + C03; float csum2 = p00 + x00 + C20 + C30; float C21 = 3.0F*p11 - 2.0F*x01 - x11 - 3.0F*csum1 - C20; float C31 = -2.0F*p11 + x01 + x11 + 2.0F*csum1 - C30; float C12 = 3.0F*p11 - 2.0F*y10 - y11 - 3.0F*csum2 - C02; float C13 = -2.0F*p11 + y10 + y11 + 2.0F*csum2 - C03; float C11 = x01 - C13 - C12 - x00; float u = p.x - cx; float u2 = u*u; float u3 = u*u2; float v = p.y - cy; float v2 = v*v; float v3 = v*v2; float density = p00 + x00*u + y00*v + C20*u2 + C02*v2 + C30*u3 + C03*v3 + C21*u2*v + C31*u3*v + C12* u*v2 + C13*u*v3 + C11*u*v; //DENSITY MULTIPLIER: //Higher value = lower density density = density*1.0f; float pressure = density - 1.0F; if (pressure > 2.0F) { pressure = 2.0F; } float fx = 0.0F; float fy = 0.0F; if (p.x < 4.0F) fx += water.m*(4.0F - p.x); else if (p.x > gsizeX - 5) { fx += water.m*(gsizeX - 5 - p.x); } if (p.y < 4.0F) fy += water.m*(4.0F - p.y); else if (p.y > gsizeY - 5) { fy += water.m*(gsizeY - 5 - p.y); } for (int i = 0; i < 3; i++) { for (int j = 0; j < 3; j++) { Node n = grid[(p.cx + i)][(p.cy + j)]; float phi = p.px*p.py[j]; float gx = p.gx*p.py[j]; float gy = p.px*p.gy[j]; // H-L Poorly labled Force: n.ax += (-(gx*pressure) + fx*phi)*1.4f; //Some interesting parameter of the fluid n.ay += (-(gy*pressure) + fy*phi)*1.4f; } } } // H-L: Nothing magical... finding acceleration. for (Node n : active) { //What happens in this loop? if (n.m > 0.0F) { n.ax /= n.m; n.ay /= n.m; n.ay += 0.03F; //Gravity? >0.06f is negative gravity? } } // H-L: Computes new velocity from acceleration. for (Particle p : particles) { //What happens in this loop? for (int i = 0; i < 3; i++) { for (int j = 0; j < 3; j++) { Node n = grid[(p.cx + i)][(p.cy + j)]; float phi = p.px*p.py[j]; p.u += phi*n.ax; p.v += phi*n.ay; } } p.v += -0.06f; //-0.06f float mu = water.m*p.u; float mv = water.m*p.v; for (int i = 0; i < 3; i++) { for (int j = 0; j < 3; j++) { Node n = grid[(p.cx + i)][(p.cy + j)]; float phi = p.px*p.py[j]; //Viscosity? n.u += phi*mu*1.0f; //Try multipliers between 0.8 and 1.01 n.v += phi*mv*1.0f; } } } // Computes velocity by dividing momentum by mass. for (Node n : active) { //What happens in this loop? if (n.m > 0.0F) { n.u /= n.m; n.v /= n.m; } } // Positions are calculated then sent to rendering // pipeline. for (Particle p : particles) { //What happens in this loop? float gu = 0.0F; float gv = 0.0F; for (int i = 0; i < 3; i++) { for (int j = 0; j < 3; j++) { Node n = grid[(p.cx + i)][(p.cy + j)]; float phi = p.px*p.py[j]; gu += phi*n.u; gv += phi*n.v; //Mouse interaction if (n.affectedByMouse &amp;&amp; mousePressed) { gu +=(0.07f*mVec.x); gv -=(0.07f*mVec.y); } } } p.x += gu; p.y += gv; p.u += 1.0f*(gu - p.u); //Surface Tension? p.v += 1.0f*(gv - p.v); if (p.x < 1.0F) { p.x = (1.0F + (float) Math.random()*0.01F); p.u = 0.0F; } else if (p.x > gsizeX - 2) { p.x = (gsizeX - 2 - (float) Math.random()*0.01F); p.u = 0.0F; } if (p.y < 1.0F) { p.y = (1.0F + (float) Math.random()*0.01F); p.v = 0.0F; } else if (p.y > gsizeY - 2) { p.y = (gsizeY - 2 - (float) Math.random()*0.01F); p.v = 0.0F; } } //end of simulation } } 

You're brave for taking this on. That should be respected.

This is one of the best books that I've ever read on this subject. Liu is one of clearest, most concise authors I've ever read.

##### Share on other sites
That code has no comments because it was decompiled from my fluid simulator. My own code is in fact very well commented. He also decided to place a license in it which has been confusing for a lot of people who have assumed it was open source. That being said, I am now working with somebody on a more official and well documented version that will be legally open source. The open source part aside, the more deadly part has actually been that he has introduced a few bugs and magic numbers into the simulation, for whatever reason, such as applying gravity twice first with .03f and then .06f. Or that multiplying of grid force by 1.4f. Most importantly, it has led to a string of buggy implementations that tarnish the awesomeness of MPM and keep it second to SPH.

MPM should not be regarded as elf magic and in reality is simpler to implement and faster than SPH.

##### Share on other sites
I see. I apologize, I didn't realize that it was decompiled. I was under the impression that I was look at the original source code. My comments were made based upon what I had in front of me.

It's not elf magic; you're right, it's merely a useful way to handle PDEs.

• ### Game Developer Survey

We are looking for qualified game developers to participate in a 10-minute online survey. Qualified participants will be offered a \$15 incentive for your time and insights. Click here to start!

• 19
• 36
• 9
• 16
• 75