Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
Browse files

SPH not working yet, debugging densities and forces

  • Loading branch information...
commit 149002ec6b492cee62c4e9e309598453359a0ee9 1 parent 765d256
@hrvthzs authored
View
2  Makefile
@@ -46,6 +46,8 @@ CUFILES := \
# CUDA dependency files
CU_DEPS := \
+ src/boundary.cu \
+ src/boundary_walls.cu \
src/grid_utils.cu \
src/grid_kernel.cu \
src/sph_density.cu \
View
19 shaders/shader.fs
@@ -1,9 +1,26 @@
#version 150
in vec4 color;
+in vec2 vTextCoord;
out vec4 fragColor;
void main()
{
- fragColor = vec4(0.09, 0.31, 0.98, 1.0);
+ vec4 c = vec4(0.09, 0.31, 0.98, 1.0);
+
+ const vec3 lightDir = vec3(0.577, 0.577, 0.577);
+
+ // calculate normal from texture coordinates
+ //vec3 N;
+ //N.xy = vTextCoord.xy*vec2(2.0, -2.0) + vec2(-1.0, 1.0);
+ //float mag = dot(N.xy, N.xy);
+ //if (mag > 1.0) discard; // kill pixels outside circle
+ //N.z = sqrt(1.0-mag);
+
+ // calculate lighting
+ //float diffuse = max(1.0, dot(lightDir, N));
+
+
+ //fragColor = c * diffuse;
+ fragColor = c;
}
View
12 shaders/shader.vs
@@ -1,12 +1,18 @@
-#version 150
-
-
+#version 130
in vec4 position;
uniform mat4 mvp;
+uniform mat4 mv;
+uniform float pointRadius;
+uniform float pointScale;
+out vec2 vTexCoord;
void main()
{
+ vec3 posEye = vec3(mv * position);
+ float dist = length(posEye);
+ //vTexCoord = gl_MultiTexCoord0.xy;
+ gl_PointSize = pointScale / (pointRadius * dist);
gl_Position = mvp * vec4(position);
}
View
83 src/boundary.cu
@@ -0,0 +1,83 @@
+#ifndef __BOUNDARY_CU__
+#define __BOUNDARY_CU__
+
+namespace Boundary {
+
+ //for collision detection
+ #define EPSILON 0.00001f
+
+
+ __device__ float3 calculateRepulsionForce(
+ float3 const& vel,
+ float3 const& normal,
+ float const& boundary_distance,
+ float const& boundary_dampening,
+ float const& boundary_stiffness
+ ) {
+
+ // from ama06
+ return (boundary_stiffness * boundary_distance - boundary_dampening * dot(normal, vel)) * normal;
+ }
+
+
+ /*
+ COLLISION RESPONSE
+ SIMPLE REFLECTION
+ ................
+ ....^...........
+ Vn..|.../.V.....
+ ....|../........
+ ....|./.........
+ ....|/________>.
+ ........Vt......
+ ................
+
+ Vn = (Vbc * n)Vbc
+ Vt = Vbc –vn
+ Vbc = velocity before collision
+ Vn = normal component of velocity
+ Vt == tangential component of velocity
+ V = (1-u)Vt – eVn
+ u = dynamic friction (affects tangent velocity)
+ e = resilience (affects normal velocity)
+ */
+
+ __device__ float3 calculateFrictionForce(
+ float3 const& vel,
+ float3 const& force,
+ float3 const& normal,
+ float const& friction_kinetic,
+ float const& friction_static_limit
+ )
+ {
+ float3 friction_force = make_float3(0,0,0);
+
+ // the normal part of the force vector (ie, the part that is going "towards" the boundary
+ float3 f_n = force * dot(normal, force);
+ // tangent on the terrain along the force direction (unit vector of tangential force)
+ float3 f_t = force - f_n;
+
+ // the normal part of the velocity vector (ie, the part that is going "towards" the boundary
+ float3 v_n = vel * dot(normal, vel);
+ // tangent on the terrain along the velocity direction (unit vector of tangential velocity)
+ float3 v_t = vel - v_n;
+
+ if((v_t.x + v_t.y + v_t.z)/3.0f > friction_static_limit)
+ friction_force = -v_t;
+ else
+ friction_force = friction_kinetic * -v_t;
+
+ // above static friction limit?
+ // friction_force.x = f_t.x > friction_static_limit ? friction_kinetic * -v_t.x : -v_t.x;
+ // friction_force.y = f_t.y > friction_static_limit ? friction_kinetic * -v_t.y : -v_t.y;
+ // friction_force.z = f_t.z > friction_static_limit ? friction_kinetic * -v_t.z : -v_t.z;
+
+ //TODO; friction should cause energy/heat in contact particles!
+ friction_force = friction_kinetic * -v_t;
+
+ return friction_force;
+ }
+
+};
+
+#endif // __BOUNDARY_CU__
View
132 src/boundary_walls.cu
@@ -0,0 +1,132 @@
+#ifndef __BOUNDARY_WALLS_CU__
+#define __BOUNDARY_WALLS_CU__
+
+#include "boundary.cu"
+
+namespace Boundary {
+
+ namespace Walls {
+
+ __device__ float3 calculateWallsNoPenetrationForce(
+ float3 const& pos,
+ float3 const& vel,
+ float3 const& grid_min,
+ float3 const& grid_max,
+ float const& boundary_distance,
+ float const& boundary_stiffness,
+ float const& boundary_dampening,
+ float const& scale_to_simulation)
+ {
+ float3 repulsion_force = make_float3(0,0,0);
+ float diff;
+
+ // simple limit for "wall" in Y direction (min of simulated volume)
+ diff = boundary_distance - ((pos.y - grid_min.y ) * scale_to_simulation);
+ if (diff > EPSILON) {
+ float3 normal = make_float3(0,1,0);
+ repulsion_force += calculateRepulsionForce(vel, normal, diff, boundary_dampening, boundary_stiffness);
+ }
+
+ // simple limit for "wall" in Y direction (max of simulated volume)
+ diff = boundary_distance - ((grid_max.y - pos.y ) * scale_to_simulation);
+ if (diff > EPSILON) {
+ float3 normal = make_float3(0,-1,0);
+ repulsion_force += calculateRepulsionForce(vel, normal, diff, boundary_dampening, boundary_stiffness);
+ }
+
+ // simple limit for "wall" in Z direction (min of simulated volume)
+ diff = boundary_distance - ((pos.z - grid_min.z ) * scale_to_simulation);
+ if (diff > EPSILON ) {
+ float3 normal = make_float3(0,0,1);
+ repulsion_force += calculateRepulsionForce(vel, normal, diff, boundary_dampening, boundary_stiffness);
+ }
+
+ // simple limit for "wall" in Z direction (max of simulated volume)
+ diff = boundary_distance - ((grid_max.z - pos.z ) * scale_to_simulation);
+ if (diff > EPSILON) {
+ float3 normal = make_float3(0,0,-1);
+ float adj = boundary_stiffness * diff - boundary_dampening * dot(normal, vel);
+ repulsion_force += adj * normal;
+ }
+
+ // simple limit for "wall" in X direction (min of simulated volume)
+ diff = boundary_distance - ((pos.x - grid_min.x ) * scale_to_simulation);
+ if (diff > EPSILON ) {
+ float3 normal = make_float3(1,0,0);
+ repulsion_force += calculateRepulsionForce(vel, normal, diff, boundary_dampening, boundary_stiffness);
+ }
+
+ // simple limit for "wall" in X direction (max of simulated volume)
+ diff = boundary_distance - ((grid_max.x - pos.x ) * scale_to_simulation);
+ if (diff > EPSILON) {
+ float3 normal = make_float3(-1,0,0);
+ repulsion_force += calculateRepulsionForce(vel, normal, diff, boundary_dampening, boundary_stiffness);
+ }
+
+ return repulsion_force;
+ }
+
+
+ __device__ float3 calculateWallsNoSlipForce(
+ float3 const& pos,
+ float3 const& vel,
+ float3 const& force,
+ float3 const& grid_min,
+ float3 const& grid_max,
+ float const& boundary_distance,
+ float const& friction_kinetic,
+ float const& friction_static_limit,
+ float const& scale_to_simulation)
+ {
+ float3 friction_force = make_float3(0,0,0);
+ float diff;
+
+ // simple limit for "wall" in Y direction (min of simulated volume)
+ diff = boundary_distance - ((pos.y - grid_min.y ) * scale_to_simulation);
+ if (diff > EPSILON) {
+ float3 normal = make_float3(0,1,0);
+ friction_force += calculateFrictionForce(vel, force, normal, friction_kinetic, friction_static_limit);
+ }
+
+ // simple limit for "wall" in Y direction (max of simulated volume)
+ diff = boundary_distance - ((grid_max.y - pos.y ) * scale_to_simulation);
+ if (diff > EPSILON) {
+ float3 normal = make_float3(0,-1,0);
+ friction_force += calculateFrictionForce(vel, force, normal, friction_kinetic, friction_static_limit);
+ }
+
+ // simple limit for "wall" in Z direction (min of simulated volume)
+ diff = boundary_distance - ((pos.z - grid_min.z ) * scale_to_simulation);
+ if (diff > EPSILON ) {
+ float3 normal = make_float3(0,0,1);
+ friction_force += calculateFrictionForce(vel, force, normal, friction_kinetic, friction_static_limit);
+ }
+
+ // simple limit for "wall" in Z direction (max of simulated volume)
+ diff = boundary_distance - ((grid_max.z - pos.z ) * scale_to_simulation);
+ if (diff > EPSILON) {
+ float3 normal = make_float3(0,0,-1);
+ friction_force += calculateFrictionForce(vel, force, normal, friction_kinetic, friction_static_limit);
+ }
+
+ // simple limit for "wall" in X direction (min of simulated volume)
+ diff = boundary_distance - ((pos.x - grid_min.x ) * scale_to_simulation);
+ if (diff > EPSILON ) {
+ float3 normal = make_float3(1,0,0);
+ friction_force += calculateFrictionForce(vel, force, normal, friction_kinetic, friction_static_limit);
+ }
+
+ // simple limit for "wall" in X direction (max of simulated volume)
+ diff = boundary_distance - ((grid_max.x - pos.x ) * scale_to_simulation);
+ if (diff > EPSILON) {
+ float3 normal = make_float3(-1,0,0);
+ friction_force += calculateFrictionForce(vel, force, normal, friction_kinetic, friction_static_limit);
+ }
+
+ return friction_force;
+ }
+
+ };
+};
+
+#endif // __BOUNDARY_WALLS_CU__
View
3  src/grid.cuh
@@ -28,12 +28,11 @@
typedef struct sGridData GridData;
typedef struct sGridParams GridParams;
- __constant__ GridParams cudaGridParams;
-
/**
* If cell is empty in grid this value idicates it
*/
#define EMPTY_CELL_VALUE 0xffffffff
+ __constant__ GridParams cudaGridParams;
#endif // __GRID_H__
View
13 src/grid_uniform.cu
@@ -185,7 +185,7 @@ namespace Grid {
////////////////////////////////////////////////////////////////////////////
void Uniform::_calcParams(float cellSize, float gridSize) {
- this->_params.min = make_float3(0.0, 0.0, 0.0);
+ this->_params.min = make_float3(-1.0f, -1.0f, -1.0f);
this->_params.max =
make_float3(
this->_params.min.x + gridSize,
@@ -266,9 +266,14 @@ namespace Grid {
std::cout
<< "Delta["
- << _params.delta.x << ", "
- << _params.delta.y << ", "
- << _params.delta.z << "]"
+ << this->_params.delta.x << ", "
+ << this->_params.delta.y << ", "
+ << this->_params.delta.z << "]"
+ << std::endl;
+
+ std::cout
+ << "Num cells: "
+ << this->_numCells
<< std::endl;
}
View
3  src/kernel.cu
@@ -13,7 +13,8 @@ __global__ void kernel(float4* pos, unsigned int width, unsigned int height, flo
float v = y / (float) height;
u = u*2.0f - 1.0f;
v = v*2.0f - 1.0f;
-;
+ //u *= 0.5f;
+ //v *= 0.5f;
// calculate simple sine wave pattern
float freq = 4.0f;
float w = sinf(u*freq + time) * cosf(v*freq + time) * 0.5f;
View
2  src/particle_system.cu
@@ -69,7 +69,7 @@ struct integrate_functor
pos = normalize(pos) * 0.99;
float velL = length(vel);
vel -= (pos / (len / velL));
- vel *= boundaryDamping;
+ vel *= -boundaryDamping;
}
// store new position and velocity
View
33 src/particles_renderer.cpp
@@ -12,8 +12,8 @@ namespace Particles {
Renderer::Renderer(Simulator* simulator) {
this->_simulator = simulator;
- this->_meshWidth = 8;
- this->_meshHeight = 8;
+ this->_meshWidth = 128;
+ this->_meshHeight = 128;
this->_animate = false;
this->_deltaTime = 0.0;
@@ -46,7 +46,7 @@ namespace Particles {
atexit(SDL_Quit);
this->_initSDL(24, 0);
- this->_simulator->init();
+ this->_simulator->init(this->_meshHeight * this->_meshWidth);
this->_onInit();
this->_runCuda();
this->_render(10);
@@ -237,9 +237,19 @@ namespace Particles {
this->_positionAttribute =
this->_shaderProgram->getAttributeLocation("position");
+ this->_pointScale =
+ this->_shaderProgram->getUniformLocation("pointScale");
+ this->_pointRadius =
+ this->_shaderProgram->getUniformLocation("pointRadius");
+
+ cout << this->_positionAttribute << endl;
+ cout << this->_pointScale << endl;
+ cout << this->_pointRadius << endl;
this->_mvpUniform =
this->_shaderProgram->getUniformLocation("mvp");
+ this->_mvUniform =
+ this->_shaderProgram->getUniformLocation("mv");
/*uint3 gridSize;
gridSize.x = gridSize.y = gridSize.z = 10;
@@ -281,12 +291,16 @@ namespace Particles {
glCullFace(GL_BACK);
glEnable(GL_CULL_FACE);
+ glEnable(GL_POINT_SPRITE_ARB);
+ glTexEnvi(GL_POINT_SPRITE_ARB, GL_COORD_REPLACE_ARB, GL_TRUE);
+ glEnable(GL_VERTEX_PROGRAM_POINT_SIZE_NV);
+
// Calculate ModelViewProjection matrix
glm::mat4 projection =
glm::perspective(
45.0f,
(float) this->_windowWidth / (float) this->_windowHeight,
- 1.0f,
+ 0.001f,
1000.0f
);
@@ -309,8 +323,16 @@ namespace Particles {
this->_shaderProgram->use();
// Set matrices
+ glUniformMatrix4fv(this->_mvUniform, 1, GL_FALSE, glm::value_ptr(mv));
glUniformMatrix4fv(this->_mvpUniform, 1, GL_FALSE, glm::value_ptr(mvp));
glEnableVertexAttribArray(this->_positionAttribute);
+ glUniform1f(
+ this->_pointScale,
+ this->_windowWidth / tanf(45.0f*0.5f*(float)M_PI/180.0f)
+ );
+
+ //cout << (this->_windowWidth / tanf(45.0f*(float)M_PI/180.0f)) << endl;
+ glUniform1f(this->_pointRadius, 1000.f);
// Draw data
//glBindBuffer(GL_ARRAY_BUFFER, this->_vbo);
@@ -325,7 +347,8 @@ namespace Particles {
);
glDrawArrays(GL_POINTS, 0, this->_meshWidth * this->_meshHeight);
-
+ glDisable(GL_POINT_SPRITE_ARB);
+ glDisable(GL_VERTEX_PROGRAM_POINT_SIZE_NV);
SDL_GL_SwapBuffers();
this->_deltaTime += 0.01f;
View
3  src/particles_renderer.h
@@ -72,6 +72,9 @@ namespace Particles {
GLuint _positionAttribute;
GLuint _mvpUniform;
+ GLuint _mvUniform;
+ GLuint _pointScale;
+ GLuint _pointRadius;
int _windowWidth;
int _windowHeight;
View
2  src/particles_simulator.h
@@ -14,7 +14,7 @@ namespace Particles{
Simulator();
virtual ~Simulator();
- virtual void init() = 0;
+ virtual void init(unsigned int numParticles) = 0;
virtual void stop() = 0;
virtual void update() = 0;
virtual float* getPositions() = 0;
View
4 src/shader_program.cpp
@@ -40,6 +40,10 @@ GLuint ShaderProgram::getUniformLocation (string uniform) {
return glGetUniformLocation(this->program, uniform.c_str());
}
+GLuint ShaderProgram::getReference() {
+ return this->program;
+}
+
void ShaderProgram::add(const GLenum type, const char *source) {
GLuint shader = this->compile(type, this->load(source).c_str());
View
3  src/shader_program.h
@@ -80,6 +80,9 @@ class ShaderProgram {
*/
GLuint getUniformLocation (string uniform);
+
+ GLuint getReference();
+
protected:
/**
* Add shader
View
2  src/sph.h
@@ -11,12 +11,14 @@ namespace SPH {
Forces,
Positions,
Pressures,
+ Velevals,
Velocities,
SortedColors,
SortedDensities,
SortedForces,
SortedPositions,
SortedPressures,
+ SortedVelevals,
SortedVelocities
};
View
3  src/sph_density.cu
@@ -37,7 +37,8 @@ namespace SPH {
) {
data.density +=
SPH::Kernels::Poly6::getVariable(
- cudaPrecalcParams.smoothLenSq, r, rLenSq);
+ cudaPrecalcParams.smoothLenSq, r, rLenSq
+ );
}
////////////////////////////////////////////////////////////////////
View
54 src/sph_force.cu
@@ -10,6 +10,7 @@ namespace SPH {
public:
////////////////////////////////////////////////////////////////////
+
struct Data {
float3 veleval;
float density;
@@ -49,9 +50,31 @@ namespace SPH {
float const &rLen,
float const &rLenSq
) {
- data.density +=
- SPH::Kernels::Poly6::getVariable(
- cudaPrecalcParams.smoothLenSq, r, rLenSq);
+ data.velevalN = make_float3(data.sorted.veleval[indexN]);
+ data.densityN = data.sorted.density[indexN];
+ data.pressureN = data.sorted.pressure[indexN];
+
+ data.pressureForce +=
+ (
+ (data.pressure + data.pressureN) /
+ (data.densityN * data.density)
+ ) *
+ SPH::Kernels::Spiky::getGradientVariable(
+ cudaFluidParams.smoothingLength,
+ r,
+ rLen
+ );
+
+ data.viscosityForce +=
+ (
+ (data.velevalN - data. velevalN ) /
+ (data.densityN * data.densityN)
+ ) *
+ SPH::Kernels::Viscosity::getLaplacianVariable(
+ cudaFluidParams.smoothingLength,
+ r,
+ rLen
+ );
}
////////////////////////////////////////////////////////////////////
@@ -60,18 +83,19 @@ namespace SPH {
Data &data,
uint const &index
) {
- float density =
- cudaFluidParams.particleMass *
- cudaPrecalcParams.poly6Coeff *
- data.density;
-
- density = max(1.0, density);
-
- data.sorted.density[index] = density;
- data.sorted.pressure[index] =
- cudaFluidParams.restPressure +
- cudaFluidParams.gasStiffness *
- (density - cudaFluidParams.restDensity);
+
+ float3 sumForce = (
+ cudaPrecalcParams.pressurePrecalc *
+ data.pressureForce +
+ cudaPrecalcParams.viscosityPrecalc *
+ data.viscosityForce
+ );
+
+ // Calculate the force, the particle_mass is added here because
+ // it is constant and thus there is no need for it to
+ // be inside the sum loop.
+ data.sorted.force[index] =
+ make_float4(sumForce * cudaFluidParams.particleMass, 1.0f);
}
////////////////////////////////////////////////////////////////////
View
54 src/sph_kernel.cu
@@ -4,6 +4,7 @@
__constant__ SPH::FluidParams cudaFluidParams;
__constant__ SPH::PrecalcParams cudaPrecalcParams;
+#include "boundary_walls.cu"
#include "grid.cuh"
#include "grid_utils.cu"
#include "sph_density.cu"
@@ -41,26 +42,73 @@ namespace SPH {
float3 externalForce = make_float3(0.0f, 0.0f, 0.0f);
externalForce.y -= 9.8f;
+ // add no-penetration force due to "walls"
+ /*externalForce += Boundary::Walls::calculateWallsNoPenetrationForce(
+ position, veleval,
+ cudaGridParams.min,
+ cudaGridParams.max,
+ cudaFluidParams.boundaryDistance,
+ cudaFluidParams.boundaryStiffness,
+ cudaFluidParams.boundaryDampening,
+ cudaFluidParams.scaleToSimulation);
+
+ // add no-slip force due to "walls"
+ externalForce += Boundary::Walls::calculateWallsNoSlipForce(
+ position, veleval, force + externalForce,
+ cudaGridParams.min,
+ cudaGridParams.max,
+ cudaFluidParams.boundaryDistance,
+ cudaFluidParams.frictionKinetic/deltaTime,
+ cudaFluidParams.frictionStaticLimit,
+ cudaFluidParams.scaleToSimulation);
+ */
+
+
+
float3 f = force + externalForce;
float speed = length(force);
if (speed > cudaFluidParams.velocityLimit) {
- f -= cudaFluidParams.velocityLimit / speed;
+ f *= cudaFluidParams.velocityLimit / speed;
}
float3 vnext = velocity + f * deltaTime;
veleval = (velocity + vnext) * 0.5;
velocity = veleval;
- //position += vnext * (deltaTime);
position += vnext * (deltaTime / cudaFluidParams.scaleToSimulation);
uint sortedIndex = gridData.index[index];
+ if (position.x < cudaGridParams.min.x) {
+ position.x = cudaGridParams.min.x;
+ //velocity *= -cudaFluidParams.boundaryDampening;
+ }
+ if (position.x > cudaGridParams.max.x) {
+ position.x = cudaGridParams.max.x;
+ //velocity *= -cudaFluidParams.boundaryDampening;
+ }
+ if (position.y < cudaGridParams.min.y) {
+ position.y = cudaGridParams.min.y;
+ //velocity *= -cudaFluidParams.boundaryDampening;
+ }
+ if (position.y > cudaGridParams.max.y) {
+ position.y = cudaGridParams.max.y;
+ //velocity *= -cudaFluidParams.boundaryDampening;
+ }
+ if (position.z < cudaGridParams.min.z) {
+ position.z = cudaGridParams.min.z;
+ //velocity *= -cudaFluidParams.boundaryDampening;
+ }
+ if (position.z > cudaGridParams.max.z) {
+ position.z = cudaGridParams.max.z;
+ //velocity *= -cudaFluidParams.boundaryDampening;
+ }
+
data.position[sortedIndex] = make_float4(position, 1.0);
data.velocity[sortedIndex] = make_float4(velocity, 1.0);
- //data.veleval[sortedIndex] = make_float4(veleval, 1.0);
+ data.veleval[sortedIndex] = make_float4(veleval, 1.0);
}
View
2  src/sph_neighbours.cu
@@ -28,7 +28,7 @@ namespace SPH {
float const &rLen,
float const &rLenSq
) {
- C::process(index, indexN, r, rLen, rLenSq);
+ C::process(data, index, indexN, r, rLen, rLenSq);
}
////////////////////////////////////////////////////////////////////
View
92 src/sph_simulator.cu
@@ -34,31 +34,27 @@ namespace SPH {
////////////////////////////////////////////////////////////////////////////
- void Simulator::init() {
- this->_numParticles = 8*8;
+ void Simulator::init(uint numParticles) {
+ this->_numParticles = numParticles;
this->_createBuffers();
- this->_grid = new Grid::Uniform();
- this->_grid->allocate(this->_numParticles, 1.0f/128.0f, 1.0f);
- this->_grid->printParams();
-
// DATABASE
this->_database = new Database();
this->_database->addUpdateCallback(this);
this->_database
->insert(ParticleNumber, "Particles", this->_numParticles)
- ->insert(GridSize, "Grid size", 256.0)
- ->insert(Timestep, "Timestep", 0.0f, 1.0f, 0.002f)
+ ->insert(GridSize, "Grid size", 2.0f)
+ ->insert(Timestep, "Timestep", 0.0f, 1.0f, 0.01f)
->insert(RestDensity, "Rest density", 0.0f, 10000.0f, 1000.0f)
->insert(RestPressure, "Rest pressure", 0.0f, 10000.0f, 0.0f)
->insert(GasStiffness, "Gas Stiffness", 0.001f, 10.0f, 1.0f)
->insert(Viscosity, "Viscosity", 0.0f, 100.0f, 1.0f)
- ->insert(BoundaryDampening, "Bound. damp.", 0.0f, 10000.0f, 256.0f)
+ ->insert(BoundaryDampening, "Bound. damp.", 0.0f, 10000.0f, 0.5f)
->insert(BoundaryStiffness, "Bound. stiff.", 0.0f, 100000.0f, 20000.0f)
->insert(VelocityLimit, "Veloc. limit", 0.0f, 10000.0f, 600.0f)
- ->insert(SimulationScale, "Sim. scale", 0.0f, 1.0f, 0.01f)
+ ->insert(SimulationScale, "Sim. scale", 0.0f, 1.0f, 1.0f)
->insert(KineticFriction, "Kinet. fric.", 0.0f, 10000.0f, 0.0f)
->insert(StaticFrictionLimit, "Stat. Fric. Lim.", 0.0f, 10000.0f, 0.0f);
@@ -72,8 +68,8 @@ namespace SPH {
);
float boundaryDist = 0.5 * particleRestDist;
float smoothingLength = 2.0 * particleRestDist;
- float cellSize =
- smoothingLength * this->_database->selectValue(SimulationScale);
+ float cellSize = 1.0f/64.0f;
+ //smoothingLength * this->_database->selectValue(SimulationScale);
this->_database
->insert(ParticleMass, "Particle mass", particleMass)
@@ -84,6 +80,15 @@ namespace SPH {
//this->_database->print();
+ this->_grid = new Grid::Uniform();
+ this->_grid->allocate(
+ this->_numParticles,
+ cellSize,
+ this->_database->selectValue(GridSize)
+ );
+ this->_grid->printParams();
+ this->_gridParams = this->_grid->getParams();
+
this->_updateParams();
@@ -150,9 +155,6 @@ namespace SPH {
minBlockSize = 416;
Utils::computeGridSize(numParticles, minBlockSize, numBlocks, numThreads);
-
- //this->_particleData.position = (float4*)this->getPositions();
-
Kernel::integrate<Data><<<numBlocks, numThreads>>>(
numParticles,
deltaTime,
@@ -170,22 +172,22 @@ namespace SPH {
this->_grid->sort();
this->_orderData();
- Buffer::Memory<uint>* buffer =
- new Buffer::Memory<uint>(new Buffer::Allocator(), Buffer::Host);
+ /*Buffer::Memory<float>* buffer =
+ new Buffer::Memory<float>(new Buffer::Allocator(), Buffer::Host);
buffer->allocate(this->_numParticles);
GridData gridData = this->_grid->getData();
- cudaMemcpy(buffer->get(), gridData.index, this->_numParticles * sizeof(uint), cudaMemcpyDeviceToHost);
+ cudaMemcpy(buffer->get(), this->_sortedData.density, this->_numParticles * sizeof(float), cudaMemcpyDeviceToHost);
- uint* e = buffer->get();
+ float* e = buffer->get();
for(uint i=0;i< this->_numParticles; i++) {
- //cout << e[i] << endl;
- }
+ cout << e[i] << endl;
+ }*/
- Buffer::Memory<float4>* posBuffer =
+ /*Buffer::Memory<float4>* posBuffer =
new Buffer::Memory<float4>(new Buffer::Allocator(), Buffer::Host);
posBuffer->allocate(this->_numParticles);
@@ -199,9 +201,12 @@ namespace SPH {
std::cout << pos[i].x << " " << pos[i].y << " " << pos[i].z << endl;
}
- //this->_step1();
- //this->_step2();
+ std::cout << "____________________" << std::endl;
+ */
+ this->_step1();
+ this->_step2();
this->integrate(this->_numParticles, this->_database->selectValue(Timestep));
+
//cutilSafeCall(cutilDeviceSynchronize());
}
@@ -223,6 +228,7 @@ namespace SPH {
Buffer::Memory<float4>* force = new Buffer::Memory<float4>();
Buffer::Vertex<float4>* position = new Buffer::Vertex<float4>();
Buffer::Memory<float>* pressure = new Buffer::Memory<float>();
+ Buffer::Memory<float4>* veleval = new Buffer::Memory<float4>();
Buffer::Memory<float4>* velocity = new Buffer::Memory<float4>();
Buffer::Memory<float4>* sColor = new Buffer::Memory<float4>();
@@ -230,6 +236,7 @@ namespace SPH {
Buffer::Memory<float4>* sForce = new Buffer::Memory<float4>();
Buffer::Vertex<float4>* sPosition = new Buffer::Vertex<float4>();
Buffer::Memory<float>* sPressure = new Buffer::Memory<float>();
+ Buffer::Memory<float4>* sVeleval = new Buffer::Memory<float4>();
Buffer::Memory<float4>* sVelocity = new Buffer::Memory<float4>();
this->_positionsVBO = position->getVBO();
@@ -240,12 +247,14 @@ namespace SPH {
->addBuffer(Forces, (Buffer::Abstract<void>*) force)
->addBuffer(Positions, (Buffer::Abstract<void>*) position)
->addBuffer(Pressures, (Buffer::Abstract<void>*) pressure)
+ ->addBuffer(Velevals, (Buffer::Abstract<void>*) veleval)
->addBuffer(Velocities, (Buffer::Abstract<void>*) velocity)
->addBuffer(SortedColors, (Buffer::Abstract<void>*) sColor)
->addBuffer(SortedDensities, (Buffer::Abstract<void>*) sDensity)
->addBuffer(SortedForces, (Buffer::Abstract<void>*) sForce)
->addBuffer(SortedPositions, (Buffer::Abstract<void>*) sPosition)
->addBuffer(SortedPressures, (Buffer::Abstract<void>*) sPressure)
+ ->addBuffer(SortedVelevals, (Buffer::Abstract<void>*) sVeleval)
->addBuffer(SortedVelocities, (Buffer::Abstract<void>*) sVelocity);
this->_bufferManager->allocateBuffers(this->_numParticles);
@@ -270,6 +279,7 @@ namespace SPH {
this->_particleData.force = force->get();
this->_particleData.position = position->get();
this->_particleData.pressure = pressure->get();
+ this->_particleData.veleval = veleval->get();
this->_particleData.velocity = velocity->get();
this->_sortedData.color = sColor->get();
@@ -277,6 +287,7 @@ namespace SPH {
this->_sortedData.force = sForce->get();
this->_sortedData.position = sPosition->get();
this->_sortedData.pressure = sPressure->get();
+ this->_sortedData.veleval = sVeleval->get();
this->_sortedData.velocity = sVelocity->get();
this->_bufferManager->unbindBuffers();
@@ -365,7 +376,7 @@ namespace SPH {
Kernels::Spiky::getGradientConstant(smoothLen);
this->_precalcParams.viscosityLapCoeff =
- Kernels::Spiky::getLaplacianConstant(smoothLen);
+ Kernels::Viscosity::getLaplacianConstant(smoothLen);
this->_precalcParams.pressurePrecalc =
-0.5 * this->_precalcParams.spikyGradCoeff;
@@ -374,6 +385,29 @@ namespace SPH {
this->_fluidParams.viscosity *
this->_precalcParams.viscosityLapCoeff;
+ // DEBUG
+ /*
+ cout
+ << "SmoothLen: "
+ << this->_precalcParams.smoothLenSq
+ << endl
+ << "Poly6Coeff: "
+ << this->_precalcParams.poly6Coeff
+ << endl
+ << "SpikyGradCoeff: "
+ << this->_precalcParams.spikyGradCoeff
+ << endl
+ << "ViscosityLapCoeff: "
+ << this->_precalcParams.viscosityLapCoeff
+ << endl
+ << "PressurePrecalc: "
+ << this->_precalcParams.pressurePrecalc
+ << endl
+ << "ViscosityPrecalc: "
+ << this->_precalcParams.viscosityPrecalc
+ << endl;
+ */
+
// Copy parameters to GPU's constant memory
// declarations of symbols are in sph_kernel.cu
CUDA_SAFE_CALL(
@@ -392,6 +426,14 @@ namespace SPH {
)
);
+ CUDA_SAFE_CALL(
+ cudaMemcpyToSymbol(
+ cudaGridParams,
+ &this->_gridParams,
+ sizeof(GridParams)
+ )
+ );
+
CUDA_SAFE_CALL(cudaThreadSynchronize());
}
View
3  src/sph_simulator.cuh
@@ -26,7 +26,7 @@ namespace SPH {
* else can cause segmentation fault
* Must be called as first method
*/
- void init();
+ void init(uint numParticles);
void stop();
/**
*
@@ -55,6 +55,7 @@ namespace SPH {
FluidParams _fluidParams;
PrecalcParams _precalcParams;
+ GridParams _gridParams;
private:
void _createBuffers();
Please sign in to comment.
Something went wrong with that request. Please try again.