tadded Zdisplacement component - sphere - GPU-based 3D discrete element method … | |
git clone git://src.adamsgaard.dk/sphere | |
Log | |
Files | |
Refs | |
LICENSE | |
--- | |
commit 45015224a9aa5cdd2491c36a4d671a7d2972a877 | |
parent af3360fc4686a3ebc3423e44b148a4bcebdc9ec1 | |
Author: Anders Damsgaard <[email protected]> | |
Date: Sun, 22 Jun 2014 09:11:05 +0200 | |
added Zdisplacement component | |
Diffstat: | |
M doc/html/python_api.html | 2 +- | |
M doc/html/searchindex.js | 4 ++-- | |
M doc/pdf/sphere.pdf | 0 | |
M python/sphere.py | 47 ++++++++++++++++++-----------… | |
M src/datatypes.h | 2 +- | |
M src/device.cu | 14 +++++++------- | |
M src/file_io.cpp | 12 +++++++----- | |
M src/integration.cuh | 21 ++++++++++++--------- | |
M src/raytracer.cuh | 4 ++-- | |
M src/sphere.cpp | 2 +- | |
M src/sphere.h | 2 +- | |
M tests/io_tests_fluid.py | 2 +- | |
12 files changed, 63 insertions(+), 49 deletions(-) | |
--- | |
diff --git a/doc/html/python_api.html b/doc/html/python_api.html | |
t@@ -328,7 +328,7 @@ are returned.</p> | |
<dl class="method"> | |
<dt id="sphere.sim.addParticle"> | |
-<tt class="descname">addParticle</tt><big>(</big><em>x</em>, <em>radius</em>, … | |
+<tt class="descname">addParticle</tt><big>(</big><em>x</em>, <em>radius</em>, … | |
<dd><p>Add a single particle to the simulation object. The only required | |
parameters are the position (x) and the radius (radius).</p> | |
<table class="docutils field-list" frame="void" rules="none"> | |
diff --git a/doc/html/searchindex.js b/doc/html/searchindex.js | |
t@@ -1 +1 @@ | |
-Search.setIndex({objects:{"":{sphere:[5,0,1,""]},sphere:{status:[5,1,1,""],con… | |
-\ No newline at end of file | |
+Search.setIndex({objects:{"":{sphere:[5,0,1,""]},sphere:{status:[5,1,1,""],con… | |
+\ No newline at end of file | |
diff --git a/doc/pdf/sphere.pdf b/doc/pdf/sphere.pdf | |
Binary files differ. | |
diff --git a/python/sphere.py b/python/sphere.py | |
t@@ -86,7 +86,7 @@ class sim: | |
self.radius = numpy.ones(self.np, dtype=numpy.float64) | |
# The sums of x and y movement [m] | |
- self.xysum = numpy.zeros((self.np, 2), dtype=numpy.float64) | |
+ self.xyzsum = numpy.zeros((self.np, 3), dtype=numpy.float64) | |
# The linear velocities [m/s] | |
self.vel = numpy.zeros((self.np, self.nd), dtype=numpy.float64) | |
t@@ -375,7 +375,7 @@ class sim: | |
elif ((self.radius != other.radius).any()): | |
print(15) | |
return 15 | |
- elif ((self.xysum != other.xysum).any()): | |
+ elif ((self.xyzsum != other.xyzsum).any()): | |
print(16) | |
return 16 | |
elif ((self.vel != other.vel).any()): | |
t@@ -592,7 +592,7 @@ class sim: | |
def addParticle(self, | |
x, | |
radius, | |
- xysum = numpy.zeros(2), | |
+ xyzsum = numpy.zeros(3), | |
vel = numpy.zeros(3), | |
fixvel = numpy.zeros(1), | |
force = numpy.zeros(3), | |
t@@ -642,7 +642,7 @@ class sim: | |
self.x = numpy.append(self.x, [x], axis=0) | |
self.radius = numpy.append(self.radius, radius) | |
self.vel = numpy.append(self.vel, [vel], axis=0) | |
- self.xysum = numpy.append(self.xysum, [xysum], axis=0) | |
+ self.xyzsum = numpy.append(self.xyzsum, [xyzsum], axis=0) | |
self.fixvel = numpy.append(self.fixvel, fixvel) | |
self.force = numpy.append(self.force, [force], axis=0) | |
self.angpos = numpy.append(self.angpos, [angpos], axis=0) | |
t@@ -668,7 +668,7 @@ class sim: | |
self.x = numpy.delete(self.x, i, axis=0) | |
self.radius = numpy.delete(self.radius, i) | |
self.vel = numpy.delete(self.vel, i, axis=0) | |
- self.xysum = numpy.delete(self.xysum, i, axis=0) | |
+ self.xyzsum = numpy.delete(self.xyzsum, i, axis=0) | |
self.fixvel = numpy.delete(self.fixvel, i) | |
self.force = numpy.delete(self.force, i, axis=0) | |
self.angpos = numpy.delete(self.angpos, i, axis=0) | |
t@@ -688,7 +688,7 @@ class sim: | |
self.np[0] = 0 | |
self.x = numpy.zeros((self.np, self.nd), dtype=numpy.float64) | |
self.radius = numpy.ones(self.np, dtype=numpy.float64) | |
- self.xysum = numpy.zeros((self.np, 2), dtype=numpy.float64) | |
+ self.xyzsum = numpy.zeros((self.np, 3), dtype=numpy.float64) | |
self.vel = numpy.zeros((self.np, self.nd), dtype=numpy.float64) | |
self.fixvel = numpy.zeros(self.np, dtype=numpy.float64) | |
self.force = numpy.zeros((self.np, self.nd), dtype=numpy.float64) | |
t@@ -755,7 +755,7 @@ class sim: | |
# Allocate array memory for particles | |
self.x = numpy.empty((self.np, self.nd), dtype=numpy.float64) | |
self.radius = numpy.empty(self.np, dtype=numpy.float64) | |
- self.xysum = numpy.empty((self.np, 2), dtype=numpy.float64) | |
+ self.xyzsum = numpy.empty((self.np, 3), dtype=numpy.float64) | |
self.vel = numpy.empty((self.np, self.nd), dtype=numpy.float64) | |
self.fixvel = numpy.empty(self.np, dtype=numpy.float64) | |
self.es_dot = numpy.empty(self.np, dtype=numpy.float64) | |
t@@ -777,8 +777,8 @@ class sim: | |
self.radius[i] =\ | |
numpy.fromfile(fh, dtype=numpy.float64, count=1) | |
- self.xysum = numpy.fromfile(fh, dtype=numpy.float64,\ | |
- count=self.np*2).reshape(self.np,2) | |
+ self.xyzsum = numpy.fromfile(fh, dtype=numpy.float64,\ | |
+ count=self.np*3).reshape(self.np,3) | |
for i in range(self.np): | |
self.vel[i,:] =\ | |
t@@ -1001,7 +1001,7 @@ class sim: | |
fh.write(self.radius[i].astype(numpy.float64)) | |
if (self.np[0] > 0): | |
- fh.write(self.xysum.astype(numpy.float64)) | |
+ fh.write(self.xyzsum.astype(numpy.float64)) | |
for i in range(self.np): | |
fh.write(self.vel[i,:].astype(numpy.float64)) | |
t@@ -1256,21 +1256,30 @@ class sim: | |
fh.write('\n') | |
fh.write(' </DataArray>\n') | |
- # xysum.x | |
+ # xyzsum.x | |
fh.write(' <DataArray type="Float32" Name="Xdisplacement" ' | |
+ 'format="ascii">\n') | |
fh.write(' ') | |
for i in range(self.np): | |
- fh.write('{} '.format(self.xysum[i,0])) | |
+ fh.write('{} '.format(self.xyzsum[i,0])) | |
fh.write('\n') | |
fh.write(' </DataArray>\n') | |
- # xysum.y | |
+ # xyzsum.y | |
fh.write(' <DataArray type="Float32" Name="Ydisplacement" ' | |
+ 'format="ascii">\n') | |
fh.write(' ') | |
for i in range(self.np): | |
- fh.write('{} '.format(self.xysum[i,1])) | |
+ fh.write('{} '.format(self.xyzsum[i,1])) | |
+ fh.write('\n') | |
+ fh.write(' </DataArray>\n') | |
+ | |
+ # xyzsum.z | |
+ fh.write(' <DataArray type="Float32" Name="Zdisplacement" ' | |
+ + 'format="ascii">\n') | |
+ fh.write(' ') | |
+ for i in range(self.np): | |
+ fh.write('{} '.format(self.xyzsum[i,2])) | |
fh.write('\n') | |
fh.write(' </DataArray>\n') | |
t@@ -2128,8 +2137,8 @@ class sim: | |
.reshape(self.np, self.nd) | |
self.es = numpy.zeros(self.np, dtype=numpy.float64) | |
self.ev = numpy.zeros(self.np, dtype=numpy.float64) | |
- self.xysum = numpy.zeros(self.np*2, dtype=numpy.float64)\ | |
- .reshape(self.np, 2) | |
+ self.xyzsum = numpy.zeros(self.np*3, dtype=numpy.float64)\ | |
+ .reshape(self.np, 3) | |
def adjustUpperWall(self, z_adjust = 1.1): | |
''' | |
t@@ -3520,14 +3529,14 @@ class sim: | |
I = numpy.nonzero((self.x[:,2] > zlower) & (self.x[:,2] < zupper)) | |
# Save mean x displacement | |
- xdisp[iz] = numpy.mean(self.xysum[I,0]) | |
+ xdisp[iz] = numpy.mean(self.xyzsum[I,0]) | |
# Save x displacement standard deviation | |
- err[iz] = numpy.std(self.xysum[I,0]) | |
+ err[iz] = numpy.std(self.xyzsum[I,0]) | |
plt.figure(figsize=[4, 4]) | |
ax = plt.subplot(111) | |
- ax.scatter(self.xysum[:,0], self.x[:,2], c='gray', marker='+') | |
+ ax.scatter(self.xyzsum[:,0], self.x[:,2], c='gray', marker='+') | |
ax.errorbar(xdisp, zpos, xerr=err, | |
c='black', linestyle='-', linewidth=1.4) | |
ax.set_xlabel("Horizontal particle displacement, [m]") | |
diff --git a/src/datatypes.h b/src/datatypes.h | |
t@@ -15,7 +15,7 @@ | |
// Structure containing kinematic particle values | |
struct Kinematics { | |
Float4 *x; // Positions + radii (w) | |
- Float2 *xysum; // Horizontal distance traveled | |
+ Float4 *xyzsum; // Horizontal distance traveled | |
Float4 *vel; // Translational velocities + fixvels (w) | |
Float4 *acc; // Translational accelerations | |
Float4 *force; // Sums of forces | |
diff --git a/src/device.cu b/src/device.cu | |
t@@ -244,7 +244,7 @@ __host__ void DEM::allocateGlobalDeviceMemory(void) | |
// Kinematics arrays | |
cudaMalloc((void**)&dev_x, memSizeF4); | |
- cudaMalloc((void**)&dev_xysum, memSizeF4); | |
+ cudaMalloc((void**)&dev_xyzsum, memSizeF4); | |
cudaMalloc((void**)&dev_vel, memSizeF4); | |
cudaMalloc((void**)&dev_vel0, memSizeF4); | |
cudaMalloc((void**)&dev_acc, memSizeF4); | |
t@@ -312,7 +312,7 @@ __host__ void DEM::freeGlobalDeviceMemory() | |
printf("\nFreeing device memory: "); | |
// Particle arrays | |
cudaFree(dev_x); | |
- cudaFree(dev_xysum); | |
+ cudaFree(dev_xyzsum); | |
cudaFree(dev_vel); | |
cudaFree(dev_vel0); | |
cudaFree(dev_acc); | |
t@@ -379,8 +379,8 @@ __host__ void DEM::transferToGlobalDeviceMemory(int status… | |
// Kinematic particle values | |
cudaMemcpy( dev_x, k.x, | |
memSizeF4, cudaMemcpyHostToDevice); | |
- cudaMemcpy( dev_xysum, k.xysum, | |
- sizeof(Float2)*np, cudaMemcpyHostToDevice); | |
+ cudaMemcpy( dev_xyzsum, k.xyzsum, | |
+ memSizeF4, cudaMemcpyHostToDevice); | |
cudaMemcpy( dev_vel, k.vel, | |
memSizeF4, cudaMemcpyHostToDevice); | |
cudaMemcpy( dev_vel0, k.vel, | |
t@@ -459,8 +459,8 @@ __host__ void DEM::transferFromGlobalDeviceMemory() | |
// Kinematic particle values | |
cudaMemcpy( k.x, dev_x, | |
memSizeF4, cudaMemcpyDeviceToHost); | |
- cudaMemcpy( k.xysum, dev_xysum, | |
- sizeof(Float2)*np, cudaMemcpyDeviceToHost); | |
+ cudaMemcpy( k.xyzsum, dev_xyzsum, | |
+ memSizeF4, cudaMemcpyDeviceToHost); | |
cudaMemcpy( k.vel, dev_vel, | |
memSizeF4, cudaMemcpyDeviceToHost); | |
cudaMemcpy( k.acc, dev_acc, | |
t@@ -1505,7 +1505,7 @@ __host__ void DEM::startTime() | |
dev_angacc, | |
dev_vel0, | |
dev_angvel0, | |
- dev_xysum, | |
+ dev_xyzsum, | |
dev_gridParticleIndex, | |
iter); | |
cudaThreadSynchronize(); | |
diff --git a/src/file_io.cpp b/src/file_io.cpp | |
t@@ -86,7 +86,7 @@ void DEM::readbin(const char *target) | |
cout << " Allocating host memory: "; | |
// Allocate more host arrays | |
k.x = new Float4[np]; | |
- k.xysum = new Float2[np]; | |
+ k.xyzsum = new Float4[np]; | |
k.vel = new Float4[np]; | |
k.force = new Float4[np]; | |
k.angpos = new Float4[np]; | |
t@@ -120,8 +120,9 @@ void DEM::readbin(const char *target) | |
ifs.read(as_bytes(k.x[i].w), sizeof(Float)); | |
} | |
for (i = 0; i<np; ++i) { | |
- ifs.read(as_bytes(k.xysum[i].x), sizeof(Float)); | |
- ifs.read(as_bytes(k.xysum[i].y), sizeof(Float)); | |
+ ifs.read(as_bytes(k.xyzsum[i].x), sizeof(Float)); | |
+ ifs.read(as_bytes(k.xyzsum[i].y), sizeof(Float)); | |
+ ifs.read(as_bytes(k.xyzsum[i].z), sizeof(Float)); | |
} | |
for (i = 0; i<np; ++i) { | |
ifs.read(as_bytes(k.vel[i].x), sizeof(Float)); | |
t@@ -348,8 +349,9 @@ void DEM::writebin(const char *target) | |
ofs.write(as_bytes(k.x[i].w), sizeof(Float)); | |
} | |
for (i = 0; i<np; ++i) { | |
- ofs.write(as_bytes(k.xysum[i].x), sizeof(Float)); | |
- ofs.write(as_bytes(k.xysum[i].y), sizeof(Float)); | |
+ ofs.write(as_bytes(k.xyzsum[i].x), sizeof(Float)); | |
+ ofs.write(as_bytes(k.xyzsum[i].y), sizeof(Float)); | |
+ ofs.write(as_bytes(k.xyzsum[i].z), sizeof(Float)); | |
} | |
for (i = 0; i<np; ++i) { | |
ofs.write(as_bytes(k.vel[i].x), sizeof(Float)); | |
diff --git a/src/integration.cuh b/src/integration.cuh | |
t@@ -17,7 +17,7 @@ __global__ void integrate(Float4* dev_x_sorted, Float4* dev_… | |
Float4* dev_force, Float4* dev_torque, Float4* dev_angpos, // Input | |
Float4* dev_acc, Float4* dev_angacc, | |
Float4* dev_vel0, Float4* dev_angvel0, | |
- Float2* dev_xysum, | |
+ Float4* dev_xyzsum, | |
unsigned int* dev_gridParticleIndex, // Input: Sorted-Unsorted key | |
unsigned int iter) | |
{ | |
t@@ -35,7 +35,7 @@ __global__ void integrate(Float4* dev_x_sorted, Float4* dev_… | |
const Float4 x = dev_x_sorted[idx]; | |
const Float4 vel = dev_vel_sorted[idx]; | |
const Float4 angvel = dev_angvel_sorted[idx]; | |
- Float2 xysum = dev_xysum[orig_idx]; | |
+ Float4 xyzsum = dev_xyzsum[orig_idx]; | |
// Get old accelerations for three-term Taylor expansion. These values | |
// don't exist in the first time step | |
t@@ -139,8 +139,9 @@ __global__ void integrate(Float4* dev_x_sorted, Float4* de… | |
// Add horizontal-displacement for this time step to the sum of | |
// horizontal displacements | |
- xysum.x += vel.x*dt; | |
- xysum.y += vel.y*dt; | |
+ xyzsum.x += vel.x*dt; | |
+ xyzsum.y += vel.y*dt; | |
+ xyzsum.z += vel.z*dt; | |
#endif | |
#ifdef TY2 | |
t@@ -168,8 +169,9 @@ __global__ void integrate(Float4* dev_x_sorted, Float4* de… | |
// Add horizontal-displacement for this time step to the sum of | |
// horizontal displacements | |
- xysum.x += vel.x*dt + 0.5*acc.x*dt*dt; | |
- xysum.y += vel.y*dt + 0.5*acc.y*dt*dt; | |
+ xyzsum.x += vel.x*dt + 0.5*acc.x*dt*dt; | |
+ xyzsum.y += vel.y*dt + 0.5*acc.y*dt*dt; | |
+ xyzsum.z += vel.z*dt + 0.5*acc.z*dt*dt; | |
#endif | |
#ifdef TY3 | |
t@@ -211,8 +213,9 @@ __global__ void integrate(Float4* dev_x_sorted, Float4* de… | |
// Add horizontal-displacement for this time step to the sum of | |
// horizontal displacements | |
- xysum.x += vel.x*dt + 0.5*acc.x*dt*dt + 1.0/6.0*dacc_dt.x*dt*dt*dt; | |
- xysum.y += vel.y*dt + 0.5*acc.y*dt*dt + 1.0/6.0*dacc_dt.y*dt*dt*dt; | |
+ xyzsum.x += vel.x*dt + 0.5*acc.x*dt*dt + 1.0/6.0*dacc_dt.x*dt*dt*dt; | |
+ xyzsum.y += vel.y*dt + 0.5*acc.y*dt*dt + 1.0/6.0*dacc_dt.y*dt*dt*dt; | |
+ xyzsum.z += vel.z*dt + 0.5*acc.z*dt*dt + 1.0/6.0*dacc_dt.z*dt*dt*dt; | |
#endif | |
// Move particles outside the domain across periodic boundaries | |
t@@ -238,7 +241,7 @@ __global__ void integrate(Float4* dev_x_sorted, Float4* de… | |
__syncthreads(); | |
// Store data in global memory at original, pre-sort positions | |
- dev_xysum[orig_idx] = xysum; | |
+ dev_xyzsum[orig_idx] = xyzsum; | |
dev_acc[orig_idx] = acc; | |
dev_angacc[orig_idx] = angacc; | |
dev_angvel[orig_idx] = angvel_new; | |
diff --git a/src/raytracer.cuh b/src/raytracer.cuh | |
t@@ -476,11 +476,11 @@ __host__ void DEM::render( | |
unit = "rad/s"; | |
} else if (method == 4) { // Visualize xdisp | |
- // Convert xysum to xsum | |
+ // Convert xyzsum to xsum | |
linarr = new Float[np]; | |
#pragma omp parallel for if(np>100) | |
for (i = 0; i<np; ++i) { | |
- linarr[i] = k.xysum[i].x; | |
+ linarr[i] = k.xyzsum[i].x; | |
} | |
transfer = 1; | |
desc = "X-axis displacement"; | |
diff --git a/src/sphere.cpp b/src/sphere.cpp | |
t@@ -82,7 +82,7 @@ DEM::~DEM(void) | |
if (verbose == 1) | |
std::cout << "Freeing host memory: "; | |
delete[] k.x; | |
- delete[] k.xysum; | |
+ delete[] k.xyzsum; | |
delete[] k.vel; | |
delete[] k.force; | |
delete[] k.angpos; | |
diff --git a/src/sphere.h b/src/sphere.h | |
t@@ -58,7 +58,7 @@ class DEM { | |
// Particle kinematics arrays | |
Float4 *dev_x; | |
- Float2 *dev_xysum; | |
+ Float4 *dev_xyzsum; | |
Float4 *dev_vel; | |
Float4 *dev_vel0; | |
Float4 *dev_acc; | |
diff --git a/tests/io_tests_fluid.py b/tests/io_tests_fluid.py | |
t@@ -42,7 +42,7 @@ if numpy.allclose(orig.x, cuda.x, 0.01): | |
cuda.x = orig.x # ignore small changes | |
if numpy.max(numpy.abs(cuda.vel - orig.vel)) < 1.0e-5: | |
cuda.vel = orig.vel # ignore small changes | |
- cuda.xysum = orig.xysum | |
+ cuda.xyzsum = orig.xyzsum | |
cuda.force = orig.force | |
compare(orig, cuda, "CUDA IO: ") | |