Howdy, Stranger!

It looks like you're new here. If you want to get involved, click one of these buttons!

Compute Shader for Cloth Simulation (PBD)

NixonNixon Member
edited May 20 in Programming

Hi everyone! Is there anyone making cloth simulations using compute shader? I am currently on journey to use GPU to simulate cloth and soft body physics. Does anyone has any solutions or guides on how this can be achieved? I have searched online for alot of solutions but just can't come up with a decent solution. Alot of people used mass spring cloth simulation but it is actually really expensive to simulate as you need to dispatch a kernel a few hundred times per FixedUpdate() to get a decent simulation.

I wanna achieve sth like this but unfortunately the person who made this, Jim Hugunin, did not post any links or open source it...


  • jtok4jjtok4j Member


    Regarding, mass spring cloth simulation, perhaps this might be a job for the Burst Compiler with unity ESC/JOBS, only run on the GPU.


    Burst and Synchronouscompilation:

     Keep on Creating! 

    Justin of JustinTime Studios ( 

  • NixonNixon Member
    edited May 22

    Thanks @jtok4j but isn't burst compiler for platform specific only? And i thought mass spring cloth simulation is a very old technique already, is there any more recent techniques that can just runs on the compute shader?

  • NixonNixon Member

    Hi everyone hope you’re safe and healthy. Recently, I am making a cloth simulation inside unity using compute shader. Through my weeks of research, I came to a simple solution called mass spring solution, where every vertices of the cloth is a particle and each particle is held together with an invisible spring. This solution is great but it is very expensive to calculate because it requires a lot of iterations to achieve a realistic result. Think about it, if our cloth is a 5x5 vertex cloth and the top row of the vertices is fixed. When the simulation starts, all the vertices that is not fixed will immediate fall due to gravity. The first to notice the fixed vertices will be the second row which will stop falling by the second simulation step and the third at the third simulation step. What I want to bring here is that if we use this model and our cloth has about 1k rows of vertices, then we need 1k simulation step to simulate a non - stretchy cloth (which is what real world cloth looks like. ). So I’m here to ask if there are any different approach to this problem so that it is not expensive to compute and behaves realistically.  Like this video: 

    Thanks for reading guys!

  • I know nothing about cloth physics but i did see a video about someone using machine learning or neural networks to predict or emulate cloth physics.i know its a long shot but i hear that it was very performant..sorry wish i could help more

  • NixonNixon Member

    @TumiTheSoleMaker Thanks for the reply! Is it from 2 minutes paper video? I have seen it too, but before training AI simulated cloth we need to have the right simulation first... I have read alot of research paper on PBD for cloth simulation but unfortunately I dun fully understand them as I don't have the knowledge especially when they write out the math equation...

  • NixonNixon Member
    edited May 23

    Hi guys, I think I should at least post my most recent solution for cloth simulation to give you guys an idea on where I am right now. This is my mass spring cloth simulation: (this is the compute shader code)// Each #kernel tells which function to compile; you can have many kernels

    #pragma kernel Cloth
    float3 gravity;
    uint totalParticles;
    float particleMass;
    float invertParticleMass;
    uint springK;
    float maxStretch;
    float damping;
    float deltaT;
    struct Neighbor
      uint n1;
      uint n2;
      uint n3;
      uint n4;
      uint n5;
      uint n6;
      uint n7;
      uint n8;
    struct Dist
      float d1;
      float d2;
      float d3;
      float d4;
      float d5;
      float d6;
      float d7;
      float d8;
    struct NeighborVert
      uint origin;
      Neighbor neighbor;
      Dist dist;
    RWStructuredBuffer<NeighborVert> neighborVerts;
    RWStructuredBuffer<uint> controlled;
    RWStructuredBuffer<float3> pos;
    RWStructuredBuffer<float3> vel;
    float3 SpingForce(uint origin, uint neighbor, float dist)
      float3 originPos = pos[origin];
      float3 neighborPos = pos[neighbor];
      float3 r = neighborPos - originPos;
      float3 force = normalize(r) * springK * (length(r) - dist);
      return force;
    [numthreads(8, 1, 1)]
    void Cloth (uint3 id : SV_DispatchThreadID)
      uint idx = id.x;
      NeighborVert nv = neighborVerts[idx];
      float3 f = gravity * particleMass;
      float3 a;
      float3 v = vel[idx];
      float3 p = pos[idx];
      if (controlled[idx] == 0)
        if (nv.neighbor.n1 != -1) f += SpingForce(nv.origin, nv.neighbor.n1, nv.dist.d1);
        if (nv.neighbor.n2 != -1) f += SpingForce(nv.origin, nv.neighbor.n2, nv.dist.d2);
        if (nv.neighbor.n3 != -1) f += SpingForce(nv.origin, nv.neighbor.n3, nv.dist.d3);
        if (nv.neighbor.n4 != -1) f += SpingForce(nv.origin, nv.neighbor.n4, nv.dist.d4);
        if (nv.neighbor.n5 != -1) f += SpingForce(nv.origin, nv.neighbor.n5, nv.dist.d5);
        if (nv.neighbor.n6 != -1) f += SpingForce(nv.origin, nv.neighbor.n6, nv.dist.d6);
        if (nv.neighbor.n7 != -1) f += SpingForce(nv.origin, nv.neighbor.n7, nv.dist.d7);
        if (nv.neighbor.n8 != -1) f += SpingForce(nv.origin, nv.neighbor.n8, nv.dist.d8);
        // apply damping force
        f += -damping * v;
        // calculate acceleration using Newton's second law (f=ma)
        a = f * invertParticleMass;
        v += a * deltaT;
        p += v * deltaT;
        pos[idx] = p;
        vel[idx] = v;

    As you can see, all particles is controlled by a spring connected to each neighbor particles. The stiffness of the spring is defined by a spring constant springK

    To achieve a stable result for a high resolution mesh, you will need a very small deltaT and you will need to Dispatch this kernel alot of times.

  • I will try to make your problem my problem so i am going to do some research..and let you know any useful info or techniques to optimize

  • NixonNixon Member

    Hi @TumiTheSoleMaker , thanks for engaging in my problem, I've done more research after I posted my source code, here are some resources that will help you:

    1. Basics of Position Based Dynamics:
    2. A better solver for PBD, Hierarchical Position Based Dynamics research paper by Nvidia:
    3. A much better solver, XPBD which solves the "well known limitation is that PBD’s behavior, time step and iteration count of the simulation changes the stiffness of the mesh" (also by NVIDIA)

    During your research you will come into two different popular solver method, Gauss-Seidel and Jacobi. Generally, the Gauss-Seidel method will be a much more better solver but it will bring bias to our solver as it depends on where you start to solve it first. But the good news is, for the sake of our use case, we're going to use Jacobi solver method because we are solving it parallel on a GPU.

    1. Comparison between both method:

    I hope these help you alot in your progress of research! Cheers for reading!

  • NixonNixon Member
    edited May 25

    Hey thanks for replying! In the meantime, I also have been research for quite awhile already so I have come up with a solution proposed on a NVIDIA paper about Position Based Dynamics, let me briefly explain it here and once I got the code up and running I might post it here or at least post the base code XD.

    In PBD, we directly manipulate the position of the particles when solving constraints. The easiest example is the distance constraint:

    let m = mass of particle, d = rest length, w = 1/m, = position of particle

    Imagine we have 2 particles constraint based on a distance d, then we can formulate 2 simple equation:

    1. delta x1 = 1/2 * (distance(x1, x2) - d) * -normalize(x1 - x2)
    2. delta x2 = 1/2 * (distance(x1, x2) - d) * normalize(x1 - x2)

    we divide into half by multiplying 1/2 because the constraint needs to be distributed to both particles, next we multiply the difference between particle distance and rest length to calculate how much displacement is needed in order to satisfy the distance constraint. Next we multiply the normalized difference of particles to get the direction of where the particle should go in 3D space. Hope this make sense!

Sign In or Register to comment.