tAdded ev and ev_dot arrays for energy lost to viscous damping. Values not written - sphere - GPU-based 3D discrete element method algorithm with optional fluid coupling
(HTM) git clone git://src.adamsgaard.dk/sphere
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
(DIR) commit 9f513367e78001788203a6e76e56f94894e2a3af
(DIR) parent 507464b4e3dc102985c1490e416d5ee0da7925f0
(HTM) Author: Anders Damsgaard <adc@geo.au.dk>
Date: Thu, 4 Oct 2012 20:44:51 +0200
Added ev and ev_dot arrays for energy lost to viscous damping. Values not written
Diffstat:
M python/sphere.py | 8 ++++++++
M src/contactmodels.cuh | 8 +++++---
M src/contactsearch.cuh | 40 +++++++++++++++++++-------------
M src/datatypes.h | 2 ++
M src/device.cu | 13 ++++++++++++-
M src/file_io.cpp | 6 ++++++
M src/main.cpp | 8 ++++++++
7 files changed, 65 insertions(+), 20 deletions(-)
---
(DIR) diff --git a/python/sphere.py b/python/sphere.py
t@@ -53,7 +53,9 @@ class Spherebin:
self.mu_d = numpy.zeros(self.np, dtype=numpy.float64)
self.mu_r = numpy.zeros(self.np, dtype=numpy.float64)
self.es_dot = numpy.zeros(self.np, dtype=numpy.float64)
+ self.ev_dot = numpy.zeros(self.np, dtype=numpy.float64)
self.es = numpy.zeros(self.np, dtype=numpy.float64)
+ self.ev = numpy.zeros(self.np, dtype=numpy.float64)
self.p = numpy.zeros(self.np, dtype=numpy.float64)
self.bonds = numpy.ones(self.np*4, dtype=numpy.uint32).reshape(self.np,4) * self.np
t@@ -124,7 +126,9 @@ class Spherebin:
self.mu_d = numpy.zeros(self.np, dtype=numpy.float64)
self.mu_r = numpy.zeros(self.np, dtype=numpy.float64)
self.es_dot = numpy.zeros(self.np, dtype=numpy.float64)
+ self.ev_dot = numpy.zeros(self.np, dtype=numpy.float64)
self.es = numpy.zeros(self.np, dtype=numpy.float64)
+ self.ev = numpy.zeros(self.np, dtype=numpy.float64)
self.p = numpy.zeros(self.np, dtype=numpy.float64)
self.bonds = numpy.zeros(self.np, dtype=numpy.float64)
t@@ -158,7 +162,9 @@ class Spherebin:
self.mu_d[j] = numpy.fromfile(fh, dtype=numpy.float64, count=1)
self.mu_r[j] = numpy.fromfile(fh, dtype=numpy.float64, count=1)
self.es_dot[j] = numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ self.ev_dot[j] = numpy.fromfile(fh, dtype=numpy.float64, count=1)
self.es[j] = numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ self.ev[j] = numpy.fromfile(fh, dtype=numpy.float64, count=1)
self.p[j] = numpy.fromfile(fh, dtype=numpy.float64, count=1)
# Constant, single-value physical parameters
t@@ -258,7 +264,9 @@ class Spherebin:
fh.write(self.mu_d[j].astype(numpy.float64))
fh.write(self.mu_r[j].astype(numpy.float64))
fh.write(self.es_dot[j].astype(numpy.float64))
+ fh.write(self.ev_dot[j].astype(numpy.float64))
fh.write(self.es[j].astype(numpy.float64))
+ fh.write(self.ev[j].astype(numpy.float64))
fh.write(self.p[j].astype(numpy.float64))
# Constant, single-value physical parameters
(DIR) diff --git a/src/contactmodels.cuh b/src/contactmodels.cuh
t@@ -7,7 +7,8 @@
// Linear viscoelastic contact model for particle-wall interactions
// with tangential friction and rolling resistance
-__device__ Float contactLinear_wall(Float3* F, Float3* T, Float* es_dot, Float* p,
+__device__ Float contactLinear_wall(Float3* F, Float3* T, Float* es_dot,
+ Float* ev_dot, Float* p,
unsigned int idx_a, Float radius_a,
Float4* dev_vel_sorted, Float4* dev_angvel_sorted,
Float3 n, Float delta, Float wvel)
t@@ -103,7 +104,8 @@ __device__ Float contactLinear_wall(Float3* F, Float3* T, Float* es_dot, Float*
// Linear vicoelastic contact model for particle-particle interactions
// with tangential friction and rolling resistance
-__device__ void contactLinearViscous(Float3* F, Float3* T, Float* es_dot, Float* p,
+__device__ void contactLinearViscous(Float3* F, Float3* T,
+ Float* es_dot, Float* ev_dot, Float* p,
unsigned int idx_a, unsigned int idx_b,
Float4* dev_vel_sorted,
Float4* dev_angvel_sorted,
t@@ -246,7 +248,7 @@ __device__ void contactLinearViscous(Float3* F, Float3* T, Float* es_dot, Float*
// Linear elastic contact model for particle-particle interactions
__device__ void contactLinear(Float3* F, Float3* T,
- Float* es_dot, Float* p,
+ Float* es_dot, Float* ev_dot, Float* p,
unsigned int idx_a_orig,
unsigned int idx_b_orig,
Float4 vel_a,
(DIR) diff --git a/src/contactsearch.cuh b/src/contactsearch.cuh
t@@ -84,7 +84,8 @@ __device__ void findAndProcessContactsInCell(int3 targetCell,
unsigned int idx_a,
Float4 x_a, Float radius_a,
Float3* F, Float3* T,
- Float* es_dot, Float* p,
+ Float* es_dot, Float* ev_dot,
+ Float* p,
Float4* dev_x_sorted,
Float* dev_radius_sorted,
Float4* dev_vel_sorted,
t@@ -142,7 +143,7 @@ __device__ void findAndProcessContactsInCell(int3 targetCell,
// Check for particle overlap
if (delta_ab < 0.0f) {
- contactLinearViscous(F, T, es_dot, p,
+ contactLinearViscous(F, T, es_dot, ev_dot, p,
idx_a, idx_b,
dev_vel_sorted,
dev_angvel_sorted,
t@@ -158,7 +159,7 @@ __device__ void findAndProcessContactsInCell(int3 targetCell,
// Check wether particles are bonded together
/*if (bonds.x == idx_b || bonds.y == idx_b ||
bonds.z == idx_b || bonds.w == idx_b) {
- bondLinear(F, T, es_dot, p,
+ bondLinear(F, T, es_dot, p, % ev_dot missing
idx_a, idx_b,
dev_x_sorted, dev_vel_sorted,
dev_angvel_sorted,
t@@ -282,7 +283,7 @@ __device__ void findContactsInCell(int3 targetCell,
// Check wether particles are bonded together
/*if (bonds.x == idx_b || bonds.y == idx_b ||
bonds.z == idx_b || bonds.w == idx_b) {
- bondLinear(F, T, es_dot, p,
+ bondLinear(F, T, es_dot, p, % ev_dot missing
idx_a, idx_b,
dev_x_sorted, dev_vel_sorted,
dev_angvel_sorted,
t@@ -364,7 +365,11 @@ __global__ void interact(unsigned int* dev_gridParticleIndex, // Input: Unsorted
Float4* dev_vel_sorted, Float4* dev_angvel_sorted,
Float4* dev_vel, Float4* dev_angvel,
Float4* dev_force, Float4* dev_torque,
- Float* dev_es_dot, Float* dev_es, Float* dev_p,
+ Float* dev_es_dot,
+ Float* dev_ev_dot,
+ Float* dev_es,
+ Float* dev_ev,
+ Float* dev_p,
Float4* dev_w_nx, Float4* dev_w_mvfd,
Float* dev_w_force, //uint4* dev_bonds_sorted,
unsigned int* dev_contacts,
t@@ -400,6 +405,7 @@ __global__ void interact(unsigned int* dev_gridParticleIndex, // Input: Unsorted
// Initiate shear friction loss rate at 0.0
Float es_dot = 0.0f;
+ Float ev_dot = 0.0f;
// Initiate pressure on particle at 0.0
Float p = 0.0f;
t@@ -442,14 +448,14 @@ __global__ void interact(unsigned int* dev_gridParticleIndex, // Input: Unsorted
// Process collision if the particles are overlapping
if (delta_n < 0.0f) {
- /*contactLinearViscous(&F, &T, &es_dot, &p,
+ /*contactLinearViscous(&F, &T, &es_dot, &ev_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);*/
- contactLinear(&F, &T, &es_dot, &p,
+ contactLinear(&F, &T, &es_dot, &ev_dot, &p,
idx_a_orig,
idx_b_orig,
vel_a,
t@@ -495,7 +501,7 @@ __global__ void interact(unsigned int* dev_gridParticleIndex, // Input: Unsorted
for (int x_dim=-1; x_dim<2; ++x_dim) { // x-axis
targetPos = gridPos + make_int3(x_dim, y_dim, z_dim);
findAndProcessContactsInCell(targetPos, idx_a, x_a, radius_a,
- &F, &T, &es_dot, &p,
+ &F, &T, &es_dot, &ev_dot, &p,
dev_x_sorted, dev_radius_sorted,
dev_vel_sorted, dev_angvel_sorted,
dev_cellStart, dev_cellEnd,
t@@ -515,7 +521,7 @@ __global__ void interact(unsigned int* dev_gridParticleIndex, // Input: Unsorted
delta_w = w_up_nx.w - (x_a.z + radius_a);
w_n = MAKE_FLOAT3(0.0f, 0.0f, -1.0f);
if (delta_w < 0.0f) {
- w_force = contactLinear_wall(&F, &T, &es_dot, &p, idx_a, radius_a,
+ w_force = contactLinear_wall(&F, &T, &es_dot, &ev_dot, &p, idx_a, radius_a,
dev_vel_sorted, dev_angvel_sorted,
w_n, delta_w, w_up_mvfd.y);
}
t@@ -524,7 +530,7 @@ __global__ void interact(unsigned int* dev_gridParticleIndex, // Input: Unsorted
delta_w = x_a.z - radius_a - origo.z;
w_n = MAKE_FLOAT3(0.0f, 0.0f, 1.0f);
if (delta_w < 0.0f) {
- (void)contactLinear_wall(&F, &T, &es_dot, &p, idx_a, radius_a,
+ (void)contactLinear_wall(&F, &T, &es_dot, &ev_dot, &p, idx_a, radius_a,
dev_vel_sorted, dev_angvel_sorted,
w_n, delta_w, 0.0f);
}
t@@ -536,7 +542,7 @@ __global__ void interact(unsigned int* dev_gridParticleIndex, // Input: Unsorted
delta_w = x_a.x - radius_a - origo.x;
w_n = MAKE_FLOAT3(1.0f, 0.0f, 0.0f);
if (delta_w < 0.0f) {
- (void)contactLinear_wall(&F, &T, &es_dot, &p, idx_a, radius_a,
+ (void)contactLinear_wall(&F, &T, &es_dot, &ev_dot, &p, idx_a, radius_a,
dev_vel_sorted, dev_angvel_sorted,
w_n, delta_w, 0.0f);
}
t@@ -545,7 +551,7 @@ __global__ void interact(unsigned int* dev_gridParticleIndex, // Input: Unsorted
delta_w = L.x - (x_a.x + radius_a);
w_n = MAKE_FLOAT3(-1.0f, 0.0f, 0.0f);
if (delta_w < 0.0f) {
- (void)contactLinear_wall(&F, &T, &es_dot, &p, idx_a, radius_a,
+ (void)contactLinear_wall(&F, &T, &es_dot, &ev_dot, &p, idx_a, radius_a,
dev_vel_sorted, dev_angvel_sorted,
w_n, delta_w, 0.0f);
}
t@@ -554,7 +560,7 @@ __global__ void interact(unsigned int* dev_gridParticleIndex, // Input: Unsorted
delta_w = x_a.y - radius_a - origo.y;
w_n = MAKE_FLOAT3(0.0f, 1.0f, 0.0f);
if (delta_w < 0.0f) {
- (void)contactLinear_wall(&F, &T, &es_dot, &p, idx_a, radius_a,
+ (void)contactLinear_wall(&F, &T, &es_dot, &ev_dot, &p, idx_a, radius_a,
dev_vel_sorted, dev_angvel_sorted,
w_n, delta_w, 0.0f);
}
t@@ -563,7 +569,7 @@ __global__ void interact(unsigned int* dev_gridParticleIndex, // Input: Unsorted
delta_w = L.y - (x_a.y + radius_a);
w_n = MAKE_FLOAT3(0.0f, -1.0f, 0.0f);
if (delta_w < 0.0f) {
- (void)contactLinear_wall(&F, &T, &es_dot, &p, idx_a, radius_a,
+ (void)contactLinear_wall(&F, &T, &es_dot, &ev_dot, &p, idx_a, radius_a,
dev_vel_sorted, dev_angvel_sorted,
w_n, delta_w, 0.0f);
}
t@@ -574,7 +580,7 @@ __global__ void interact(unsigned int* dev_gridParticleIndex, // Input: Unsorted
delta_w = x_a.y - radius_a - origo.y;
w_n = MAKE_FLOAT3(0.0f, 1.0f, 0.0f);
if (delta_w < 0.0f) {
- (void)contactLinear_wall(&F, &T, &es_dot, &p, idx_a, radius_a,
+ (void)contactLinear_wall(&F, &T, &es_dot, &ev_dot, &p, idx_a, radius_a,
dev_vel_sorted, dev_angvel_sorted,
w_n, delta_w, 0.0f);
}
t@@ -583,7 +589,7 @@ __global__ void interact(unsigned int* dev_gridParticleIndex, // Input: Unsorted
delta_w = L.y - (x_a.y + radius_a);
w_n = MAKE_FLOAT3(0.0f, -1.0f, 0.0f);
if (delta_w < 0.0f) {
- (void)contactLinear_wall(&F, &T, &es_dot, &p, idx_a, radius_a,
+ (void)contactLinear_wall(&F, &T, &es_dot, &ev_dot, &p, idx_a, radius_a,
dev_vel_sorted, dev_angvel_sorted,
w_n, delta_w, 0.0f);
}
t@@ -598,7 +604,9 @@ __global__ void interact(unsigned int* dev_gridParticleIndex, // Input: Unsorted
dev_force[orig_idx] = MAKE_FLOAT4(F.x, F.y, F.z, 0.0f);
dev_torque[orig_idx] = MAKE_FLOAT4(T.x, T.y, T.z, 0.0f);
dev_es_dot[orig_idx] = es_dot;
+ dev_ev_dot[orig_idx] = ev_dot;
dev_es[orig_idx] += es_dot * devC_dt;
+ dev_ev[orig_idx] += ev_dot * devC_dt;
dev_p[orig_idx] = p;
dev_w_force[orig_idx] = w_force;
}
(DIR) diff --git a/src/datatypes.h b/src/datatypes.h
t@@ -82,7 +82,9 @@ struct Particles {
Float *mu_r;
Float *rho;
Float *es_dot;
+ Float *ev_dot;
Float *es;
+ Float *ev;
Float *p;
Float *m;
Float *I;
(DIR) diff --git a/src/device.cu b/src/device.cu
t@@ -187,7 +187,9 @@ __host__ void gpuMain(Float4* host_x,
Float4* dev_torque; // Sum of torques
Float* dev_radius; // Particle radius
Float* dev_es_dot; // Current shear energy producion rate
+ Float* dev_ev_dot; // Current viscous energy producion rate
Float* dev_es; // Total shear energy excerted on particle
+ Float* dev_ev; // Total viscous energy excerted on particle
Float* dev_p; // Pressure excerted onto particle
//uint4* dev_bonds; // Particle bond pairs
t@@ -241,7 +243,9 @@ __host__ void gpuMain(Float4* host_x,
cudaMalloc((void**)&dev_radius, memSizeF);
cudaMalloc((void**)&dev_radius_sorted, memSizeF);
cudaMalloc((void**)&dev_es_dot, memSizeF);
+ cudaMalloc((void**)&dev_ev_dot, memSizeF);
cudaMalloc((void**)&dev_es, memSizeF);
+ cudaMalloc((void**)&dev_ev, memSizeF);
cudaMalloc((void**)&dev_p, memSizeF);
//cudaMalloc((void**)&dev_bonds, sizeof(uint4) * p->np);
//cudaMalloc((void**)&dev_bonds_sorted, sizeof(uint4) * p->np);
t@@ -279,7 +283,9 @@ __host__ void gpuMain(Float4* host_x,
//cudaMemcpy(dev_bonds, host_bonds, sizeof(uint4) * p->np, cudaMemcpyHostToDevice);
cudaMemcpy(dev_radius, p->radius, memSizeF, cudaMemcpyHostToDevice);
cudaMemcpy(dev_es_dot, p->es_dot, memSizeF, cudaMemcpyHostToDevice);
+ cudaMemcpy(dev_ev_dot, p->ev_dot, memSizeF, cudaMemcpyHostToDevice);
cudaMemcpy(dev_es, p->es, memSizeF, cudaMemcpyHostToDevice);
+ cudaMemcpy(dev_ev, p->ev, memSizeF, cudaMemcpyHostToDevice);
cudaMemcpy(dev_p, p->p, memSizeF, cudaMemcpyHostToDevice);
// Wall data (wall mass and number in constant memory)
t@@ -509,7 +515,8 @@ __host__ void gpuMain(Float4* host_x,
dev_vel_sorted, dev_angvel_sorted,
dev_vel, dev_angvel,
dev_force, dev_torque,
- dev_es_dot, dev_es, dev_p,
+ dev_es_dot, dev_ev_dot,
+ dev_es, dev_ev, dev_p,
dev_w_nx, dev_w_mvfd, dev_w_force,
//dev_bonds_sorted,
dev_contacts,
t@@ -593,7 +600,9 @@ __host__ void gpuMain(Float4* host_x,
cudaMemcpy(host_torque, dev_torque, memSizeF4, cudaMemcpyDeviceToHost);
//cudaMemcpy(host_bonds, dev_bonds, sizeof(uint4) * p->np, cudaMemcpyDeviceToHost);
cudaMemcpy(p->es_dot, dev_es_dot, memSizeF, cudaMemcpyDeviceToHost);
+ cudaMemcpy(p->ev_dot, dev_ev_dot, memSizeF, cudaMemcpyDeviceToHost);
cudaMemcpy(p->es, dev_es, memSizeF, cudaMemcpyDeviceToHost);
+ cudaMemcpy(p->ev, dev_ev, memSizeF, cudaMemcpyDeviceToHost);
cudaMemcpy(p->p, dev_p, memSizeF, cudaMemcpyDeviceToHost);
// Wall data
t@@ -689,7 +698,9 @@ __host__ void gpuMain(Float4* host_x,
cudaFree(dev_radius);
cudaFree(dev_radius_sorted);
cudaFree(dev_es_dot);
+ cudaFree(dev_ev_dot);
cudaFree(dev_es);
+ cudaFree(dev_ev);
cudaFree(dev_p);
//cudaFree(dev_bonds);
//cudaFree(dev_bonds_sorted);
(DIR) diff --git a/src/file_io.cpp b/src/file_io.cpp
t@@ -105,7 +105,9 @@ int fwritebin(char *target,
fwrite(&p->mu_d[j], sizeof(p->mu_d[j]), 1, fp);
fwrite(&p->mu_r[j], sizeof(p->mu_r[j]), 1, fp);
fwrite(&p->es_dot[j], sizeof(p->es_dot[j]), 1, fp);
+ fwrite(&p->ev_dot[j], sizeof(p->ev_dot[j]), 1, fp);
fwrite(&p->es[j], sizeof(p->es[j]), 1, fp);
+ fwrite(&p->ev[j], sizeof(p->ev[j]), 1, fp);
fwrite(&p->p[j], sizeof(p->p[j]), 1, fp);
}
t@@ -254,8 +256,12 @@ int fwritebin(char *target,
fwrite(&d, sizeof(d), 1, fp);
d = (double)p->es_dot[j];
fwrite(&d, sizeof(d), 1, fp);
+ d = (double)p->ev_dot[j];
+ fwrite(&d, sizeof(d), 1, fp);
d = (double)p->es[j];
fwrite(&d, sizeof(d), 1, fp);
+ d = (double)p->ev[j];
+ fwrite(&d, sizeof(d), 1, fp);
d = (double)p->p[j];
fwrite(&d, sizeof(d), 1, fp);
}
(DIR) diff --git a/src/main.cpp b/src/main.cpp
t@@ -176,7 +176,9 @@ int main(int argc, char *argv[])
p.mu_d = new Float[p.np]; // Inter-particle dynamic shear contact friction coefficients
p.mu_r = new Float[p.np]; // Inter-particle rolling contact friction coefficients
p.es_dot = new Float[p.np]; // Rate of shear energy dissipation
+ p.ev_dot = new Float[p.np]; // Rate of viscous energy dissipation
p.es = new Float[p.np]; // Total shear energy dissipation
+ p.ev = new Float[p.np]; // Total viscous energy dissipation
p.p = new Float[p.np]; // Pressure excerted onto particle
params.wmode = new int[MAXWALLS]; // Wall BC's, 0: devs, 1: vel
t@@ -274,8 +276,12 @@ int main(int argc, char *argv[])
exit(1);
if (fread(&p.es_dot[j], sizeof(p.es_dot[j]), 1, fp) != 1)
exit(1);
+ if (fread(&p.ev_dot[j], sizeof(p.ev_dot[j]), 1, fp) != 1)
+ exit(1);
if (fread(&p.es[j], sizeof(p.es[j]), 1, fp) != 1)
exit(1);
+ if (fread(&p.ev[j], sizeof(p.ev[j]), 1, fp) != 1)
+ exit(1);
if (fread(&p.p[j], sizeof(p.p[j]), 1, fp) != 1)
exit(1);
}
t@@ -472,7 +478,9 @@ int main(int argc, char *argv[])
delete[] p.mu_r;
delete[] p.rho;
delete[] p.es_dot;
+ delete[] p.ev_dot;
delete[] p.es;
+ delete[] p.ev;
delete[] p.p;
// Wall arrays