Solution for compressibility problem in Eulerian method

As I mentioned before, FLIP solver suffers from compressibility problem. In fact, all the grid-based method might be suffering from comressibility problem.

You might consider me to be naive, but I'll try to convince you:

The grid based solver calculate the velocity field based on two formula:

1st is the Navier Stokes equation, 2nd is the incompressible constraint.

For grid solver, all the calculations are done based on one simple assumption:
all the fluid cells have the same density as the rest density.

And this is the key why these solvers are compressible. Though, it's just a small difference, which is hard to tell visually. This problem would become apparent in FLIP solver: with the same particle input, using different grid resolution would lead to different result, some shrink the volume(high resolution), some increase the volume(low resolution).

I tried to associate a density coefficient for each cell for the pressure solve, but it does not contribute. If you combine the rhow and p term in the 1st equation as a whole, you'll find out although the density coefficient would affect the value for pressure, when calculating back to velocity field, this term would be eliminated.

So I realize the problem lies in the second formula. There should be something more than just velocity field to conserve the volume.

It turns out I'm right. Well it's a pity that I found out someone has already done research into this before, or maybe I could be the first one.
With a simple search I found 2 papers dealing with this problem:
This one by Nuttapong Chentanez and Matthias Muller,
and this one by Michael Lentine, Mridul Aanjaneya and Ronald Fedkiw.

What is interesting is that what I've been doing always match these guys research. The papers I've been referring are always done by these names.

So next step: reading these papers and integrate into the FLIP solver. This might be useful for SPH as well, but I need to take a further look.