From a83d395885ec82f26fb99fd2bcdd13e57852f153 Mon Sep 17 00:00:00 2001
From: adany292 <adany292@student.liu.se>
Date: Sun, 29 Dec 2024 19:56:57 +0100
Subject: [PATCH] Itterative SESPH/PCISPH improvement, precalulateed stifness
cosnt
---
src/FluidSimulationApp.cpp | 5 --
src/cuda/include/cuda_device_vars.h | 1 +
src/cuda/include/init_cuda.h | 9 +++
src/cuda/init_cuda.cu | 85 +++++++++++++++++++++++++++++
src/cuda/pcisph_kernels.cu | 50 ++++++++---------
src/cuda/pcisph_update.cu | 4 +-
src/settings/constants.h | 6 +-
7 files changed, 126 insertions(+), 34 deletions(-)
diff --git a/src/FluidSimulationApp.cpp b/src/FluidSimulationApp.cpp
index db4d400..45773eb 100644
--- a/src/FluidSimulationApp.cpp
+++ b/src/FluidSimulationApp.cpp
@@ -245,28 +245,23 @@ namespace FluidSimulation {
break;
case RenderStage::ThicknessBuffer:
- printf("Render Thickness \n");
- printf("Thickness stage active \n");
m_fluidRenderer->renderThicknessFrame(particleCount, boundaryCount);
m_fluidRenderer->visualizeBuffer(RenderStage::ThicknessBuffer);
break;
case RenderStage::FilteredThickness:
- printf("Render Filtered Thickness \n");
m_fluidRenderer->renderThicknessFrame(particleCount, boundaryCount);
m_fluidRenderer->renderThicknessFilterFrame();
m_fluidRenderer->visualizeBuffer(RenderStage::FilteredThickness);
break;
case RenderStage::LightAttenuation:
- printf("Render Filtered Thickness \n");
m_fluidRenderer->renderThicknessFrame(particleCount, boundaryCount);
m_fluidRenderer->renderThicknessFilterFrame();
m_fluidRenderer->renderThicknessAttenFrame();
break;
case RenderStage::ScreenSpacedRender:
- printf("Render Screen Spaced \n");
m_fluidRenderer->renderDepthFrame(particleCount, boundaryCount);
m_fluidRenderer->renderDepthFilterFrame();
m_fluidRenderer->renderDepthNormalsFrame();
diff --git a/src/cuda/include/cuda_device_vars.h b/src/cuda/include/cuda_device_vars.h
index b4eeb12..69a8f20 100644
--- a/src/cuda/include/cuda_device_vars.h
+++ b/src/cuda/include/cuda_device_vars.h
@@ -14,6 +14,7 @@ extern cudaGraphicsResource_t cuda_vbo_boundary_resource;
extern __constant__ float PI;
extern __constant__ float H;
extern __constant__ float H2;
+extern __constant__ float K_PCI;
extern __constant__ float POLY6_CONST;
extern __constant__ float SPIKY_GRAD_CONST;
extern __constant__ float VISC_LAP_CONST;
diff --git a/src/cuda/include/init_cuda.h b/src/cuda/include/init_cuda.h
index a554776..ac02e55 100644
--- a/src/cuda/include/init_cuda.h
+++ b/src/cuda/include/init_cuda.h
@@ -37,6 +37,15 @@ void updateConstantsOnDevice(
float smoothingRadius
);
+void ComputePCISPHStiffness(
+ float timeStep,
+ float particleMass,
+ float restDensity,
+ float kernelRadius,
+ float restSpacing,
+ float spiky_d1_normalization
+);
+
} // namespace FluidSimulationGPU
#endif // INIT_CUDA_H
diff --git a/src/cuda/init_cuda.cu b/src/cuda/init_cuda.cu
index b5b629e..fa7bc1a 100644
--- a/src/cuda/init_cuda.cu
+++ b/src/cuda/init_cuda.cu
@@ -4,6 +4,7 @@
#include "../settings/constants.h"
#include "include/pcisph_kernels.h"
+#include "include/float3_helper_math.h"
#include <cuda_runtime.h>
#include <cuda_gl_interop.h>
@@ -20,6 +21,7 @@ cudaGraphicsResource_t cuda_vbo_boundary_resource;
__constant__ float PI;
__constant__ float H;
__constant__ float H2;
+__constant__ float K_PCI;
__constant__ float POLY6_CONST;
__constant__ float SPIKY_GRAD_CONST;
__constant__ float VISC_LAP_CONST;
@@ -71,6 +73,78 @@ unsigned int* d_boundaryCellEnd = nullptr;
unsigned int* d_boundaryCellIndices = nullptr;
unsigned int* d_boundaryParticleIndices = nullptr;
+void ComputePCISPHStiffness(
+ float timeStep,
+ float particleMass,
+ float restDensity,
+ float kernelRadius,
+ 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);
+ printf("Rest Density: %f \n", restDensity);
+
+ //cudaMemcpyToSymbol(PARTICLE_MASS, &mass, sizeof(float));
+ //cudaMemcpyToSymbol(REST_DENSITY, &restDensity, sizeof(float));
+
+ float dtm_rho = (timeStep * mass) / restDensity;
+ 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
+ for (float z = -kernelRadius; z <= kernelRadius + 1e-5f; z += restSpacing)
+ {
+ for (float y = -kernelRadius; y <= kernelRadius + 1e-5f; y += restSpacing)
+ {
+ for (float x = -kernelRadius; x <= kernelRadius + 1e-5f; x += restSpacing)
+ {
+ 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
+ grad_sum += grad;
+ sum_dot_grad += dot(grad, grad);
+ }
+ }
+ }
+ }
+
+ // 4) 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:
+ float grad_term = - (dot_sum_grad + sum_dot_grad);
+
+ // 6) 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);
+ cudaMemcpyToSymbol(K_PCI, &k_PCI, sizeof(float));
+}
+
void initializeConstants()
{
// Host constants
@@ -98,6 +172,15 @@ void initializeConstants()
cudaMemcpyToSymbol(BOUND_Y_MAX, &h_BOUND_Y_MAX, sizeof(float));
cudaMemcpyToSymbol(BOUND_Z_MIN, &h_BOUND_Z_MIN, sizeof(float));
cudaMemcpyToSymbol(BOUND_Z_MAX, &h_BOUND_Z_MAX, sizeof(float));
+
+ ComputePCISPHStiffness(
+ 0.006,
+ h_PARTICLE_MASS,
+ h_REST_DENSITY,
+ h_H,
+ h_H/96.0f,
+ h_SPIKY_GRAD_CONST
+ );
}
// Register VBO with CUDA
@@ -166,6 +249,8 @@ void initializeFluidDeviceArrays(
// Initialize positions_pred to initial positions
cudaMemcpy(d_positions_pred, h_positions, numParticles * sizeof(float3), cudaMemcpyHostToDevice);
+
+
}
void initializeBoundryDeviceArrays(
diff --git a/src/cuda/pcisph_kernels.cu b/src/cuda/pcisph_kernels.cu
index 0357fd8..d5c995e 100644
--- a/src/cuda/pcisph_kernels.cu
+++ b/src/cuda/pcisph_kernels.cu
@@ -139,7 +139,7 @@ __global__ void updatePressures(
int idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < numParticles) {
- float deltaP = -density_errors[idx] / A;
+ float deltaP = -density_errors[idx] / K_PCI;
pressure_increments[idx] = deltaP;
}
}
@@ -211,29 +211,29 @@ __global__ void computePressureForces(
// Boundary neighbors
// ------------
- /*
- unsigned int b_startIdx = boundaryCellStart[neighborCellIndex];
- unsigned int b_endIdx = boundaryCellEnd[neighborCellIndex];
-
- if (b_startIdx != 0xFFFFFFFF) {
- for (unsigned int j = b_startIdx; j < b_endIdx; ++j) {
- float3 pos_j = boundaryPositions[j];
- // Mirror pressure: p_b = p_i, and use REST_DENSITY for boundary
- float deltaP_j = deltaP_i; // Pressure mirroring
-
- float3 r = pos_i - pos_j;
- float distSq = dot(r, r);
- if (distSq < H2 && distSq > 0.0001f) {
- float dist = sqrtf(distSq);
- float gradW = spikyKernelGrad(dist);
- float3 gradWij = gradW * (r / dist);
- // Using Equation (85): p_i and p_b = p_i; densities: ρ_i and ρ_0
- float factor = - (deltaP_i + deltaP_j) / (REST_DENSITY * REST_DENSITY);
- fPressure += factor * gradWij * PARTICLE_MASS * PARTICLE_MASS;
- }
- }
- }
- */
+
+ // unsigned int b_startIdx = boundaryCellStart[neighborCellIndex];
+ // unsigned int b_endIdx = boundaryCellEnd[neighborCellIndex];
+
+ // if (b_startIdx != 0xFFFFFFFF) {
+ // for (unsigned int j = b_startIdx; j < b_endIdx; ++j) {
+ // float3 pos_j = boundaryPositions[j];
+ // // Mirror pressure: p_b = p_i, and use REST_DENSITY for boundary
+ // float deltaP_j = deltaP_i; // Pressure mirroring
+
+ // float3 r = pos_i - pos_j;
+ // float distSq = dot(r, r);
+ // if (distSq < H2 && distSq > 0.0001f) {
+ // float dist = sqrtf(distSq);
+ // float gradW = spikyKernelGrad(dist);
+ // float3 gradWij = gradW * (r / dist);
+ // // Using Equation (85): p_i and p_b = p_i; densities: ρ_i and ρ_0
+ // float factor = - (deltaP_i + deltaP_j) / (REST_DENSITY * REST_DENSITY);
+ // fPressure += factor * gradWij * PARTICLE_MASS * PARTICLE_MASS;
+ // }
+ // }
+ // }
+
}
}
}
@@ -369,7 +369,7 @@ __global__ void computeViscosityAndExternalForces(
float3 vel_i = velocities[idx];
float3 fViscosity = make_float3(0.0f, 0.0f, 0.0f);
- float3 fGravity = make_float3(0.0f, -20.0f, 0.0f) * PARTICLE_MASS;
+ float3 fGravity = make_float3(0.0f, -10.0f, 0.0f) * PARTICLE_MASS;
unsigned int cellIndex = cellIndices[idx];
diff --git a/src/cuda/pcisph_update.cu b/src/cuda/pcisph_update.cu
index 757f723..4b296be 100644
--- a/src/cuda/pcisph_update.cu
+++ b/src/cuda/pcisph_update.cu
@@ -22,7 +22,7 @@ void updateParticles(float deltaTime)
int blocksPerGrid = (NUM_PARTICLES + threadsPerBlock - 1) / threadsPerBlock;
// Time step should be small but we can run updateParticles multiple times during display
- float timeStep = 0.004f;
+ float timeStep = 0.006f;
// -------------------- Neighbor Search --------------------
int gridSizeX, gridSizeY, gridSizeZ;
@@ -99,6 +99,8 @@ void updateParticles(float deltaTime)
float rho0 = h_REST_DENSITY;
float A = 2.0f * timeStep * timeStep * m * m / (rho0 * rho0) * totalSumGradW2;
+ //printf("K_PCI: %f \n", A);
+
// ------------ Iterative Pressure Correction Loop --------------------
const int maxIterations = 3;
diff --git a/src/settings/constants.h b/src/settings/constants.h
index d34d1eb..de8dea1 100644
--- a/src/settings/constants.h
+++ b/src/settings/constants.h
@@ -12,11 +12,11 @@ const float h_BOUND_Z_MAX = 4.0f;
const float h_GAS_CONSTANT = 0.02f;
const float h_VISCOSITY = 0.04f;
-const float h_PARTICLE_MASS = 0.01f;
-const float h_REST_DENSITY = 800.0f;
+const float h_PARTICLE_MASS = 0.005f;
+const float h_REST_DENSITY = 600.0f;
const float h_PI = 3.14159265358979323846f;
-const float h_H = 0.25f; // Smoothing radius
+const float h_H = 0.20f; // Smoothing radius
const int NUM_PARTICLES = 50000;
--
GitLab