GopherProxy

	tBroad commit - sphere - GPU-based 3D discrete element method algorithm with op…
	git clone git://src.adamsgaard.dk/sphere
	Log
	Files
	Refs
	LICENSE
	---
	commit 10800d34259399b383022bcd6aa50811023445c5
	parent 1736587a16f7500de97a3cf3ff04eb47558a5546
	Author: Anders Damsgaard <[email protected]>
	Date: Thu, 30 Aug 2012 14:04:54 +0200

	Broad commit

	Diffstat:
	M src/cohesion.cuh \| 2 +-
	M src/contactsearch.cuh \| 52 ++++++++++++++++-------------…
	M src/device.cu \| 1 +
	M src/integration.cuh \| 4 ++--

	4 files changed, 30 insertions(+), 29 deletions(-)
	---
	diff --git a/src/cohesion.cuh b/src/cohesion.cuh
	t@@ -66,7 +66,7 @@ __device__ void bondLinear(Float3* N, Float3* T, Float* es_d…

	if (length(vel_t_ab) > 0.f) {
	// Shear force component: Viscous
	- f_t = -1.0f * devC_gamma_s * vel_t_ab;
	+ f_t = -1.0f * devC_gamma_t * vel_t_ab;

	// Shear friction production rate [W]
	//*es_dot += -dot(vel_t_ab, f_t);
	diff --git a/src/contactsearch.cuh b/src/contactsearch.cuh
	t@@ -8,7 +8,7 @@
	// Calculate the distance modifier between a contact pair. The modifier
	// accounts for periodic boundaries. If the target cell lies outside
	// the grid, it returns -1.
	-__device__ int findDistMod(int3 targetCell, Float3* distmod)
	+__device__ int findDistMod(int3* targetCell, Float3* distmod)
	{
	// Check whether x- and y boundaries are to be treated as periodic
	// 1: x- and y boundaries periodic
	t@@ -16,22 +16,22 @@ __device__ int findDistMod(int3 targetCell, Float3* distmo…
	if (devC_periodic == 1) {

	// Periodic x-boundary
	- if (targetCell.x < 0) {
	- targetCell.x = devC_num[0] - 1;
	+ if (targetCell->x < 0) {
	+ targetCell->x = devC_num[0] - 1;
	*distmod += MAKE_FLOAT3(devC_L[0], 0.0f, 0.0f);
	}
	- if (targetCell.x == devC_num[0]) {
	- targetCell.x = 0;
	+ if (targetCell->x == devC_num[0]) {
	+ targetCell->x = 0;
	*distmod -= MAKE_FLOAT3(devC_L[0], 0.0f, 0.0f);
	}

	// Periodic y-boundary
	- if (targetCell.y < 0) {
	- targetCell.y = devC_num[1] - 1;
	+ if (targetCell->y < 0) {
	+ targetCell->y = devC_num[1] - 1;
	*distmod += MAKE_FLOAT3(0.0f, devC_L[1], 0.0f);
	}
	- if (targetCell.y == devC_num[1]) {
	- targetCell.y = 0;
	+ if (targetCell->y == devC_num[1]) {
	+ targetCell->y = 0;
	*distmod -= MAKE_FLOAT3(0.0f, devC_L[1], 0.0f);
	}

	t@@ -40,17 +40,17 @@ __device__ int findDistMod(int3 targetCell, Float3* distmo…
	} else if (devC_periodic == 2) {

	// Periodic x-boundary
	- if (targetCell.x < 0) {
	- targetCell.x = devC_num[0] - 1;
	+ if (targetCell->x < 0) {
	+ targetCell->x = devC_num[0] - 1;
	*distmod += MAKE_FLOAT3(devC_L[0], 0.0f, 0.0f);
	}
	- if (targetCell.x == devC_num[0]) {
	- targetCell.x = 0;
	+ if (targetCell->x == devC_num[0]) {
	+ targetCell->x = 0;
	*distmod -= MAKE_FLOAT3(devC_L[0], 0.0f, 0.0f);
	}

	// Hande out-of grid cases on y-axis
	- if (targetCell.y < 0 \|\| targetCell.y == devC_num[1])
	+ if (targetCell->y < 0 \|\| targetCell->y == devC_num[1])
	return -1;


	t@@ -58,14 +58,14 @@ __device__ int findDistMod(int3 targetCell, Float3* distmo…
	} else {

	// Hande out-of grid cases on x- and y-axes
	- if (targetCell.x < 0 \|\| targetCell.x == devC_num[0])
	+ if (targetCell->x < 0 \|\| targetCell->x == devC_num[0])
	return -1;
	- if (targetCell.y < 0 \|\| targetCell.y == devC_num[1])
	+ if (targetCell->y < 0 \|\| targetCell->y == devC_num[1])
	return -1;
	}

	// Handle out-of-grid cases on z-axis
	- if (targetCell.z < 0 \|\| targetCell.z == devC_num[2])
	+ if (targetCell->z < 0 \|\| targetCell->z == devC_num[2])
	return -1;

	// Return successfully
	t@@ -95,7 +95,7 @@ __device__ void overlapsInCell(int3 targetCell,
	// Get distance modifier for interparticle
	// vector, if it crosses a periodic boundary
	Float3 distmod = MAKE_FLOAT3(0.0f, 0.0f, 0.0f);
	- if (findDistMod(targetCell, &distmod) == -1)
	+ if (findDistMod(&targetCell, &distmod) == -1)
	return; // Target cell lies outside the grid


	t@@ -189,7 +189,7 @@ __device__ void findContactsInCell(int3 targetCell,
	// Get distance modifier for interparticle
	// vector, if it crosses a periodic boundary
	Float3 distmod = MAKE_FLOAT3(0.0f, 0.0f, 0.0f);
	- if (findDistMod(targetCell, &distmod) == -1)
	+ if (findDistMod(&targetCell, &distmod) == -1)
	return; // Target cell lies outside the grid


	t@@ -403,9 +403,9 @@ __global__ void interact(unsigned int* dev_gridParticleInd…
	Float delta_n, x_ab_length, radius_b;
	Float3 x_ab;
	Float4 x_b, distmod;
	- //Float4 vel_a = dev_vel_sorted[idx_a];
	- //Float4 angvel4_a = dev_angvel_sorted[idx_a];
	- //Float3 angvel_a = MAKE_FLOAT3(angvel4_a.x, angvel4_a.y, angvel4_a.z);
	+ Float4 vel_a = dev_vel_sorted[idx_a];
	+ Float4 angvel4_a = dev_angvel_sorted[idx_a];
	+ Float3 angvel_a = MAKE_FLOAT3(angvel4_a.x, angvel4_a.y, angvel4_a.z);

	// Loop over all possible contacts, and remove contacts
	// that no longer are valid (delta_n > 0.0)
	t@@ -433,15 +433,15 @@ __global__ void interact(unsigned int* dev_gridParticleI…
	if (delta_n < 0.0f) {
	//cuPrintf("\nProcessing contact, idx_a_orig = %u, idx_b_orig = %u…
	// idx_a_orig, idx_b_orig, i, delta_n);
	- contactLinearViscous(&N, &T, &es_dot, &p,
	+ /*contactLinearViscous(&N, &T, &es_dot, &p,
	idx_a_orig, idx_b_orig,
	dev_vel,
	dev_angvel,
	radius_a, radius_b,
	x_ab, x_ab_length,
	delta_n, devC_kappa);
	- dev_delta_t[mempos] = MAKE_FLOAT4(0.0f, 0.0f, 0.0f, 0.0f);
	- /*contactLinear(&N, &T, &es_dot, &p,
	+ dev_delta_t[mempos] = MAKE_FLOAT4(0.0f, 0.0f, 0.0f, 0.0f);*/
	+ contactLinear(&N, &T, &es_dot, &p,
	idx_a_orig,
	idx_b_orig,
	vel_a,
	t@@ -451,7 +451,7 @@ __global__ void interact(unsigned int* dev_gridParticleInd…
	radius_a, radius_b,
	x_ab, x_ab_length,
	delta_n, dev_delta_t,
	- mempos);*/
	+ mempos);
	} else {
	__syncthreads();
	// Remove this contact (there is no particle with index=np)
	diff --git a/src/device.cu b/src/device.cu
	t@@ -136,6 +136,7 @@ __host__ void transferToConstantMemory(Particles* p,
	cudaMemcpyToSymbol("devC_db", &params->db, sizeof(Float));
	cudaMemcpyToSymbol("devC_V_b", &params->V_b, sizeof(Float));
	cudaMemcpyToSymbol("devC_shearmodel", &params->shearmodel, sizeof(unsigned…
	+ cudaMemcpyToSymbol("devC_wmode", &params->wmode, sizeof(int)*MAXWALLS);
	} else {
	//printf("(params.global == %d) ", params.global);
	// Copy params structure with individual physical values from host to glob…
	diff --git a/src/integration.cuh b/src/integration.cuh
	t@@ -193,12 +193,12 @@ __global__ void integrateWalls(Float4* dev_w_nx,

	// Calculate resulting acceleration of wall
	// (Wall mass is stored in w component of position Float4)
	- acc = (w_mvfd.z+N)/w_mvfd.x;
	+ acc = (w_mvfd.z + N)/w_mvfd.x;

	// Update linear velocity
	w_mvfd.y += acc * dt;

	- // Wall BC is controlled by velocity
	+ // Wall BC is controlled by velocity, which should not change
	} else if (wmode == 1) {
	acc = 0.0f;
	}