I'm having trouble understanding the algorithm behind liquid simulations.

I've found a source-code and therefore seen what the code looks like inside the simulation, but I still don't understand how the lines are related to each function that describes the motion of a fluid in real world. Of course, I can look at how each variable is created in the simulation, but I don't know why the calculations are what they are, as they don't have any resemblance to the real world fluid functions. And since the particles' and grid's variables have so strange abbreviations, the code is just some mumbo jumbo for me.

I've read sevelar articles and small e-books about fluid simulations, but none of them really describes how the theory is converted into a working code.

Some of the papers that I've read are:

An Informal Tutorial on Smoothed Particle Hydrodynamics for Interactive Simulation - Brandon Pelfrey

Lagrangian Fluid Dynamics Using Smoothed Particle Hydrodynamics - Micky Kelager

Particle-based Viscoelastic Fluid Simulation - Simon Clavet

Right now I'm trying to write a 2D fluid simulation using Material Point Method (MPM algorithm), which is very similar to SPH (or that's what I've heard, at least).

I have figured out that the MPM uses only two kinds of elements in the simulation: fluid particles and grid cells.

The fluid particles are the visual part of the simulation and the grid is used for storing the data about density, pressure (and other physical quantities?) at any given point in space, and that data is gathered and calculated from the nearby fluid particles. The data from the grid cells is then used by the particles later on in the simulation step to calculate their new velocity and direction.

That's basically all I know about it right now. The next step is to understand how each line of code is related to the functions that describe a fluid's behaviour in the real world. I want to know where are the code chunks that control the density, pressure, viscosity, elasticity, smoothing, gravity and stiffness of the fluid - or if some of those parts do even exist in the source-code I've found.

I find it strange that first of all, very few people are willing to share the source-code of their fluid simulations, and secondly, nobody on the Internet is willing to write a comprehensive tutorial on this subject. I haven't found any source-code that has comments on the lines.

If I ever learn the logic behind fluid simulations, I surely will do a video tutorial where I teach everything about how fluids behave and I will describe every line of code and tell why it's written the way it is.

So, to put it shortly: I'm asking you to comment the code and describe what happens in each chunk of the code.

Here's the 200-line code of the fluid simulation step that I'm studying:
http://servut.us/sna...uid/Liquid.java

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

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

Thanks, I will check it out.

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?

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

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?

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!

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.

http://books.google....id=_cwFMmEQvZQC</grid.length;></node></node></particle></particle>

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.