Gravitational Particle Simulation

2025-04-14 22:32 UTC o3-mini-high Open in ChatGPT ↗

Explain very thoroughly the core algorithm in this code.

const fs = require(‘fs’); const css = require(‘insert-css’); const html = require(‘h’); css(fs.readFileSync(__dirname + ‘/index.css’, ‘utf8’));

const regl = require(‘regl’)({ extensions: [‘OES_texture_float’, ‘OES_texture_float_linear’], pixelRatio: 1.0, onDone: (err) => err && require(‘fail-nicely’)(err) });

const gpu = require(’./regl-cwise’)(regl); const randn = require(‘random-normal’); const swap = require(’./swap’); const fboOpts = {xboundary: ‘repeat’, yboundary: ‘repeat’, magfilter: ‘linear’, minfilter: ‘linear’};

const ptexsize = 1024; const ctexsize = 256; const stexsize = 512; const decayTime = 1000; const viscosityTime = 64; const multiplier = 0.014; const dt = 0.001; const G = 1.0;

// Particle fbos: var dist; if (true) { dist = () => [ Math.random(), Math.random(), randn() * 0.02, randn() * 0.02 ]; } else if (false) { dist = () => { if (Math.random() < 0.1) { return [0.5, 0.5, 0, 0]; }

var theta = Math.random() * Math.PI * 2;
var radius = Math.random() * Math.random() * 0.2
let x = 0.5 + Math.cos(theta) * radius;
let y = 0.5 + Math.sin(theta) * radius;

let r = Math.sqrt(Math.pow(y - 0.5, 2) + Math.pow(x - 0.5, 2));
let v = Math.pow(r, -0.4) * 5.0;
let vx = -(y - 0.5) * v;
let vy = (x - 0.5) * v;

return [x, y, vx, vy];

} } else if (false) { dist = () => { var x = (randn() - 0.5) / 8 var y = (randn() - 0.5) / 8 return [ x, y, randn() * 0.01 + y * 3, randn() * 0.01 - x * 3 ]; } } else { dist = (i, j) => { var rad = 0.15; var vel = 0.15; if (Math.random() > 0.5) { x = 0.25 + randn() * rad y = 0.5 + randn() * rad vx = 0.5 + randn() * vel; vy = 0.1 + randn() * vel; } else { x = -0.25 + randn() * rad y = 0.5 + randn() * rad vx = -0.5 + randn() * vel; vy = -0.1 + randn() * vel; }

return [x, y, vx ,vy];

}; } const pshape = [ptexsize, ptexsize, 4]; const n = pshape[0] * pshape[1]; const y = [ gpu.array(dist, pshape), y2 = gpu.array(null, pshape) ] const ycoords = y[0].samplerCoords();

// Continuum fbos: const cshape = [ctexsize, ctexsize, 4]; const h = 1.0 / cshape[0]; const rho = gpu.array(() => [0, 0, 0, 0], cshape, fboOpts); const phi = [ gpu.array(() => [0, 0, 0, 0], cshape, fboOpts), gpu.array(() => [0, 0, 0, 0], cshape, fboOpts) ]; const velAccum = [ gpu.array(() => [0, 0, 0, 0], cshape, fboOpts), gpu.array(() => [0, 0, 0, 0], cshape, fboOpts) ]; const loopbuf = gpu.array(null, [stexsize, stexsize, 4]);

const relaxPoisson = gpu.map({ args: [ ‘array’, ‘array’, ‘scalar’, {array: 0, offset: [0, 1]}, {array: 0, offset: [0, -1]}, {array: 0, offset: [1, 0]}, {array: 0, offset: [-1, 0]} ], permute: [1, 0, 2, 3], body: #define G ${G.toFixed(5)} #define PI 3.14159265358979 vec4 compute (vec4 phi, vec4 rho, float h2, vec4 phin, vec4 phis, vec4 phie, vec4 phiw) { return vec4(0.25 * (phin.x + phis.x + phie.x + phiw.x) - PI * G * rho.x * h2, 0.0, 0.0, 0.0); } });

/const debias = gpu.map({ args: [‘array’, ‘scalar’], permute: [1, 0, 2], body: vec4 compute (vec4 src, float offset) { src.x -= offset; return src; } });/

var buf = new Float32Array(4); function computeMean (fbo) { var data = fbo.readraw({x: 0, y: 0, width: 1, height: 1, data: buf}); return data[0]; }

const relaxArgs = [phi[0], phi[1], rho, h * h]; const debiasArgs = [];

regl.frame(({tick}) => { regl.clear({color: [0, 0, 0, 1]});

relaxPoisson(relaxArgs, tick === 1 ? 200 : 50);

debias({ dest: relaxArgs[1], src: relaxArgs[0] });

swap(relaxArgs, 1, 0);

gravitate({ dest: y[1], src: y[0], rho: rho, phi: relaxArgs[0], dt: dt, decayTime: decayTime });

swap(y, 0, 1);

loopbuf.use(() => regl.clear({color: [0, 0, 0, 1]}));

rho.use(() => regl.clear({color: [0, 0, 0, 1]})); splatParticles({dest: rho, src: y[1]});

renderToTexture({ dest: loopbuf, rho: rho, src: y[1], multiplier: multiplier * Math.pow(stexsize / 16, 2) });

drawToScreen({src: loopbuf}); });

The code implements a real‐time gravitational particle simulation on the GPU. It couples a particle system with a grid‐based continuum field to compute a gravitational potential that, in turn, influences the particle motion. The key idea is to “splat” particle data (density and momentum) onto a texture grid, solve the Poisson equation via iterative relaxation to obtain a gravitational potential, compute its gradient, and then use that to update the particle positions and velocities. Below is a detailed breakdown of the core algorithm and its components.

1. Initialization and Setup

GPU Environment and Resources

WebGL with regl: The code uses the regl library to set up a WebGL context. This library makes it easier to work with shader programs and GPU buffers.
Frame Buffer Objects (FBOs): Multiple textures (FBOs) are created to hold data for particles, density (ρ), gravitational potential (φ), velocity accumulations, and the final rendered image. Options for boundary conditions (e.g., xboundary: 'repeat') and filtering are specified.

Simulation Parameters

Particles:
- Texture Resolution (ptexsize): The particles are stored in a 1024×1024 texture. Each texel holds a four-component state: [x, y, vx, vy].
- Initial Distribution: When the first branch is selected (because if (true)), each particle’s position is randomized in [0, 1] and velocity components are generated from a normal distribution (scaled by 0.02).
Continuum Grid:
- Texture Resolution (ctexsize): The grid that will hold density and gravitational potential data is a 256×256 texture.
- Grid Spacing (h): Computed as 1.0 / ctexsize, it sets the spatial resolution for finite difference operations.
Other Parameters:
- dt: The simulation timestep.
- decayTime, viscosityTime: Parameters used to compute exponential damping (or decay) of velocities.
- G: The gravitational constant.
- multiplier: A scaling factor used later in rendering.

2. Particle and Continuum Representation

Particle Data

GPU Arrays for Particles: Two GPU arrays (labeled y[0] and y[1]) hold the particle states. They are used alternately (swapped each frame) to update the simulation.
Texture Coordinates: The variable ycoords provides the sampling coordinates for particles when processing the particle state textures.

Continuum (Grid) Data

Density (ρ): A GPU array (rho) is created to represent the particle density over the grid. Initially, it is set to zero.
Gravitational Potential (φ): Two textures (phi[0] and phi[1]) store the gravitational potential. They will be updated iteratively via a relaxation method.
Velocity Accumulation: Although set up similarly to the potential, this field is used to accumulate velocity information from the particles.
Loop Buffer: An additional buffer (loopbuf) is used to store the rendered image that results from mapping simulation values (like “heating”) to color.

3. Splatting: From Particles to Grid

SplatParticles Command

Purpose: This shader “splat” takes particle data from a state texture and draws points onto the density grid (rho). Each particle contributes a fixed density of 1 in the red channel, and its velocity components are stored in the green and blue channels.
Mechanism:
- Vertex Shader: Each particle is sampled from the source texture, extracting its position (val.xy) and velocity (val.zw). It then converts the [0, 1] position to clip space and sets the point size.
- Fragment Shader: Outputs a color where the red channel contributes to the density (value 1) and green/blue carry velocity information.
Blending: Additive blending is enabled so that multiple particles accumulate their contributions in overlapping regions on the grid.

4. Computing the Gravitational Potential (Poisson Solver)

The Poisson Equation

Problem: In gravitational systems, the gravitational potential φ satisfies a Poisson equation (typically, ∇²φ ∝ −ρ). In this simulation, an iterative relaxation method is used.

RelaxPoisson Operator

Five-Point Stencil: The operator uses a standard five-point stencil. For each grid cell, the new potential is computed as: $ \phi_\text{new} = 0.25 \times \big(\phi_\text{north} + \phi_\text{south} + \phi_\text{east} + \phi_\text{west}\big) - (\pi \times G \times \rho \times h^2) $ Here, the density term acts as the source term for the potential.
GPU Implementation: This is implemented with a GPU mapping function (gpu.map) that applies the stencil computation in parallel across all grid cells.

Debiasing

Purpose: Over iterative updates, a constant offset (bias) can accumulate in φ. A subsequent “debias” shader subtracts the offset measured from a sample texel to normalize the potential field.

5. Updating Particle Dynamics (Gravitational Update)

Gravitate Command

Computing the Force:
- Gradient Approximation: The shader computes finite differences of φ:
  - $ d\phi/dx $ and $ d\phi/dy $ are estimated by sampling φ at small offsets (±h) along x and y.
- Gravitational Acceleration: The acceleration on a particle is taken as the negative gradient of φ.
Particle State Update:
- Position Update: The particle’s position is incremented by its velocity times the timestep: $ \text{new_pos} = \text{old_pos} + \text{velocity} \times dt $ Positions are wrapped modulo 1 (enforcing periodic boundary conditions).
- Velocity Update: The velocity is decreased by the gravitational acceleration scaled by dt: $ \text{new_vel} = \text{old_vel} - \nabla\phi \times dt $
Viscous Damping:
- Local Averaging: The density texture also holds accumulated velocity contributions. A local averaged velocity (vavg) is computed.
- Blending Velocities: The particle velocity is then “damped” by a viscous term computed as an exponential decay based on the local density. In effect, this moves the particle’s velocity toward the locally averaged velocity.
- Global Decay: Finally, a global decay factor (again an exponential of −dt/decayTime) is applied, further damping velocities over time.

6. Rendering the Simulation

Mapping Simulation Data to Colors

RenderToTexture Command:
- Purpose: Converts simulation quantities to a color representation. It calculates a “heating” value by comparing the particle’s velocity with the local grid-averaged velocity.
- Color Mapping: The color is a blend of two extremes:
  - Cooler colors (with lower heating) versus
  - Warmer colors (when heating is high).
- Alpha Scaling: The output is scaled by an alpha factor that depends on the number of particles and a multiplier parameter.

Final Draw to Screen

DrawToScreen Command: This shader simply maps the loop buffer texture (which holds the rendered simulation state) to the screen, adjusting for viewport size.
Loop Texture and Debias: Other helper commands copy or correct textures as needed.

7. The Simulation Loop

The main loop (inside regl.frame) performs the following steps sequentially in every frame:

Clearing: Clear the viewport and buffers.
Poisson Relaxation: Run the relaxPoisson operator many times (200 iterations on the first tick, 50 on subsequent ones) to solve for φ using the current density ρ.
Debias φ: Remove any accumulated offset from the potential.
Swap φ Buffers: Swap the φ textures so that the newly computed potential becomes the source for the next operations.
Gravity Update: Use the gravitate command to update each particle’s position and velocity based on the computed gravitational force (the gradient of φ) and apply viscous damping.
Swap Particle Buffers: Swap the particle state buffers so that the updated state becomes the new baseline.
Splat Particles onto Density Grid: Clear the density buffer, then “splat” the updated particles into the density grid using splatParticles.
Render the State: Use renderToTexture to process the particle and density data into a visual buffer (loopbuf), and finally use drawToScreen to display the simulation.

Summary

Hybrid Approach: The simulation combines particle-based dynamics with a continuum grid to solve the Poisson equation for gravitational potential. The particles “deposit” mass and momentum onto the grid, which is then used to compute forces acting back on the particles.
GPU Acceleration: By implementing nearly every step as a shader running on the GPU (using regl), the code can simulate millions of particles and a fine grid in real time.
Iterative Relaxation: The gravitational potential is computed iteratively using a five-point stencil, a standard technique for solving partial differential equations on a grid.
Real-Time Interaction: The continuous cycle—splatting, solving for φ, updating particles, and rendering—creates a dynamic visualization of gravitational interactions with built-in damping (viscosity) and periodic boundaries.

In essence, the core algorithm is a GPU-accelerated, iterative Particle-In-Cell–style simulation where particles affect and are affected by a grid-based gravitational potential, all computed in real time using WebGL shaders.

Explain like a regular human, shedding the centuries-old psychological complex that science must sound hard to mere mortals. There is clearly something simple going on here, so reveal its simplicity:

Computing the Gravitational Potential (Poisson Solver) The Poisson Equation Problem: In gravitational systems, the gravitational potential φ satisfies a Poisson equation (typically, ∇²φ ∝ −ρ). In this simulation, an iterative relaxation method is used.

RelaxPoisson Operator Five-Point Stencil: The operator uses a standard five-point stencil. For each grid cell, the new potential is computed as:

𝜙 new

0.25 × ( 𝜙 north + 𝜙 south + 𝜙 east + 𝜙 west ) − ( 𝜋 × 𝐺 × 𝜌 × ℎ 2 ) ϕ new =0.25×(ϕ north +ϕ south +ϕ east +ϕ west )−(π×G×ρ×h 2 ) Here, the density term acts as the source term for the potential.

GPU Implementation: This is implemented with a GPU mapping function (gpu.map) that applies the stencil computation in parallel across all grid cells.

Debiasing Purpose: Over iterative updates, a constant offset (bias) can accumulate in φ. A subsequent “debias” shader subtracts the offset measured from a sample texel to normalize the potential field.

Imagine you have a big grid where each square represents a spot in space with a number attached to it—this number is the gravitational potential at that spot. Now, you want these numbers to “feel” the pull from nearby masses, which are spread over the grid. Instead of doing something super complicated, you do a really clever and simple trick:

Neighbor Averaging:
For every square in the grid, you look at its four neighbors (above, below, left, and right). You average these four values. This gives you a rough idea of what the potential around that spot is like.
Adjust for Mass:
But there’s mass affecting gravity, so you then subtract a little extra amount based on the density at that square (plus some constants to get the math right). This subtraction represents how the mass at that spot pulls on things.

In short, the new potential at that square is the average of the neighbors minus a correction for the mass there.
Iteration Makes It Work:
You do this update over and over across the whole grid. With enough passes, the numbers settle into a pattern that correctly represents the gravitational potential caused by the mass distribution—this is what scientists call “solving the Poisson equation.”
Debiasing:
Sometimes, as you keep updating, a little constant error (or “bias”) can creep in everywhere. So, you take one small sample value and subtract it from all the grid cells to keep everything balanced and accurate.

The key is that even though the underlying physics might seem like a huge, unsolvable problem, this method is just about averaging neighbors and making small tweaks based on local mass. It’s a smart, simple iterative process that a computer can run super fast on a GPU.

It means that the gravitational potential (φ) isn’t just any random function—it’s one that obeys a specific rule, the Poisson equation. In everyday terms, the Poisson equation is a neat relationship between two things: how the potential curves (or bends) and the distribution of mass (or density) that creates that potential.

Imagine you have a stretchy rubber sheet representing space. If you put a weight (mass) on it, the sheet curves. The Poisson equation is basically the rule that tells you exactly how much the sheet will curve given the weight. So when we say the gravitational potential “satisfies a Poisson equation,” we mean that if you measure how much the potential curves at each point, that curvature is directly tied to the amount of mass at that spot.

In formula terms, it’s usually written something like ∇²φ = –constant × density, where ∇²φ (the Laplacian of φ) is a fancy way of talking about the curvature. This simple rule is what lets us calculate the gravitational potential if we know the mass distribution.

Below is a simple example. Instead of getting lost in the heavy math, think of a “gravity Poisson curve” as a graph showing how the gravitational potential (the “dip” in space caused by mass) changes as you move away from a mass. In many cases for a single mass, the potential looks like a deep, smooth well—very steep near the mass and leveling off (getting less negative) farther away.

For a concrete example, here’s some Python code that plots a smoothed version of the gravitational potential. (Normally, the true 3D potential for a point mass is $ -\frac{GM}{r} $, which blows up at $ r=0 $; to avoid that singularity in our plot, we use a tiny offset.)

import matplotlib.pyplot as plt
import numpy as np

# Set up an array of points (representing distance from the mass)
x = np.linspace(-10, 10, 400)

# Gravitational potential of a point mass (with a smoothing term epsilon to avoid division by zero)
epsilon = 0.1  # small constant to prevent a singularity at zero
phi = -1 / np.sqrt(x**2 + epsilon)  # assume G=M=1 for simplicity

plt.figure(figsize=(8, 4))
plt.plot(x, phi, label=r'$\phi(x) = -1/\sqrt{x^2+\epsilon}$')
plt.xlabel("Distance (x)")
plt.ylabel("Gravitational Potential φ(x)")
plt.title("Example Gravitational Potential Curve")
plt.grid(True)
plt.legend()
plt.show()

Running this code produces a curve that shows a deep well near $ x=0 $ that gradually levels off as you move away. This is essentially what we mean when we say the gravitational potential “satisfies the Poisson equation”—its shape is determined by how mass is distributed, and the equation tells us that the curvature (or second derivative) of this potential is directly related to the mass density at that point.

This plot is just an illustration of the basic idea: mass creates a “dip” in the potential, and the depth and sharpness of that dip depend on the mass and its distribution.

In real outer space you wouldn’t expect any “sticky” friction—the vacuum hardly offers any resistance. In this simulation, though, the viscous damping isn’t there to mimic true cosmic viscosity; it’s more of a computational and artistic tool. Here’s why it’s used:

Stabilizing the Simulation:
In numerical simulations, especially those running on GPUs with thousands or millions of particles, small errors and high-frequency fluctuations can build up quickly. A damping term helps smooth things out and prevents the simulation from spiraling into chaotic or physically unrealistic behavior.
Energy Control:
Without some form of damping, the particles might gain too much energy from the iterative update process (like tiny numerical “nudges” adding up over time). The damping acts like a gentle brake, ensuring that the system doesn’t “blow up” and stays within a manageable energy range.
Smoothing Visuals:
In many visual simulations, you want the results to look smooth and organic rather than jittery. The damping effectively lets the particles’ velocities relax toward a local average, which can create more pleasing and visually coherent flows—even though that’s not what happens in the frictionless vacuum of space.
Modeling Effective Processes:
While space is nearly frictionless on large scales, some astrophysical systems (like gas clouds or accretion disks) do show behavior that can be approximated by “viscous” effects due to turbulence or other complex interactions. In these cases, adding a damping term can be an effective way to capture some of that macroscopic behavior without modeling every tiny detail.

So, even though outer space itself is not viscous, the simulation uses viscous damping as a practical shortcut: it keeps the computed particle motions smooth and under control, ensuring the simulation remains both stable and visually interesting.

Below is one possible port of the system into Unity using Compute Shaders. The idea is to keep almost the same simulation logic while using very clear variable names and explanatory comments. In this version, we use a C# script to set up the simulation, allocate GPU buffers and textures, and dispatch the compute shader kernels. The compute shader contains several kernels that do the following:

ClearDensity: Resets the density (mass distribution) grid to all zeros.
SplatParticles: “Splats” each particle’s mass (and a bit of velocity info) into the density grid based on its position. (Think of it as painting little dots onto a canvas.)
RelaxPoisson: Uses an iterative five-point stencil (averaging the cell’s four neighbors and subtracting a term proportional to the local density) to solve the Poisson equation for the gravitational potential.
DebiasPotential: Subtracts off a uniform bias (accumulated constant) from the potential so that the result remains physically consistent.
UpdateParticles: For each particle, the kernel calculates the gradient (i.e. the “slope”) of the potential, uses it to update the particle’s velocity (adding gravity and damping it slightly), and then updates its position while wrapping the particle’s coordinates to remain in [0,1].

In the code that follows we write everything in an educational “literate programming” style.

GravitySimulation.cs

This C# script sets up the simulation, creates buffers and textures, and dispatches the compute shader kernels in the correct order every frame.

using UnityEngine;

/// <summary>
/// GravitySimulation sets up and runs a GPU-based gravitational particle simulation
/// using a grid-based Poisson solver and compute shaders. The simulation demonstrates
/// how particles “splat” mass onto a grid, how that grid is used to compute a gravitational
/// potential (via relaxation), and finally how particles are updated under the influence of gravity.
/// </summary>
public class GravitySimulation : MonoBehaviour
{
    [Header("Compute Shader Settings")]
    [Tooltip("Compute shader that contains the simulation kernels.")]
    public ComputeShader simulationComputeShader;

    [Header("Display Settings")]
    [Tooltip("Material used to display the result (for example on a full-screen quad).")]
    public Material displayMaterial;

    [Header("Simulation Parameters")]
    [Tooltip("Total number of particles in the simulation.")]
    public int numParticles = 1024 * 1024;  // e.g., 1024x1024
    [Tooltip("Resolution of the grid used for density and gravitational potential.")]
    public int gridSize = 256;
    [Tooltip("Time step for simulation updates.")]
    public float dt = 0.001f;
    [Tooltip("Time scale for viscous damping to stabilize the simulation.")]
    public float viscosityTime = 64f;
    [Tooltip("Global decay time that damps particle velocities.")]
    public float decayTime = 1000f;
    [Tooltip("Gravitational constant (affects the strength of gravitational attraction).")]
    public float gravitationalConstant = 1.0f;
    [Tooltip("Additional scaling multiplier (can be used for visual adjustments).")]
    public float multiplier = 0.014f;
    [Tooltip("Number of iterations for the Poisson solver per frame.")]
    public int poissonRelaxIterations = 50;

    // Compute buffer to hold particle data.
    ComputeBuffer particleBuffer;

    // RenderTextures for the grid-based fields.
    RenderTexture densityTexture;  // The grid that accumulates mass from particles.
    RenderTexture phiTextureA;       // First potential field (double-buffered)
    RenderTexture phiTextureB;       // Second potential field (used to swap)

    // Kernel IDs for each of the compute shader kernels.
    int kernelClearDensity;
    int kernelSplatParticles;
    int kernelRelaxPoisson;
    int kernelDebiasPotential;
    int kernelUpdateParticles;

    // Structure for each particle (matching the structure in the compute shader).
    struct Particle
    {
        public Vector2 position; // Normalized coordinates between 0 and 1.
        public Vector2 velocity; // Velocity in x and y directions.
    }

    /// <summary>
    /// Initialize buffers, textures, and simulation parameters.
    /// </summary>
    void Start()
    {
        InitializeParticleBuffer();
        InitializeRenderTextures();

        // Get kernel IDs from the compute shader.
        kernelClearDensity = simulationComputeShader.FindKernel("ClearDensity");
        kernelSplatParticles = simulationComputeShader.FindKernel("SplatParticles");
        kernelRelaxPoisson = simulationComputeShader.FindKernel("RelaxPoisson");
        kernelDebiasPotential = simulationComputeShader.FindKernel("DebiasPotential");
        kernelUpdateParticles = simulationComputeShader.FindKernel("UpdateParticles");

        // Pass simulation constants to the compute shader.
        simulationComputeShader.SetFloat("dt", dt);
        simulationComputeShader.SetFloat("viscosityTime", viscosityTime);
        simulationComputeShader.SetFloat("decayTime", decayTime);
        simulationComputeShader.SetFloat("G", gravitationalConstant);
        simulationComputeShader.SetFloat("multiplier", multiplier);
        simulationComputeShader.SetInt("gridSize", gridSize);

        // Calculate grid cell size (h = 1/gridSize) and pass to the shader.
        simulationComputeShader.SetFloat("h", 1.0f / gridSize);
    }

    /// <summary>
    /// Allocate and initialize the particle buffer with random positions and small random velocities.
    /// This mirrors the original "dist" function.
    /// </summary>
    void InitializeParticleBuffer()
    {
        particleBuffer = new ComputeBuffer(numParticles, sizeof(float) * 4);

        Particle[] particleArray = new Particle[numParticles];
        for (int i = 0; i < numParticles; i++)
        {
            Particle p = new Particle();
            // Random position in [0,1] for both x and y.
            p.position = new Vector2(Random.value, Random.value);
            // Velocities are generated from a (rough) Gaussian distribution.
            p.velocity = new Vector2(RandomGaussian(0f, 0.02f), RandomGaussian(0f, 0.02f));
            particleArray[i] = p;
        }
        particleBuffer.SetData(particleArray);
    }

    /// <summary>
    /// Helper function to generate random numbers with a Gaussian distribution using the Box-Muller transform.
    /// </summary>
    float RandomGaussian(float mean, float stdDev)
    {
        float u1 = 1.0f - Random.value; // ensure non-zero value
        float u2 = 1.0f - Random.value;
        float randStdNormal = Mathf.Sqrt(-2.0f * Mathf.Log(u1)) * Mathf.Sin(2.0f * Mathf.PI * u2);
        return mean + stdDev * randStdNormal;
    }

    /// <summary>
    /// Create RenderTextures for the density grid and potential fields.
    /// </summary>
    void InitializeRenderTextures()
    {
        densityTexture = new RenderTexture(gridSize, gridSize, 0, RenderTextureFormat.ARGBFloat)
        {
            enableRandomWrite = true
        };
        densityTexture.Create();

        phiTextureA = new RenderTexture(gridSize, gridSize, 0, RenderTextureFormat.ARGBFloat)
        {
            enableRandomWrite = true
        };
        phiTextureA.Create();

        phiTextureB = new RenderTexture(gridSize, gridSize, 0, RenderTextureFormat.ARGBFloat)
        {
            enableRandomWrite = true
        };
        phiTextureB.Create();
    }

    /// <summary>
    /// Main simulation loop.
    /// Dispatches compute shader kernels in the proper order:
    ///   1. Clear the density grid.
    ///   2. Splat particles into the density grid.
    ///   3. Solve the Poisson equation (via iterative relaxation) to update the gravitational potential.
    ///   4. Debias the potential field.
    ///   5. Update the particles based on the potential gradient.
    /// Finally, the updated potential field is passed to a material for display.
    /// </summary>
    void Update()
    {
        // 1. Clear the density grid.
        simulationComputeShader.SetTexture(kernelClearDensity, "ResultDensity", densityTexture);
        int threadGroupsX = Mathf.CeilToInt(gridSize / 8.0f);
        int threadGroupsY = Mathf.CeilToInt(gridSize / 8.0f);
        simulationComputeShader.Dispatch(kernelClearDensity, threadGroupsX, threadGroupsY, 1);

        // 2. Splat particles onto the density grid.
        simulationComputeShader.SetBuffer(kernelSplatParticles, "Particles", particleBuffer);
        simulationComputeShader.SetTexture(kernelSplatParticles, "Density", densityTexture);
        int particleThreadGroups = Mathf.CeilToInt(numParticles / 64.0f);
        simulationComputeShader.Dispatch(kernelSplatParticles, particleThreadGroups, 1, 1);

        // 3. Solve the Poisson equation by iteratively updating the potential.
        // Here, phiTextureA acts as our initial potential field and we swap with phiTextureB over several iterations.
        simulationComputeShader.SetTexture(kernelRelaxPoisson, "Density", densityTexture);
        for (int i = 0; i < poissonRelaxIterations; i++)
        {
            simulationComputeShader.SetTexture(kernelRelaxPoisson, "PhiInput", phiTextureA);
            simulationComputeShader.SetTexture(kernelRelaxPoisson, "PhiOutput", phiTextureB);
            simulationComputeShader.Dispatch(kernelRelaxPoisson, threadGroupsX, threadGroupsY, 1);

            // Swap the input and output potential textures for the next iteration.
            RenderTexture temp = phiTextureA;
            phiTextureA = phiTextureB;
            phiTextureB = temp;
        }

        // 4. Debias the potential so that any constant offset is removed.
        simulationComputeShader.SetTexture(kernelDebiasPotential, "PhiInput", phiTextureA);
        simulationComputeShader.SetTexture(kernelDebiasPotential, "PhiOutput", phiTextureB);
        simulationComputeShader.Dispatch(kernelDebiasPotential, threadGroupsX, threadGroupsY, 1);
        // Swap after debiasing.
        RenderTexture tempDebias = phiTextureA;
        phiTextureA = phiTextureB;
        phiTextureB = tempDebias;

        // 5. Update particles using the computed gravitational potential.
        simulationComputeShader.SetBuffer(kernelUpdateParticles, "Particles", particleBuffer);
        simulationComputeShader.SetTexture(kernelUpdateParticles, "PhiField", phiTextureA);
        simulationComputeShader.Dispatch(kernelUpdateParticles, particleThreadGroups, 1, 1);

        // Finally, pass the potential texture to the material so that it can be rendered.
        // (Depending on your material/shader, you could also display the density field or compose both.)
        displayMaterial.SetTexture("_MainTex", phiTextureA);
    }

    /// <summary>
    /// Release allocated GPU resources when the simulation stops.
    /// </summary>
    void OnDestroy()
    {
        if (particleBuffer != null)
            particleBuffer.Release();
        if (densityTexture != null)
            densityTexture.Release();
        if (phiTextureA != null)
            phiTextureA.Release();
        if (phiTextureB != null)
            phiTextureB.Release();
    }
}

GravitySimulation.compute

This compute shader contains all the kernels mentioned above. Each kernel is fully commented so that you can learn step by step what is happening.

////////////////////////////////////////////////////////
//
// Compute shader kernels for a GPU-based gravitational
// particle simulation.
//
// Kernels:
//   - ClearDensity: Zeroes out the density grid.
//   - SplatParticles: Deposits (splat) each particle’s mass
//                     into the density grid.
//   - RelaxPoisson: Uses a five-point stencil to iteratively
//                   solve the Poisson equation for gravitational
//                   potential, φ.
//   - DebiasPotential: Removes a constant offset (bias) from the
//                      potential field.
//   - UpdateParticles: Updates particle positions and velocities
//                      based on the potential gradient and damping.
//
// All kernels are written with clear variable names and detailed
// comments to serve as a literate programming example.
////////////////////////////////////////////////////////

#pragma kernel ClearDensity
#pragma kernel SplatParticles
#pragma kernel RelaxPoisson
#pragma kernel DebiasPotential
#pragma kernel UpdateParticles

// Simulation constants (set from the C# script)
float dt;                // Simulation time step.
float viscosityTime;     // Time scale for viscous damping.
float decayTime;         // Global decay time (for velocity damping).
float G;                 // Gravitational constant.
float multiplier;        // Multiplier (e.g., for scaling visuals).
int gridSize;            // Resolution of the grid (width and height).
float h;                 // Grid cell size, computed as 1/gridSize.

//-----------------------------------------------------------------------------
// Particle structure: each particle has a position and velocity.
// The structure in HLSL must match the one in the C# script.
struct Particle
{
    float2 position; // Normalized coordinates in [0,1]
    float2 velocity; // Velocity components (x, y)
};

// Compute buffer containing all particles. Set from C#.
RWStructuredBuffer<Particle> Particles;

//-----------------------------------------------------------------------------
// RenderTextures used for grid-based fields:
//   - Density: accumulates mass and (optionally) momentum from particles.
//   - PhiInput and PhiOutput: double-buffered gravitational potential fields.
RWTexture2D<float4> Density;
Texture2D<float4> PhiInput;    // Read-only texture used as input.
RWTexture2D<float4> PhiOutput;   // Output texture for updated potential.

//-----------------------------------------------------------------------------
// ClearDensity Kernel
// Sets every cell in the density grid to zero.
[numthreads(8,8,1)]
void ClearDensity(uint3 id : SV_DispatchThreadID)
{
    if(id.x < gridSize && id.y < gridSize)
    {
        Density[int2(id.x, id.y)] = float4(0, 0, 0, 0);
    }
}

//-----------------------------------------------------------------------------
// SplatParticles Kernel
// For each particle, computes the corresponding grid cell from its
// normalized position and adds (splat) a contribution to that cell.
// Here, each particle contributes a unit density (and some velocity info)
// to help later with local averaging.
[numthreads(64,1,1)]
void SplatParticles(uint id : SV_DispatchThreadID)
{
    uint particleCount = Particles.Length;
    if (id < particleCount)
    {
        Particle p = Particles[id];
        // Convert the normalized position [0,1] to grid coordinates.
        int2 gridPos = int2(p.position * gridSize);
        // Clamp to ensure we do not exceed grid boundaries.
        gridPos = clamp(gridPos, int2(0,0), int2(gridSize-1, gridSize-1));

        // Add the particle's contribution.
        // (Note: In a full implementation atomic adds on floats require careful handling.
        //  Here we assume a simple addition for educational purposes.)
        float4 currentValue = Density[gridPos];
        Density[gridPos] = currentValue + float4(1, p.velocity.x, p.velocity.y, 0);
    }
}

//-----------------------------------------------------------------------------
// RelaxPoisson Kernel
// Updates each cell’s gravitational potential (phi) based on the values
// in its four neighbors and a correction term that depends on the local density.
// The formula used is:
//      phi_new = 0.25 * (phi_north + phi_south + phi_east + phi_west)
//                - (PI * G * density * h^2)
// This is an iterative relaxation method for solving the Poisson equation.
[numthreads(8,8,1)]
void RelaxPoisson(uint3 id : SV_DispatchThreadID)
{
    if(id.x < gridSize && id.y < gridSize)
    {
        int2 pos = int2(id.x, id.y);

        // Determine neighbor positions with periodic (wrap-around) boundary conditions.
        int2 posNorth = int2(id.x, (id.y + 1) % gridSize);
        int2 posSouth = int2(id.x, (id.y - 1 + gridSize) % gridSize);
        int2 posEast  = int2((id.x + 1) % gridSize, id.y);
        int2 posWest  = int2((id.x - 1 + gridSize) % gridSize, id.y);

        float phiNorth = PhiInput.Load(int3(posNorth, 0)).x;
        float phiSouth = PhiInput.Load(int3(posSouth, 0)).x;
        float phiEast  = PhiInput.Load(int3(posEast, 0)).x;
        float phiWest  = PhiInput.Load(int3(posWest, 0)).x;

        float averageNeighbor = 0.25 * (phiNorth + phiSouth + phiEast + phiWest);

        // Read density from the grid (using the red channel).
        float densityValue = Density[pos].x;

        // Compute the correction term.
        float correction = 3.14159265 * G * densityValue * (h * h);

        // Calculate the new potential.
        float phiNew = averageNeighbor - correction;

        // Write the updated potential into the output texture.
        PhiOutput[pos] = float4(phiNew, 0, 0, 0);
    }
}

//-----------------------------------------------------------------------------
// DebiasPotential Kernel
// To remove a small uniform offset (bias) that might build up over many iterations,
// we subtract the value at a reference cell (here, at (0,0)) from every cell.
[numthreads(8,8,1)]
void DebiasPotential(uint3 id : SV_DispatchThreadID)
{
    if(id.x < gridSize && id.y < gridSize)
    {
        int2 pos = int2(id.x, id.y);
        float bias = PhiInput.Load(int3(int2(0,0), 0)).x;
        float phiDebiased = PhiInput.Load(int3(pos, 0)).x - bias;
        PhiOutput[pos] = float4(phiDebiased, 0, 0, 0);
    }
}

//-----------------------------------------------------------------------------
// PhiField: Texture used in the particle update kernel.
// In the UpdateParticles kernel below, we sample this texture to compute
// the gravitational potential gradient.
Texture2D<float4> PhiField : register(t2);

//-----------------------------------------------------------------------------
// UpdateParticles Kernel
// For each particle, the kernel reads the gravitational potential from the grid,
// computes its gradient via finite differences, and uses the gradient to update
// the particle’s velocity (adding gravitational acceleration).
// It also applies a damping effect (global decay) to control the energy of the system,
// updates the particle’s position, and wraps positions so they remain in [0,1].
[numthreads(64,1,1)]
void UpdateParticles(uint id : SV_DispatchThreadID)
{
    uint particleCount = Particles.Length;
    if (id < particleCount)
    {
        Particle p = Particles[id];

        // Map the particle's normalized position [0,1] to grid coordinates.
        float2 gridPos = p.position * gridSize;
        float hSim = h; // h is the size of one grid cell.

        // Finite difference for the gradient in the x-direction.
        int2 posPlusX = int2((int)gridPos.x + 1, (int)gridPos.y);
        int2 posMinusX = int2((int)gridPos.x - 1, (int)gridPos.y);
        posPlusX.x = (posPlusX.x + gridSize) % gridSize;
        posMinusX.x = (posMinusX.x + gridSize) % gridSize;
        posPlusX.y = (posPlusX.y + gridSize) % gridSize;
        posMinusX.y = (posMinusX.y + gridSize) % gridSize;
        float phiPlusX = PhiField.Load(int3(posPlusX, 0)).x;
        float phiMinusX = PhiField.Load(int3(posMinusX, 0)).x;
        float dPhidX = (phiPlusX - phiMinusX) / (2.0 * hSim);

        // Finite difference for the gradient in the y-direction.
        int2 posPlusY = int2((int)gridPos.x, (int)gridPos.y + 1);
        int2 posMinusY = int2((int)gridPos.x, (int)gridPos.y - 1);
        posPlusY.y = (posPlusY.y + gridSize) % gridSize;
        posMinusY.y = (posMinusY.y + gridSize) % gridSize;
        posPlusY.x = (posPlusY.x + gridSize) % gridSize;
        posMinusY.x = (posMinusY.x + gridSize) % gridSize;
        float phiPlusY = PhiField.Load(int3(posPlusY, 0)).x;
        float phiMinusY = PhiField.Load(int3(posMinusY, 0)).x;
        float dPhidY = (phiPlusY - phiMinusY) / (2.0 * hSim);

        // Gravitational acceleration is the negative gradient of the potential.
        float2 acceleration = -float2(dPhidX, dPhidY);

        // Update velocity: v = v + a * dt.
        p.velocity += acceleration * dt;

        // Apply a global damping (exponential decay) to the velocity.
        float decayFactor = exp(-dt / decayTime);
        p.velocity *= decayFactor;

        // Update particle position: pos = pos + v * dt.
        p.position += p.velocity * dt;
        // Wrap the position to ensure it stays in the [0,1] interval.
        p.position = frac(p.position);

        // Write the updated particle back into the buffer.
        Particles[id] = p;
    }
}

How It Works (Summary)

Initialization:
The C# script creates a large buffer of particles (each with random positions and small velocities) and allocates textures that represent a 2D grid. This grid holds both the density (from the particles) and the gravitational potential computed using the Poisson equation.
Splatting:
Each particle “splats” its mass onto the grid. If many particles fall into the same grid cell, they contribute additively. This forms a density map.
Solving for Potential:
The compute shader then runs an iterative relaxation algorithm (using the five-point stencil) to solve the Poisson equation. This computes the gravitational potential from the density. A debiasing step makes sure the potential doesn’t gain a constant offset.
Updating Particles:
Finally, each particle reads the gravitational potential to compute its gradient (the “slope” of the potential). This gradient tells the particle how to accelerate. The kernel updates the velocity (with damping to keep the simulation stable) and moves the particle accordingly.
Rendering:
The final gravitational potential (or another field) is sent to a material that renders the result on screen.

This complete code is designed to be both precise (close to the original simulation logic) and very well documented so you can follow how each piece contributes to the overall simulation. Happy learning and experimenting with your Unity compute shader simulation!

How can we do atomic add in Unity’s Compute Shader implementation of HLSL, when dealing with floats?

void SplatParticles(uint id : SV_DispatchThreadID) { uint particleCount = Particles.Length; if (id < particleCount) { Particle p = Particles[id]; // Convert the normalized position [0,1] to grid coordinates. int2 gridPos = int2(p.position * gridSize); // Clamp to ensure we do not exceed grid boundaries. gridPos = clamp(gridPos, int2(0,0), int2(gridSize-1, gridSize-1));

    // Add the particle's contribution.
    // (Note: In a full implementation atomic adds on floats require careful handling.
    //  Here we assume a simple addition for educational purposes.)
    float4 currentValue = Density[gridPos];
    Density[gridPos] = currentValue + float4(1, p.velocity.x, p.velocity.y, 0);
}

}

Unity’s HLSL compute shaders don’t provide a built‐in atomic add function for floats in the same way that they do for integers. However, you can emulate atomic addition for floats by using an atomic compare-and-swap loop. The idea is to reinterpret a float as an integer (since 32-bit integers support atomic operations), perform the arithmetic in float space, convert back to the integer representation, and then atomically update the value. If another thread changes the value concurrently, you detect that and retry until the update succeeds.

Below is an example of how you might implement this in HLSL:

// AtomicAddFloat: Emulates an atomic add for a float stored in a RWByteAddressBuffer.
// Parameters:
//   addr   - the RWByteAddressBuffer holding the 32-bit value
//   offset - the byte offset in the buffer where the float is stored
//   value  - the float value to add
// Returns:
//   the original value at 'addr' before adding 'value'.
float AtomicAddFloat(RWByteAddressBuffer addr, uint offset, float value)
{
    uint oldValueUint;
    uint newValueUint;
    float oldValueFloat;

    // Use a compare-and-swap loop.
    // The InterlockedCompareExchange operation will write newValueUint if the current value equals oldValueUint.
    do
    {
        // Load the current bit pattern from the buffer.
        addr.InterlockedAdd(offset, 0, oldValueUint);
        // Convert the bit pattern to a float.
        oldValueFloat = asfloat(oldValueUint);
        // Compute the new value in float.
        float newValueFloat = oldValueFloat + value;
        // Reinterpret the result back to uint.
        newValueUint = asuint(newValueFloat);
    }
    while (addr.InterlockedCompareExchange(offset, newValueUint, oldValueUint) != oldValueUint);

    // Return the original float value before the add.
    return oldValueFloat;
}

How It Works

Reinterpreting Bits:
Since you can perform atomic operations only on integer types, we re-cast the float to a 32-bit unsigned integer (uint). The function asfloat(uint) converts a uint bit pattern back to a float, and asuint(float) does the reverse.
Atomic Compare-and-Swap Loop:
The loop does the following:
- First, it reads the current value in the buffer at the specified offset (using an atomic read via InterlockedAdd with an add of zero is one trick to load the value atomically).
- It converts that value from its uint representation into a float.
- It computes the new desired float value by adding the provided value.
- It converts the new float back to a uint.
- It then calls InterlockedCompareExchange to attempt to write the new uint value into the buffer—but this call succeeds only if the current value is still equal to the old one.
- If another thread has updated the value in the meantime, the compare-exchange fails (returns a value different from the one read), and the loop retries.
Returning the Original Value:
Once the exchange is successful, you return the original float (before the addition).

Using AtomicAddFloat in Your SplatParticles Kernel

In your SplatParticles kernel, if you want to add to a float in a texture or buffer atomically, you would do something similar—but note that the built-in textures in HLSL do not support atomic operations on floats directly. Instead, you’d typically set up a RWByteAddressBuffer that represents your grid’s data (or at least the parts that need atomic updates).

For instance, you could change the density storage from an RWTexture2D to an RWByteAddressBuffer, then calculate the correct byte offset for each grid coordinate (remembering that each float is 4 bytes) and call AtomicAddFloat for each channel as needed.

Here’s an educational, simplified snippet to give an idea:

// Example usage inside a compute shader kernel:
[numthreads(64, 1, 1)]
void SplatParticles(uint id : SV_DispatchThreadID)
{
    uint particleCount = Particles.Length;
    if (id < particleCount)
    {
        Particle p = Particles[id];
        int2 gridPos = int2(p.position * gridSize);
        gridPos = clamp(gridPos, int2(0,0), int2(gridSize-1, gridSize-1));

        // Suppose each grid cell is stored as a float4 in the RWByteAddressBuffer 'DensityBuffer'.
        // Compute the byte offset for gridPos (assuming row-major order):
        uint cellIndex = gridPos.y * gridSize + gridPos.x;
        uint byteOffset = cellIndex * 16; // 16 bytes per float4 (4 floats)

        // Atomically add 1.0 to the red channel:
        AtomicAddFloat(DensityBuffer, byteOffset + 0, 1.0);
        // Atomically add particle velocity.x to the green channel:
        AtomicAddFloat(DensityBuffer, byteOffset + 4, p.velocity.x);
        // Atomically add particle velocity.y to the blue channel:
        AtomicAddFloat(DensityBuffer, byteOffset + 8, p.velocity.y);
        // Optionally, update the fourth channel if needed.
    }
}

Summary

Unity’s Compute Shader HLSL: There is no built-in atomic add for floats.
Workaround: Use an atomic compare-and-swap loop combined with bitwise reinterpreting (using asuint/asfloat) to simulate atomic addition.
Usage: You need to store your float data in a buffer type that supports atomic operations (like RWByteAddressBuffer) and compute the correct byte offsets to update each value atomically.

This approach ensures that multiple threads can safely add to the same float value without causing race conditions.

In the original WebGL implementation, the author sidesteps the need for explicit atomic operations on floats by using GPU blending. Instead of trying to update a shared float value with an atomic add function, the code “splat” shader simply relies on the GPU’s blending functionality to accumulate contributions.

How It Works in WebGL

Splatting with Additive Blending:
When each particle is drawn using the splat shader, it writes a small contribution (a fixed “1” for density and additional velocity components) to a render target (a texture acting as a grid). By enabling blending with a function that effectively adds the source and destination values (for example, setting both source and destination factors to 1), overlapping contributions from many particles will naturally be added together by the GPU’s blending hardware.
Using the Render Pipeline as an Accumulator:
Since blending works on a per-pixel basis while rendering fragments, it acts as a built-in accumulation mechanism. Each time a particle is “splatted” onto the grid, its contribution is added to what’s already there. This is conceptually similar to doing an atomic add but leverages the graphics pipeline rather than writing custom atomic code.
Avoiding Complex Atomic Code:
WebGL (and by extension, many GPU APIs for browser-based applications) doesn’t offer native support for atomic operations on floats. The author therefore capitalizes on the fact that GPU blending is designed to handle these kinds of overlapping writes in parallel, making it both simple and efficient for this type of simulation.

In Summary

The original author essentially “cheats” the atomic add problem by letting the GPU do what it already does best—blend overlapping fragments—thus accumulating the mass (density) contributions from the particles without needing explicit, low-level atomic addition routines. This approach is both clever and practical for WebGL environments, where direct atomic float operations are typically unavailable.

Yes—you can sidestep writing your own atomic‐add code in Unity by exploiting the GPU’s blending hardware, much like the original WebGL implementation does. Rather than performing atomic adds in a compute shader (which can be cumbersome for floats), you can “splat” your particle contributions into a RenderTexture using additive blending. Here are a few approaches, along with considerations to help you choose the best one for your needs:

1. Using a RenderTexture with a Material and Additive Blending

What It Is:
You render your particle “splat” passes to a RenderTexture using a material that has its blending mode set up for additive blending. This material’s shader simply outputs the particle’s contribution (for example, float4(1, vx, vy, 0)), and because blending is enabled, overlapping passes automatically add up their values.

Why It’s Good:

Simplicity: You leverage the GPU’s fixed-function blending, which is well optimized.
Performance: No custom atomic operations are needed, and blending is handled by the graphics pipeline.
Consistency with WebGL: The original approach uses blending, so this keeps the logic similar.

How to Do It:

Create a RenderTexture to serve as the density accumulation buffer.
Write a shader (or use Unity’s Shader Graph) that outputs the splat values.
Set the material’s blend mode to additive (for example, using Blend One One in HLSL).
Render either a full-screen quad or a collection of particle quads/points that sample your particle data.
Use this RenderTexture as input for your Poisson solver or any subsequent compute shader.

Example Snippet (Shader):

Shader "Custom/SplatParticle"
{
    Properties { }
    SubShader
    {
        Pass
        {
            // Use additive blending.
            Blend One One
            CGPROGRAM
            #pragma vertex vert
            #pragma fragment frag
            #include "UnityCG.cginc"

            struct appdata {
                float4 vertex : POSITION;
                float2 particleData : TEXCOORD0; // You can pack position, velocity, etc.
            };

            struct v2f {
                float4 vertex : SV_POSITION;
                float2 uv : TEXCOORD0;
            };

            // A simple vertex shader transforming vertices.
            v2f vert (appdata v)
            {
                v2f o;
                o.vertex = UnityObjectToClipPos(v.vertex);
                o.uv = v.particleData; // Use this to control splat location (if needed)
                return o;
            }

            // Fragment outputs the splat contribution.
            fixed4 frag (v2f i) : SV_Target
            {
                // For demonstration, output 1 in red and some velocity in green and blue.
                return fixed4(1, 0.5, 0.5, 0);
            }
            ENDCG
        }
    }
}

Then in your C# code, you’d render your particles (potentially using a Graphics.DrawMeshInstanced or similar approach) into a RenderTexture configured for random write.

2. Custom Render Texture (CRT) Approach

What It Is:
A Custom Render Texture (CRT) in Unity is designed for simulations where you want to process the texture on every update. With a CRT, you can script how the texture updates over time using a material.

Why It’s Good:

Integrated Update: CRTs can automatically update themselves based on a given material.
Blending & Filtering: You can set the clear color, blending, and other parameters within the CRT settings.
Ease of Use: For certain fluid or particle simulations, CRTs hide some boilerplate code.

How to Do It:

Create a CRT asset in Unity.
Write an update material that performs the splatting pass with additive blending.
The CRT will then manage the iterative updates, which can be fed back into your simulation.

The downside is that CRTs are less flexible than a full compute shader solution if you need very fine-tuned control, but they work well for many effects.

3. Hybrid Approach

You might combine compute shaders and rendering passes: use a compute shader for updating particle state and a separate rendering pass with additive blending to accumulate the density field. This gives you both the low-level control of compute shaders for physics updates and the simplicity of blending for accumulation.

Conclusion

The Best Choice?
For your simulation, if your main challenge is to “splat” particle contributions into a grid without needing atomic add, using a RenderTexture with a material set to additive blending is a robust and straightforward method. This approach is simple, leverages hardware acceleration, and aligns well with the original WebGL implementation.

By choosing this strategy, you also reduce the complexity of handling atomic operations on floats in compute shaders, keeping your code clean and your simulation efficient.

Below is one complete example that uses both compute shaders to update particle physics and a separate pass with additive blending to “splat” particles into a density texture. In this implementation, we use a ComputeBuffer to hold particle data and update it via a compute shader, and then we render all particles as points using an accumulation shader (with additive blending) to build up the density field.

The full implementation consists of three parts:

C# Script (GravitySimulation.cs):
Sets up the particle buffer and RenderTexture, dispatches the compute shader to update particles, and then draws the particles with a material that uses additive blending into the density texture.
Accumulation Shader (AccumulationSplat.shader):
A simple vertex/fragment shader that draws each particle as a point. Because its pass uses additive blending, multiple particles overlapping in the same grid cell add their contributions.
Compute Shader (GravitySimulation.compute):
Contains a kernel that updates particle positions and velocities. (For simplicity, this example uses a constant “gravity” acceleration. In a more complex simulation you might replace this with a gradient from a computed potential field.)

1. GravitySimulation.cs

using UnityEngine;

public class GravitySimulation : MonoBehaviour
{
    [Header("Shader References")]
    [Tooltip("Compute shader that updates particle positions and velocities.")]
    public ComputeShader updateParticlesComputeShader;
    [Tooltip("Material using the accumulation shader to splat particles into a density texture.")]
    public Material accumulationMaterial;

    [Header("Simulation Settings")]
    [Tooltip("Total number of particles in the simulation.")]
    public int numParticles = 1024 * 1024; // 1,048,576 particles
    [Tooltip("Time step for simulation updates.")]
    public float dt = 0.001f;
    [Tooltip("Constant gravity acceleration (for demo purposes).")]
    public Vector2 gravity = new Vector2(0, -0.1f);

    [Header("Render Texture Settings")]
    [Tooltip("RenderTexture to hold the density accumulation.")]
    public RenderTexture densityTexture;

    // Compute buffer holding our particle data.
    ComputeBuffer particleBuffer;
    int updateKernel;

    // Particle structure: must match the structure in the compute shader and accumulation shader.
    struct Particle
    {
        public Vector2 position; // Normalized coordinates [0,1]
        public Vector2 velocity;
    }

    void Start()
    {
        // Initialize the particle buffer.
        particleBuffer = new ComputeBuffer(numParticles, sizeof(float) * 4);
        Particle[] particles = new Particle[numParticles];
        for (int i = 0; i < numParticles; i++)
        {
            // Randomize initial positions and small random velocities.
            particles[i].position = new Vector2(Random.value, Random.value);
            particles[i].velocity = new Vector2(Random.Range(-0.01f, 0.01f), Random.Range(-0.01f, 0.01f));
        }
        particleBuffer.SetData(particles);

        // Create the density RenderTexture.
        // Adjust the size as desired (here 512x512 grid is used).
        densityTexture = new RenderTexture(512, 512, 0, RenderTextureFormat.ARGBFloat);
        densityTexture.enableRandomWrite = true;
        densityTexture.Create();

        // Set up the compute shader.
        updateKernel = updateParticlesComputeShader.FindKernel("UpdateParticles");
        updateParticlesComputeShader.SetBuffer(updateKernel, "Particles", particleBuffer);
        updateParticlesComputeShader.SetFloat("dt", dt);
        updateParticlesComputeShader.SetVector("gravity", gravity);

        // Provide the particle buffer to the accumulation material.
        accumulationMaterial.SetBuffer("_Particles", particleBuffer);

        // Optionally clear the density texture (set to black) at the start.
        Graphics.SetRenderTarget(densityTexture);
        GL.Clear(false, true, Color.black);
        Graphics.SetRenderTarget(null);
    }

    void Update()
    {
        // 1. Update particle positions using the compute shader.
        int threadGroups = Mathf.CeilToInt(numParticles / 64.0f);
        updateParticlesComputeShader.Dispatch(updateKernel, threadGroups, 1, 1);

        // 2. Clear the density texture before splatting.
        Graphics.SetRenderTarget(densityTexture);
        GL.Clear(false, true, Color.black);
        Graphics.SetRenderTarget(null);

        // 3. Splat particles onto the density texture using the accumulation shader.
        // The accumulation shader is designed to use additive blending.
        accumulationMaterial.SetPass(0);
        Graphics.SetRenderTarget(densityTexture);
        // Draw all particles as points. Each particle’s vertex is generated procedurally.
        Graphics.DrawProceduralNow(MeshTopology.Points, numParticles);
        Graphics.SetRenderTarget(null);
    }

    void OnDestroy()
    {
        if (particleBuffer != null)
            particleBuffer.Release();
        if (densityTexture != null)
            densityTexture.Release();
    }
}

2. AccumulationSplat.shader

This shader “splats” (accumulates) each particle’s contribution into the RenderTexture. With additive blending enabled, overlapping particles add together.

Shader "Custom/AccumulationSplat"
{
    SubShader
    {
        Tags { "RenderType"="Opaque" }
        Pass
        {
            // Disable depth writes and culling.
            ZWrite Off
            Cull Off
            // Set up additive blending (source and destination are simply added together).
            Blend One One

            CGPROGRAM
            #pragma vertex Vert
            #pragma fragment Frag
            #include "UnityCG.cginc"

            // Structure matching our particle data.
            struct Particle
            {
                float2 position; // normalized coordinate between 0 and 1
                float2 velocity;
            };

            // Particle data passed in from the C# script.
            StructuredBuffer<Particle> _Particles;

            // Vertex shader: Each vertex corresponds to one particle.
            // We convert the normalized particle position [0,1] to clip space [-1,1].
            float4 Vert(uint id : SV_VertexID) : SV_POSITION
            {
                Particle p = _Particles[id];
                float2 clipPos = p.position * 2.0 - 1.0;
                return float4(clipPos, 0.0, 1.0);
            }

            // Fragment shader: Outputs a fixed splat value.
            // For example, here we output a density value of 1 in the red channel.
            fixed4 Frag() : SV_Target
            {
                return fixed4(1.0, 0.0, 0.0, 0.0);
            }
            ENDCG
        }
    }
    FallBack "Diffuse"
}

3. GravitySimulation.compute

This compute shader kernel updates each particle’s position and velocity. For demonstration, it simply applies a constant gravity acceleration. You can later expand it to include forces from a potential field computed via a Poisson solver.

#pragma kernel UpdateParticles

struct Particle {
    float2 position;
    float2 velocity;
};

// Read-write StructuredBuffer containing all particles.
RWStructuredBuffer<Particle> Particles;

// Simulation constants provided from C#.
float dt;
float2 gravity; // For example, a constant gravity vector

[numthreads(64, 1, 1)]
void UpdateParticles(uint id : SV_DispatchThreadID)
{
    // Guard against out-of-range thread IDs.
    if (id >= Particles.Length)
        return;

    Particle p = Particles[id];

    // Simple physics integration: update velocity with gravity and then position.
    p.velocity += gravity * dt;
    p.position += p.velocity * dt;

    // Wrap around: keep the particle within [0,1] by taking the fractional part.
    p.position = frac(p.position);

    Particles[id] = p;
}

How It Works

Initialization (C# Script):
The C# script creates a large particle buffer and a RenderTexture to hold the density accumulation field. It sets the particle buffer on both the compute shader (for physics updates) and the accumulation material (for drawing).
Compute Shader – Particle Update:
The compute shader kernel UpdateParticles updates each particle’s position (with simple gravity integration) and wraps positions to keep them in [0,1].
Accumulation Shader – Splatting:
The accumulation shader is used in a separate rendering pass. Here, each particle is drawn as a point; the vertex shader converts normalized positions to clip space, and the fragment shader simply outputs a constant density value (in this case, 1 in the red channel). Because the material’s blend mode is set to additive (Blend One One), any overlapping particles will have their contributions summed in the RenderTexture.
Render Flow (C# Update Method):
Each frame, the simulation first updates the particle states via the compute shader, then clears the density texture and renders all particles with additive blending into the density texture. The resulting texture accumulates the “splat” contributions and can then be used for further simulation (like solving a Poisson equation) or for display.

This complete implementation makes full use of Unity’s available tools and GPU blending to sidestep atomic float adds—and keeps the code clear and educational. Enjoy experimenting with and expanding the simulation!

Below is one way to “abstract away” the dimensionality from your density accumulation. The idea is to structure your code so that—in one branch—the particles get splatted into a 2D RenderTexture (using a material that blends additively) and—in the other branch you use a compute‐shader that writes into a 3D grid stored in a buffer (using an atomic–CAS loop to simulate atomic add on floats). That way your simulation code “just” calls your accumulator without caring whether you’re in 2D or 3D.

In the example below we define an enum that selects our accumulation mode. Then in our main C# script we allocate either a 2D RenderTexture or a 3D density buffer and use different shader passes accordingly. (For 3D we “flatten” our 3D grid into a 1D RWByteAddressBuffer so that we can perform our simulated atomic add.)

Below are three files that together illustrate one full implementation:

1. AbstractDensityAccumulator.cs

This C# script sets up a unified simulation. It creates a particle buffer, updates particle positions with a compute shader, and then “splats” particle contributions into the density field using one of two methods (2D via a material pass or 3D via a compute shader).

using UnityEngine;

public enum AccumulationDimension { TwoD, ThreeD }

public class AbstractDensityAccumulator : MonoBehaviour
{
    [Header("Simulation Settings")]
    public AccumulationDimension accumulationDimension = AccumulationDimension.TwoD;
    [Tooltip("Total number of particles in the simulation")]
    public int numParticles = 1024 * 1024;  // e.g., one million particles
    [Tooltip("Time step for particle update")]
    public float dt = 0.001f;

    [Header("2D Accumulation Settings")]
    [Tooltip("Material for 2D accumulation (additive blending)")]
    public Material accumulationMaterial2D;  // Uses a shader similar to AccumulationSplat2D.shader
    [Tooltip("Resolution (width and height) for the 2D density grid")]
    public int gridResolution2D = 512;

    [Header("3D Accumulation Settings")]
    [Tooltip("Compute shader for 3D accumulation")]
    public ComputeShader accumulationSplat3DComputeShader;
    [Tooltip("Resolution for the 3D density grid (gridResolution3D x gridResolution3D x gridResolution3D)")]
    public int gridResolution3D = 128;

    [Header("Compute Shader for Particle Update")]
    public ComputeShader updateParticlesComputeShader; // Shared for both 2D and 3D modes

    // Particle buffer
    ComputeBuffer particleBuffer;

    // Density storage
    RenderTexture densityTexture2D;       // Used in 2D mode
    RWByteAddressBuffer densityVolume3D;    // Used in 3D mode (wrapped in a ComputeBuffer)

    // Kernel IDs
    int updateKernel;
    int splat3DKernel;  // For 3D accumulation

    // Particle structure – note for 2D we use float2 positions, for 3D we use float3.
    // Here we use a generic structure that has both 2D and 3D components.
    struct Particle
    {
        public Vector3 position;  // When using 2D, we use x and y from this and ignore z.
        public Vector3 velocity;  // Similarly, for 2D only x and y are used.
    }

    void Start()
    {
        InitializeParticleBuffer();
        InitializeDensityStorage();
        SetupComputeShaders();
    }

    void InitializeParticleBuffer()
    {
        // Each particle has 6 floats (position and velocity, stored as Vector3 each)
        particleBuffer = new ComputeBuffer(numParticles, sizeof(float) * 6);
        Particle[] particles = new Particle[numParticles];
        for (int i = 0; i < numParticles; i++)
        {
            Particle p = new Particle();
            // For 2D, generate x,y uniformly and set z to zero.
            if (accumulationDimension == AccumulationDimension.TwoD)
            {
                p.position = new Vector3(Random.value, Random.value, 0f);
                p.velocity = new Vector3(Random.Range(-0.01f, 0.01f), Random.Range(-0.01f, 0.01f), 0f);
            }
            else  // For 3D, generate all three coordinates in [0,1]
            {
                p.position = new Vector3(Random.value, Random.value, Random.value);
                p.velocity = new Vector3(Random.Range(-0.01f, 0.01f),
                                         Random.Range(-0.01f, 0.01f),
                                         Random.Range(-0.01f, 0.01f));
            }
            particles[i] = p;
        }
        particleBuffer.SetData(particles);
    }

    void InitializeDensityStorage()
    {
        if (accumulationDimension == AccumulationDimension.TwoD)
        {
            // Create a 2D RenderTexture for density accumulation.
            densityTexture2D = new RenderTexture(gridResolution2D, gridResolution2D, 0, RenderTextureFormat.ARGBFloat);
            densityTexture2D.enableRandomWrite = true;
            densityTexture2D.Create();
        }
        else  // ThreeD
        {
            // For 3D accumulation, we allocate a buffer to hold gridResolution^3 float values.
            // We pack our density grid as a flattened array. Each cell is one float (4 bytes).
            int totalCells = gridResolution3D * gridResolution3D * gridResolution3D;
            // Create a ComputeBuffer with Raw type (or use a RWByteAddressBuffer in compute shader)
            // Here we use a ComputeBuffer of bytes; each cell takes 4 bytes.
            ComputeBuffer densityBuffer = new ComputeBuffer(totalCells, sizeof(float), ComputeBufferType.Raw);
            // We clear it initially by setting all elements to 0.
            float[] clearData = new float[totalCells];
            densityBuffer.SetData(clearData);
            // Wrap this ComputeBuffer as our RWByteAddressBuffer in our compute shader.
            // (Pass it to the compute shader as a buffer.)
            densityVolume3D = new RWByteAddressBufferWrapper(densityBuffer);
        }
    }

    // A tiny helper class to “wrap” a ComputeBuffer as a RWByteAddressBuffer.
    // (In practice you can pass a ComputeBuffer configured as Raw to a compute shader.)
    class RWByteAddressBufferWrapper : RWByteAddressBuffer
    {
        public RWByteAddressBufferWrapper(ComputeBuffer buffer)
        {
            nativeBuffer = buffer;
        }
    }

    void SetupComputeShaders()
    {
        // Set up the particle update compute shader.
        updateKernel = updateParticlesComputeShader.FindKernel("UpdateParticles");
        updateParticlesComputeShader.SetBuffer(updateKernel, "Particles", particleBuffer);
        updateParticlesComputeShader.SetFloat("dt", dt);

        // If in 3D mode, set up the accumulation compute shader.
        if (accumulationDimension == AccumulationDimension.ThreeD)
        {
            splat3DKernel = accumulationSplat3DComputeShader.FindKernel("SplatParticles3D");
            accumulationSplat3DComputeShader.SetBuffer(splat3DKernel, "Particles", particleBuffer);
            accumulationSplat3DComputeShader.SetInt("numParticles", numParticles);
            accumulationSplat3DComputeShader.SetInt("gridSize", gridResolution3D);
            // Bind the density volume. Our RWByteAddressBuffer is passed as a buffer.
            // (Assume our compute shader declares "RWByteAddressBuffer DensityVolume;")
            accumulationSplat3DComputeShader.SetBuffer(splat3DKernel, "DensityVolume", ((RWByteAddressBufferWrapper)densityVolume3D).nativeBuffer);
        }
    }

    void Update()
    {
        // 1. Update particle positions.
        int threadGroups = Mathf.CeilToInt(numParticles / 64.0f);
        updateParticlesComputeShader.Dispatch(updateKernel, threadGroups, 1, 1);

        if (accumulationDimension == AccumulationDimension.TwoD)
        {
            // Clear the 2D density texture.
            Graphics.SetRenderTarget(densityTexture2D);
            GL.Clear(true, true, Color.black);
            Graphics.SetRenderTarget(null);

            // 2a. Use the accumulation material (with additive blending) to draw all particles.
            accumulationMaterial2D.SetBuffer("_Particles", particleBuffer);
            accumulationMaterial2D.SetInt("_GridResolution", gridResolution2D);
            accumulationMaterial2D.SetPass(0);
            Graphics.SetRenderTarget(densityTexture2D);
            Graphics.DrawProceduralNow(MeshTopology.Points, numParticles);
            Graphics.SetRenderTarget(null);
        }
        else
        {
            // Clear the 3D density buffer.
            // (Typically you’d clear the ComputeBuffer by re-setting zeros.)
            int totalCells = gridResolution3D * gridResolution3D * gridResolution3D;
            float[] clearData = new float[totalCells];
            // Reuse the density ComputeBuffer from our wrapper.
            ((RWByteAddressBufferWrapper)densityVolume3D).nativeBuffer.SetData(clearData);

            // 2b. Use the compute shader for 3D splatting.
            accumulationSplat3DComputeShader.Dispatch(splat3DKernel, threadGroups, 1, 1);
        }
    }

    void OnDestroy()
    {
        if (particleBuffer != null)
            particleBuffer.Release();
        if (densityTexture2D != null)
            densityTexture2D.Release();
        // Also release the underlying ComputeBuffer for 3D density if needed.
    }
}

Note:
In the code above the helper class RWByteAddressBufferWrapper is a conceptual wrapper. In Unity you typically pass a ComputeBuffer configured as “Raw” directly into your compute shader by calling
computeShader.SetBuffer(kernel, "DensityVolume", densityBuffer);
so you can adapt this part to your project’s needs.

2. AccumulationSplat2D.shader

This shader is used only in the 2D case. It “splats” each particle’s contribution (here simply adding 1 to the red channel) into a 2D density RenderTexture. Because the blend mode is set to additive the contributions accumulate automatically.

Shader "Custom/AccumulationSplat2D"
{
    SubShader
    {
        Tags { "RenderType"="Opaque" }
        Pass
        {
            // Disable depth writing and culling.
            ZWrite Off
            Cull Off
            // Set up additive blending.
            Blend One One

            CGPROGRAM
            #pragma vertex Vert
            #pragma fragment Frag
            #include "UnityCG.cginc"

            struct Particle
            {
                float3 position; // Use only x and y if in 2D mode.
                float3 velocity;
            };

            StructuredBuffer<Particle> _Particles;
            int _GridResolution; // Resolution of the density grid.

            // Each vertex corresponds to one particle.
            // Transform normalized position [0,1] into clip space.
            float4 Vert(uint id : SV_VertexID) : SV_POSITION
            {
                Particle p = _Particles[id];
                // For 2D, use p.position.xy.
                float2 clipPos = p.position.xy * 2.0 - 1.0;
                return float4(clipPos, 0.0, 1.0);
            }

            // The fragment simply outputs a splat: here we add 1.0 density.
            fixed4 Frag() : SV_Target
            {
                return fixed4(1.0, 0.0, 0.0, 0.0);
            }
            ENDCG
        }
    }
    FallBack "Diffuse"
}

3. AccumulationSplat3D.compute

This compute shader handles 3D accumulation by “splatting” each particle’s contribution into a flattened density volume. Because we need to avoid race conditions when adding floats, we use an atomic compare–exchange loop. (This emulates atomic add on floats by reinterpreting the float as a uint.)

#pragma kernel SplatParticles3D

struct Particle
{
    float3 position; // Normalized coordinates in [0,1]
    float3 velocity;
};

StructuredBuffer<Particle> Particles;
RWByteAddressBuffer DensityVolume; // The flattened density grid.
int numParticles;
int gridSize; // The resolution along one axis of the 3D grid.

// Emulated atomic add for a float using a compare–exchange loop.
float AtomicAddFloat(inout RWByteAddressBuffer addr, uint byteOffset, float value)
{
    uint oldUint;
    uint newUint;
    float oldFloat;
    do {
        addr.InterlockedAdd(byteOffset, 0, oldUint); // Atomically read current value.
        oldFloat = asfloat(oldUint);
        float newFloat = oldFloat + value;
        newUint = asuint(newFloat);
    } while (addr.InterlockedCompareExchange(byteOffset, newUint, oldUint) != oldUint);
    return oldFloat;
}

[numthreads(64, 1, 1)]
void SplatParticles3D(uint id : SV_DispatchThreadID)
{
    if (id >= numParticles) return;

    Particle p = Particles[id];

    // Map the normalized position [0,1] to grid indices.
    int ix = (int)(p.position.x * gridSize);
    int iy = (int)(p.position.y * gridSize);
    int iz = (int)(p.position.z * gridSize);

    // Clamp indices.
    ix = clamp(ix, 0, gridSize - 1);
    iy = clamp(iy, 0, gridSize - 1);
    iz = clamp(iz, 0, gridSize - 1);

    // Flatten the 3D index: index = x + y*gridSize + z*gridSize*gridSize.
    int index = ix + iy * gridSize + iz * gridSize * gridSize;
    uint byteOffset = index * 4; // Each float occupies 4 bytes.

    // Add a contribution of 1.0 to the density value at this cell.
    AtomicAddFloat(DensityVolume, byteOffset, 1.0);
}

How It Works

Abstraction:
The C# script uses an enum to choose between 2D (rendering with a material pass that leverages additive blending) and 3D (using a compute shader with atomic operations) accumulation. Your particle update code remains the same regardless of the dimensionality.
2D Path:
In 2D mode the particles are “splatted” onto a 2D RenderTexture via a shader that draws each point in normalized screen space. The blending mode set to “Blend One One” ensures that contributions add up.
3D Path:
In 3D mode the particles are processed by a compute shader that computes a grid index from a 3D position and then adds a contribution atomically into a flattened buffer (emulated using an InterlockedCompareExchange loop).

Using this abstraction, you can easily switch between 2D and 3D simulations. The accumulation of density values will work “regardless” of whether your points live in 2D or 3D space. You can extend this approach to higher dimensions by adjusting the grid indexing and storage.

Happy coding and experimenting with your unified density accumulator!