Flip solver update

for the past few days I've been working on this solver.

I imported the ghost sph sampler and anisotropic kernel part for surface reconstruction. This time everything was done with tbb(Thread Building Blocks), parallelizing all the particles. The performance is really satisfying: for 150K particles, the mesh reconstruction took less than 2 seconds per frame, including neighbor search, iterative matrix decomposition and color field gathering.

Compare to my previous implementation of this part, due to everything is done on the fly, the memory deallocation cost reduced drastically, and that is one of the major reason for boosting the performance.


However, I discovered a huge problem for this solver.
It is COMPRESSIBLE!
The error is introduced by particle advection. Within each single frame, the projection is ensured to be incompressible, however, after particles being advected, the number for fluid cells would change, and that is the reason why it is compressible.

With the 150k particle configuration, using 50*50*50 grid would lead to volume increasing, while 100*100*100 would shrink volume into a thin sheet.


In fact, setting the divergence to be free would not be sufficient for FLIP solver. Because the density(or mass) is carried with the particles, while for the grid solving part, there's no way in ensuring the density of each cell to be constant.

This problem could be alleviated with a certain grid size for a certain particle distribution. Yet only alleviation is possible, cause the particle and grid are kind of coupled in terms of density.

Right now I had an idea which might be useful for decoupling the two, but I need to take a further look into the physics and math.

I'll discuss with the author of FLIP paper, and hopefully I can find a better way for it soon.

FLIP solver done!

Finally finished the FLIP solver.

Right now it's a first version, so still, lots of details need to be improved.

A quick demo for the new solver is here.

Well I'll keep on working on this project and come up with way better polished demo.

Fast sweeping

Still working on the FLIP solver.
The fast sweeping algorithm is significant for FLIP solver, in generating levelset data and extending velocity field.

Based on Hongkai Zhao's paper I implemented the fast sweeping algorithm for the FLIP solver.
However proving part for the paper was using 2d example. So I extended the proving part to 3D, and here's the scanned image of the proving(also the pseudo code for the core part is included).

Compared to the implementation I did during the summer, the current one would be more abstract because of I'm using the result of my own proving directly. As a result of which, the operation would be faster, because the number of comparison and numerical calculation are minimized.

As I always says, math and physics are real science compared to computer science. In my biased personal opinion, a huge difference between science and engineering is: science is continuous and engineering is discretized.

Thanks math and physics for always providing guidance for everything.
I really enjoy this kind of life, surrounded by science. Maybe I should go for a physics phd later. Seriously I'm not kidding.

New solver for the fluid simulation project

SPH suffers from severe compressible problem. A combination of WCSPH( http://cg.informatik.uni-freiburg.de/publications/2007_SCA_SPH.pdf ) and PCISPH( http://dl.acm.org/citation.cfm?id=1531346 ) could be a good improvement of SPH. 

However, it's far from enough. It could be a good solution for the inner particles. while for the boundary particles(either in contact with air or solid), the pressure gathering is incorrect, which will lead to an artifact.

Ghost sph(http://www.cs.ubc.ca/~rbridson/docs/schechter-siggraph2012-ghostsph.pdf) to me is the best solution for improving the quality of SPH. But the problem is ghost sph is kind of expensive. The number of sample particles for the solids is usually much more than the fluid particles. And re-sampling ghost particles is not trivial task.

In order to solve the compressibility problem and maintain the details that particles bring about, I decided to turn to FLIP(http://www.cs.ubc.ca/~rbridson/docs/zhu-siggraph05-sandfluid.pdf) for help.

FLIP is a combination of Lagrangian method and Eulerian method. Particles are used for advection, which eliminates the numerical dissipation caused by Eulerian methods, and MACGrid is used for solving the Poisson equation, which solve the pressure distribution problem perfectly.

For the past few days I've been reading lots of papers related to Eulerian fluid simulation, including Robert Bridson's book: Fluid Simulation for Computer Graphics. And finally I got a good understanding of each steps that is necessary for FLIP solver.

I'll list my understanding of the core steps in FLIP solver here:

1. transfer velocity from particles to grid. The grid is only used for solving the pressure distribution, so only the velocity field is needed. The way used for transferring is splatting each particle's velocity onto the grid using tri-linear weighting.

2. generate a level set from the particles. The level set is necessary for later use. 

3. extend the velocity field to the whole grid. Before this step, only the grid cells that intersect the particles have a non-zero velocity. In order to get the correct velocity distribution, we have to extrapolate the velocity to the whole grid based on one simple principal: the dot product of the gradient of velocity and the gradient of level set should be zero. Which means the extrapolated velocity should not change in the normal direction of fluid surface.

4. solve for Poisson equation. I haven't start coding with this part. Yet from the previous smoke simulation project, this should be similar. The key is setting boundary condition, using ghost pressure and apply pre-conditioner.

5. extend the velocity again. In the previous step, only the velocity field of the fluid cells have been updated, so in order to get a correct value in interpolating back to particles, we need to update the rest part of the grid.

6. transfer back to the particles. Update particles' velocity based on the new velocity field.

That's basically my understanding for the FLIP solver from reading during the past few days. I'm starting to write this solver part and integrate that for my fluid simulation project.