Hi,

I've been trying to get a basic fluid simulation up and running but I've run into a dead end: for several days now I've been trying to find where I've gone wrong and I'm still looking. I started with a 3D Eulerian grid based simulation but debugging that was a nightmare so I've now coded up a 2D version and I can see the issue clearer now (at least I think so). The simulation seems to run fine for the first couple of frames but after that it just seems to stop: I mean the velocity doesn't spread through the fluid (you can see small changes within the area that was initially affected). From what I can tell, the divergence I calculate to find pressure goes to 0 after a few frames (~10). Thus the pressure will also be 0 and then the projection in the end just copies previous frames velocity.

Here are a couple of screenshots of different properties: (not scale-corrected colors, e.g. there are negative values for velocity)

velocity: http://postimg.org/image/th0fraz5d/

divergence: http://postimg.org/image/am7wseu6r/

pressure is just a mirrored image of divergence (at least so it looks like, has a bigger area and disappears a bit slower though)

color property (or "dye" if you will): http://postimg.org/image/yun01ufzd/

The color/velocity values were written using a simple gaussian splat, with white as color value and (10.0f, 10.0f, 0.0f) as velocity values. Thus the direction seems about correct to me.

Since I'm doing the simulation on the GPU, there's a lot of code involved with ping-ponging render targets etc but none of my debugging tools show any issues with the pipeline and all data buffers seem to contain the correct data at each step. So, I'm left wondering if my maths is correct: I've consulted both GPU gems books that had articles about this and checked my versions against theirs and I simply can't find any differences, yet the simulation isn't working as intended.

For now, I'll just post my advection for color/velocity and divergence routines as (to me) it seems to be going wrong already at one of these steps.

input.position is the pixel position SV_Position. Also, helper func used in both cases (texDim is simulation grid width/height which currently matches the viewport):

float2 PositionToTexCoord(float2 position) {
return float2(position.x / texDim.x, position.y / texDim.y);
}

Advection:

float4 PSMain(PSInput input) : SV_TARGET {
float2 center = PositionToTexCoord(input.position.xy);
float2 pos = input.position.xy - dt * velocity.Sample(pointSampler, center).rg;


float2 samplingPoint = PositionToTexCoord(pos);


return qtyToAdvect.Sample(linearSampler, samplingPoint);
}

float PSMain(PSInput input) : SV_TARGET { 
float2 left = PositionToTexCoord(input.position.xy - float2(1.0f, 0.0f));
float2 right = PositionToTexCoord(input.position.xy + float2(1.0f, 0.0f));
float2 bottom = PositionToTexCoord(input.position.xy + float2(0.0f, 1.0f));
float2 top = PositionToTexCoord(input.position.xy - float2(0.0f, 1.0f));
 
float4 l = velocity.Sample(pointSampler, left);
float4 r = velocity.Sample(pointSampler, right);
float4 b = velocity.Sample(pointSampler, bottom);
float4 t = velocity.Sample(pointSampler, top);
 
float div = 0.5f * ((r.x - l.x) + (t.y - b.y));
 
return div;
}

float2 bottom = PositionToTexCoord(input.position.xy + float2(0.0f, 1.0f));
float2 top = PositionToTexCoord(input.position.xy - float2(0.0f, 1.0f));

and you are doing

float div = 0.5f * ((r.x - l.x) + (t.y - b.y));

try swapping the +/- where you define bottom and top

float2 bottom = PositionToTexCoord(input.position.xy + float2(0.0f, 1.0f));
float2 top = PositionToTexCoord(input.position.xy - float2(0.0f, 1.0f));

and you are doing

float div = 0.5f * ((r.x - l.x) + (t.y - b.y));

try swapping the +/- where you define bottom and top

YES! Thank you so much! I guess debugging something for several hours has a tendency to turn out like this >.<

The result: http://postimg.org/image/v9kgv9m2b/

To be honest though, doesn't the SV_Position semantic go from [0, screenWidth] x [0, screenHeight] from the top left corner? If so, if there's a pixel at (10,10), the bottom one would be (10,11) and top (10,9), wouldn't it?

Thanks again!

To be honest though, doesn't the SV_Position semantic go from [0, screenWidth] x [0, screenHeight] from the top left corner? If so, if there's a pixel at (10,10), the bottom one would be (10,11) and top (10,9), wouldn't it?

I know that OpenGL is flipped upside down (like TGA) where the bottom left is (0,0), I'm not sure about D3D though. In reality, it doesn't matter if you call it top/bottom head/toe, etc.. What matters is that the 'greater' value has the 'lesser' value subtracted - regardless of the orientation (semantics). In this case, you have the coordinate with (0,1) added, and (0,-1) added, to create two coordinates, and all that was happening here was you were subtracting the larger coordinate from the lesser coordinate, instead of the other way around.

To be honest though, doesn't the SV_Position semantic go from [0, screenWidth] x [0, screenHeight] from the top left corner? If so, if there's a pixel at (10,10), the bottom one would be (10,11) and top (10,9), wouldn't it?

I know that OpenGL is flipped upside down (like TGA) where the bottom left is (0,0), I'm not sure about D3D though. In reality, it doesn't matter if you call it top/bottom head/toe, etc.. What matters is that the 'greater' value has the 'lesser' value subtracted - regardless of the orientation (semantics). In this case, you have the coordinate with (0,1) added, and (0,-1) added, to create two coordinates, and all that was happening here was you were subtracting the larger coordinate from the lesser coordinate, instead of the other way around.

Ah, I see now. Thanks!

<snip>

The problem I had in 3D was related to the way I was setting up the grid... Funnily enough, I created a topic in the D3D section about rendering a quad with a single pixel border and I thought everything was fine. While the primitives were rendering the way they were supposed to, ie a quad with a 1px border and lines being 1px thick, I think the pixel to texel relationship was wrong. As soon as I copied the grid coordinate creation code from GPU Gems 3, I got a similar result to what happened in 2D with just the advection, ie all seems good. Could someone explain to me what on Earth is going on with the coordinate creation maths?

For the main quad:

const float w = static_cast<float>(width_);
const float h = static_cast<float>(height_); 
 
FluidFx::VS_INPUT tempVertices[4];
 
tempVertices[0].position = D3DXVECTOR3(1.0f * 2.0f / w - 1.0f, -1.0f * 2.0f / h + 1.0f, 0.0f);
tempVertices[0].cellCoord = D3DXVECTOR3(1.0f, 1.0f, z);
 
tempVertices[1].position = D3DXVECTOR3((w - 1.0f) * 2.0f / w - 1.0f, -1.0f * 2.0f / h + 1.0f, 0.0f);
tempVertices[1].cellCoord = D3DXVECTOR3(w - 1.0f, 1.0f, z);
 
tempVertices[2].position = D3DXVECTOR3((w - 1.0f) * 2.0f / w - 1.0f, -(h - 1.0f) * 2.0f / h + 1.0f, 0.0f);
tempVertices[2].cellCoord = D3DXVECTOR3(w - 1.0f, h - 1.0f, z);
 
tempVertices[3].position = D3DXVECTOR3(1.0f * 2.0f / w - 1.0f, -(h - 1.0f) * 2.0f / h + 1.0f, 0.0f);
tempVertices[3].cellCoord = D3DXVECTOR3(1.0f, h - 1.0f, z);
 
// construct a quad
vertices[vertIndex++] = tempVertices[0];
vertices[vertIndex++] = tempVertices[1];
vertices[vertIndex++] = tempVertices[2];
vertices[vertIndex++] = tempVertices[0];
vertices[vertIndex++] = tempVertices[2];
vertices[vertIndex++] = tempVertices[3];

For the lines:

int vertIndex = 0;
for(float z = 0.0f; z < d; z += 1.0f) {
SetupBoundaryLine(0.0f, 1.0f, w, 1.0f, z, vertices, vertIndex);
SetupBoundaryLine(0.0f, h, w, h, z, vertices, vertIndex);
SetupBoundaryLine(1.0f, 0.0f, 1.0f, h, z, vertices, vertIndex);
SetupBoundaryLine(w, 0.0f, w, h, z, vertices, vertIndex);
}

float w = static_cast<float>(width_);
float h = static_cast<float>(height_);
 
FluidFx::VS_INPUT tempVertex;
 
tempVertex.position = D3DXVECTOR3(x1 * 2.0f / w - 1.0f, -y1 * 2.0f / h + 1.0f, 0.5f);
tempVertex.cellCoord = D3DXVECTOR3(0.0f, 0.0f, z);
vertices[vertIndex++] = tempVertex;
 
tempVertex.position = D3DXVECTOR3(x2 * 2.0f / w - 1.0f, -y2 * 2.0f / h + 1.0f, 0.5f);
tempVertex.cellCoord = D3DXVECTOR3(0.0f, 0.0f, z);
vertices[vertIndex++] = tempVertex;