pressure_solve.comp

// Pressure (PPE) solving iteration
// we're solving a pseudo pressure which is pseudopressure = (density / deltaTime * pressure)

#version 450

#include "../per_frame_resources.glsl"
#include "../utilities.glsl"
#include "hybrid_fluid.glsl"

// layout(set = 2, binding = 0) uniform texture3D VelocityVolumeX;
// layout(set = 2, binding = 1) uniform texture3D VelocityVolumeY;
// layout(set = 2, binding = 2) uniform texture3D VelocityVolumeZ;
layout(set = 2, binding = 3) uniform utexture3D MarkerVolume;
layout(set = 2, binding = 4) uniform texture3D DivergenceVolume;
layout(set = 2, binding = 5) uniform texture3D PressureVolumeRead;
layout(set = 2, binding = 6, r32f) uniform restrict image3D PressureVolumeWrite;

COMPUTE_PASS_VOLUME

void main() {
    ivec3 gridCoord = ivec3(gl_GlobalInvocationID);

    // Inside...
    // * Air cells: Pressure is zero (could be any constant value really)
    // * Solid cells: Not defined!
    //      But for easier derivation we can just set imagine to set it to a value that would cause the velocity at the fluid/solid interface to
    //      fulfill the boundary condition. Note that this means though that we have different pressure values for every boundary face. For details
    //      check out chapter 4.1 in Bridson's SIGGRAPH2007 course notes https://www.cs.ubc.ca/~rbridson/fluidsimulation/fluids_notes.pdf

    // Boundary conditions:
    // Air: Dirichlet (-> pressure is set to zero)
    // Solid: Neumann (-> velocity of fluid is boundary velocity -> pressure gradient is negative boundary velocity)

    // Is this a fluid cell?
    if (texelFetch(MarkerVolume, gridCoord, 0).x != CELL_FLUID) {
        return;
    }

    uvec3 maxGridCoord = textureSize(MarkerVolume, 0) - uvec3(1);

    // jacobi iteration
    float divergence = texelFetch(DivergenceVolume, gridCoord, 0).x;
    uint markerX0 = texelFetch(MarkerVolume, gridCoord - ivec3(1, 0, 0), 0).x;
    uint markerX1 = texelFetch(MarkerVolume, gridCoord + ivec3(1, 0, 0), 0).x;
    uint markerY0 = texelFetch(MarkerVolume, gridCoord - ivec3(0, 1, 0), 0).x;
    uint markerY1 = texelFetch(MarkerVolume, gridCoord + ivec3(0, 1, 0), 0).x;
    uint markerZ0 = texelFetch(MarkerVolume, gridCoord - ivec3(0, 0, 1), 0).x;
    uint markerZ1 = texelFetch(MarkerVolume, gridCoord + ivec3(0, 0, 1), 0).x;

    float numNonSolid = 0.0;
    numNonSolid += float(markerX0 != CELL_SOLID);
    numNonSolid += float(markerX1 != CELL_SOLID);
    numNonSolid += float(markerY0 != CELL_SOLID);
    numNonSolid += float(markerY1 != CELL_SOLID);
    numNonSolid += float(markerZ0 != CELL_SOLID);
    numNonSolid += float(markerZ1 != CELL_SOLID);
    if (numNonSolid == 0.0)
        return;

    float pressureX0 = markerX0 == CELL_FLUID ? texelFetch(PressureVolumeRead, gridCoord - ivec3(1, 0, 0), 0).x : 0.0;
    float pressureX1 = markerX1 == CELL_FLUID ? texelFetch(PressureVolumeRead, gridCoord + ivec3(1, 0, 0), 0).x : 0.0;
    float pressureY0 = markerY0 == CELL_FLUID ? texelFetch(PressureVolumeRead, gridCoord - ivec3(0, 1, 0), 0).x : 0.0;
    float pressureY1 = markerY1 == CELL_FLUID ? texelFetch(PressureVolumeRead, gridCoord + ivec3(0, 1, 0), 0).x : 0.0;
    float pressureZ0 = markerZ0 == CELL_FLUID ? texelFetch(PressureVolumeRead, gridCoord - ivec3(0, 0, 1), 0).x : 0.0;
    float pressureZ1 = markerZ1 == CELL_FLUID ? texelFetch(PressureVolumeRead, gridCoord + ivec3(0, 0, 1), 0).x : 0.0;

    float newPressure = (pressureX0 + pressureX1 + pressureY0 + pressureY1 + pressureZ0 + pressureZ1 - divergence) / numNonSolid;

    imageStore(PressureVolumeWrite, gridCoord, vec4(newPressure));
}