c# – Cheap 2D fluid simulation at low resolution

I have a 64×64 display and an accelerometer.

I want to run a 2D fluid simulation and take a reading from the accelerometer as the gravity vector.

I’m going to be doing this on fairly low-end hardware (64-bit quad-core Arm Cortex-A76 @ 2.4GHz and ~7GB usable RAM), without the benefit of GPU acceleration.

I’ve attempted to modify a Smoothed-particle hydrodynamics implementation, but the result I get is unstable and seems to expand forever.

Even then, I’m not sure it’s the best approach.

At a minimum, I’d like:

• A [mostly] incompressible liquid
• Have momentum
• Handle pressure/Be self-levelling around obstacles
• The ability to define obstacles
• Fixed volume of liquid in a closed square container

Sources and sinks are desirable, and vorticity is a stretch goal.

I considered using a grid-based approach but that means at a minimum iterating 64x64x30 FPS
= 122,880 cells/second which seems excessive if I’m doing any heaving lifting in that loop on an ARM processor.

Through an odd toolchain, I’m using C# for this project, and my SPH attempt is below

using System;

namespace FluidSharp {
    public class FluidSimulation {
        private readonly int gridSize;
        private readonly float diffusionRate;
        private readonly float viscosity;
        private float[] sourceDensity;
        private float[] density;

        private float[] velocityX;
        private float[] velocityY;
        private float[] previousVelocityX;
        private float[] previousVelocityY;

        public FluidSimulation(int size, float diffusion, float viscosity) {
            gridSize = size + 2; // Add some space for borders
            diffusionRate = diffusion;
            this.viscosity = viscosity;

            sourceDensity = new float[gridSize * gridSize];
            density = new float[gridSize * gridSize];

            velocityX = new float[gridSize * gridSize];
            velocityY = new float[gridSize * gridSize];
            previousVelocityX = new float[gridSize * gridSize];
            previousVelocityY = new float[gridSize * gridSize];
        }


        public void SetInitialState(float fractionOfTotalVolume) {
            fractionOfTotalVolume = Math.Clamp(fractionOfTotalVolume, 0f, 1f);

            int fillHeight = (int)((gridSize - 2) * fractionOfTotalVolume);

            float liquidDensity = 1.0f;

            for (int j = 1; j <= fillHeight; j++) {
                for (int i = 1; i < gridSize - 1; i++) {
                    density[i + (gridSize - j) * gridSize] = liquidDensity;
                }
            }
        }

        public void Update(float gravityX, float gravityY, float deltaTime) {
            AddGravity(gravityX, gravityY, deltaTime);
            Diffuse(1, previousVelocityX, velocityX, viscosity, deltaTime);
            Diffuse(2, previousVelocityY, velocityY, viscosity, deltaTime);
            Project(previousVelocityX, previousVelocityY, velocityX, velocityY);
            Advect(1, velocityX, previousVelocityX, previousVelocityX, previousVelocityY, deltaTime);
            Advect(2, velocityY, previousVelocityY, previousVelocityX, previousVelocityY, deltaTime);
            Project(velocityX, velocityY, previousVelocityX, previousVelocityY);
            Diffuse(0, sourceDensity, density, diffusionRate, deltaTime);
            Advect(0, density, sourceDensity, velocityX, velocityY, deltaTime);
        }

        private void AddGravity(float gravityX, float gravityY, float deltaTime) {
            for (int i = 0; i < gridSize * gridSize; i++) {
                velocityX[i] += gravityX * deltaTime;
                velocityY[i] += gravityY * deltaTime;
            }
        }

        private void Diffuse(int boundary, float[] current, float[] previous, float rate, float dt) {
            float a = dt * rate * (gridSize - 2) * (gridSize - 2);
            LinearSolve(boundary, current, previous, a, 1 + 4 * a);
        }

        private void Project(float[] velocX, float[] velocY, float[] p, float[] div) {
            for (int j = 1; j < gridSize - 1; j++) {
                for (int i = 1; i < gridSize - 1; i++) {
                    div[i + j * gridSize] = -0.5f * (
                        velocX[(i + 1) + j * gridSize] - velocX[(i - 1) + j * gridSize] +
                        velocY[i + (j + 1) * gridSize] - velocY[i + (j - 1) * gridSize]
                    ) / gridSize;
                    p[i + j * gridSize] = 0;
                }
            }
            SetBoundaries(0, div);
            SetBoundaries(0, p);
            LinearSolve(0, p, div, 1, 4);

            for (int j = 1; j < gridSize - 1; j++) {
                for (int i = 1; i < gridSize - 1; i++) {
                    velocX[i + j * gridSize] -= 0.5f * (
                        p[(i + 1) + j * gridSize] - p[(i - 1) + j * gridSize]
                    ) * gridSize;
                    velocY[i + j * gridSize] -= 0.5f * (
                        p[i + (j + 1) * gridSize] - p[i + (j - 1) * gridSize]
                    ) * gridSize;
                }
            }
            SetBoundaries(1, velocX);
            SetBoundaries(2, velocY);
        }

        private void Advect(int boundary, float[] d, float[] d0, float[] velocX, float[] velocY, float dt) {
            float i0, i1, j0, j1;
            float dtx = dt * (gridSize - 2);
            float dty = dt * (gridSize - 2);
            float s0, s1, t0, t1;
            float tempX, tempY, x, y;

            for (int j = 1; j < gridSize - 1; j++) {
                for (int i = 1; i < gridSize - 1; i++) {
                    tempX = dtx * velocX[i + j * gridSize];
                    tempY = dty * velocY[i + j * gridSize];
                    x = i - tempX;
                    y = j - tempY;

                    if (x < 0.5f) x = 0.5f;
                    if (x > gridSize - 1.5f) x = gridSize - 1.5f;
                    i0 = (float)Math.Floor(x);
                    i1 = i0 + 1.0f;
                    if (y < 0.5f) y = 0.5f;
                    if (y > gridSize - 1.5f) y = gridSize - 1.5f;
                    j0 = (float)Math.Floor(y);
                    j1 = j0 + 1.0f;

                    s1 = x - i0;
                    s0 = 1.0f - s1;
                    t1 = y - j0;
                    t0 = 1.0f - t1;

                    int i0i = (int)i0;
                    int i1i = (int)i1;
                    int j0i = (int)j0;
                    int j1i = (int)j1;

                    d[i + j * gridSize] =
                      s0 * (t0 * d0[i0i + j0i * gridSize] + t1 * d0[i0i + j1i * gridSize]) +
                      s1 * (t0 * d0[i1i + j0i * gridSize] + t1 * d0[i1i + j1i * gridSize]);
                }
            }
            SetBoundaries(boundary, d);
        }

        private void SetBoundaries(int boundary, float[] values) {
            for (int i = 1; i < gridSize - 1; i++) {
                values[i] = (boundary == 1) ? -values[i + gridSize] : values[i + gridSize];
                values[i + (gridSize - 1) * gridSize] = (boundary == 1) ? -values[i + (gridSize - 2) * gridSize] : values[i + (gridSize - 2) * gridSize];
                values[i * gridSize] = (boundary == 2) ? -values[1 + i * gridSize] : values[1 + i * gridSize];
                values[(gridSize - 1) + i * gridSize] = (boundary == 2) ? -values[(gridSize - 2) + i * gridSize] : values[(gridSize - 2) + i * gridSize];
            }
            // Corner cases
            values[0] = 0.5f * (values[1] + values[gridSize]);
            values[(gridSize - 1) + (gridSize - 1) * gridSize] = 0.5f * (values[(gridSize - 2) + (gridSize - 1) * gridSize] + values[(gridSize - 1) + (gridSize - 2) * gridSize]);
            values[gridSize - 1] = 0.5f * (values[gridSize - 2] + values[(gridSize - 1) + gridSize]);
            values[(gridSize - 1) * gridSize] = 0.5f * (values[(gridSize - 1) * gridSize + 1] + values[(gridSize - 2) * gridSize]);
        }

        private void LinearSolve(int boundary, float[] current, float[] previous, float a, float c) {
            float cRecip = 1.0f / c;
            for (int k = 0; k < 20; k++) {
                for (int j = 1; j < gridSize - 1; j++) {
                    for (int i = 1; i < gridSize - 1; i++) {
                        current[i + j * gridSize] =
                            (previous[i + j * gridSize] + a * (current[i + 1 + j * gridSize] + current[i - 1 + j * gridSize] + current[i + (j + 1) * gridSize] + current[i + (j - 1) * gridSize])) * cRecip;
                    }
                }
                SetBoundaries(boundary, current);
            }
        }

        public void PrintDensity() {
            for (int j = 1; j < gridSize - 1; j++) {
                for (int i = 1; i < gridSize - 1; i++) {
                    // Print "X" for water presence, space for none. Adjust the threshold as needed.
                    Console.Write(density[i + j * gridSize] > 0.1f ? "XX" : "  ");
                }
                Console.WriteLine();
            }
            Console.WriteLine();
        }

    }
}

(For debugging purposes, I’m just calling PrintDensity every step and watching the output in a console)

Is there a well known approach I should be adopting? Or is SPH really the best option on the table? (If so, can you see an obvious mistake in my implementation?)

I’m invoking the sim with the following…

Stopwatch stopwatch = new Stopwatch();
stopwatch.Restart();

var sim = new FluidSimulation(64, .001f, .001f);
sim.SetInitialState(.1f);

var shouldAbort = false;

while (!shouldAbort) {
    stopwatch.Stop();
    TimeSpan timeBetweenFrames = stopwatch.Elapsed;
    Thread.Sleep(15);
    stopwatch.Restart();
    sim.Update(0, -.098f, (float)timeBetweenFrames.TotalMilliseconds / 1000);

    Console.WriteLine(stopwatch.Elapsed);

    // Print the density grid to the console
    sim.PrintDensity();
    if (Console.KeyAvailable && Console.ReadKey(true).Key == ConsoleKey.Escape) {
        shouldAbort = true;
    }
}
```