Master Thesis about Real-time Liquid Simulation with Unity DOTS

This master thesis project investigates whether Position-Based Material Point Method can be implemented as a practical real-time liquid solver in Unity DOTS. The other part of this research is to determine the feasability in modern games to use a solution like this for liquid simulations during run time.

The work follows Lewin's PB-MPM formulation and his WebGPU prototype closely at the algorithmic level, but the Unity version changes two important constraints: it runs on the CPU rather than the GPU, and it extends the method from 2D to 3D. The thesis implementation is intentionally scoped to the liquid branch of PB-MPM so the evaluation can focus on solver architecture, particle-grid transfers, collision handling, and measured runtime behavior.

Important Links

Software Architecture

The central runtime system isPBMPMSolverSystem. It owns the fixed-size solver step scheduling, runs the PB-MPM stages as separate jobs, rebuilds the collider cache used inside those jobs, and applies a final smoothing pass for presentation. Rather than attaching the simulation to Unity's default fixed physics loop, the system accumulates frame time manually and executes the required number of solver substeps at the configured solver frequency.

A solver step advances the simulation by a fixedDeltaTime. Inside one solver step, the PB-MPM loop can be executed multiple times before the particles are finally integrated. That separation between solver frequency and iteration count is important in practice, because it makes the runtime cost of stability tuning visible and measurable instead of hiding it inside one monolithic update.

private void ScheduleParticleJobs(
    ref SystemState state,
    int iterationCount,
    GridInterpolationMode interpolationMode,
    bool useGridVolumePreservation)
{
    var gridEntities = _gridQuery.ToEntityArray(Allocator.TempJob);

    for (int substepIndex = 0; substepIndex < _solverSubsteps; substepIndex++)
    {
        for (int iterationIndex = 0; iterationIndex < iterationCount; iterationIndex++)
        {
            _gridIterationId++;
            RunSolverIteration(
                ref state,
                gridEntities,
                interpolationMode,
                useGridVolumePreservation,
                applyParticleGravity: iterationIndex == iterationCount - 1);
        }

        ScheduleIntegrateParticles(ref state);
    }

    ScheduleSmoothing(ref state);
    state.Dependency = gridEntities.Dispose(state.Dependency);
}

private void RunSolverIteration(
    ref SystemState state,
    NativeArray<Entity> gridEntities,
    GridInterpolationMode interpolationMode,
    bool useGridVolumePreservation,
    bool applyParticleGravity)
{
    ScheduleSolveConstraints(ref state, useGridVolumePreservation);
    ScheduleParticleToGrid(ref state, gridEntities, interpolationMode);
    ScheduleUpdateGrid(ref state, gridEntities, interpolationMode);
    ScheduleGridToParticle(
        ref state,
        gridEntities,
        interpolationMode,
        useGridVolumePreservation,
        applyParticleGravity);
}

Core runtime data

The architecture stays close to the PB-MPM paper, but the data layout is adapted to Unity's ECS model. Instead of storing the whole simulation in one object graph, the runtime is organized around a small set of explicit data containers:

• Particles are ECS entities with a ParticleComponent storing physical state such as position, displacement, deformation gradient, deformation displacement, liquid density, and velocity.
• The background grid is a separate ECS entity with a GridComponent and a DynamicBuffer for GridCells containing the per-node mass, volume, displacement, and last-touched iteration.
• Collision data is copied from Unity Physics box colliders into a compact native collider cache before the solver jobs run, so the jobs do not need to query scene objects directly.
• Particle initialization is handled by a dedicated spawn pipeline that turns authored spawn volumes into runtime ECS particle entities once at startup.

One detail that mattered a lot here is that the grid topology is baked once and then reused. The grid entity stores its global dimensions in GridComponent, while the dynamic buffer holds one GridCell entry per grid node. That keeps the memory layout dense and makes it straightforward for jobs to iterate over the grid without constantly rebuilding topology information.

public override void Bake(GridAuthoring authoring)
{
    var entity = GetEntity(authoring, TransformUsageFlags.Dynamic);
    AddComponent(entity, authoring.CreateGridComponent());
    AddComponent(entity, new GridAuthoringReference
    {
        AuthoringInstanceId = authoring.GetInstanceID()
    });

    var gridCells = AddBuffer<GridCell>(entity);
    authoring.RebuildGridCells(gridCells);
}

public void RebuildGridCells(DynamicBuffer<GridCell> gridCells)
{
    GridComponent grid = CreateGridComponent();
    int3 cellCounts = GridUtilities.GetCellCounts(grid);
    int3 nodeCounts = cellCounts + 1;

    gridCells.Clear();
    gridCells.EnsureCapacity(nodeCounts.x * nodeCounts.y * nodeCounts.z);

    for (int x = 0; x < nodeCounts.x; x++)
    {
        for (int y = 0; y < nodeCounts.y; y++)
        {
            for (int z = 0; z < nodeCounts.z; z++)
            {
                gridCells.Add(new GridCell
                {
                    Coordinates = new int3(x, y, z),
                    WeightedDisplacement = float3.zero,
                    Displacement = float3.zero,
                    Mass = 0f,
                    Volume = 0f,
                    LastTouchedIteration = 0
                });
            }
        }
    }

    gridCells.TrimExcess();
}

Authoring And Data Preparation

The simulation can still be authored from regular Unity scene objects. Bakers convert those authoring objects into compact ECS data before runtime. For the grid, that means turning a configured simulation volume into one grid entity plus its node buffer. For particle sources, the baker stores a compact spawn description rather than every particle position explicitly.

The spawn pipeline only supports a box volume in the current version, but the data representation is already set up in a way that scales cleanly: a spawn component stores counts per axis, a local start position, and a local step size. TheParticleSpawnSystem reconstructs the actual particle positions from that compact description, writes the initial transform and ParticleComponent values, computes the particle volume from the local step size, and then disables itself after the one-time initialization pass.

Solver Stages

Each solver iteration follows the same stage order as the PB-MPM paper: solve constraints, transfer particle state to the grid, update grid nodes, transfer the corrected state back to the particles, and finally integrate particle state after all iterations of the current solver step are complete. Splitting the pipeline into separate jobs made it possible to profile each stage in isolation and adjust the data flow where Unity's threading model demanded it.

Liquid constraint projection

The constraint stage implements the liquid branch of PB-MPM. For every particle, the solver reads the current deformation displacement tensor, splits it into volumetric and shear parts, and applies hydrostatic and viscous corrections directly to that tensor. The implementation can optionally use a grid-measured density value instead of relying only on the particle's local deformation history, which improves volume preservation when the iterative solve has not fully converged.

private void Execute(ref ParticleComponent particle)
{
    float3x3 D = particle.DeformationDisplacement;
    float c = D.c0.x + D.c1.y + D.c2.z;
    float deformationVolume =
        math.max(math.determinant(particle.DeformationGradient), 1e-4f);
    float particleLocalLiquidDensity = 1f / deformationVolume;
    float liquidDensity = UseGridVolumePreservation
        ? math.max(particle.GridMeasuredLiquidDensity, 1e-4f)
        : math.max(particleLocalLiquidDensity, 1e-4f);

    float3x3 dHydro = float3x3.identity * (liquidDensity - 1f - c);
    float3x3 dVisc = -(D - c * float3x3.identity);

    D += particle.LiquidHydroFactor * dHydro;
    D += particle.LiquidViscosityFactor * dVisc;

    particle.DeformationDisplacement = D;
    particle.LiquidDensity = liquidDensity;
}

This is the point where the project differs most clearly from a force-accumulation solver. Material behavior is expressed through local constraint corrections, while the particle-grid transfers are responsible for propagating those corrections across the simulation domain.

Two-stage particle-to-grid transfer

The particle-to-grid transfer is the most important concurrency problem in the whole solver. Every particle influences a quadratic B-spline support of 3 x 3 x 3 grid nodes, so naive parallel writes would immediately collide. The solution used in the thesis is a tiled two-stage transfer.

The first job does not write to the grid at all. It determines which grid tiles a particle touches and stores a compact particle record for each of those tiles. A second job then processes the active tiles independently, clears the cells inside the current tile, and accumulates only the particle records assigned to that tile. This is where mass, volume, and the affine displacement-related term are written into the grid.

public void Execute(int index)
{
    int tileIndex = ActiveTiles[index];
    int3 nodeCounts = GridUtilities.GetNodeCounts(Grid);
    int3 tileCounts = GridUtilities.GetTileCounts(nodeCounts, TileSize);
    int3 tileCoordinates = GridUtilities.GetTileCoordinates(tileCounts, tileIndex);
    int3 tileMin = tileCoordinates * TileSize;
    int3 tileMax = math.min(tileMin + TileSize, nodeCounts);

    ClearTile(tileMin, tileMax, nodeCounts);

    NativeParallelMultiHashMapIterator<int> iterator;
    GridTileParticleRecord particleRecord;
    if (!TileParticles.TryGetFirstValue(tileIndex, out particleRecord, out iterator))
    {
        return;
    }

    do
    {
        AccumulateParticle(ref particleRecord, tileMin, tileMax, nodeCounts);
    } while (TileParticles.TryGetNextValue(out particleRecord, ref iterator));
}

private void AccumulateParticle(
    ref GridTileParticleRecord particle,
    int3 tileMin,
    int3 tileMax,
    int3 nodeCounts)
{
    if (!GridUtilities.TryGetQuadraticWeights(
            Grid,
            particle.Position,
            out QuadraticWeights3D weights))
    {
        return;
    }

    for (int xOffset = 0; xOffset < 3; xOffset++)
    {
        for (int yOffset = 0; yOffset < 3; yOffset++)
        {
            for (int zOffset = 0; zOffset < 3; zOffset++)
            {
                int3 offset = new int3(xOffset, yOffset, zOffset);
                int3 supportCoordinates = weights.BaseCoordinate + offset;
                if (!GridUtilities.IsInsideNodeCounts(nodeCounts, supportCoordinates) ||
                    !GridUtilities.IsInsideTile(tileMin, tileMax, supportCoordinates))
                {
                    continue;
                }

                float weight = weights.GetWeight(offset);
                if (weight <= 0f)
                {
                    continue;
                }

                int gridIndex = GridUtilities.GetGridIndex(
                    nodeCounts,
                    supportCoordinates.x,
                    supportCoordinates.y,
                    supportCoordinates.z);
                float3 supportPosition = GridUtilities.GetGridPosition(
                    Grid,
                    supportCoordinates,
                    InterpolationMode);
                float3 relativePosition = supportPosition - particle.Position;

                GridCell cell = GridCells[gridIndex];
                cell.Mass += weight * particle.Mass;
                cell.WeightedDisplacement +=
                    weight *
                    particle.Mass *
                    (particle.Displacement +
                     math.mul(particle.DeformationDisplacement, relativePosition));
                cell.Volume += weight * particle.Volume;
                GridCells[gridIndex] = cell;
            }
        }
    }
}

That design keeps the implementation close to the PB-MPM data flow while fitting Unity's multithreaded execution model. Importantly, it also means inactive tiles and untouched cells can be skipped later in the pipeline instead of paying a full-grid cost on every iteration.

Grid update and collisions

After accumulation, each active grid node holds mass, volume, and mass-weighted displacement rather than a final nodal state.GridUpdateJob normalizes the accumulated displacement by nodal mass, reconstructs the node's world-space support position, and resolves collisions directly on the grid.

Collision handling is implemented against a cached set of box colliders built from Unity Physics data before the jobs start. For each collider, the displaced support position is tested and the penetrating normal component is projected away. TheLastTouchedIteration marker is used here to skip stale cells, which keeps the update tied to the actually active parts of the grid.

public void Execute(int index)
{
    GridCell cell = GridCells[index];

    if (cell.LastTouchedIteration != CurrentGridIteration || cell.Mass <= 0f)
    {
        return;
    }

    float3 supportPosition = GridUtilities.GetGridPosition(
        Grid,
        cell.Coordinates,
        InterpolationMode);

    cell.Displacement = cell.WeightedDisplacement / cell.Mass;

    foreach (var collider in Colliders)
    {
        float3 displacedSupportPosition = supportPosition + cell.Displacement;
        GridBoxCollider.CollisionResult collision =
            collider.Collide(displacedSupportPosition);
        if (collision.Collides)
        {
            float gap = math.min(
                0f,
                math.dot(collision.Normal,
                         collision.PointOnCollider - supportPosition));
            float penetration =
                math.dot(collision.Normal, cell.Displacement) - gap;
            float radialImpulse = math.max(penetration, 0f);
            cell.Displacement -= radialImpulse * collision.Normal;
        }
    }

    GridCells[index] = cell;
}

Grid-to-particle reconstruction

Once the grid contains corrected nodal displacements, the solver transfers that state back to the particles. The same quadratic B-spline neighborhood is sampled again, but now it is used to reconstruct three different quantities: the particle displacement, an affine local deformation term, and a grid-measured liquid density value.

The affine term is written back in displacement form asDeformationDisplacement, so the later integration step can update the deformation gradient. Gravity is added only on the final iteration of the current solver step. Applying it during every iteration would effectively multiply gravity by the iteration count.

for (int xOffset = 0; xOffset < 3; xOffset++)
{
    for (int yOffset = 0; yOffset < 3; yOffset++)
    {
        for (int zOffset = 0; zOffset < 3; zOffset++)
        {
            int3 offset = new int3(xOffset, yOffset, zOffset);
            int3 supportCoordinates = weights.BaseCoordinate + offset;
            if (!GridUtilities.IsInsideNodeCounts(nodeCounts, supportCoordinates))
            {
                continue;
            }

            float weight = weights.GetWeight(offset);
            if (weight <= 0f)
            {
                continue;
            }

            int cellIndex = GridUtilities.GetGridIndex(
                nodeCounts,
                supportCoordinates.x,
                supportCoordinates.y,
                supportCoordinates.z);
            GridCell cell = gridCells[cellIndex];
            if (cell.LastTouchedIteration != CurrentGridIteration)
            {
                continue;
            }

            float3 supportPosition = GridUtilities.GetGridPosition(
                grid,
                supportCoordinates,
                InterpolationMode);
            float3 relativePosition = supportPosition - particle.Position;

            displacement += weight * cell.Displacement;
            c += 4f * GridUtilities.OuterProduct(cell.Displacement, relativePosition) *
                (weight * inverseCellSizeSquared);
            gridMeasuredLiquidDensity +=
                weight * cell.Volume * inverseCellVolume;
        }
    }
}

particle.Displacement = displacement;
if (ApplyParticleGravity)
{
    particle.Displacement += GravityDisplacement;
}

particle.DeformationDisplacement = c * DeltaTime;
particle.GridMeasuredLiquidDensity = math.max(gridMeasuredLiquidDensity, 1e-4f);
if (UseGridVolumePreservation)
{
    particle.LiquidDensity = particle.GridMeasuredLiquidDensity;
}

Particle integration and visual smoothing

After the iteration loop finishes for the current solver step, particle state is integrated. The implementation advances particle position by the reconstructed displacement, updates the deformation gradient with I + D, performs a final particle-level collision correction, reconstructs velocity from displacement and time step, and writes the result back toLocalTransform.

private void Execute(
    ref ParticleComponent particle,
    ref LocalTransform transform)
{
    particle.Position += particle.Displacement;

    particle.DeformationGradient = math.mul(
        float3x3.identity + particle.DeformationDisplacement,
        particle.DeformationGradient);

    foreach (var collider in Colliders)
    {
        GridBoxCollider.CollisionResult collision =
            collider.Collide(particle.Position);
        if (collision.Collides)
        {
            particle.Displacement -= collision.Penetration * collision.Normal;
            particle.Position -= collision.Penetration * collision.Normal;
        }
    }

    particle.Velocity = particle.Displacement / DeltaTime;
    transform.Position = particle.Position;
}

The final smoothing pass is intentionally separated from the physical solve. When the solver frequency is low, for example 15 Hz, the physical state advances in visible steps. The smoothing system updates only the rendered transform every frame using the current particle velocity so the simulation remains visually smoother without altering the actual solver state.

Runtime Controls And Profiling

A large part of the thesis work was building the prototype in a way that supports profiling rather than only demonstration. The Unity debug overlay exposes solver frequency, iteration count, interpolation mode, grid volume preservation, and several scene parameters at runtime. That made it practical to move back and forth between implementation changes and controlled benchmark runs.

Key runtime measurements from the thesis:

• Baseline: 31,250 particles, 1 solver iteration, 137,700 grid cells, and a 30 Hz solver frequency produced a mean PBMPMSolverSystem update time of 1.08 ms.
• Running the same baseline configuration at 60 Hz raised the mean update time to 1.98 ms.
• At 60 Hz, increasing the particle count from 31,250 to 54,000 raised the solver time to 4.57 ms.
• At 30 Hz, increasing the particle count to 85,750 raised the solver time to 5.71 ms.
• At 30 Hz, increasing solver iterations from 1 to 2, 4, and 8 produced 1.90 ms, 4.33 ms, and 22.88 ms respectively.

The results show that particle count is the dominant scaling factor and stays close to linear within the tested range, while solver iterations are the most expensive stability control. Grid resolution also matters, but less strongly than particle count. Very coarse grids were not automatically faster either, because more particles ended up interacting through the same cells and increased the transfer cost.

Current Limits And Next Steps

The current prototype is intentionally narrow. It implements only the liquid branch of PB-MPM, resolves collisions only against box colliders, and renders the fluid as individual particles rather than a reconstructed surface. The most direct follow-up work is support for elastic materials and sand, broader collider types, GPU acceleration for the particle-to-grid stages, and a surface reconstruction and shading pass that presents the liquid as a continuous volume instead of debug spheres.

Within that scope, the thesis shows that a PB-MPM liquid solver can be mapped onto Unity DOTS in a way that is explicit, profileable, and still compatible with real-time budgets when fluid behavior is meant to support gameplay rather than act as a purely cosmetic effect.