Skip to content

Commit

Permalink
Clean up g2p2g implementation and also fuse the grid clearing op into…
Browse files Browse the repository at this point in the history 
… g2p2g
  • Loading branch information
Chris Lewin committed Aug 8, 2024
1 parent 70a9233 commit 54e7dae
Show file tree
Hide file tree
Showing 6 changed files with 269 additions and 269 deletions.
217 changes: 179 additions & 38 deletions shaders/g2p2g.wgsl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,9 +12,11 @@
@group(0) @binding(1) var<storage, read_write> g_particles : array<Particle>;
@group(0) @binding(2) var<storage> g_gridSrc : array<i32>;
@group(0) @binding(3) var<storage, read_write> g_gridDst : array<atomic<i32>>;
@group(0) @binding(4) var<storage> g_bukkitThreadData : array<BukkitThreadData>;
@group(0) @binding(5) var<storage> g_bukkitParticleData : array<u32>;
@group(0) @binding(6) var<storage> g_shapes : array<SimShape>;
@group(0) @binding(4) var<storage, read_write> g_gridToBeCleared : array<i32>;
@group(0) @binding(5) var<storage> g_bukkitThreadData : array<BukkitThreadData>;
@group(0) @binding(6) var<storage> g_bukkitParticleData : array<u32>;
@group(0) @binding(7) var<storage> g_shapes : array<SimShape>;
@group(0) @binding(8) var<storage, read_write> g_freeIndices : array<atomic<i32>>;

const TotalBukkitEdgeLength = BukkitSize + BukkitHaloSize*2;
const TileDataSizePerEdge = TotalBukkitEdgeLength * 4;
Expand Down Expand Up @@ -47,8 +49,6 @@ fn csMain( @builtin(local_invocation_index) indexInGroup: u32, @builtin(workgrou

var gridVertexIsValid = all(gridVertex >= vec2i(0)) && all(gridVertex <= vec2i(g_simConstants.gridSize));

var appliedGridCorrection = false;

if(gridVertexIsValid)
{
let gridVertexAddress = gridVertexIndex(vec2u(gridVertex), g_simConstants.gridSize);
Expand All @@ -72,42 +72,36 @@ fn csMain( @builtin(local_invocation_index) indexInGroup: u32, @builtin(workgrou
dy = dy / w;
}

//dx = dx * 0.995;
//dy = dy * 0.995;




var gridDisplacement = vec2f(dx, dy);

// Collision detection against collider shapes
for(var shapeIndex = 0u; shapeIndex < g_simConstants.shapeCount; shapeIndex++)
{
let shape = g_shapes[shapeIndex];

// if(shape.functionality != ShapeFunctionCollider)
// {
// continue;
// }
if(shape.functionality != ShapeFunctionCollider)
{
continue;
}

// let displacedGridPosition = gridPosition + gridDisplacement;
let displacedGridPosition = gridPosition + gridDisplacement;

// let collideResult = collide(shape, displacedGridPosition);
let collideResult = collide(shape, displacedGridPosition);

// if(collideResult.collides)
// {
// let gap = min(0, dot(collideResult.normal, collideResult.pointOnCollider - gridPosition));
// let penetration = dot(collideResult.normal, gridDisplacement) - gap;
if(collideResult.collides)
{
let gap = min(0, dot(collideResult.normal, collideResult.pointOnCollider - gridPosition));
let penetration = dot(collideResult.normal, gridDisplacement) - gap;

// // Prevent any further penetration in radial direction
// let radialImpulse = max(penetration, 0);
// gridDisplacement -= radialImpulse*collideResult.normal;
// }
// Prevent any further penetration in radial direction
let radialImpulse = max(penetration, 0);
gridDisplacement -= radialImpulse*collideResult.normal;
}
}

// Collision detection against guardian shape

// Grid vertices near or inside the guardian region should have their displacenent values
// Grid vertices near or inside the guardian region should have their displacement values
// corrected in order to prevent particles moving into the guardian.
// We do this by finding whether a grid vertex would be inside the guardian region after displacement
// with the current velocity and, if it is, setting the displacement so that no further penetration can occur.
Expand All @@ -118,15 +112,11 @@ fn csMain( @builtin(local_invocation_index) indexInGroup: u32, @builtin(workgrou
if(projectedDifference.x != 0)
{
gridDisplacement.x = 0;
gridDisplacement.y = 0;
appliedGridCorrection = true;
}

if(projectedDifference.y != 0)
{
gridDisplacement.x = 0;
gridDisplacement.y = 0;
appliedGridCorrection = true;
}

dx = gridDisplacement.x;
Expand All @@ -140,11 +130,6 @@ fn csMain( @builtin(local_invocation_index) indexInGroup: u32, @builtin(workgrou
atomicStore(&s_tileData[tileDataIndex+2], encodeFixedPoint(w, g_simConstants.fixedPointMultiplier));
atomicStore(&s_tileData[tileDataIndex+3], encodeFixedPoint(v, g_simConstants.fixedPointMultiplier));

atomicStore(&s_tileDataDst[tileDataIndex], 0);
atomicStore(&s_tileDataDst[tileDataIndex+1], 0);
atomicStore(&s_tileDataDst[tileDataIndex+2], 0);
atomicStore(&s_tileDataDst[tileDataIndex+3], 0);

workgroupBarrier();

if(indexInGroup < threadData.rangeCount)
Expand Down Expand Up @@ -219,12 +204,163 @@ fn csMain( @builtin(local_invocation_index) indexInGroup: u32, @builtin(workgrou
particle.deformationDisplacement = B * 4.0;
particle.displacement = d;

// Save particle
g_particles[myParticleIndex] = particle;
}
// Do integration
if(g_simConstants.iteration == g_simConstants.iterationCount -1)
{
if(particle.material == MaterialLiquid)
{
// The liquid material only cares about the determinant of the deformation gradient.
// We can use the regular MPM integration below to evolve the deformation gradient, but
// this approximation actually conserves more volume.
// This is based on det(F^n+1) = det((I + D)*F^n) = det(I+D)*det(F^n)
// and noticing that D is a small velocity, we can use the identity det(I + D) ≈ 1 + tr(D) to first order
// ending up with det(F^n+1) = (1+tr(D))*det(F^n)
// Then we directly set particle.liquidDensity to reflect the approximately integrated volume.
// The liquid material does not actually use the deformation gradient matrix.
particle.liquidDensity *= (tr(particle.deformationDisplacement) + 1.0);

// Safety clamp to avoid instability with very small densities.
particle.liquidDensity = max(particle.liquidDensity, 0.05);
}
else
{
// Integrate transform using standard MPM formula
particle.deformationGradient = (Identity + particle.deformationDisplacement) * particle.deformationGradient;
}

if(particle.material != MaterialLiquid)
{
// SVD is necessary at least for safety clamp
var svdResult = svd(particle.deformationGradient);

// Clamp the lower bound of scale to prevent situations where large forces are generated leading to explosions
svdResult.Sigma = clamp(svdResult.Sigma, vec2f(0.1), vec2f(10000.0));

// Plasticity implementations
if(particle.material == MaterialSand)
{
// Drucker-Prager sand based on:
// Gergely Klár, Theodore Gast, Andre Pradhana, Chuyuan Fu, Craig Schroeder, Chenfanfu Jiang, and Joseph Teran. 2016.
// Drucker-prager elastoplasticity for sand animation. ACM Trans. Graph. 35, 4, Article 103 (July 2016), 12 pages.
// https://doi.org/10.1145/2897824.2925906
let sinPhi = sin(g_simConstants.frictionAngle/180.0 * 3.14159);
let alpha = sqrt(2.0/3.0)*2.0*sinPhi/(3.0 - sinPhi);
let beta = 0.5;

let eDiag = log(max(abs(svdResult.Sigma), vec2f(1e-6)));

let eps = diag(eDiag);
let trace = tr(eps) + particle.logJp;

let eHat = eps - (trace / 2) * Identity;
let frobNrm = length(vec2f(eHat[0][0], eHat[1][1]));

if(trace >= 0.0)
{
// In this case the motion is expansionary and we should not resist it at all.
// This means we should forget about any deformation that has occurred, which we do by setting Sigma to 1.
svdResult.Sigma = vec2f(1);
particle.logJp = beta * trace;
}
else
{
particle.logJp = 0;
let deltaGammaI = frobNrm + (g_simConstants.elasticityRatio + 1.0) * trace * alpha;
if(deltaGammaI > 0)
{
// Project to cone surface.
// This means we have to forget some deformation that the particle has undergone.
let h = eDiag - deltaGammaI/frobNrm * (eDiag - (trace*0.5)) ;
svdResult.Sigma = exp(h);
}
}
}
else if(particle.material == MaterialVisco)
{
// Very simple plasticity with volume preservation
let yieldSurface = exp(1-g_simConstants.plasticity);

// Record the volume before plasticity calculation
let J = svdResult.Sigma.x*svdResult.Sigma.y;

// Forget any deformation beyond the yield surface
svdResult.Sigma = clamp(svdResult.Sigma, vec2f(1.0/yieldSurface), vec2f(yieldSurface));

// Re-scale to original volume
let newJ = svdResult.Sigma.x*svdResult.Sigma.y;
svdResult.Sigma *= sqrt(J/newJ);
}

particle.deformationGradient = svdResult.U * diag(svdResult.Sigma) * svdResult.Vt;
}

// Integrate position
particle.position += particle.displacement;

// Mouse interaction
if(g_simConstants.mouseActivation > 0)
{
let offset = particle.position - g_simConstants.mousePosition;
let lenOffset = max(length(offset), 0.0001);
if(lenOffset < g_simConstants.mouseRadius)
{
let normOffset = offset/lenOffset;

if(g_simConstants.mouseFunction == MouseFunctionPush)
{
particle.displacement += normOffset * g_simConstants.mouseActivation;
}
else if(g_simConstants.mouseFunction == MouseFunctionGrab)
{
particle.displacement = 0.7*g_simConstants.mouseVelocity*g_simConstants.deltaTime;
}
}
}

// Gravity acceleration is normalized to the vertical size of the window.
particle.displacement.y -= f32(g_simConstants.gridSize.y)*g_simConstants.gravityStrength*g_simConstants.deltaTime*g_simConstants.deltaTime;

// Free count may be negative because of emission. So make sure it is at last zero before incrementing.
atomicMax(&g_freeIndices[0], 0i);

for(var shapeIndex = 0u; shapeIndex < g_simConstants.shapeCount; shapeIndex++)
{
let shape = g_shapes[shapeIndex];

// Push particles out of colliders. Most of the work should have been done at the grid level already.
if(shape.functionality == ShapeFunctionCollider)
{
let collideResult = collide(shape, particle.position);

if(collideResult.collides)
{
particle.displacement -= collideResult.penetration*collideResult.normal;
}
}

// Delete particles if they are inside a drain shape
if(shape.functionality == ShapeFunctionDrain)
{
if(collide(shape, particle.position).collides)
{
particle.enabled = 0;

// Add index of this particle to free list
let freeIndex = atomicAdd(&g_freeIndices[0], 1i);
atomicStore(&g_freeIndices[1 + u32(freeIndex)], i32(myParticleIndex));
}
}
}

// Ensure particles are inside the simulation limits
particle.position = projectInsideGuardian(particle.position, g_simConstants.gridSize, GuardianSize);
}



// Save particle
g_particles[myParticleIndex] = particle;
}


//if(g_simConstants.iteration != g_simConstants.iterationCount-1)
Expand Down Expand Up @@ -345,5 +481,10 @@ fn csMain( @builtin(local_invocation_index) indexInGroup: u32, @builtin(workgrou
atomicAdd(&g_gridDst[gridVertexAddress + 1], dyi);
atomicAdd(&g_gridDst[gridVertexAddress + 2], wi);
atomicAdd(&g_gridDst[gridVertexAddress + 3], vi);

g_gridToBeCleared[gridVertexAddress + 0] = 0;
g_gridToBeCleared[gridVertexAddress + 1] = 0;
g_gridToBeCleared[gridVertexAddress + 2] = 0;
g_gridToBeCleared[gridVertexAddress + 3] = 0;
}
}
7 changes: 3 additions & 4 deletions shaders/particleEmit.wgsl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,8 +13,7 @@
@group(0) @binding(1) var<storage, read_write> g_particleCount : array<atomic<u32>>;
@group(0) @binding(2) var<storage, read_write> g_particles : array<Particle>;
@group(0) @binding(3) var<storage> g_shapes : array<SimShape>;
@group(0) @binding(4) var<storage, read_write> g_freeCount : array<atomic<i32>>;
@group(0) @binding(5) var<storage, read_write> g_freeIndices : array<u32>;
@group(0) @binding(4) var<storage, read_write> g_freeIndices : array<atomic<i32>>;

fn createParticle(position: vec2f, material: f32, mass: f32, volume: f32, color: vec3f) -> Particle
{
Expand All @@ -38,10 +37,10 @@ fn addParticle(position: vec2f, material: f32, volume: f32, density: f32, jitter
{
var particleIndex = 0u;
// First check the free list to see if we can reuse a particle slot
let freeIndexSlot = atomicSub(&g_freeCount[0], 1i) - 1i;
let freeIndexSlot = atomicSub(&g_freeIndices[0], 1i) - 1i;
if(freeIndexSlot >= 0)
{
particleIndex = g_freeIndices[u32(freeIndexSlot)];
particleIndex = u32(atomicLoad(&g_freeIndices[u32(freeIndexSlot) + 1]));
}
else // If free list is empty then grow the particle count
{
Expand Down
Loading

0 comments on commit 54e7dae

Please sign in to comment.