tbonds working, tests passed - 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 b67be3fda4d26ef813cb6af47cbae812adf68e7a
(DIR) parent c530064dfebbdf97dafcdd137bfcdb90c63d954d
(HTM) Author: Anders Damsgaard <adc@geo.au.dk>
Date: Mon, 11 Mar 2013 21:54:36 +0100
bonds working, tests passed
Diffstat:
M src/cohesion.cuh | 105 +++++++++++++++++++-------------
M src/device.cu | 7 +++++++
M src/sphere.cpp | 5 ++++-
M tests/bond_tests.py | 63 +++++++++++++++++--------------
M tests/pytestutils.py | 12 ++++++++++++
5 files changed, 120 insertions(+), 72 deletions(-)
---
(DIR) diff --git a/src/cohesion.cuh b/src/cohesion.cuh
t@@ -32,21 +32,23 @@ __global__ void bondsLinear(
if (bond.x >= devC_np)
return;
- const Float4 delta_t0_4 = dev_bonds_delta[idx];
- const Float4 omega_t0_4 = dev_bonds_omega[idx];
+ const Float4 delta0_4 = dev_bonds_delta[idx];
+ const Float4 omega0_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;
+ // Uncorrected tangential component of displacement
+ Float3 delta0_t = MAKE_FLOAT3(
+ delta0_4.x,
+ delta0_4.y,
+ delta0_4.z);
+ const Float delta0_n = delta0_4.w;
+
+ // Uncorrected tangential component of rotation
+ Float3 omega0_t = MAKE_FLOAT3(
+ omega0_4.x,
+ omega0_4.y,
+ omega0_4.z);
+ const Float omega0_n = omega0_4.w;
// Read particle data
const Float4 x_i = dev_x[bond.x];
t@@ -59,9 +61,6 @@ __global__ void bondsLinear(
const Float3 angvel_i = MAKE_FLOAT3(angvel4_i.x, angvel4_i.y, angvel4_i.z);
const Float3 angvel_j = MAKE_FLOAT3(angvel4_j.x, angvel4_j.y, angvel4_j.z);
- // Initialize force- and torque vectors
- Float3 f, t, f_n, f_t, t_n, t_t;
-
//// Bond geometry and inertia
t@@ -89,16 +88,16 @@ __global__ void bondsLinear(
const Float overlap = fmin(0.0, x_length - (x_i.w + x_j.w));
// Normal vector of contact (points from i to j)
- const Float3 n = x/x_length;
+ 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));
+ delta0_t = delta0_t - (n * dot(n, delta0_t));
- // Contact displacement
+ // Contact displacement, Luding 2008 eq. 10
const Float3 ddelta = (
MAKE_FLOAT3(
vel_i.x - vel_j.x,
t@@ -113,68 +112,84 @@ __global__ void bondsLinear(
const Float ddelta_n = -dot(ddelta, n);
// Normal component of the total displacement
- const Float delta_n = delta_t0_n + ddelta_n;
+ const Float delta_n = delta0_n + ddelta_n;
// Tangential component of the displacement increment
- //const Float3 ddelta_t = ddelta - dot(ddelta, n);
+ // Luding 2008, eq. 9
const Float3 ddelta_t = ddelta - n * dot(n, ddelta);
// Tangential component of the total displacement
- const Float3 delta_t = delta_t0 + ddelta_t;
+ const Float3 delta_t = delta0_t + ddelta_t;
// Normal force: Visco-elastic contact model
// The elastic component caused by the overlap is subtracted.
//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;
+ const Float3 f_n = (devC_params.k_n * A * delta_n + devC_params.gamma_n * ddelta_n/devC_dt) * n;
//f_n += devC_params.k_n * overlap * n;
// Tangential force: Visco-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;
+ const Float3 f_t = -devC_params.k_t * A * delta_t - devC_params.gamma_t * ddelta_t/devC_dt;
// Force vector
- f = f_n + f_t;
+ const Float3 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));
+ omega0_t = omega0_t - (-n) * dot(omega0_t, -n);
+ //omega0_t = omega0_t - (n * dot(n, omega0_t));
// 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;
+ /*const Float3 domega = MAKE_FLOAT3(
+ angvel_i.x - angvel_j.x,
+ angvel_i.y - angvel_j.y,
+ angvel_i.z - angvel_j.z) * devC_dt;*/
// Normal component of the rotational increment
- //const Float domega_n = dot(domega, n);
- const Float domega_n = -dot(n, domega);
+ const Float domega_n = dot(domega, -n);
+ //const Float domega_n = dot(-n, domega);
+ //const Float domega_n = -dot(n, domega);
- // Normal component of the total displacement
- const Float omega_n = omega_t0_n + domega_n;
+ // Normal component of the total rotation
+ const Float omega_n = omega0_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 rotational increment
+ //const Float3 domega_t = domega - (-n) * dot(domega, -n);
+ const Float3 domega_t = domega - domega_n * (-n);
+ //const Float3 domega_t = domega - n * dot(n, domega);
- // Tangential component of the total displacement
- const Float3 omega_t = omega_t0 + domega_t;
+ // Tangential component of the total rotation
+ const Float3 omega_t = omega0_t + domega_t;
// Twisting torque: Visco-elastic contact model
- t_n = -devC_params.k_t * J * omega_n * n;
+ //const Float3 t_n = -devC_params.k_t * J * omega_n * n;
+ const Float3 t_n = -devC_params.k_t * J * omega_n * -n;
+ //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: Visco-elastic contact model
//t_t = -devC_params.k_n * I * omega_t;
+ //const Float3 t_t = devC_params.k_n * I * omega_t;
+ const Float3 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;
+ //t_t = devC_params.k_n * I * omega_t - devC_params.gamma_n * domega_t/devC_dt;
// Torque vector
- t = t_n + t_t;
-
+ //t = t_n + t_t;
+ //Float3 t_i = t_n + cross(-(x_i.w + overlap*0.5) * n, t_t);
+ //Float3 t_j = t_n + cross(-(x_j.w + overlap*0.5) * n, t_t);
+ //const Float3 t_i = t_n + cross(-(x_i.w + overlap*0.5) * n, f_t + t_t);
+ //const Float3 t_j = t_n + cross(-(x_j.w + overlap*0.5) * n, f_t - t_t);
+ const Float3 t_j = t_n + t_t;
+ //const Float3 t_i = t_n - t_t; //t_n - t_t;
+ //const Float3 t_j = t_n + t_t;
//// Save values
__syncthreads();
t@@ -186,8 +201,14 @@ __global__ void bondsLinear(
// 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);
+ //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);
+ //dev_torque[bond.x] += MAKE_FLOAT4(t_i.x, t_i.y, t_i.z, 0.0);
+ dev_torque[bond.x] -= MAKE_FLOAT4(t_j.x, t_j.y, t_j.z, 0.0);
+ dev_torque[bond.y] += MAKE_FLOAT4(t_j.x, t_j.y, t_j.z, 0.0);
+ //dev_torque[bond.y] += MAKE_FLOAT4(t_j.x, t_j.y, t_j.z, 0.0);
+ //dev_torque[bond.y] -= MAKE_FLOAT4(t_j.x, t_j.y, t_j.z, 0.0);
+ //dev_torque[bond.y] -= MAKE_FLOAT4(t.x, t.y, t.z, 0.0);
// make sure to remove write conflicts
}
(DIR) diff --git a/src/device.cu b/src/device.cu
t@@ -800,6 +800,13 @@ __host__ void DEM::startTime()
stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed, &t_integrateWalls);
checkForCudaErrors("Post integrateWalls");
+ /*for (int a=0; a<params.nb0; ++a)
+ std::cout << "bond " << a << ":\n"
+ << k.bonds_delta[a].x << ", "
+ << k.bonds_delta[a].y << ", "
+ << k.bonds_delta[a].z << ", "
+ << k.bonds_delta[a].w << std::endl;
+ break;*/
// Update timers and counters
time.current += time.dt;
(DIR) diff --git a/src/sphere.cpp b/src/sphere.cpp
t@@ -627,14 +627,17 @@ void DEM::forcechains(const std::string format, const int threedim,
cout << "set term epslatex size 8.6 cm, 5.6 cm\n";
cout << "set out 'plots/" << s << "-fc.tex'\n";
} else if (format == "epslatex-color") {
- cout << "set term epslatex color size 8.6 cm, 5.6 cm\n";
+ //cout << "set term epslatex color size 12 cm, 8.6 cm\n";
+ cout << "set term epsla color size 8.6 cm, 5.6 cm\n";
cout << "set out 'plots/" << s << "-fc.tex'\n";
}
cout << "set xlabel '\\sffamily $x_1$, [m]'\n";
if (threedim == 1) {
cout << "set ylabel '\\sffamily $x_2$, [m]'\n"
<< "set zlabel '\\sffamily $x_3$, [m]' offset 2\n";
+ //<< "set zlabel '\\sffamily $x_3$, [m]' offset 0\n";
} else
+ //cout << "set ylabel '\\sffamily $x_3$, [m]' offset 0\n";
cout << "set ylabel '\\sffamily $x_3$, [m]' offset 2\n";
cout << "set cblabel '\\sffamily $||\\boldsymbol{f}_n||$, [N]'\n"
(DIR) diff --git a/tests/bond_tests.py b/tests/bond_tests.py
t@@ -27,7 +27,7 @@ s2_1 = numpy.ones((2,1))*smallval
# Inter-particle distances to try (neg. for overlap)
distances = [0.2, 0.0, -0.2]
-#distances = [-0.2]
+#distances = [0.2]
for d in distances:
t@@ -39,7 +39,6 @@ for d in distances:
# setup particles, bond, and simulation
sb.x[0,:] = numpy.array((2.0, 2.0, 2.0))
- #sb.x[1,:] = numpy.array((2.2, 2.0, 2.0))
sb.x[1,:] = numpy.array((2.0+2.0*radii+d, 2.0, 2.0))
sb.radius = numpy.ones(sb.np)*radii
sb.initGridAndWorldsize(margin = 10, periodic = 1, contactmodel = 2, g = numpy.array([0.0, 0.0, 0.0]))
t@@ -51,6 +50,8 @@ for d in distances:
#visualize(sb.sid, "energy")
print("# Stability test")
+ sb.x[0,:] = numpy.array((2.0, 2.0, 2.0))
+ sb.x[1,:] = numpy.array((2.0+2.0*radii+d, 2.0, 2.0))
sb.zeroKinematics()
sb.initTemporal(total=0.2, file_dt=0.01)
#sb.initTemporal(total=0.01, file_dt=0.0001)
t@@ -67,9 +68,10 @@ for d in distances:
#printKinematics(sb)
print("# Normal expansion")
+ sb.x[0,:] = numpy.array((2.0, 2.0, 2.0))
+ sb.x[1,:] = numpy.array((2.0+2.0*radii+d, 2.0, 2.0))
sb.zeroKinematics()
- #sb.initTemporal(total=0.2, file_dt=0.01)
- sb.initTemporal(total=0.2, file_dt=0.001)
+ sb.initTemporal(total=0.2, file_dt=0.01)
sb.vel[1,0] = 1e-4
Ekinrot0 = sb.energy("kin") + sb.energy("rot")
sb.writebin(verbose=False)
t@@ -105,6 +107,8 @@ for d in distances:
#printKinematics(sb)
print("# Shear")
+ sb.x[0,:] = numpy.array((2.0, 2.0, 2.0))
+ sb.x[1,:] = numpy.array((2.0+2.0*radii+d, 2.0, 2.0))
sb.zeroKinematics()
sb.initTemporal(total=0.2, file_dt=0.01)
sb.vel[1,2] = 1e-4
t@@ -135,53 +139,54 @@ for d in distances:
#visualize(sb.sid, "energy")
- print("# Bend")
+ #'''
+ print("# Twist")
+ sb.x[0,:] = numpy.array((2.0, 2.0, 2.0))
+ sb.x[1,:] = numpy.array((2.0+2.0*radii+d, 2.0, 2.0))
sb.zeroKinematics()
sb.initTemporal(total=0.2, file_dt=0.01)
- sb.angvel[1,1] = 1e-4
+ #sb.initTemporal(total=0.001, file_dt=0.00001)
+ sb.angvel[1,0] = 1e-4
Ekinrot0 = sb.energy("kin") + sb.energy("rot")
sb.writebin(verbose=False)
sb.run(verbose=False)
sb.readlast(verbose=False)
compareFloats(Ekinrot0, sb.energy("kin") + sb.energy("rot") + sb.energy("bondpot"), "Energy cons.")
- print("vel[:,0]"),
- if (sb.vel[0,0] <= 0.0) or (sb.vel[1,0] >= 0.0):
- print(failed())
- else :
- print(passed())
- compareNumpyArrays(sb.vel[:,1], z2_1, "vel[:,1]")
- print("vel[:,2]"),
- if ((sb.vel[0,2] <= 0.0) or (sb.vel[1,2] >= 0.0)):
- print(failed())
- else :
- print(passed())
- compareNumpyArrays(sb.angvel[:,0:2:2], z2_2, "angvel[:,0:2:2]")
- print("angvel[:,1]"),
- if ((sb.angvel[0,1] <= 0.0) or (sb.angvel[1,1] <= 0.0)):
+ #compareNumpyArrays(sb.vel, z2_3, "vel\t")
+ print("angvel[:,0]"),
+ if ((sb.angvel[0,0] <= 0.0) or (sb.angvel[1,0] <= 0.0)):
+ raise Exception("Failed")
print(failed())
else :
print(passed())
+ compareNumpyArrays(sb.angvel[:,1:2], z2_2, "angvel[:,1:2]")
#printKinematics(sb)
#visualize(sb.sid, "energy")
+
-
- print("# Twist")
+ #'''
+ print("# Bend")
+ sb.x[0,:] = numpy.array((2.0, 2.0, 2.0))
+ sb.x[1,:] = numpy.array((2.0+2.0*radii+d, 2.0, 2.0))
sb.zeroKinematics()
sb.initTemporal(total=0.2, file_dt=0.01)
- sb.angvel[1,0] = 1e-4
+ sb.angvel[0,1] = -1e-4
+ sb.angvel[1,1] = 1e-4
Ekinrot0 = sb.energy("kin") + sb.energy("rot")
sb.writebin(verbose=False)
sb.run(verbose=False)
sb.readlast(verbose=False)
compareFloats(Ekinrot0, sb.energy("kin") + sb.energy("rot") + sb.energy("bondpot"), "Energy cons.")
- compareNumpyArrays(sb.vel, z2_3, "vel\t")
+ #compareNumpyArrays(sb.vel, z2_3, "vel\t")
+ #compareNumpyArrays(sb.angvel[:,0:2:2], z2_2, "angvel[:,0:2:2]")
print("angvel[:,1]"),
- if ((sb.angvel[0,1] <= 0.0) or (sb.angvel[0,1] <= 0.0)):
+ if ((sb.angvel[0,1] == 0.0) or (sb.angvel[1,1] == 0.0)):
+ raise Exception("Failed")
print(failed())
else :
print(passed())
- compareNumpyArrays(sb.angvel[:,0:2:2], z2_2, "angvel[:,0:2:2]")
- printKinematics(sb)
- visualize(sb.sid, "energy")
+ #printKinematics(sb)
+ #visualize(sb.sid, "energy")
+ #'''
-#cleanup(sb)
+cleanup(sb)
(DIR) diff --git a/tests/pytestutils.py b/tests/pytestutils.py
t@@ -2,11 +2,13 @@
from sphere import *
import subprocess
+import sys
def passed():
return "\tPassed"
def failed():
+ raise Exception("Failed")
return "\tFailed"
def compare(first, second, string):
t@@ -20,6 +22,16 @@ def compareFloats(first, second, string, criterion=1e-5):
print(string + passed())
else :
print(string + failed())
+ print("First: " + str(first))
+ print("Second: " + str(second))
+ print("Difference: " + str(second-first))
+
+def compareNumpyArrays(first, second, string):
+ if ((first == second).all()):
+ print(string + passed())
+ else :
+ print(string + failed())
+
def cleanup(spherebin):
'Remove temporary files'