// Real GPU-vs-CPU physics benchmark for COMPLEX/COMPOUND factory-part bodies. // Measures per-step wall-clock for: (1) per-body semi-implicit Euler integration and // (2) an all-pairs box-box contact impulse pass (the O(n) integration + O(contacts) solve // that dominate a compound-body sub-step). Same seeded scene + same math on GPU (CUDA) and // CPU (single-thread), across body counts. Reports cpuMs, gpuMs, speedup, and mean penetration // error (a correctness proxy — GPU vs CPU must match within tolerance). HONEST measured only. // // Build: nvcc -O3 -arch=native -o bench_gpu bench/bench.cu // Run: ./bench_gpu --json results/gpu_results.json #include #include #include #include #include #include #include struct Body { float px,py,pz, vx,vy,vz, hx,hy,hz, invMass; }; // A compound "factory part" = a cluster of child boxes; here we flatten each compound body's // bounding box into one collider per body for the contact pass, but give it a randomized // heavier mass + non-uniform extents (brackets/shafts) so it is NOT a uniform box. static unsigned lcg(unsigned &s){ s = s*1664525u + 1013904223u; return s; } static float frand(unsigned &s){ return (lcg(s) & 0xFFFFFF) / float(0x1000000); } void initScene(std::vector& b, int n, unsigned seed){ b.resize(n); unsigned s = seed; for(int i=0;i separated // Minimum translation axis float nx=0,ny=0,nz=0,depth; if(ox<=oy&&ox<=oz){ nx=dx<0?-1:1; depth=ox; } else if(oy<=oz){ ny=dy<0?-1:1; depth=oy; } else { nz=dz<0?-1:1; depth=oz; } float rvx=B.vx-A.vx, rvy=B.vy-A.vy, rvz=B.vz-A.vz; float vn=rvx*nx+rvy*ny+rvz*nz; float im=A.invMass+B.invMass; if(im<1e-9f) return; float bias=0.2f*fmaxf(depth-0.001f,0.0f)*60.0f; float j=-(vn - bias)/im; if(j<0)j=0; float ix=nx*j, iy=ny*j, iz=nz*j; A.vx-=ix*A.invMass; A.vy-=iy*A.invMass; A.vz-=iz*A.invMass; B.vx+=ix*B.invMass; B.vy+=iy*B.invMass; B.vz+=iz*B.invMass; } double stepCPU(std::vector& b, int iters, float dt, double* penOut){ int n=b.size(); // integrate velocities for(int i=0;i0) b[i].vy += -9.81f*dt; } // contact solve (Gauss-Seidel), all pairs, iters passes double penSum=0; int penCnt=0; for(int it=0; it0&&oy>0&&oz>0){ solvePairCPU(b[i],b[j]); if(it==iters-1){ penSum+=fminf(fminf(ox,oy),oz); penCnt++; } } } } // ground contact + integrate positions for(int i=0;i0){ x.vy += 0.2f*gpen*60.0f; if(x.vy<0)x.vy=0; } x.px+=x.vx*dt; x.py+=x.vy*dt; x.pz+=x.vz*dt; } if(penOut) *penOut = penCnt? penSum/penCnt : 0.0; return 0; } // ---- GPU kernels ---- __global__ void integrateVel(Body* b, int n, float dt){ int i=blockIdx.x*blockDim.x+threadIdx.x; if(i>=n)return; if(b[i].invMass>0) b[i].vy += -9.81f*dt; } // Jacobi-style contact solve: each thread = one body, accumulate impulses vs all others, // apply at end (order-independent — the natural GPU batching; slightly softer than Gauss-Seidel). __global__ void solveContactsJacobi(const Body* in, Body* out, int n){ int i=blockIdx.x*blockDim.x+threadIdx.x; if(i>=n)return; Body a=in[i]; float ax=0,ay=0,az=0; for(int j=0;j=n)return; Body&x=b[i]; if(x.invMass<=0)return; float gpen=x.hy-x.py; if(gpen>0){ x.vy+=0.2f*gpen*60.0f; if(x.vy<0)x.vy=0; } x.px+=x.vx*dt; x.py+=x.vy*dt; x.pz+=x.vz*dt; } double stepGPU(Body* dA, Body* dB, int n, int iters, float dt){ int T=128, G=(n+T-1)/T; integrateVel<<>>(dA,n,dt); for(int it=0;it>>(dA,dB,n); cudaMemcpy(dA,dB,n*sizeof(Body),cudaMemcpyDeviceToDevice); } integratePos<<>>(dA,n,dt); return 0; } int main(int argc,char**argv){ const char* jsonPath="results/gpu_results.json"; for(int i=1;i cpu; initScene(cpu,n,seed); std::vector gpuInit=cpu; // identical seeded scene // CPU timing double penCPU=0; for(int i=0;i(t1-t0).count()/STEPS; // GPU timing Body *dA,*dB; cudaMalloc(&dA,n*sizeof(Body)); cudaMalloc(&dB,n*sizeof(Body)); cudaMemcpy(dA,gpuInit.data(),n*sizeof(Body),cudaMemcpyHostToDevice); for(int i=0;i gback(n); cudaMemcpy(gback.data(),dA,n*sizeof(Body),cudaMemcpyDeviceToHost); double div=0; for(int i=0;i