tadded triaxial simulation possibilities - 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 a715612236ad2daa95f223638876f0b63b387605
(DIR) parent 31a7e78397358a6c7b91b5f037d61a71b928bc90
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Wed, 21 Aug 2013 14:46:40 +0200
added triaxial simulation possibilities
Diffstat:
M python/sphere.py | 65 +++++++++++++++++++++++++++----
M src/cohesion.cuh | 1 +
M src/contactsearch.cuh | 151 ++++++++++++++++++++-----------
M src/device.cu | 3 ++-
4 files changed, 156 insertions(+), 64 deletions(-)
---
(DIR) diff --git a/python/sphere.py b/python/sphere.py
t@@ -88,12 +88,25 @@ class Spherebin:
self.V_b = numpy.zeros(1, dtype=numpy.float64)
# Wall data
+ # nw: Number of dynamic walls
+ # nw = 1: Uniaxial
+ # nw = 2: Biaxial
+ # nw = 5: Triaxial
self.nw = numpy.ones(1, dtype=numpy.uint32) * nw
self.wmode = numpy.zeros(self.nw, dtype=numpy.int32)
self.w_n = numpy.zeros((self.nw, self.nd), dtype=numpy.float64)
- if (self.nw > 0):
+ if (self.nw >= 1):
self.w_n[0,2] = -1.0
+ if (self.nw >= 2):
+ self.w_n[1,0] = -1.0
+ if (self.nw >= 3):
+ self.w_n[2,0] = 1.0
+ if (self.nw >= 4):
+ self.w_n[3,1] = -1.0
+ if (self.nw >= 5):
+ self.w_n[4,1] = 1.0
+
self.w_x = numpy.ones(self.nw, dtype=numpy.float64)
self.w_m = numpy.zeros(self.nw, dtype=numpy.float64)
self.w_vel = numpy.zeros(self.nw, dtype=numpy.float64)
t@@ -875,21 +888,38 @@ class Spherebin:
.reshape(self.np, 2)
def adjustUpperWall(self, z_adjust = 1.1):
- 'Adjust grid and dynamic upper wall to max. particle height'
+ 'Included for legacy purposes, calls adjustWall with idx=0'
+ self.adjustWall(idx=0, adjust = z_adjust)
+
+ def adjustWall(self, idx, adjust):
+ 'Adjust grid and dynamic wall to max. particle position'
+
+ if (idx == 0):
+ dim = 2
+ elif (idx == 1 or idx == 2):
+ dim = 0
+ elif (idx == 3 or idx == 4):
+ dim = 1
+ else:
+ print("adjustWall: idx value not understood")
+
+ xmin = numpy.min(self.x[:,dim] - self.radius)
+ xmax = numpy.max(self.x[:,dim] + self.radius)
- # Compute new grid, scaled to fit max. and min. particle positions
- z_min = numpy.min(self.x[:,2] - self.radius)
- z_max = numpy.max(self.x[:,2] + self.radius)
cellsize = self.L[0] / self.num[0]
- self.num[2] = numpy.ceil(((z_max-z_min)*z_adjust + z_min)/cellsize)
- self.L[2] = (z_max-z_min)*z_adjust + z_min
+
+ self.num[dim] = numpy.ceil(((xmax-xmin)*adjust + xmin)/cellsize)
+ self.L[dim] = (xmax-xmin)*adjust + xmin
# Initialize upper wall
self.nw = numpy.ones(1)
self.wmode = numpy.zeros(1) # fixed BC
self.w_n = numpy.zeros(self.nw*self.nd, dtype=numpy.float64).reshape(self.nw,self.nd)
self.w_n[0,2] = -1.0
- self.w_x = numpy.array([z_max])
+ if (idx == 0 or idx == 1 or idx == 3):
+ self.w_x = numpy.array([xmax])
+ else:
+ self.w_x = numpy.array([xmin])
self.w_m = numpy.array([self.rho[0] * self.np * math.pi * (cellsize/2.0)**3])
self.w_vel = numpy.zeros(1)
self.w_force = numpy.zeros(1)
t@@ -926,6 +956,25 @@ class Spherebin:
self.wmode = numpy.array([2]) # strain rate BC
self.w_vel = numpy.array([wvel])
+ def triaxial(self, wvel = -0.001, deviatoric_stress = 10.0e3):
+ """ Setup triaxial experiment. The upper wall is moved at a fixed
+ velocity in m/s, default values is -0.001 m/s (i.e. downwards).
+ The side walls are exerting a deviatoric stress
+ """
+
+ # zero kinematics
+ self.zeroKinematics()
+
+ # Initialize walls
+ self.nw[0] = 5 # five dynamic walls
+ for i in range(5):
+ self.adjustWall(i)
+ self.w_m[:] = numpy.array([self.rho[0] * self.np * math.pi * (cellsize/2.0)**3])
+ self.wmode = numpy.array([2,1,1,1,1]) # define BCs (vel, stress, stress, ...)
+ self.w_vel = numpy.array([1,0,0,0,0]) * wvel
+ self.w_devs = numpy.array([0,1,1,1,1]) * deviatoric_stress
+
+
def shear(self,
shear_strain_rate = 1,
periodic = 1):
(DIR) diff --git a/src/cohesion.cuh b/src/cohesion.cuh
t@@ -219,6 +219,7 @@ __global__ void bondsLinear(
dev_bonds_omega[idx] = MAKE_FLOAT4(omega_t.x, omega_t.y, omega_t.z, omega_n);
// Save forces and torques to the particle pairs
+ // !!! This is probably wrong, see Obermayer et al. 2013, C & GT (49)
dev_force[bond.x] += MAKE_FLOAT4(f.x, f.y, f.z, 0.0);
dev_force[bond.y] -= MAKE_FLOAT4(f.x, f.y, f.z, 0.0);
//dev_torque[bond.x] += MAKE_FLOAT4(t.x, t.y, t.z, 0.0);
(DIR) diff --git a/src/contactsearch.cuh b/src/contactsearch.cuh
t@@ -409,14 +409,43 @@ __global__ void interact(unsigned int* dev_gridParticleIndex, // Input: Unsorted
devC_grid.L[2]);
// Fetch wall data in global read
- Float4 w_up_nx;
- Float4 w_up_mvfd;
- if (devC_nw > 0) {
- w_up_nx = dev_walls_nx[0];
- w_up_mvfd = dev_walls_mvfd[0];
- } else {
- w_up_nx = MAKE_FLOAT4(0.0f, 0.0f, -1.0f, L.z);
- w_up_mvfd = MAKE_FLOAT4(0.0f, 0.0f, 0.0f, 0.0f);
+ Float4 w_0_nx, w_1_nx, w_2_nx, w_3_nx, w_4_nx;
+ Float4 w_0_mvfd, w_1_mvfd, w_2_mvfd, w_3_mvfd, w_4_mvfd;
+
+ // default wall normals and positions
+ w_0_nx = MAKE_FLOAT4( 0.0f, 0.0f,-1.0f, L.z);
+ w_1_nx = MAKE_FLOAT4(-1.0f, 0.0f, 0.0f, L.x);
+ w_2_nx = MAKE_FLOAT4( 1.0f, 0.0f, 0.0f, 0.0f);
+ w_3_nx = MAKE_FLOAT4( 0.0f,-1.0f, 0.0f, L.y);
+ w_4_nx = MAKE_FLOAT4( 0.0f, 1.0f, 0.0f, 0.0f);
+
+ // default wall mass, vel, force, and devs
+ w_0_mvfd = MAKE_FLOAT4(0.0f, 0.0f, 0.0f, 0.0f);
+ w_1_mvfd = MAKE_FLOAT4(0.0f, 0.0f, 0.0f, 0.0f);
+ w_2_mvfd = MAKE_FLOAT4(0.0f, 0.0f, 0.0f, 0.0f);
+ w_3_mvfd = MAKE_FLOAT4(0.0f, 0.0f, 0.0f, 0.0f);
+ w_4_mvfd = MAKE_FLOAT4(0.0f, 0.0f, 0.0f, 0.0f);
+
+ // fetch data for dynamic walls
+ if (devC_nw >= 1) {
+ w_0_nx = dev_walls_nx[0];
+ w_0_mvfd = dev_walls_mvfd[0];
+ if (devC_nw >= 2) {
+ w_1_nx = dev_walls_nx[1];
+ w_1_mvfd = dev_walls_mvfd[1];
+ }
+ if (devC_nw >= 3) {
+ w_2_nx = dev_walls_nx[2];
+ w_2_mvfd = dev_walls_mvfd[2];
+ }
+ if (devC_nw >= 4) {
+ w_3_nx = dev_walls_nx[3];
+ w_3_mvfd = dev_walls_mvfd[3];
+ }
+ if (devC_nw >= 5) {
+ w_4_nx = dev_walls_nx[4];
+ w_4_mvfd = dev_walls_mvfd[4];
+ }
}
// Index of particle which is bonded to particle A.
t@@ -549,83 +578,87 @@ __global__ void interact(unsigned int* dev_gridParticleIndex, // Input: Unsorted
//// Interact with walls
Float delta_w; // Overlap distance
Float3 w_n; // Wall surface normal
- Float w_force = 0.0; // Force on wall from particle A
+ Float w_0_force = 0.0; // Force on wall 0 from particle A
+ Float w_1_force = 0.0; // Force on wall 1 from particle A
+ Float w_2_force = 0.0; // Force on wall 2 from particle A
+ Float w_3_force = 0.0; // Force on wall 3 from particle A
+ Float w_4_force = 0.0; // Force on wall 4 from particle A
// Upper wall (idx 0)
- delta_w = w_up_nx.w - (x_a.z + radius_a);
- w_n = MAKE_FLOAT3(0.0f, 0.0f, -1.0f);
+ delta_w = w_0_nx.w - (x_a.z + radius_a);
+ w_n = MAKE_FLOAT3(w_0_nx.x, w_0_nx.y, w_0_nx.z);
if (delta_w < 0.0f) {
- 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);
+ w_0_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_0_mvfd.y);
}
// Lower wall (force on wall not stored)
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, &ev_dot, &p, idx_a, radius_a,
- dev_vel_sorted, dev_angvel_sorted,
+ (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);
}
- if (devC_grid.periodic == 0) {
+ if (devC_grid.periodic == 0) { // no periodic walls
- // Left wall
- delta_w = x_a.x - radius_a - origo.x;
- w_n = MAKE_FLOAT3(1.0f, 0.0f, 0.0f);
+ // Right wall (idx 1)
+ delta_w = w_1_nx.w - (x_a.x + radius_a);
+ w_n = MAKE_FLOAT3(w_1_nx.x, w_1_nx.y, w_1_nx.z);
if (delta_w < 0.0f) {
- (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);
+ w_1_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_1_mvfd.y);
}
- // Right wall
- delta_w = L.x - (x_a.x + radius_a);
- w_n = MAKE_FLOAT3(-1.0f, 0.0f, 0.0f);
+ // Left wall (idx 2)
+ delta_w = x_a.x - radius_a - w_2_nx.w;
+ w_n = MAKE_FLOAT3(w_2_nx.x, w_2_nx.y, w_2_nx.z);
if (delta_w < 0.0f) {
- (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);
+ w_2_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_2_mvfd.y);
}
- // Front wall
- delta_w = x_a.y - radius_a - origo.y;
- w_n = MAKE_FLOAT3(0.0f, 1.0f, 0.0f);
+ // Back wall (idx 3)
+ delta_w = w_3_nx.w - (x_a.y + radius_a);
+ w_n = MAKE_FLOAT3(w_3_nx.x, w_3_nx.y, w_3_nx.z);
if (delta_w < 0.0f) {
- (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);
+ w_3_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_3_mvfd.y);
}
- // Back wall
- delta_w = L.y - (x_a.y + radius_a);
- w_n = MAKE_FLOAT3(0.0f, -1.0f, 0.0f);
+ // Front wall (idx 4)
+ delta_w = x_a.y - radius_a - w_4_nx.w;
+ w_n = MAKE_FLOAT3(w_4_nx.x, w_4_nx.y, w_4_nx.z);
if (delta_w < 0.0f) {
- (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);
+ w_4_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_4_mvfd.y);
}
- } else if (devC_grid.periodic == 2) {
+ } else if (devC_grid.periodic == 2) { // right and left walls periodic
- // Front wall
- delta_w = x_a.y - radius_a - origo.y;
- w_n = MAKE_FLOAT3(0.0f, 1.0f, 0.0f);
+ // Back wall (idx 3)
+ delta_w = w_3_nx.w - (x_a.y + radius_a);
+ w_n = MAKE_FLOAT3(w_3_nx.x, w_3_nx.y, w_3_nx.z);
if (delta_w < 0.0f) {
- (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);
+ w_3_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_3_mvfd.y);
}
- // Back wall
- delta_w = L.y - (x_a.y + radius_a);
- w_n = MAKE_FLOAT3(0.0f, -1.0f, 0.0f);
+ // Front wall (idx 4)
+ delta_w = x_a.y - radius_a - w_4_nx.w;
+ w_n = MAKE_FLOAT3(w_4_nx.x, w_4_nx.y, w_4_nx.z);
if (delta_w < 0.0f) {
- (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);
+ w_4_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_4_mvfd.y);
}
}
t@@ -643,7 +676,15 @@ __global__ void interact(unsigned int* dev_gridParticleIndex, // Input: Unsorted
dev_ev[orig_idx] += ev_dot * devC_dt;
dev_p[orig_idx] = p;
if (devC_nw > 0)
- dev_walls_force_pp[orig_idx] = w_force;
+ dev_walls_force_pp[orig_idx] = w_0_force;
+ if (devC_nw > 1)
+ dev_walls_force_pp[orig_idx+devC_np] = w_1_force;
+ if (devC_nw > 2)
+ dev_walls_force_pp[orig_idx+devC_np*2] = w_2_force;
+ if (devC_nw > 3)
+ dev_walls_force_pp[orig_idx+devC_np*3] = w_3_force;
+ if (devC_nw > 4)
+ dev_walls_force_pp[orig_idx+devC_np*4] = w_4_force;
}
} // End of interact(...)
(DIR) diff --git a/src/device.cu b/src/device.cu
t@@ -595,7 +595,8 @@ __host__ void DEM::startTime()
unsigned int smemSize = sizeof(unsigned int)*(threadsPerBlock+1);
// Pre-sum of force per wall
- cudaMalloc((void**)&dev_walls_force_partial, sizeof(Float)*dimGrid.x);
+ cudaMalloc((void**)&dev_walls_force_partial,
+ sizeof(Float)*dimGrid.x*walls.nw);
// Report to stdout
if (verbose == 1) {