Unity Physics Simulation Method Environment Hodgepodge

Comparison Overview

It was only when I was writing the title that I realized how miscellaneous my work was… Taking advantage of the need to select a paper for a school symposium, I chose PBD and XPBD, and also did a comparison between those and tradition PBA. I remembered that I had always wanted to compare the performance of C# and C++ during my internship but never did, so I wrote two comparisons of the semi-implicit method in Unity and C++ respectively (how can I miss parallelization). Here is a general comparison table first. The specific values don’t make much sense, but it’s still okay to look at the trend. (35x35 cloth grid, 2000 stiffness coefficient, 32 iterations, 64 threads, GTX1650 laptop. Parallelism is all lock-free multi-threading)

	Unity	Unity Job + Burst	C++	C++ Thread
Time (ms)	36.9	1.6	12.2	38.8

The performance of C++ Thread is a bit unexpected. I think it might be because Thread too costly? However, the performance of C++ is much better than Unity. I think I can try C++ when I want to optimize the code that is not suitable for Job + Burst parallelization the next time.

	Semi-Implicit Euler	Implicit Euler	PBD	XPBD
Main Thread (ms)	36.9	N/A	70.9	52.2
Job + Burst (ms)	1.6	N/A	2.1	2.2

Implicit crashes and not be taken into account. XPBD is better than PBD in the main thread may because I used Jacobi instead of Gauss-Seidel when writing PBD, but just look at the parallelization. What disappoint me is that Unity’s Cloth component only needs 0.2ms to simulate these, while my XPBD can only reach 0.4ms while maintaining stability. I will study how to optimize it later when I have time

Semi-Implicit Euler

I wanted to write an explicit one at the begining, but I haven’t written Euler for a long time, so after I finished I find it is a semi-implicit one. P is the grid points array, which is read in as local space. In order to be able to drag and observe the effect in Unity, it is converted to world space and then converted back after calculation.

void Semi_Implict()
{
    for (int i = 0; i < P.Length; i++) 
    { 
        P2[i] = transform.TransformPoint(P[i]); 
        V[i] += (P2[i] - P1[i]) / dt;
    }

    float idt = dt / iter, d = Mathf.Pow(damping, 1f / iter);
    for (int i = 0; i < iter; i++)
    {
        // Constrant
        for (int j = 0; j < C.Length; j++)
        {
            int x = C[j].x, y = C[j].y; float z = C[j].z;
            Vector3 f = k * (1 - z / (P2[x] - P2[y]).magnitude) * (P2[x] - P2[y]);
            F[x] -= f; F[y] += f;
        }

        // Gravity & Update & Calculate Local Position
        for (int j = 0; j < P.Length; j++)
        {
            V[j] = d * V[j] + idt * F[j];
            if (W[j] != 0) P2[j] += idt * V[j];
            P1[j] = P2[j]; F[j] = g;
        }
    }

    for (int i = 0; i < P.Length; i++) P[i] = transform.InverseTransformPoint(P2[i]);
}

Implicit Euler

There is nothing much to say, it is just gradient iteration. I tried Chebyshev acceleration but it doesn’t seem to make much difference. It is worth mentioning that in the previous version I wrote, I used a normalize and a magnitute to calculate the gradient. After changing to this version, the performance has improved a lot. Daily doubt the implementation of Unity functions.

void Implict()
{
    for (int i = 0; i < P.Length; i++)
    {
        V[i] += ((P2[i] = transform.TransformPoint(P[i])) - P1[i]) / dt;
        if (W[i] != 0) P[i] = P2[i] += (dt * (V[i] *= damping));
    }

    float invdt2 = 1 / dt / dt, factor = 1 / (invdt2 + 4 * k), rho2 = rho * rho, w = 1;
    for (int i = 0; i < iter; i++)
    {
        //Momentum and Gravity.
        for (int j = 0; j < P.Length; j++) G[j] = invdt2 * (P2[j] - P[j]) - g;

        //Spring Force.
        for (int j = 0; j < C.Length; j++)
        {
            int x = C[j].x, y = C[j].y; float z = C[j].z;
            Vector3 grad = k * (1 - z / (P2[x] - P2[y]).magnitude) * (P2[x] - P2[y]);
            G[x] += grad; G[y] -= grad;
        }

        if (i == 0) w = 1;
        else if (i == 1) w = 2 / (2 - rho2);
        else w = 4 / (4 - w * rho2);

        //Update P2 by gradient.
        for (int j = 0; j < P.Length; j++)
        {
            if (W[j] == 0) continue;
            Vector3 old = P2[j];
            P2[j] = w * (P2[j] - factor * G[j]) + (1 - w) * P1[j];
            P1[j] = old;
        }
    }

    for (int i = 0; i < P.Length; i++) 
    { 
        V[i] += (P2[i] - P[i]) / dt; 
        P[i] = transform.InverseTransformPoint(P2[i]); 
    }
}

PBD

I won’t go into the formulas and algorithm flow, but the general idea is to first move freely, then impose constraints. For a detailed paper interpretation, you can read this jeffzzz: Position Based Dynamics Notes. But PBD actually has some problems. The biggest problem is that it is not so physical. The physical properties it displays are only linked to the time step and the number of iterations. But the time step is related to performance, and the number of iterations is related to performance. It is a landmine to give different time steps and iterations to different objects in order to change the physical properties…

void PBD()
{
    for (int i = 0; i < P.Length; i++) P2[i] = transform.TransformPoint(P[i]);

    // Free Update
    for (int i = 0; i < P.Length; i++)
    {
        if (W[i] == 0) continue;
        V[i] = damping * V[i] + dt * g;
        P2[i] += dt * V[i] + damping * (P2[i] - P1[i]);
    }

    for (int i = 0; i < iter; i++)
    {
        // Strain Limit
        for (int j = 0; j < C.Length; j++)
        {
            int x = C[j].x, y = C[j].y; float z = C[j].z;
            Vector3 p = z * (P2[x] - P2[y]).normalized;
            SX[x] += 0.5f * (P2[x] + P2[y] + p); SN[x]++; SX[y] += 0.5f * (P2[x] + P2[y] - p); SN[y]++;
        }
        for (int j = 0; j < P.Length; j++)
        {
            if (W[j] == 0) continue;
            P1[j] = P2[j];
            P2[j] = (0.2f * P2[j] + SX[j]) / (0.2f + SN[j]);
            V[j] += (P2[j] - P1[j]) / dt;
            SX[j] = Vector3.zero; SN[j] = 0;
        }
    }
    
    for (int i = 0; i < P.Length; i++) P[i] = transform.InverseTransformPoint(P1[i] = P2[i]);
}

XPBD

XPBD redefines PBD constraints by introducing elastic potential energy, solving the above PBD problem with less than 5% increase in time (the paper seems to say 2-3%, but it’s roughly the same). The recommended paper interpretation here is zilch: XPBD Paper Interpretation, or you can also take a look at WebXPBDCloth written by the founder Müller. It will be faster to understand it in combination with actual engineering.

void XPBD()
{
    float idt = dt / iter, a = 1f / k / idt / idt, d = Mathf.Pow(damping, 1f / iter);
    for (int i = 0; i < P.Length; i++) { V[i] += ((P2[i] = transform.TransformPoint(P[i])) - P1[i]) / idt; }

    for (int i = 0; i < iter; i++)
    {
        // Free Update
        for (int j = 0; j < P.Length; j++)
        {
            if (W[j] == 0) continue;
            V[j] = d * V[j] + idt * g;
            P2[j] += idt * V[j];
        }

        // Solve Constraint (Gauss Seidel)
        for (int j = 0; j < C.Length; j++)
        {
            int x = C[j].x, y = C[j].y; float z = C[j].z, w = W[x] + W[y];
            Vector3 grad = -(1 - z / (P2[x] - P2[y]).magnitude) * (P2[x] - P2[y]) / (a + w);
            P2[x] += grad * W[x]; P2[y] -= grad * W[y];
        }

        // Record Velocity
        for (int j = 0; j < P.Length; j++) { V[j] = (P2[j] - P1[j]) / idt; P1[j] = P2[j]; }
    }

    for (int i = 0; i < P.Length; i++) P[i] = transform.InverseTransformPoint(P2[i]);
}