Skip to content
Snippets Groups Projects
Commit b48d9a4f authored by Adam Nyberg's avatar Adam Nyberg
Browse files

Surface tension WIP

parent 28a3ec4d
Branches
No related tags found
No related merge requests found
......@@ -19,7 +19,7 @@ void main()
discard;
// For thickness: each particle in this pixel adds some amount.
float thicknessIncrement = 0.08;
float thicknessIncrement = 0.015;
// Output that in the red channel. Using GL32F only reed chanele sotred
FragColor = vec4(thicknessIncrement, 0.0, 0.0, 1.0);
......
......@@ -16,7 +16,7 @@ void main()
// Sample center
float centerVal = texture(uThicknessTex, TexCoords).r;
if (centerVal < 0.01) discard;
if (centerVal < 0.005) discard;
// Accumulate for a 1D Gaussian blur
float sum = 0.0;
......
......@@ -637,7 +637,7 @@ void FluidRenderer::renderDepthFrame(size_t particleCount, size_t boundaryCount)
glUniformMatrix4fv(m_depthRenderProgram.viewMatrixLoc, 1, GL_TRUE, viewMatrix.m);
if (m_depthRenderProgram.sphereRadiusLoc != -1) {
glUniform1f(m_depthRenderProgram.sphereRadiusLoc, 0.25f);
glUniform1f(m_depthRenderProgram.sphereRadiusLoc, 0.05f);
}
if (m_depthRenderProgram.pointRadiusLoc != -1) {
glUniform1f(m_depthRenderProgram.pointRadiusLoc, m_pointRadius);
......
......@@ -216,67 +216,66 @@ namespace FluidSimulation {
if (m_updatePhysics)
{
m_fluidSimulation->updateSimulation();
}
size_t particleCount = m_fluidSimulation->getParticleCount();
size_t boundaryCount = m_fluidSimulation->getBoundaryParticleCount();
// Precompute particle counts to avoid multiple function calls
size_t particleCount = m_fluidSimulation->getParticleCount();
size_t boundaryCount = m_fluidSimulation->getBoundaryParticleCount();
switch (selectedStage)
{
case RenderStage::ParticleView:
m_fluidRenderer->renderNormalFrame(particleCount, boundaryCount);
break;
case RenderStage::DepthBuffer:
m_fluidRenderer->renderDepthFrame(particleCount, boundaryCount);
m_fluidRenderer->visualizeBuffer(RenderStage::DepthBuffer);
break;
case RenderStage::FilteredDepth:
m_fluidRenderer->renderDepthFrame(particleCount, boundaryCount);
m_fluidRenderer->renderDepthFilterFrame();
m_fluidRenderer->visualizeBuffer(RenderStage::FilteredDepth);
break;
case RenderStage::DepthNormals:
m_fluidRenderer->renderDepthFrame(particleCount, boundaryCount);
m_fluidRenderer->renderDepthFilterFrame();
m_fluidRenderer->renderDepthNormalsFrame();
break;
case RenderStage::ThicknessBuffer:
m_fluidRenderer->renderThicknessFrame(particleCount, boundaryCount);
m_fluidRenderer->visualizeBuffer(RenderStage::ThicknessBuffer);
break;
case RenderStage::FilteredThickness:
m_fluidRenderer->renderThicknessFrame(particleCount, boundaryCount);
m_fluidRenderer->renderThicknessFilterFrame();
m_fluidRenderer->visualizeBuffer(RenderStage::FilteredThickness);
break;
case RenderStage::LightAttenuation:
m_fluidRenderer->renderThicknessFrame(particleCount, boundaryCount);
m_fluidRenderer->renderThicknessFilterFrame();
m_fluidRenderer->renderThicknessAttenFrame();
break;
switch (selectedStage)
{
case RenderStage::ParticleView:
m_fluidRenderer->renderNormalFrame(particleCount, boundaryCount);
break;
case RenderStage::DepthBuffer:
m_fluidRenderer->renderDepthFrame(particleCount, boundaryCount);
m_fluidRenderer->visualizeBuffer(RenderStage::DepthBuffer);
break;
case RenderStage::FilteredDepth:
m_fluidRenderer->renderDepthFrame(particleCount, boundaryCount);
m_fluidRenderer->renderDepthFilterFrame();
m_fluidRenderer->visualizeBuffer(RenderStage::FilteredDepth);
break;
case RenderStage::DepthNormals:
m_fluidRenderer->renderDepthFrame(particleCount, boundaryCount);
m_fluidRenderer->renderDepthFilterFrame();
m_fluidRenderer->renderDepthNormalsFrame();
break;
case RenderStage::ThicknessBuffer:
m_fluidRenderer->renderThicknessFrame(particleCount, boundaryCount);
m_fluidRenderer->visualizeBuffer(RenderStage::ThicknessBuffer);
break;
case RenderStage::FilteredThickness:
m_fluidRenderer->renderThicknessFrame(particleCount, boundaryCount);
m_fluidRenderer->renderThicknessFilterFrame();
m_fluidRenderer->visualizeBuffer(RenderStage::FilteredThickness);
break;
case RenderStage::LightAttenuation:
m_fluidRenderer->renderThicknessFrame(particleCount, boundaryCount);
m_fluidRenderer->renderThicknessFilterFrame();
m_fluidRenderer->renderThicknessAttenFrame();
break;
case RenderStage::ScreenSpacedRender:
m_fluidRenderer->renderDepthFrame(particleCount, boundaryCount);
m_fluidRenderer->renderDepthFilterFrame();
m_fluidRenderer->renderDepthNormalsFrame();
case RenderStage::ScreenSpacedRender:
m_fluidRenderer->renderDepthFrame(particleCount, boundaryCount);
m_fluidRenderer->renderDepthFilterFrame();
m_fluidRenderer->renderDepthNormalsFrame();
m_fluidRenderer->renderThicknessFrame(particleCount, boundaryCount);
m_fluidRenderer->renderThicknessFilterFrame();
m_fluidRenderer->renderThicknessAttenFrame();
m_fluidRenderer->renderFinalCombineFrame();
break;
m_fluidRenderer->renderThicknessFrame(particleCount, boundaryCount);
m_fluidRenderer->renderThicknessFilterFrame();
m_fluidRenderer->renderThicknessAttenFrame();
default:
std::cerr << "Unknown RenderStage selected." << std::endl;
break;
}
m_fluidRenderer->renderFinalCombineFrame();
break;
default:
std::cerr << "Unknown RenderStage selected." << std::endl;
break;
}
// Display FBO texture in ImGui window
......
......@@ -19,6 +19,11 @@ extern __constant__ float POLY6_CONST;
extern __constant__ float SPIKY_GRAD_CONST;
extern __constant__ float VISC_LAP_CONST;
extern __constant__ float POLY6_D1_CONST;
extern __constant__ float SURFACE_TENS_CONST;
extern __constant__ float SURFACE_TENS_TERM;
extern __constant__ float SURFACE_TENS_COEFFICIENT;
extern __constant__ float PARTICLE_MASS;
extern __constant__ float REST_DENSITY;
extern __constant__ float VISCOSITY;
......@@ -50,6 +55,7 @@ extern float* d_density_errors;
extern float* d_pressure_increments;
extern float3* d_pressure_forces;
extern float3* d_accelerations;
extern float3* d_normals;
// Additional arrays for PCISPH
extern float3* d_velocities_pred;
......
......@@ -7,6 +7,19 @@
namespace FluidSimulationGPU {
__global__ void initialDensity(
float3* positions,
float* densities,
unsigned int* cellIndices,
unsigned int* cellStart,
unsigned int* cellEnd,
float3* boundaryPositions,
unsigned int* boundaryCellStart,
unsigned int* boundaryCellEnd,
int numParticles,
int gridSizeX, int gridSizeY, int gridSizeZ
);
// Function to compute density errors
__global__ void computeDensityError(
float3* positions,
......@@ -63,10 +76,24 @@ __global__ void correctBoundaryPenetrations(
);
// Function to compute viscosity and external forces
__global__ void computeViscosityAndExternalForces(
__global__ void computeInitialForces(
float3* positions,
float3* velocities,
float3* accelerations,
float3* normals,
float* densities,
unsigned int* cellIndices,
unsigned int* cellStart,
unsigned int* cellEnd,
int numParticles,
int gridSizeX, int gridSizeY, int gridSizeZ
);
__global__ void calculateNormals(
float3* positions,
float3* velocities,
float3* normals,
float* densities,
unsigned int* cellIndices,
unsigned int* cellStart,
unsigned int* cellEnd,
......
......@@ -6,8 +6,10 @@
namespace FluidSimulationGPU {
__device__ float poly6Kernel(float r2);
__device__ float kernel_poly6_d1(float r2);
__device__ float spikyKernelGrad(float r);
__device__ float viscosityKernelLap(float r);
__device__ float kernel_surface_tension(float r);
} // namespace FluidSimulationGPU
......
......@@ -26,6 +26,11 @@ __constant__ float POLY6_CONST;
__constant__ float SPIKY_GRAD_CONST;
__constant__ float VISC_LAP_CONST;
__constant__ float POLY6_D1_CONST;
__constant__ float SURFACE_TENS_CONST;
__constant__ float SURFACE_TENS_TERM;
__constant__ float SURFACE_TENS_COEFFICIENT;
__constant__ float PARTICLE_MASS;
__constant__ float REST_DENSITY;
__constant__ float VISCOSITY;
......@@ -57,6 +62,7 @@ float* d_density_errors = nullptr;
float* d_pressure_increments = nullptr;
float3* d_pressure_forces = nullptr;
float3* d_accelerations = nullptr;
float3* d_normals = nullptr;
// Additional arrays for PCISPH
float3* d_velocities_pred = nullptr;
......@@ -81,8 +87,6 @@ void ComputePCISPHStiffness(
float restSpacing,
float spiky_d1_normalization)
{
// 1) Compute the factor often called 'beta' or something similar
// in your code. (Matches your "2 * (dt*m / rho0)^2" idea)
float mass = restDensity * powf(kernelRadius/2.0f, 3.f);
printf("Particle mass: %f \n", mass);
......@@ -95,13 +99,12 @@ void ComputePCISPHStiffness(
float beta = 2.0f * dtm_rho * dtm_rho;
// 2) Variables to accumulate gradient sums
float3 grad_sum = make_float3(0.f, 0.f, 0.f);
float sum_dot_grad = 0.f;
float kernelRadius2 = kernelRadius * kernelRadius;
// 3) Sample points in [-kernelRadius, +kernelRadius] in increments of restSpacing
// Sample points in [-kernelRadius, +kernelRadius] in increments of restSpacing
for (float z = -kernelRadius; z <= kernelRadius + 1e-5f; z += restSpacing)
{
for (float y = -kernelRadius; y <= kernelRadius + 1e-5f; y += restSpacing)
......@@ -111,15 +114,10 @@ void ComputePCISPHStiffness(
float3 r = make_float3(x, y, z);
float r2 = dot(r, r);
// Only compute if within kernel support
if (r2 < kernelRadius2 && r2 > 0.0f) // avoid r=0 to prevent division by 0
{
float r_norm = sqrtf(r2);
// Spiky kernel gradient in 3D usually:
// grad W_spiky(r) = -C * ( (h - r)^2 / r ) * r̂
//
// If you're using a negative sign in spiky_d1_normalization,
// you might want to omit the minus here, etc. Just be consistent.
float3 grad = (r / r_norm) * spiky_d1_normalization * powf(kernelRadius - r_norm, 2.f);
// Accumulate sums
......@@ -130,15 +128,14 @@ void ComputePCISPHStiffness(
}
}
// 4) Compute dot of the sum of gradients
// Compute dot of the sum of gradients
float dot_sum_grad = dot(grad_sum, grad_sum);
// 5) Combine to get grad_term:
// grad_term ~ - [ (Σ gradW)^2 + Σ (gradW·gradW) ]
// so we define:
// Combine to get grad_term:
// grad_term ~ - [ (Σ gradW)^2 + Σ (gradW·gradW) ]
float grad_term = - (dot_sum_grad + sum_dot_grad);
// 6) Finally, PCISPH stiffness factor k_PCI ~ -1 / (beta * grad_term)
// Finally, PCISPH stiffness factor k_PCI ~ -1 / (beta * grad_term)
float k_PCI = -1.f / (beta * grad_term);
printf("K_PCi precalc: %f \n", k_PCI);
......@@ -149,9 +146,15 @@ void initializeConstants()
{
// Host constants
const float h_H2 = h_H * h_H;
const float h_POLY6_CONST = 315.0f / (64.0f * h_PI * powf(h_H, 9));
const float h_SPIKY_GRAD_CONST = -45.0f / (h_PI * powf(h_H, 6));
const float h_VISC_LAP_CONST = 45.0f / (h_PI * powf(h_H, 6));
const float h_POLY6_CONST = 315.0f / (64.0f * M_PI * powf(h_H, 9));
const float h_SPIKY_GRAD_CONST = -45.0f / (M_PI * powf(h_H, 6));
const float h_VISC_LAP_CONST = 45.0f / (M_PI * powf(h_H, 6));
const float poly6_d1 = -945.f / (32.f * M_PI * powf(h_H, 9.f));
const float surface_tens_normalization = 32.f / (M_PI * powf(h_H, 9.f));
const float surface_tens_term = powf(h_H, 6.f) / 64.f;
// Copy constants to device
cudaMemcpyToSymbol(PI, &h_PI, sizeof(float));
......@@ -161,6 +164,13 @@ void initializeConstants()
cudaMemcpyToSymbol(SPIKY_GRAD_CONST, &h_SPIKY_GRAD_CONST, sizeof(float));
cudaMemcpyToSymbol(VISC_LAP_CONST, &h_VISC_LAP_CONST, sizeof(float));
cudaMemcpyToSymbol(POLY6_D1_CONST, &poly6_d1, sizeof(float));
cudaMemcpyToSymbol(SURFACE_TENS_CONST, &surface_tens_normalization, sizeof(float));
cudaMemcpyToSymbol(SURFACE_TENS_TERM, &surface_tens_term, sizeof(float));
cudaMemcpyToSymbol(SURFACE_TENS_COEFFICIENT, &h_surface_tenns_coef, sizeof(float));
cudaMemcpyToSymbol(PARTICLE_MASS, &h_PARTICLE_MASS, sizeof(float));
cudaMemcpyToSymbol(REST_DENSITY, &h_REST_DENSITY, sizeof(float));
cudaMemcpyToSymbol(VISCOSITY, &h_VISCOSITY, sizeof(float));
......@@ -243,14 +253,13 @@ void initializeFluidDeviceArrays(
cudaMalloc(&d_pressure_increments, numParticles * sizeof(float));
cudaMalloc(&d_pressure_forces, numParticles * sizeof(float3));
cudaMalloc(&d_accelerations, numParticles * sizeof(float3));
cudaMalloc(&d_normals, numParticles * sizeof(float3));
// Allocate device memory for sum_gradW2
cudaMalloc(&d_sum_gradW2, numParticles * sizeof(float));
// Initialize positions_pred to initial positions
cudaMemcpy(d_positions_pred, h_positions, numParticles * sizeof(float3), cudaMemcpyHostToDevice);
}
void initializeBoundryDeviceArrays(
......
......@@ -26,6 +26,73 @@ namespace FluidSimulationGPU {
// Compute density errors (used in PCISPH)
__global__ void initialDensity(
float3* positions,
float* densities,
unsigned int* cellIndices,
unsigned int* cellStart,
unsigned int* cellEnd,
float3* boundaryPositions,
unsigned int* boundaryCellStart,
unsigned int* boundaryCellEnd,
int numParticles,
int gridSizeX, int gridSizeY, int gridSizeZ
)
{
int idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx >= numParticles) return;
float3 pos_i = positions[idx];
float density = 0.0f;
unsigned int cellIndex = cellIndices[idx];
unsigned int neighborCellIdx;
FOREACH_NEIGHBOR_CELL_BEGIN(cellIndex, gridSizeX, gridSizeY, gridSizeZ, neighborCellIdx)
{
unsigned int startIdx = cellStart[neighborCellIdx];
unsigned int endIdx = cellEnd[neighborCellIdx];
if (startIdx != 0xFFFFFFFF)
{
for (unsigned int j = startIdx; j < endIdx; ++j)
{
float3 pos_j = positions[j];
float3 r = pos_i - pos_j;
float r2 = dot(r, r);
// Check if in range
if (r2 < H2 && r2 > 0.0001f) {
float W = poly6Kernel(r2);
density += PARTICLE_MASS * W;
}
}
}
// Boundary neighbors
unsigned int b_startIdx = boundaryCellStart[neighborCellIdx];
unsigned int b_endIdx = boundaryCellEnd[neighborCellIdx];
if (b_startIdx != 0xFFFFFFFF)
{
for (unsigned int j = b_startIdx; j < b_endIdx; ++j)
{
float3 pos_j = boundaryPositions[j];
float3 r = pos_i - pos_j;
float r2 = dot(r, r);
if (r2 < H2 && r2 > 0.0001f) {
float W = poly6Kernel(r2);
density += PARTICLE_MASS * W;
}
}
}
}
FOREACH_NEIGHBOR_CELL_END();
densities[idx] = density;
}
__global__ void computeDensityError(
float3* positions,
float* densities,
......@@ -96,6 +163,57 @@ __global__ void computeDensityError(
density_errors[idx] = density - REST_DENSITY;
}
__global__ void calculateNormals(
float3* positions,
float3* velocities,
float3* normals,
float* densities,
unsigned int* cellIndices,
unsigned int* cellStart,
unsigned int* cellEnd,
int numParticles,
int gridSizeX, int gridSizeY, int gridSizeZ
)
{
int idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx >= numParticles) return;
float3 pos_i = positions[idx];
float density = 0.0f;
float3 normal = make_float3(0.f, 0.f, 0.f);
unsigned int cellIndex = cellIndices[idx];
unsigned int neighborCellIdx;
FOREACH_NEIGHBOR_CELL_BEGIN(cellIndex, gridSizeX, gridSizeY, gridSizeZ, neighborCellIdx)
{
unsigned int startIdx = cellStart[neighborCellIdx];
unsigned int endIdx = cellEnd[neighborCellIdx];
if (startIdx != 0xFFFFFFFF)
{
for (unsigned int j = startIdx; j < endIdx; ++j)
{
float3 pos_j = positions[j];
float dens_j = densities[j];
float3 r = pos_i - pos_j;
float r2 = dot(r, r);
normal += kernel_poly6_d1(r2) * r / dens_j;
// poly6GradKernel
// normal += poly6grad(r) / other_density;
}
}
}
FOREACH_NEIGHBOR_CELL_END();
normal *= H * PARTICLE_MASS * POLY6_D1_CONST;
normals[idx] = normal;
}
// Update pressures based on density errors
__global__ void updatePressures(
......@@ -285,10 +403,12 @@ __global__ void correctBoundaryPenetrations(
}
// Compute viscosity and external forces (excluding pressure forces)
__global__ void computeViscosityAndExternalForces(
__global__ void computeInitialForces(
float3* positions,
float3* velocities,
float3* accelerations,
float3* normals,
float* densities,
unsigned int* cellIndices,
unsigned int* cellStart,
unsigned int* cellEnd,
......@@ -301,10 +421,14 @@ __global__ void computeViscosityAndExternalForces(
float3 pos_i = positions[idx];
float3 vel_i = velocities[idx];
float3 normal_i = normals[idx];
float dens_i = densities[idx];
float3 fViscosity = make_float3(0.0f, 0.0f, 0.0f);
float3 fGravity = make_float3(0.0f, -10.0f, 0.0f) * PARTICLE_MASS;
float3 st_cohesion = make_float3(0.0f, 0.0f, 0.0f);
float3 st_curvature = make_float3(0.0f, 0.0f, 0.0f);
// Current particle's cell
unsigned int cellIndex = cellIndices[idx];
// We'll store neighbor cell index here
......@@ -326,23 +450,46 @@ __global__ void computeViscosityAndExternalForces(
float3 pos_j = positions[j];
float3 vel_j = velocities[j];
float3 normal_j = normals[j];
float dens_j = densities[j];
float3 r = pos_i - pos_j;
float distSq = dot(r, r);
if (distSq < H2 && distSq > 0.0001f)
float dist = sqrtf(distSq);
if (dist < H && dist > 0.0001f)
{
float dist = sqrtf(distSq);
// Viscosity force
float3 velDiff = vel_j - vel_i;
float viscLap = viscosityKernelLap(dist);
fViscosity += VISCOSITY * velDiff * viscLap * PARTICLE_MASS / REST_DENSITY;
// Implement surface tension here
float st_correction_factor = 2.f * REST_DENSITY / (dens_i + dens_j);
// -> surface tension (cohesion)
float st_kernel = kernel_surface_tension(dist);
float3 direction = r/dist;
st_cohesion += st_correction_factor * st_kernel * direction;
// -> surface tension (curvature)
st_curvature += st_correction_factor * (normal_i - normal_j);
}
}
}
}
FOREACH_NEIGHBOR_CELL_END();
// finish surface tensiont calculations
//st_cohesion *= -SURFACE_TENS_COEFFICIENT * PARTICLE_MASS * PARTICLE_MASS * SURFACE_TENS_CONST;
//st_curvature *= -SURFACE_TENS_COEFFICIENT * PARTICLE_MASS;
st_cohesion *= 0.f;
st_curvature *= 0.f;
float3 fSurfaceTension = (st_cohesion + st_curvature);
// Total acceleration (excluding pressure)
accelerations[idx] = (fViscosity + fGravity) / PARTICLE_MASS;
}
......
......@@ -43,11 +43,37 @@ void updateParticles(float deltaTime)
gridSizeZ);
// -------------------- Prediction Step --------------------
initialDensity<<<blocksPerGrid, threadsPerBlock>>>(
d_positions_pred,
d_densities,
d_cellIndices,
d_cellStart,
d_cellEnd,
d_boundary_positions,
d_boundaryCellStart,
d_boundaryCellEnd,
NUM_PARTICLES,
gridSizeX, gridSizeY, gridSizeZ);
cudaDeviceSynchronize();
calculateNormals<<<blocksPerGrid, threadsPerBlock>>>(
d_positions,
d_velocities,
d_normals,
d_densities,
d_cellIndices,
d_cellStart,
d_cellEnd,
NUM_PARTICLES,
gridSizeX, gridSizeY, gridSizeZ);
cudaDeviceSynchronize();
computeViscosityAndExternalForces<<<blocksPerGrid, threadsPerBlock>>>(
computeInitialForces<<<blocksPerGrid, threadsPerBlock>>>(
d_positions,
d_velocities,
d_accelerations,
d_normals,
d_densities,
d_cellIndices,
d_cellStart,
d_cellEnd,
......
......@@ -14,6 +14,17 @@ __device__ float poly6Kernel(float r2)
}
}
__device__ float kernel_poly6_d1(float r2) {
if(r2 > H2) {
return 0.f;
} else {
float term = r2 - H2;
return term * term; // * r need to multiply my diff in calulation
}
}
__device__ float spikyKernelGrad(float r)
{
if (r > 0.0f && r <= H) {
......@@ -33,4 +44,19 @@ __device__ float viscosityKernelLap(float r)
}
}
__device__ float kernel_surface_tension(float r) {
if( r <= 0.5 * H ) {
float t1 = (H - r);
float t2 = (t1 * t1 * t1) * (r * r * r);
return 2.f * t2 - SURFACE_TENS_TERM;
}
else if( r <= H ) {
float t1 = (H - r);
return (t1 * t1 * t1) * (r * r * r);
}
else {
return 0.f;
}
}
} // namespace FluidSimulationGPU
......@@ -11,12 +11,13 @@ const float h_BOUND_Z_MAX = 4.0f;
const float h_GAS_CONSTANT = 0.02f;
const float h_VISCOSITY = 0.04f;
const float h_surface_tenns_coef = 0.0f;
const float h_PARTICLE_MASS = 0.005f;
const float h_PARTICLE_MASS = 0.01f;
const float h_REST_DENSITY = 600.0f;
const float h_PI = 3.14159265358979323846f;
const float h_H = 0.18f; // Smoothing radius
const float h_H = 0.20f; // Smoothing radius
const int NUM_PARTICLES = 200000;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment