tAdded parallel bonds, still untested - 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 11c8ccc15c6f107407771020cf27c7bd34233df4
(DIR) parent a2dc863676d9ddd731528d2c330be10b1aec48f2
(HTM) Author: Anders Damsgaard <adc@geo.au.dk>
Date: Sat, 23 Feb 2013 21:09:09 +0100
Added parallel bonds, still untested
Diffstat:
A python/bondtest.py | 44 +++++++++++++++++++++++++++++++
M python/sphere.py | 157 +++++++++++++++++++++++++------
M src/boost-unit-tests.cpp | 3 ++-
M src/cohesion.cuh | 174 ++++++++++++++++++++++++++++++-
M src/constants.cuh | 1 +
M src/datatypes.h | 63 +++++++++++++++++--------------
M src/device.cu | 47 ++++++++++++++++++++++++++++++-
M src/file_io.cpp | 54 +++++++++++++++++++++++++++++--
M src/forcechains.cpp | 17 ++++++++++++++---
M src/sphere.cpp | 21 +++++++++++++++------
M src/sphere.h | 7 ++++++-
11 files changed, 514 insertions(+), 74 deletions(-)
---
(DIR) diff --git a/python/bondtest.py b/python/bondtest.py
t@@ -0,0 +1,44 @@
+#!/usr/bin/env python
+
+from sphere import *
+
+sb = Spherebin(np=2, sid='bondtest')
+
+cleanup(sb)
+
+sb.x[0,:] = numpy.array((2,2,2))
+sb.x[1,:] = numpy.array((3.5,2,2))
+sb.radius = numpy.ones(sb.np)*0.5
+
+#sb.vel[1,2] = 1
+sb.angvel[1,1] = 0.01
+
+
+sb.initGridAndWorldsize(margin = 10, periodic = 1, contactmodel = 2, g = numpy.array([0.0, 0.0, 0.0]))
+
+sb.bond(0, 1)
+
+sb.defaultParams()
+#sb.gamma_n[0] = 10000
+#sb.initTemporal(total=4.0)
+sb.initTemporal(total=0.01, file_dt=2e-4)
+
+sb.writebin()
+
+sb.run(dry=True)
+sb.run()
+
+sb.render(verbose=False)
+
+visualize(sb.sid, "energy")
+
+sb.readlast()
+print(sb.bonds_delta_n)
+print(sb.bonds_delta_t)
+print(sb.bonds_omega_n)
+print(sb.bonds_omega_t)
+print()
+print(sb.force)
+print(sb.torque)
+print(sb.vel)
+print(sb.angvel)
(DIR) diff --git a/python/sphere.py b/python/sphere.py
t@@ -50,15 +50,15 @@ class Spherebin:
self.periodic = numpy.zeros(1, dtype=numpy.uint32)
# Particle data
- self.x = numpy.zeros(self.np*self.nd, dtype=numpy.float64).reshape(self.np,self.nd)
+ 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).reshape(self.np,2)
- self.vel = numpy.zeros(self.np*self.nd, dtype=numpy.float64).reshape(self.np,self.nd)
+ self.xysum = numpy.zeros((self.np, 2), 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).reshape(self.np,self.nd)
- self.angpos = numpy.zeros(self.np*self.nd, dtype=numpy.float64).reshape(self.np,self.nd)
- self.angvel = numpy.zeros(self.np*self.nd, dtype=numpy.float64).reshape(self.np,self.nd)
- self.torque = numpy.zeros(self.np*self.nd, dtype=numpy.float64).reshape(self.np,self.nd)
+ self.force = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
+ self.angpos = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
+ self.angvel = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
+ self.torque = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
self.es_dot = numpy.zeros(self.np, dtype=numpy.float64)
self.es = numpy.zeros(self.np, dtype=numpy.float64)
t@@ -90,7 +90,7 @@ class Spherebin:
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).reshape(self.nw,self.nd)
+ self.w_n = numpy.zeros((self.nw, self.nd), dtype=numpy.float64)
if (self.nw > 0):
self.w_n[0,2] = -1.0
self.w_x = numpy.ones(self.nw, dtype=numpy.float64)
t@@ -99,7 +99,13 @@ class Spherebin:
self.w_force = numpy.zeros(self.nw, dtype=numpy.float64)
self.w_devs = numpy.zeros(self.nw, dtype=numpy.float64)
+ self.lambda_bar = numpy.ones(1, dtype=numpy.float64)
self.nb0 = numpy.zeros(1, dtype=numpy.uint32)
+ self.bonds = numpy.zeros((self.nb0, 2), dtype=numpy.uint32)
+ self.bonds_delta_n = numpy.zeros(self.nb0, dtype=numpy.float64)
+ self.bonds_delta_t = numpy.zeros((self.nb0, self.nd), dtype=numpy.float64)
+ self.bonds_omega_n = numpy.zeros(self.nb0, dtype=numpy.float64)
+ self.bonds_omega_t = numpy.zeros((self.nb0, self.nd), dtype=numpy.float64)
def __cmp__(self, other):
""" Called when to Spherebin objects are compared.
t@@ -155,7 +161,13 @@ class Spherebin:
(self.w_devs == other.w_devs).all() and\
self.gamma_wn == other.gamma_wn and\
self.gamma_wt == other.gamma_wt and\
- self.nb0 == other.nb0\
+ self.lambda_bar == other.lambda_bar and\
+ self.nb0 == other.nb0 and\
+ self.bonds == other.bonds and\
+ self.bonds_delta_n == other.bonds_delta_n and\
+ self.bonds_delta_t == other.bonds_delta_t and\
+ self.bonds_omega_n == other.bonds_omega_n and\
+ self.bonds_omega_t == other.bonds_omega_t\
).all() == True):
return 0 # All equal
else :
t@@ -184,15 +196,15 @@ class Spherebin:
self.time_step_count = numpy.fromfile(fh, dtype=numpy.uint32, count=1)
# Allocate array memory for particles
- self.x = numpy.zeros(self.np*self.nd, dtype=numpy.float64).reshape(self.np,self.nd)
+ self.x = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
self.radius = numpy.zeros(self.np, dtype=numpy.float64)
- self.xysum = numpy.zeros(self.np*2, dtype=numpy.float64).reshape(self.np, 2)
- self.vel = numpy.zeros(self.np*self.nd, dtype=numpy.float64).reshape(self.np,self.nd)
+ self.xysum = numpy.zeros((self.np, 2), 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).reshape(self.np,self.nd)
- self.angpos = numpy.zeros(self.np*self.nd, dtype=numpy.float64).reshape(self.np,self.nd)
- self.angvel = numpy.zeros(self.np*self.nd, dtype=numpy.float64).reshape(self.np,self.nd)
- self.torque = numpy.zeros(self.np*self.nd, dtype=numpy.float64).reshape(self.np,self.nd)
+ #self.force = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
+ #self.angpos = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
+ #self.angvel = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
+ #self.torque = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
self.es_dot = numpy.zeros(self.np, dtype=numpy.float64)
self.es = numpy.zeros(self.np, dtype=numpy.float64)
self.ev_dot = numpy.zeros(self.np, dtype=numpy.float64)
t@@ -271,11 +283,20 @@ class Spherebin:
self.w_devs[i] = numpy.fromfile(fh, dtype=numpy.float64, count=1)
# Inter-particle bonds
+ self.lambda_bar = numpy.fromfile(fh, dtype=numpy.float64, count=1)
self.nb0 = numpy.fromfile(fh, dtype=numpy.uint32, count=1)
self.bonds = numpy.zeros((self.nb0, 2), dtype=numpy.uint32)
for i in range(self.nb0):
self.bonds[i,0] = numpy.fromfile(fh, dtype=numpy.uint32, count=1)
self.bonds[i,1] = numpy.fromfile(fh, dtype=numpy.uint32, count=1)
+ #self.bonds_delta_n = numpy.zeros(self.nb0, dtype=numpy.float64)
+ #self.bonds_delta_t = numpy.zeros((self.nb0, seld.nd), dtype=numpy.float64)
+ #self.bonds_omega_n = numpy.zeros(self.nb0, dtype=numpy.float64)
+ #self.bonds_omega_t = numpy.zeros((self.nb0, seld.nd), dtype=numpy.float64)
+ self.bonds_delta_n = numpy.fromfile(fh, dtype=numpy.float64, count=self.nb0)
+ self.bonds_delta_t = numpy.fromfile(fh, dtype=numpy.float64, count=self.nb0*self.nd).reshape(self.nb0*self.nd)
+ self.bonds_omega_n = numpy.fromfile(fh, dtype=numpy.float64, count=self.nb0)
+ self.bonds_omega_t = numpy.fromfile(fh, dtype=numpy.float64, count=self.nb0*self.nd).reshape(self.nb0*self.nd)
finally:
if fh is not None:
t@@ -366,15 +387,26 @@ class Spherebin:
fh.write(self.w_force[i].astype(numpy.float64))
fh.write(self.w_devs[i].astype(numpy.float64))
+ fh.write(self.lambda_bar.astype(numpy.float64))
fh.write(self.nb0.astype(numpy.uint32))
for i in range(self.nb0):
fh.write(self.bonds[i,0].astype(numpy.uint32))
fh.write(self.bonds[i,1].astype(numpy.uint32))
+ fh.write(self.bonds_delta_n.astype(numpy.float64))
+ fh.write(self.bonds_delta_t.astype(numpy.float64))
+ fh.write(self.bonds_omega_n.astype(numpy.float64))
+ fh.write(self.bonds_omega_t.astype(numpy.float64))
+
finally:
if fh is not None:
fh.close()
+ def readlast(self, verbose=True):
+ lastfile = status(self.sid)
+ fn = "../output/{0}.output{1:0=5}.bin".format(self.sid, lastfile)
+ self.readbin(fn, verbose)
+
def generateRadii(self, psd = 'logn',
radius_mean = 440e-6,
radius_variance = 8.8e-9,
t@@ -870,11 +902,35 @@ class Spherebin:
def bond(self, i, j):
""" Create a bond between particles i and j """
+
+ self.lambda_bar[0] = 1.0 # Radius multiplier to parallel-bond radii
if (hasattr(self, 'bonds') == False):
self.bonds = numpy.array([[i,j]], dtype=numpy.uint32)
else :
- self.bonds = numpy.vstack((bonds, [i,j]))
+ self.bonds = numpy.vstack((self.bonds, [i,j]))
+
+ if (hasattr(self, 'bonds_delta_n') == False):
+ self.bonds_delta_n = numpy.array([0.0], dtype=numpy.uint32)
+ else :
+ #self.bonds_delta_n = numpy.vstack((self.bonds_delta_n, [0.0]))
+ self.bonds_delta_n = numpy.append(self.bonds_delta_n, [0.0])
+
+ if (hasattr(self, 'bonds_delta_t') == False):
+ self.bonds_delta_t = numpy.array([[0.0, 0.0, 0.0]], dtype=numpy.uint32)
+ else :
+ self.bonds_delta_t = numpy.vstack((self.bonds_delta_t, [0.0, 0.0, 0.0]))
+
+ if (hasattr(self, 'bonds_omega_n') == False):
+ self.bonds_omega_n = numpy.array([0.0], dtype=numpy.uint32)
+ else :
+ #self.bonds_omega_n = numpy.vstack((self.bonds_omega_n, [0.0]))
+ self.bonds_omega_n = numpy.append(self.bonds_omega_n, [0.0])
+
+ if (hasattr(self, 'bonds_omega_t') == False):
+ self.bonds_omega_t = numpy.array([[0.0, 0.0, 0.0]], dtype=numpy.uint32)
+ else :
+ self.bonds_omega_t = numpy.vstack((self.bonds_omega_t, [0.0, 0.0, 0.0]))
# Increment the number of bonds with one
self.nb0 += 1
t@@ -983,20 +1039,27 @@ class Spherebin:
return numpy.array(porosity), numpy.array(depth)
- def run(self, verbose=True, hideinputfile=False, dry=False):
+ def run(self, verbose=True, hideinputfile=False, dry=False, valgrind=False, cudamemcheck=False):
'Execute sphere with target project'
quiet = ""
stdout = ""
dryarg = ""
+ valgrindbin = ""
+ cudamemchk = ""
if (verbose == False):
quiet = "-q "
if (hideinputfile == True):
stdout = " > /dev/null"
if (dry == True):
dryarg = "--dry "
+ if (valgrind == True):
+ valgrindbin = "valgrind -q "
+ if (cudamemcheck == True):
+ cudamemchk = "cuda-memcheck "
- subprocess.call("cd ..; ./sphere " + quiet + dryarg + "input/" + self.sid + ".bin " + stdout, shell=True)
+
+ subprocess.call("cd ..; " + valgrindbin + cudamemchk + "./sphere " + quiet + dryarg + "input/" + self.sid + ".bin " + stdout, shell=True)
def torqueScript(self,
t@@ -1030,11 +1093,16 @@ class Spherebin:
quiet = "-q"
# Render images using sphere raytracer
- subprocess.call("cd ..; for F in `ls output/" + self.sid + "*.bin`; do ./sphere " + quiet \
- + " --method " + method + " {}".format(max_val) \
- + " -l {}".format(lower_cutoff) \
- + " --render $F; done" \
- , shell=True)
+ if (method == "normal"):
+ subprocess.call("cd ..; for F in `ls output/" + self.sid + "*.bin`; do ./sphere " + quiet \
+ + " --render $F; done" \
+ , shell=True)
+ else :
+ subprocess.call("cd ..; for F in `ls output/" + self.sid + "*.bin`; do ./sphere " + quiet \
+ + " --method " + method + " {}".format(max_val) \
+ + " -l {}".format(lower_cutoff) \
+ + " --render $F; done" \
+ , shell=True)
# Convert images to compressed format
convert()
t@@ -1066,6 +1134,7 @@ class Spherebin:
graphicsformat = 'png',
cbmax = None,
arrowscale = 0.01,
+ velarrowscale = 1.0,
slipscale = 1.0,
verbose = False):
''' Produce a 2D image of particles on a x1,x3 plane, intersecting the second axis at x2.
t@@ -1094,6 +1163,8 @@ class Spherebin:
aylist = []
daxlist = []
daylist = []
+ dvxlist = []
+ dvylist = []
# Loop over all particles, find intersections
for i in range(self.np):
t@@ -1133,6 +1204,10 @@ class Spherebin:
daylist.append(0.0) # delta y for arrow end point
daylist.append(0.0) # delta y for arrow end point
+ # Store linear velocity data
+ dvxlist.append(self.vel[i,0]*velarrowscale) # delta x for arrow end point
+ dvylist.append(self.vel[i,2]*velarrowscale) # delta y for arrow end point
+
if (r_circ > self.radius[i]):
raise Exception("Error, circle radius is larger than the particle radius")
if (self.p[i] > pmax):
t@@ -1171,6 +1246,22 @@ class Spherebin:
if fh is not None:
fh.close()
+ # Save linear velocity data
+ # Output format: x, y, deltax, deltay
+ # gnuplot> plot '-' using 1:2:3:4 with vectors head filled lt 2
+ filename = '../gnuplot/data/' + self.sid + '-ts-x1x3-velarrows.txt'
+ fh = None
+ try :
+ fh = open(filename, 'w')
+
+ for (x, y, dvx, dvy) in zip(xlist, ylist, dvxlist, dvylist):
+ fh.write("{}\t{}\t{}\t{}\n".format(x, y, dvx, dvy))
+
+ finally :
+ if fh is not None:
+ fh.close()
+
+
# Check whether there are slips between the particles intersecting the plane
sxlist = []
t@@ -1303,11 +1394,14 @@ def render(binary,
quiet = "-q"
# Render images using sphere raytracer
- subprocess.call("cd .. ; ./sphere " + quiet + \
- " --method " + method + " {}".format(max_val) + \
- " -l {}".format(lower_cutoff) + \
- " --render " + binary, shell=True)
-
+ if (method == "normal"):
+ subprocess.call("cd .. ; ./sphere " + quiet + \
+ " --render " + binary, shell=True)
+ else :
+ subprocess.call("cd .. ; ./sphere " + quiet + \
+ " --method " + method + " {}".format(max_val) + \
+ " -l {}".format(lower_cutoff) + \
+ " --render " + binary, shell=True)
# Convert images to compressed format
convert()
t@@ -1670,6 +1764,11 @@ def status(project):
if fh is not None:
fh.close()
+def cleanup(spherebin):
+ 'Remove input/output files and images from simulation'
+ subprocess.call("rm -f ../input/" + spherebin.sid + ".bin", shell=True)
+ subprocess.call("rm -f ../output/" + spherebin.sid + ".*.bin", shell=True)
+ subprocess.call("rm -f ../img_out/" + spherebin.sid + ".*", shell=True)
def V_sphere(r):
""" Returns the volume of a sphere with radius r
(DIR) diff --git a/src/boost-unit-tests.cpp b/src/boost-unit-tests.cpp
t@@ -7,10 +7,11 @@
// define the name of the test suite
BOOST_AUTO_TEST_SUITE (spheretest)
-BOOST_AUTO_TEST_SUITE (inittest)
+BOOST_AUTO_TEST_CASE (inittest)
{
DEM testDEM();
}
BOOST_AUTO_TEST_SUITE_END()
+
(DIR) diff --git a/src/cohesion.cuh b/src/cohesion.cuh
t@@ -4,10 +4,182 @@
// cohesion.cuh
// Functions governing attractive forces between contacts
+// Check bond pair list, apply linear contact model to pairs
+__global__ void bondsLinear(
+ uint2* dev_bonds,
+ Float4* dev_bonds_delta, // Contact displacement
+ Float4* dev_bonds_omega, // Contact rotational displacement
+ Float4* dev_x,
+ Float4* dev_vel,
+ Float4* dev_angvel,
+ Float4* dev_force,
+ Float4* dev_torque)
+{
+ // Find thread index
+ unsigned int idx = blockIdx.x * blockDim.x + threadIdx.x;
+ if (idx >= devC_params.nb0)
+ return;
+
+
+
+ //// Read values
+
+ // Read bond data
+ __syncthreads();
+ const uint2 bond = dev_bonds[idx]; // Particle indexes in bond pair
+
+ // Check if the bond has been erased
+ if (bond.x >= devC_np)
+ return;
+
+ const Float4 delta_t0_4 = dev_bonds_delta[idx];
+ const Float4 omega_t0_4 = dev_bonds_omega[idx];
+
+ // Convert tangential vectors to Float3's
+ const Float3 delta_t0_uncor = MAKE_FLOAT3(
+ delta_t0_4.x,
+ delta_t0_4.y,
+ delta_t0_4.z);
+ const Float delta_t0_n = delta_t0_4.w;
+
+ const Float3 omega_t0_uncor = MAKE_FLOAT3(
+ omega_t0_4.x,
+ omega_t0_4.y,
+ omega_t0_4.z);
+ const Float omega_t0_n = omega_t0_4.w;
+
+ // Read particle data
+ const Float4 x_i = dev_x[bond.x];
+ const Float4 x_j = dev_x[bond.y];
+ const Float4 vel_i = dev_vel[bond.x];
+ const Float4 vel_j = dev_vel[bond.y];
+ const Float4 angvel_i = dev_angvel[bond.x];
+ const Float4 angvel_j = dev_angvel[bond.y];
+
+ // Initialize force- and torque vectors
+ Float3 f, t, f_n, f_t, t_n, t_t;
+
+
+ //// Bond geometry and inertia
+
+ // Parallel-bond radius (Potyondy and Cundall 2004, eq. 12)
+ const Float R_bar = devC_params.lambda_bar * fmin(x_i.w, x_j.w);
+
+ // Bond cross section area (Potyondy and Cundall 2004, eq. 15)
+ const Float A = PI * R_bar*R_bar;
+
+ // Bond moment of inertia (Potyondy and Cundall 2004, eq. 15)
+ const Float I = 0.25 * PI * R_bar*R_bar*R_bar*R_bar;
+
+ // Bond polar moment of inertia (Potyondy and Cundall 2004, eq. 15)
+ const Float J = 0.50 * PI * R_bar*R_bar*R_bar*R_bar;
+
+ // Inter-particle vector
+ const Float3 x = MAKE_FLOAT3(
+ x_i.x - x_j.x,
+ x_i.y - x_j.y,
+ x_i.z - x_j.z);
+ const Float x_length = length(x);
+
+ // Normal vector of contact (points from i to j)
+ const Float3 n = x/x_length;
+
+
+ //// Force
+
+ // Correct tangential displacement vector for rotation of the contact plane
+ //const Float3 delta_t0 = delta_t0_uncor - dot(delta_t0_uncor, n);
+ const Float3 delta_t0 = delta_t0_uncor - (n * dot(n, delta_t0_uncor));
+
+ // Contact displacement (should this include rolling?)
+ const Float3 ddelta = MAKE_FLOAT3(
+ vel_i.x - vel_j.x,
+ vel_i.y - vel_j.y,
+ vel_i.z - vel_j.z) * devC_dt;
+
+ // Normal component of the displacement increment
+ //const Float ddelta_n = dot(ddelta, n);
+ const Float ddelta_n = -dot(ddelta, n);
+
+ // Normal component of the total displacement
+ const Float delta_n = delta_t0_n + ddelta_n;
+
+ // Tangential component of the displacement increment
+ //const Float3 ddelta_t = ddelta - dot(ddelta, n);
+ const Float3 ddelta_t = ddelta - n * dot(n, ddelta);
+
+ // Tangential component of the total displacement
+ const Float3 delta_t = delta_t0 + ddelta_t;
+
+ // Normal force: Elastic contact model
+ //f_n = devC_params.k_n * A * delta_n * n;
+ f_n = (devC_params.k_n * A * delta_n + devC_params.gamma_n * ddelta_n/devC_dt) * n;
+
+ // Tangential force: Elastic contact model
+ //f_t = -devC_params.k_t * A * delta_t;
+ f_t = -devC_params.k_t * A * delta_t - devC_params.gamma_t * ddelta_t/devC_dt;
+
+ // Force vector
+ f = f_n + f_t;
+
+
+ //// Torque
+
+ // Correct tangential rotational vector for rotation of the contact plane
+ //Float3 omega_t0 = omega_t0_uncor - dot(omega_t0_uncor, n);
+ Float3 omega_t0 = omega_t0_uncor - (n * dot(n, omega_t0_uncor));
+
+ // Contact rotational velocity
+ Float3 domega = MAKE_FLOAT3(
+ angvel_j.x - angvel_i.x,
+ angvel_j.y - angvel_i.y,
+ angvel_j.z - angvel_i.z) * devC_dt;
+
+ // Normal component of the rotational increment
+ //const Float domega_n = dot(domega, n);
+ const Float domega_n = -dot(n, domega);
+
+ // Normal component of the total displacement
+ const Float omega_n = omega_t0_n + domega_n;
+
+ // Tangential component of the displacement increment
+ //const Float3 domega_t = domega - dot(domega, n);
+ const Float3 domega_t = domega - n * dot(n, domega);
+
+ // Tangential component of the total displacement
+ const Float3 omega_t = omega_t0 + domega_t;
+
+ // Twisting torque: Elastic contact model
+ //t_n = -devC_params.k_t * J * omega_n * n;
+ t_n = (devC_params.k_t * J * omega_n + devC_params.gamma_t * domega_n/devC_dt) * n;
+
+ // Bending torque: Elastic contact model
+ //t_t = -devC_params.k_n * I * omega_t;
+ //t_t = -devC_params.k_n * I * omega_t - devC_params.gamma_n * domega_t/devC_dt;
+ t_t = devC_params.k_n * I * omega_t - devC_params.gamma_n * domega_t/devC_dt;
+
+ // Torque vector
+ t = t_n + t_t;
+
+
+ //// Save values
+ __syncthreads();
+
+ // Save updated displacements in global memory
+ dev_bonds_delta[idx] = MAKE_FLOAT4(delta_t.x, delta_t.y, delta_t.z, delta_n);
+ 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
+ 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);
+ dev_torque[bond.y] -= MAKE_FLOAT4(t.x, t.y, t.z, 0.0);
+ // make sure to remove write conflicts
+}
// Linear-elastic bond: Attractive force with normal- and shear components
// acting upon particle A in a bonded particle pair
-__device__ void bondLinear(Float3* N, Float3* T, Float* es_dot, Float* p,
+__device__ void bondLinear_old(Float3* N, Float3* T, Float* es_dot, Float* p,
unsigned int idx_a, unsigned int idx_b,
Float4* dev_x_sorted, Float4* dev_vel_sorted,
Float4* dev_angvel_sorted,
(DIR) diff --git a/src/constants.cuh b/src/constants.cuh
t@@ -12,6 +12,7 @@ __constant__ unsigned int devC_np; // Number of particles
__constant__ unsigned int devC_nw; // Number of dynamic walls
__constant__ int devC_nc; // Max. number of contacts a particle can have
__constant__ Float devC_dt; // Computational time step length
+__constant__ unsigned int devC_nb0; // Number of inter-particle bonds at t=0
// Device constant memory structures
__constant__ Params devC_params;
(DIR) diff --git a/src/datatypes.h b/src/datatypes.h
t@@ -15,18 +15,21 @@
// Structure containing kinematic particle values
struct Kinematics {
- Float4 *x; // Positions + radii (w)
- Float2 *xysum; // Horizontal distance traveled
- Float4 *vel; // Translational velocities + fixvels (w)
- Float4 *acc; // Translational accelerations
- Float4 *force; // Sums of forces
- Float4 *angpos; // Angular positions
- Float4 *angvel; // Angular velocities
- Float4 *angacc; // Angular accelerations
- Float4 *torque; // Sums of torques
+ Float4 *x; // Positions + radii (w)
+ Float2 *xysum; // Horizontal distance traveled
+ Float4 *vel; // Translational velocities + fixvels (w)
+ Float4 *acc; // Translational accelerations
+ Float4 *force; // Sums of forces
+ Float4 *angpos; // Angular positions
+ Float4 *angvel; // Angular velocities
+ Float4 *angacc; // Angular accelerations
+ Float4 *torque; // Sums of torques
unsigned int *contacts; // List of contacts per particle
- Float4 *distmod; // Distance modifiers for contacts across periodic boundaries
- Float4 *delta_t; // Accumulated shear distance of contacts
+ Float4 *distmod; // Distance modifiers for contacts across periodic boundaries
+ Float4 *delta_t; // Accumulated shear distance of contacts
+ uint2 *bonds; // Particle bond pairs
+ Float4 *bonds_delta; // Particle bond displacement
+ Float4 *bonds_omega; // Particle bond rotation
};
// Structure containing individual particle energies
t@@ -68,25 +71,27 @@ struct Time {
// Structure containing constant, global physical parameters
struct Params {
- Float g[ND]; // Gravitational acceleration
- Float k_n; // Normal stiffness
- Float k_t; // Tangential stiffness
- Float k_r; // Rotational stiffness
- Float gamma_n; // Normal viscosity
- Float gamma_t; // Tangential viscosity
- Float gamma_r; // Rotational viscosity
- Float mu_s; // Static friction coefficient
- Float mu_d; // Dynamic friction coefficient
- Float mu_r; // Rotational friction coefficient
- Float gamma_wn; // Wall normal viscosity
- Float gamma_wt; // Wall tangential viscosity
- Float mu_ws; // Wall static friction coefficient
- Float mu_wd; // Wall dynamic friction coefficient
- Float rho; // Material density
+ Float g[ND]; // Gravitational acceleration
+ Float k_n; // Normal stiffness
+ Float k_t; // Tangential stiffness
+ Float k_r; // Rotational stiffness
+ Float gamma_n; // Normal viscosity
+ Float gamma_t; // Tangential viscosity
+ Float gamma_r; // Rotational viscosity
+ Float mu_s; // Static friction coefficient
+ Float mu_d; // Dynamic friction coefficient
+ Float mu_r; // Rotational friction coefficient
+ Float gamma_wn; // Wall normal viscosity
+ Float gamma_wt; // Wall tangential viscosity
+ Float mu_ws; // Wall static friction coefficient
+ Float mu_wd; // Wall dynamic friction coefficient
+ Float rho; // Material density
unsigned int contactmodel; // Inter-particle contact model
- Float kappa; // Capillary bond prefactor
- Float db; // Capillary bond debonding distance
- Float V_b; // Volume of fluid in capillary bond
+ Float kappa; // Capillary bond prefactor
+ Float db; // Capillary bond debonding distance
+ Float V_b; // Volume of fluid in capillary bond
+ Float lambda_bar; // Radius multiplier to parallel-bond radii
+ unsigned int nb0; // Number of inter-particle bonds at t=0
};
// Structure containing wall parameters
(DIR) diff --git a/src/device.cu b/src/device.cu
t@@ -137,7 +137,9 @@ __global__ void checkConstantValues(int* dev_equal,
dev_params->contactmodel != devC_params.contactmodel ||
dev_params->kappa != devC_params.kappa ||
dev_params->db != devC_params.db ||
- dev_params->V_b != devC_params.V_b)
+ dev_params->V_b != devC_params.V_b ||
+ dev_params->lambda_bar != devC_params.lambda_bar ||
+ dev_params->nb0 != devC_params.nb0)
*dev_equal = 2; // Not ok
}
t@@ -250,6 +252,9 @@ __host__ void DEM::allocateGlobalDeviceMemory(void)
cudaMalloc((void**)&dev_contacts, sizeof(unsigned int)*np*NC); // Max NC contacts per particle
cudaMalloc((void**)&dev_distmod, memSizeF4*NC);
cudaMalloc((void**)&dev_delta_t, memSizeF4*NC);
+ cudaMalloc((void**)&dev_bonds, sizeof(uint2)*params.nb0);
+ cudaMalloc((void**)&dev_bonds_delta, sizeof(Float4)*params.nb0);
+ cudaMalloc((void**)&dev_bonds_omega, sizeof(Float4)*params.nb0);
// Sorted arrays
cudaMalloc((void**)&dev_x_sorted, memSizeF4);
t@@ -311,6 +316,9 @@ __host__ void DEM::freeGlobalDeviceMemory()
cudaFree(dev_contacts);
cudaFree(dev_distmod);
cudaFree(dev_delta_t);
+ cudaFree(dev_bonds);
+ cudaFree(dev_bonds_delta);
+ cudaFree(dev_bonds_omega);
cudaFree(dev_es_dot);
cudaFree(dev_es);
t@@ -381,6 +389,12 @@ __host__ void DEM::transferToGlobalDeviceMemory()
memSizeF4*NC, cudaMemcpyHostToDevice);
cudaMemcpy( dev_delta_t, k.delta_t,
memSizeF4*NC, cudaMemcpyHostToDevice);
+ cudaMemcpy( dev_bonds, k.bonds,
+ sizeof(uint2)*params.nb0, cudaMemcpyHostToDevice);
+ cudaMemcpy( dev_bonds_delta, k.bonds_delta,
+ sizeof(Float4)*params.nb0, cudaMemcpyHostToDevice);
+ cudaMemcpy( dev_bonds_omega, k.bonds_omega,
+ sizeof(Float4)*params.nb0, cudaMemcpyHostToDevice);
// Individual particle energy values
cudaMemcpy( dev_es_dot, e.es_dot,
t@@ -447,6 +461,12 @@ __host__ void DEM::transferFromGlobalDeviceMemory()
memSizeF4*NC, cudaMemcpyDeviceToHost);
cudaMemcpy( k.delta_t, dev_delta_t,
memSizeF4*NC, cudaMemcpyDeviceToHost);
+ cudaMemcpy( k.bonds, dev_bonds,
+ sizeof(uint2)*params.nb0, cudaMemcpyDeviceToHost);
+ cudaMemcpy( k.bonds_delta, dev_bonds_delta,
+ sizeof(Float4)*params.nb0, cudaMemcpyDeviceToHost);
+ cudaMemcpy( k.bonds_omega, dev_bonds_omega,
+ sizeof(Float4)*params.nb0, cudaMemcpyDeviceToHost);
// Individual particle energy values
cudaMemcpy( e.es_dot, dev_es_dot,
t@@ -468,6 +488,9 @@ __host__ void DEM::transferFromGlobalDeviceMemory()
cudaMemcpy( walls.mvfd, dev_walls_mvfd,
sizeof(Float4)*walls.nw, cudaMemcpyDeviceToHost);
+ // Bond parameters
+
+
checkForCudaErrors("End of transferFromGlobalDeviceMemory");
}
t@@ -557,6 +580,7 @@ __host__ void DEM::startTime()
double t_reorderArrays = 0.0;
double t_topology = 0.0;
double t_interact = 0.0;
+ double t_bondsLinear = 0.0;
double t_integrate = 0.0;
double t_summation = 0.0;
double t_integrateWalls = 0.0;
t@@ -698,6 +722,25 @@ __host__ void DEM::startTime()
stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed, &t_interact);
checkForCudaErrors("Post interact - often caused if particles move outside the grid", iter);
+ // Process particle pairs
+ if (params.nb0 > 0) {
+ bondsLinear<<< 1, params.nb0 >>>(
+ dev_bonds,
+ dev_bonds_delta,
+ dev_bonds_omega,
+ dev_x,
+ dev_vel,
+ dev_angvel,
+ dev_force,
+ dev_torque);
+ // Synchronization point
+ cudaThreadSynchronize();
+ //cudaPrintfDisplay(stdout, true);
+ if (PROFILING == 1)
+ stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed, &t_bondsLinear);
+ checkForCudaErrors("Post bondsLinear", iter);
+ }
+
// Update particle kinematics
if (PROFILING == 1)
t@@ -875,6 +918,8 @@ __host__ void DEM::startTime()
<< "\t(" << 100.0*t_topology/t_sum << " %)\n"
<< " - interact:\t\t" << t_interact/1000.0 << " s"
<< "\t(" << 100.0*t_interact/t_sum << " %)\n"
+ << " - bondsLinear:\t" << t_bondsLinear/1000.0 << " s"
+ << "\t(" << 100.0*t_bondsLinear/t_sum << " %)\n"
<< " - integrate:\t\t" << t_integrate/1000.0 << " s"
<< "\t(" << 100.0*t_integrate/t_sum << " %)\n"
<< " - summation:\t\t" << t_summation/1000.0 << " s"
(DIR) diff --git a/src/file_io.cpp b/src/file_io.cpp
t@@ -74,9 +74,9 @@ void DEM::readbin(const char *target)
if (verbose == 1)
cout << " Allocating host memory: ";
// Allocate more host arrays
- k.x = new Float4[np];
+ k.x = new Float4[np];
k.xysum = new Float2[np];
- k.vel = new Float4[np];
+ k.vel = new Float4[np];
k.force = new Float4[np];
k.angpos = new Float4[np];
k.angvel = new Float4[np];
t@@ -86,7 +86,7 @@ void DEM::readbin(const char *target)
e.es = new Float[np];
e.ev_dot = new Float[np];
e.ev = new Float[np];
- e.p = new Float[np];
+ e.p = new Float[np];
if (verbose == 1)
cout << "Done\n";
t@@ -205,6 +205,31 @@ void DEM::readbin(const char *target)
ifs.read(as_bytes(walls.mvfd[i].w), sizeof(Float));
}
+ // Read bond parameters
+ ifs.read(as_bytes(params.lambda_bar), sizeof(params.lambda_bar));
+ ifs.read(as_bytes(params.nb0), sizeof(params.nb0));
+ k.bonds = new uint2[params.nb0];
+ k.bonds_delta = new Float4[np];
+ k.bonds_omega = new Float4[np];
+ for (i = 0; i<params.nb0; ++i) {
+ ifs.read(as_bytes(k.bonds[i].x), sizeof(unsigned int));
+ ifs.read(as_bytes(k.bonds[i].y), sizeof(unsigned int));
+ }
+ for (i = 0; i<params.nb0; ++i) // Normal component
+ ifs.read(as_bytes(k.bonds_delta[i].w), sizeof(Float));
+ for (i = 0; i<params.nb0; ++i) { // Tangential component
+ ifs.read(as_bytes(k.bonds_delta[i].x), sizeof(Float));
+ ifs.read(as_bytes(k.bonds_delta[i].y), sizeof(Float));
+ ifs.read(as_bytes(k.bonds_delta[i].z), sizeof(Float));
+ }
+ for (i = 0; i<params.nb0; ++i) // Normal component
+ ifs.read(as_bytes(k.bonds_omega[i].w), sizeof(Float));
+ for (i = 0; i<params.nb0; ++i) { // Tangential component
+ ifs.read(as_bytes(k.bonds_omega[i].x), sizeof(Float));
+ ifs.read(as_bytes(k.bonds_omega[i].y), sizeof(Float));
+ ifs.read(as_bytes(k.bonds_omega[i].z), sizeof(Float));
+ }
+
// Close file if it is still open
if (ifs.is_open())
ifs.close();
t@@ -337,12 +362,35 @@ void DEM::writebin(const char *target)
ofs.write(as_bytes(walls.mvfd[i].w), sizeof(Float));
}
+ // Write bond parameters
+ ofs.write(as_bytes(params.lambda_bar), sizeof(params.lambda_bar));
+ ofs.write(as_bytes(params.nb0), sizeof(params.nb0));
+ for (i = 0; i<params.nb0; ++i) {
+ ofs.write(as_bytes(k.bonds[i].x), sizeof(unsigned int));
+ ofs.write(as_bytes(k.bonds[i].y), sizeof(unsigned int));
+ }
+ for (i = 0; i<params.nb0; ++i) // Normal component
+ ofs.write(as_bytes(k.bonds_delta[i].w), sizeof(Float));
+ for (i = 0; i<params.nb0; ++i) { // Tangential component
+ ofs.write(as_bytes(k.bonds_delta[i].x), sizeof(Float));
+ ofs.write(as_bytes(k.bonds_delta[i].y), sizeof(Float));
+ ofs.write(as_bytes(k.bonds_delta[i].z), sizeof(Float));
+ }
+ for (i = 0; i<params.nb0; ++i) // Normal component
+ ofs.write(as_bytes(k.bonds_omega[i].w), sizeof(Float));
+ for (i = 0; i<params.nb0; ++i) { // Tangential component
+ ofs.write(as_bytes(k.bonds_omega[i].x), sizeof(Float));
+ ofs.write(as_bytes(k.bonds_omega[i].y), sizeof(Float));
+ ofs.write(as_bytes(k.bonds_omega[i].z), sizeof(Float));
+ }
+
// Close file if it is still open
if (ofs.is_open())
ofs.close();
} else {
std::cerr << "Can't write output when in single precision mode.\n";
+ exit(1);
}
}
(DIR) diff --git a/src/forcechains.cpp b/src/forcechains.cpp
t@@ -30,6 +30,8 @@ int main(const int argc, const char *argv[])
int verbose = 0;
int nfiles = 0; // number of input files
int threedim = 1; // 0 if 2d, 1 if 3d
+ double lowercutoff = 0.0;
+ double uppercutoff = 1.0e9;
std::string format = "interactive"; // gnuplot terminal type
// Process input parameters
t@@ -45,9 +47,12 @@ int main(const int argc, const char *argv[])
<< "-h, --help\t\tPrint help\n"
<< "-V, --version\t\tPrint version information and exit\n"
<< "-v, --verbose\t\tDisplay in-/output file names\n"
+ << "-lc <val>, --lower-cutoff <val>\t\tOnly show contacts where the force value is greater\n"
+ << "-uc <val>, --upper-cutoff <val>\t\tOnly show contacts where the force value is greater\n"
<< "-f, --format\t\tOutput format to stdout, interactive default. Possible values:\n"
- << "\t\t\tinteractive, png, epslatex, epslatex-color"
- << "-2d\t\t\twrite output as 2d coordinates (3d default)"
+ << "\t\t\tinteractive, png, epslatex, epslatex-color\n"
+ << "-2d\t\t\twrite output as 2d coordinates (3d default)\n"
+ << "The values below the cutoff are not visualized, the values above are truncated to the upper limit\n"
<< std::endl;
return 0; // Exit with success
}
t@@ -65,6 +70,12 @@ int main(const int argc, const char *argv[])
else if (argvi == "-f" || argvi == "--format")
format = argv[++i];
+ else if (argvi == "-lc" || argvi == "--lower-cutoff")
+ lowercutoff = atof(argv[++i]);
+
+ else if (argvi == "-uc" || argvi == "--upper-cutoff")
+ uppercutoff = atof(argv[++i]);
+
else if (argvi == "-2d")
threedim = 0;
t@@ -79,7 +90,7 @@ int main(const int argc, const char *argv[])
DEM dem(argvi, verbose, 0, 0, 0);
// Calculate porosity and save as file
- dem.forcechains(format, threedim);
+ dem.forcechains(format, threedim, lowercutoff, uppercutoff);
}
}
(DIR) diff --git a/src/sphere.cpp b/src/sphere.cpp
t@@ -304,6 +304,8 @@ void DEM::reportValues()
<< grid.num[1] << " * "
<< grid.num[2];
cout << " cells\n";
+
+ cout << " - No. of particle bonds: " << params.nb0 << endl;
}
// Returns the volume of a spherical cap
t@@ -569,7 +571,9 @@ void DEM::findOverlaps(
}
// Calculate force chains and generate visualization script
-void DEM::forcechains(const std::string format, const int threedim)
+void DEM::forcechains(const std::string format, const int threedim,
+ const double lower_cutoff,
+ const double upper_cutoff)
{
using std::cout;
using std::endl;
t@@ -589,8 +593,8 @@ void DEM::forcechains(const std::string format, const int threedim)
// Define limits of visualization [0;1]
//Float lim_low = 0.1;
- Float lim_low = 0.05;
- Float lim_high = 0.25;
+ //Float lim_low = 0.15;
+ //Float lim_high = 0.25;
if (format == "txt") {
// Write text header
t@@ -640,7 +644,8 @@ void DEM::forcechains(const std::string format, const int threedim)
//<< "set palette defined (0 'gray', 0.5 'blue', 1 'red')\n"
//<< "set palette defined (0 'white', 0.5 'gray', 1 'red')\n"
<< "set palette defined ( 1 '#000fff', 2 '#0090ff', 3 '#0fffee', 4 '#90ff70', 5 '#ffee00', 6 '#ff7000', 7 '#ee0000', 8 '#7f0000')\n"
- << "set cbrange [" << f_n_max*lim_low << ':' << f_n_max*lim_high << "]\n"
+ //<< "set cbrange [" << f_n_max*lim_low << ':' << f_n_max*lim_high << "]\n"
+ << "set cbrange [" << lower_cutoff << ':' << upper_cutoff << "]\n"
<< endl;
}
t@@ -660,6 +665,9 @@ void DEM::forcechains(const std::string format, const int threedim)
// Normal force on contact
f_n = -params.k_n * delta_n;
+ if (f_n < lower_cutoff)
+ continue; // skip the rest of this iteration
+
// Line weight
ratio = f_n/f_n_max;
t@@ -674,13 +682,14 @@ void DEM::forcechains(const std::string format, const int threedim)
if (threedim == 1)
cout << k.x[j].y << '\t';
cout << k.x[j].z << '\t';
- cout << -delta_n*params.k_n << endl;
+ cout << f_n << endl;
} else {
// Gnuplot output
// Save contact pairs if they are above the lower limit
// and not fixed at their horizontal velocity
- if (ratio > lim_low && (k.vel[i].w + k.vel[j].w) == 0.0) {
+ //if (ratio > lim_low && (k.vel[i].w + k.vel[j].w) == 0.0) {
+ if (f_n > lower_cutoff && (k.vel[i].w + k.vel[j].w) == 0.0) {
// Plot contact as arrow without tip
cout << "set arrow " << n+1 << " from "
(DIR) diff --git a/src/sphere.h b/src/sphere.h
t@@ -85,6 +85,9 @@ class DEM {
Float *dev_walls_force_partial; // Pre-sum per wall
Float *dev_walls_force_pp; // Force per particle per wall
Float *dev_walls_vel0; // Half-step velocity
+ uint2 *dev_bonds; // Particle bond pairs
+ Float4 *dev_bonds_delta; // Particle bond displacement
+ Float4 *dev_bonds_omega; // Particle bond rotation
unsigned char *dev_img;
float4 *dev_ray_origo; // Ray data always single precision
float4 *dev_ray_direction;
t@@ -192,7 +195,9 @@ class DEM {
// Calculate force chains and save as Gnuplot script
void forcechains(
const std::string format = "interactive",
- const int threedim = 1);
+ const int threedim = 1,
+ const double lower_cutoff = 0.0,
+ const double upper_cutoff = 1.0e9);
};