Odd behavior in fluid simulation when converting from my working shadertoy to unity compute shaders

Question

I am trying to replicate my working 2d fluid shadertoy in a compute shader in Unity, with the hopes of moving it to 3D soon. When I replicate the algorithm the same way, I get some very weird behavior (seen in this video I took). I have tried to debug everything I can think of, but I can't figure out why they aren't the same. I am visualizing the vector matrix in this capture (the same as pressing Space while viewing my shadertoy).

I created a pastebin with the code that I am using to perform the Navier-Stokes equations that drive the velocity matrix. The meat of the simulation boils down to:

float4 S(RWTexture2D<float4> target, uint2 id)
{
    return target[(id.xy + resolution)%resolution];
}

void Fluid(RWTexture2D<float4> target, uint2 id, float2 offset, float4 values, inout float2 velocity, inout float pressure, inout float divergence, inout float neighbors)
{
    float2 v = S(target, id.xy + offset);
    float4 s = S(target, id.xy + offset.xy - v);
    
    float2 o= normalize(offset);
    
    velocity += o * (s.w - values.w);
    
    pressure += s.w;
    
    divergence += dot(o, s.xy);
    
    ++neighbors;
}

void StepVelocity(uint3 id, RWTexture2D<float4> write, RWTexture2D<float4> read, bool addJet)
{
    //sample our current values, then sample the values from the cell our velocity is coming from
    float4 values = S(read, id.xy);
    values = S(read, id.xy - values.xy);
    
    float2 velocity = float2(0,0);
    float pressure = 0;
    float neighbors = 0;
    float divergence = 0;
    
    //sample neighboring cells
    Fluid(read, id.xy, float2(0, 1), values, velocity, pressure, divergence, neighbors);
    Fluid(read, id.xy, float2(0, -1), values, velocity, pressure, divergence, neighbors);
    Fluid(read, id.xy, float2(1, 0), values, velocity, pressure, divergence, neighbors);
    Fluid(read, id.xy, float2(-1, 0), values, velocity, pressure, divergence, neighbors);

    velocity = velocity / neighbors;
    divergence = divergence / neighbors;
    pressure = pressure/neighbors;
    
    values.w = pressure-divergence;
    values.xy -= velocity;
    
    if (addJet && distance(id.xy, float2(resolution / 2.0, resolution / 2.0)) < 10)
        values = float4(0, .25, 0, values.w);
    
    write[id.xy] = values;
}

It should be pretty straightforward, i did my best to comment the shadertoy excessively to make it easy to understand (they are the same code, setup for different environments). Here is the C# code that drives drives the simulation.

I know its an inconvenient request to ask someone to dive into my code, but I am not at all sure what I am doing wrong and its driving me crazy. I have the exact same algorithm working on shadertoy but it behaves very strangely in Unity compute shaders and I can't figure out why. Every time I "step" i make sure to switch the read/write textures as to cleanly move forward in the simulation and not interfere.

Any ideas/tips towards fixing my issue would be greatly appreciated.

Answer 1

There are a couple of things I changed to make this work:

First, the pixel coordinates need to be offset by 0.5. In the compute shader the coordinates are ints, but in ShaderToy the coordinates that are passed in come from a SV_Position input, which already has that offset applied.

SV_Position describes the pixel location. Available in all shaders to get the pixel center with a 0.5 offset.

The same is true for gl_FragCoord. By default it is offset by 0.5. This matters when you add velocities to positions. Without the offset you wouldn't have the same behaviour in all directions.

Second, I used samplers for the inputs, point filtering doesn't seem to work with this. Otherwise the code is basically the same and should be fairly easy to follow.

Here is how it looks:

FluidSimulation.cs:

public class FluidSimulation : MonoBehaviour
{
    private RenderTexture In;
    private RenderTexture Out;
    private RenderTexture DrawIn;
    private RenderTexture DrawOut;
    
    private int nthreads = 8;
    private int threadresolution => (resolution / nthreads);
    private int stepKernel;

    [Range(8, 1024)] public int resolution = 800;
    [Range(0, 50)] public int stepsPerFrame = 8;
    public ComputeShader Compute;
    public Material OutputMaterial;

    void Start() => Reset();   

    private void Reset()
    {
        In = CreateTexture(RenderTextureFormat.ARGBHalf);
        Out = CreateTexture(RenderTextureFormat.ARGBHalf);
        DrawIn = CreateTexture(RenderTextureFormat.ARGBHalf);
        DrawOut = CreateTexture(RenderTextureFormat.ARGBHalf);
        
        stepKernel = Compute.FindKernel("StepKernel");
        Compute.SetFloat("resolution", resolution);        
    }

    void Update()
    {
        for(int i = 0; i<stepsPerFrame; i++) Step();
    }

    void Step()
    {
        Compute.SetTexture(stepKernel, "In", In);
        Compute.SetTexture(stepKernel, "Out", Out);
        Compute.SetTexture(stepKernel, "DrawIn", DrawIn);
        Compute.SetTexture(stepKernel, "DrawOut", DrawOut);     

        Compute.Dispatch(stepKernel, threadresolution, threadresolution, 1);
        
        OutputMaterial.SetTexture("_MainTex", DrawOut);

        SwapTex(ref In, ref Out);
        SwapTex(ref DrawIn, ref DrawOut);        
    }

    protected RenderTexture CreateTexture(RenderTextureFormat format)
    {
        RenderTexture tex = new RenderTexture(resolution, resolution, 0, format);

        //IMPORTANT FOR GPU SHADERS, allows random access (like gpus will do)
        tex.enableRandomWrite = true;
        tex.dimension = UnityEngine.Rendering.TextureDimension.Tex2D;        
        tex.filterMode = FilterMode.Bilinear;
        tex.wrapMode = TextureWrapMode.Clamp;
        tex.useMipMap = false;
        tex.Create();

        return tex;
    }

    void SwapTex(ref RenderTexture In, ref RenderTexture Out)
    {
        RenderTexture tmp = In;
        In = Out;
        Out = tmp;       
    }
}

Compute Shader:

#pragma kernel StepKernel

float resolution;

Texture2D<float4> In;
SamplerState samplerIn;
RWTexture2D<float4> Out;

Texture2D<float4> DrawIn;
SamplerState samplerDrawIn;
RWTexture2D<float4> DrawOut;

float4 Sample(Texture2D<float4> t, SamplerState s, float2 coords) {
    return t.SampleLevel(samplerIn, coords / resolution, 0);
}

void Fluid(float2 coord, float2 offset, inout float2 velocity, inout float pressure, inout float divergence, inout float neighbors)
{
    // Sample buffer C, which samples B, which samples A, making our feedback loop    
    float4 s = Sample(In, samplerIn, coord + offset - Sample(In, samplerIn, coord + offset).xy);

    // gradient of pressure from the neighboring cell to ours
    float sampledPressure = s.w;    

    //add the velocity scaled by the pressure that its exerting
    velocity += offset * sampledPressure;

    // add pressure
    pressure += sampledPressure;

    // divergence of velocity
    divergence += dot(offset, s.xy);

    //increase number of neighbors sampled
    neighbors++;
}

float4 StepVelocity(float2 id) {

    //sample from the previous state    
    float4 values = Sample(In, samplerIn, id - Sample(In, samplerIn, id).xy);
    
    float2 velocity = float2(0, 0);
    float divergence = 0.;
    float pressure = 0., neighbors = 0.;

    Fluid(id.xy, float2( 0., 1.), velocity, pressure, divergence, neighbors);
    Fluid(id.xy, float2( 0.,-1.), velocity, pressure, divergence, neighbors);
    Fluid(id.xy, float2( 1., 0.), velocity, pressure, divergence, neighbors);
    Fluid(id.xy, float2(-1., 0.), velocity, pressure, divergence, neighbors);

    //average the samples
    velocity /= neighbors;
    divergence /= neighbors;
    pressure /= neighbors;

    //output pressure in w, velocity in xy
    values.w = pressure - divergence;
    values.xy -= velocity;
    
    float2 p1 = float2(.47, .2);
    if (length(id.xy - resolution * p1) < 10.) {
        values.xy = float2(0, .5);
    }
    float2 p2 = float2(.53, .8);
    if (length(id.xy - resolution * p2) < 10.) {
        values.xy = float2(0, -.5);
    }
    return values;
}

float4 StepFluid(float2 id) {

    for (int i = 0; i < 4; i++)
        id -= Sample(In, samplerIn, id).xy;

    float4 color = Sample(DrawIn, samplerDrawIn, id);

    float2 p1 = float2(.47, .2);
    if (length(id.xy - resolution * p1) < 10.) {
        color = float4(0, 1, 0, 1);
    }
    float2 p2 = float2(.53, .8);
    if (length(id.xy - resolution * p2) < 10.) {
        color = float4(1, 0, 0, 1);
    }
    color *= .999;
    return color;
}

[numthreads(8, 8, 1)]
void StepKernel (uint3 id : SV_DispatchThreadID)
{     
    float2 coord = float2(id.x + .5, id.y + .5);
    Out[id.xy] = StepVelocity(coord);
    DrawOut[id.xy] = StepFluid(coord);
}

One more thing I wanted to mention, You can look at the translated HLSL code of your ShaderToy shaders by clicking the analyze button on the bottom of ShaderToys editor window. I didn't know this existed, but it can be quite useful when trying to convert shaders from ShaderToy to Unity.