I'm using smoothed particle hydrodynamics for a real-time simulation intended for computer graphics application (rather than accuracy), which runs at a good frame rate with as much as 20,000 particles on a quad core Intel (I spent some time getting things to parallelize nicely). My problem is stability in a tall column of fluid--I'm using a weakly compressible fluid (pressure computed from density by Tait's equation) and if restricted to a tall and thin vessel resulting into many layers of particles, under gravity, the liquid bounces a lot. If the gravity is high, the simulation blows up. One solution is applying a stability criterion to determine the maximum allowable timestep, but reducing the timestep artificially makes the simulation non-realtime, so this is not an acceptable solution for this application.

I've read that a fully incompressible SPH simulation can be done by using a pressure projection step, which needs solving a Poission equation. This is not something I'm familiar with and I'm wondering which solving algorithm to use, and how to specify boundary conditions (the fluid is free-surface).