GopherProxy

	ttesting pressure gradient force - sphere - GPU-based 3D discrete element metho…
	git clone git://src.adamsgaard.dk/sphere
	Log
	Files
	Refs
	LICENSE
	---
	commit ea7d8aa29788bcbca20557265e50c4247dfc4a26
	parent 48692b733e6e121ddf0f6bd984db99ee677ba2cb
	Author: Anders Damsgaard <[email protected]>
	Date: Fri, 2 May 2014 13:14:46 +0200

	ttesting pressure gradient force

	Diffstat:
	M src/device.cu \| 1 +
	M src/navierstokes.cuh \| 23 ++++++++++++++++++-----
	M tests/stokes_law.py \| 2 +-

	3 files changed, 20 insertions(+), 6 deletions(-)
	---
	diff --git a/src/device.cu b/src/device.cu
	t@@ -884,6 +884,7 @@ __host__ void DEM::startTime()
	applyParticleInteractionForce<<<dimGridFluid, dimBlockFluid>>>(
	dev_ns_fi,
	dev_ns_phi,
	+ dev_ns_p,
	dev_gridParticleIndex,
	dev_cellStart,
	dev_cellEnd,
	diff --git a/src/navierstokes.cuh b/src/navierstokes.cuh
	t@@ -2632,6 +2632,7 @@ __global__ void findInteractionForce(
	__global__ void applyParticleInteractionForce(
	Float3* dev_ns_fi, // in
	Float* dev_ns_phi, // in
	+ Float* dev_ns_p, // in
	unsigned int* dev_gridParticleIndex, // in
	unsigned int* dev_cellStart, // in
	unsigned int* dev_cellEnd, // in
	t@@ -2643,13 +2644,24 @@ __global__ void applyParticleInteractionForce(
	const unsigned int y = blockDim.y * blockIdx.y + threadIdx.y;
	const unsigned int z = blockDim.z * blockIdx.z + threadIdx.z;

	+ // Grid dimensions
	+ const unsigned int nx = devC_grid.num[0];
	+ const unsigned int ny = devC_grid.num[1];
	+ const unsigned int nz = devC_grid.num[2];
	+
	+ // Cell sizes
	+ const Float dx = devC_grid.L[0]/nx;
	+ const Float dy = devC_grid.L[1]/ny;
	+ const Float dz = devC_grid.L[2]/nz;
	+
	// Check that we are not outside the fluid grid
	- if (x < devC_grid.num[0] && y < devC_grid.num[1] && z < devC_grid.num[2]) {
	+ if (x < nx && y < ny && z < nz) {

	const unsigned int cellidx = idx(x,y,z);

	__syncthreads();
	const Float3 fi = dev_ns_fi[cellidx];
	+ const Float3 grad_p = gradient(dev_ns_p, x, y, z, dx, dy, dz);

	// apply to all particle in the cell
	// Calculate linear cell ID
	t@@ -2675,10 +2687,11 @@ __global__ void applyParticleInteractionForce(
	r = dev_x_sorted[i].w; // radius
	//phi = dev_ns_phi[idx(x,y,z)];

	- // this term could include the pressure gradient
	- //fd = (-grad_p + fi/(1.0 - phi))(4.0/3.0M_PIrr*r);
	- //fd = (fi/(1.0 - phi))(4.0/3.0M_PIrr*r);
	- fd = fi(4.0/3.0M_PIrr*r);
	+ // stokes drag force
	+ //fd = fi(4.0/3.0M_PIrr*r);
	+
	+ // pressure gradient force + stokes drag force
	+ fd = (-1.0grad_p + fi)(4.0/3.0M_PIrrr);

	__syncthreads();
	dev_force[origidx] += MAKE_FLOAT4(fd.x, fd.y, fd.z, 0.0);
	diff --git a/tests/stokes_law.py b/tests/stokes_law.py
	t@@ -31,7 +31,7 @@ orig.vel[0,2] = -0.1
	#orig.vel[0,2] = -0.001
	#orig.setBeta(0.5)
	orig.setTolerance(1.0e-4)
	-#orig.setDEMstepsPerCFDstep(100)
	+orig.setDEMstepsPerCFDstep(100)
	orig.run(dry=True)
	orig.run(verbose=True)
	py = sphere.sim(sid = orig.sid, fluid = True)