tParameters k_s and gamma_s renamed to k_t and gamma_t to refer to geometry - 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 b318e0cd9d99bc8a3b15e0372bb989c6ca917538
(DIR) parent f5c826ebd0189733df3104e874d91082aeeb97ef
(HTM) Author: Anders Damsgaard <adc@geo.au.dk>
Date: Mon, 27 Aug 2012 16:51:45 +0200
Parameters k_s and gamma_s renamed to k_t and gamma_t to refer to geometry
Diffstat:
M python/sphere.py | 30 +++++++++++++++---------------
M src/constants.cuh | 4 ++--
M src/contactmodels.cuh | 10 +++++-----
M src/datatypes.h | 8 ++++----
M src/device.cu | 8 ++++----
M src/file_io.cpp | 8 ++++----
M src/main.cpp | 12 ++++++------
7 files changed, 40 insertions(+), 40 deletions(-)
---
(DIR) diff --git a/python/sphere.py b/python/sphere.py
t@@ -44,10 +44,10 @@ class Spherebin:
self.radius = numpy.zeros(self.np, dtype=numpy.float64)
self.rho = numpy.zeros(self.np, dtype=numpy.float64)
self.k_n = numpy.zeros(self.np, dtype=numpy.float64)
- self.k_s = numpy.zeros(self.np, dtype=numpy.float64)
+ self.k_t = numpy.zeros(self.np, dtype=numpy.float64)
self.k_r = numpy.zeros(self.np, dtype=numpy.float64)
self.gamma_n = numpy.zeros(self.np, dtype=numpy.float64)
- self.gamma_s = numpy.zeros(self.np, dtype=numpy.float64)
+ self.gamma_t = numpy.zeros(self.np, dtype=numpy.float64)
self.gamma_r = numpy.zeros(self.np, dtype=numpy.float64)
self.mu_s = numpy.zeros(self.np, dtype=numpy.float64)
self.mu_d = numpy.zeros(self.np, dtype=numpy.float64)
t@@ -114,10 +114,10 @@ class Spherebin:
self.radius = numpy.zeros(self.np, dtype=numpy.float64)
self.rho = numpy.zeros(self.np, dtype=numpy.float64)
self.k_n = numpy.zeros(self.np, dtype=numpy.float64)
- self.k_s = numpy.zeros(self.np, dtype=numpy.float64)
+ self.k_t = numpy.zeros(self.np, dtype=numpy.float64)
self.k_r = numpy.zeros(self.np, dtype=numpy.float64)
self.gamma_n = numpy.zeros(self.np, dtype=numpy.float64)
- self.gamma_s = numpy.zeros(self.np, dtype=numpy.float64)
+ self.gamma_t = numpy.zeros(self.np, dtype=numpy.float64)
self.gamma_r = numpy.zeros(self.np, dtype=numpy.float64)
self.mu_s = numpy.zeros(self.np, dtype=numpy.float64)
self.mu_d = numpy.zeros(self.np, dtype=numpy.float64)
t@@ -148,10 +148,10 @@ class Spherebin:
self.radius[j] = numpy.fromfile(fh, dtype=numpy.float64, count=1)
self.rho[j] = numpy.fromfile(fh, dtype=numpy.float64, count=1)
self.k_n[j] = numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.k_s[j] = numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ self.k_t[j] = numpy.fromfile(fh, dtype=numpy.float64, count=1)
self.k_r[j] = numpy.fromfile(fh, dtype=numpy.float64, count=1)
self.gamma_n[j] = numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.gamma_s[j] = numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ self.gamma_t[j] = numpy.fromfile(fh, dtype=numpy.float64, count=1)
self.gamma_r[j] = numpy.fromfile(fh, dtype=numpy.float64, count=1)
self.mu_s[j] = numpy.fromfile(fh, dtype=numpy.float64, count=1)
self.mu_d[j] = numpy.fromfile(fh, dtype=numpy.float64, count=1)
t@@ -247,10 +247,10 @@ class Spherebin:
fh.write(self.radius[j].astype(numpy.float64))
fh.write(self.rho[j].astype(numpy.float64))
fh.write(self.k_n[j].astype(numpy.float64))
- fh.write(self.k_s[j].astype(numpy.float64))
+ fh.write(self.k_t[j].astype(numpy.float64))
fh.write(self.k_r[j].astype(numpy.float64))
fh.write(self.gamma_n[j].astype(numpy.float64))
- fh.write(self.gamma_s[j].astype(numpy.float64))
+ fh.write(self.gamma_t[j].astype(numpy.float64))
fh.write(self.gamma_r[j].astype(numpy.float64))
fh.write(self.mu_s[j].astype(numpy.float64))
fh.write(self.mu_d[j].astype(numpy.float64))
t@@ -614,8 +614,8 @@ class Spherebin:
r_min = numpy.amin(self.radius)
# Computational time step (O'Sullivan et al, 2003)
- #self.time_dt[0] = 0.17 * math.sqrt((4.0/3.0 * math.pi * r_min**3 * self.rho[0]) / numpy.amin([self.k_n[:], self.k_s[:]]) )
- self.time_dt[0] = 0.12 * math.sqrt((4.0/3.0 * math.pi * r_min**3 * self.rho[0]) / numpy.amin([self.k_n[:], self.k_s[:]]) )
+ #self.time_dt[0] = 0.17 * math.sqrt((4.0/3.0 * math.pi * r_min**3 * self.rho[0]) / numpy.amax([self.k_n[:], self.k_t[:]]) )
+ self.time_dt[0] = 0.12 * math.sqrt((4.0/3.0 * math.pi * r_min**3 * self.rho[0]) / numpy.amax([self.k_n[:], self.k_t[:]]) )
# Time at start
self.time_current[0] = current
t@@ -628,10 +628,10 @@ class Spherebin:
ang_r = 0,
rho = 2660,
k_n = 1e9,
- k_s = 1e9,
+ k_t = 1e9,
k_r = 4e6,
gamma_n = 1e3,
- gamma_s = 1e3,
+ gamma_t = 1e3,
gamma_r = 2e3,
gamma_wn = 1e3,
gamma_ws = 1e3,
t@@ -650,7 +650,7 @@ class Spherebin:
self.k_n = numpy.ones(self.np, dtype=numpy.float64) * k_n
# Contact shear elastic stiffness (for shearmodel = 2), N/m
- self.k_s = numpy.ones(self.np, dtype=numpy.float64) * k_s
+ self.k_t = numpy.ones(self.np, dtype=numpy.float64) * k_t
# Contact rolling elastic stiffness (for shearmodel = 2), N/m
self.k_r = numpy.ones(self.np, dtype=numpy.float64) * k_r
t@@ -663,7 +663,7 @@ class Spherebin:
self.gamma_n = numpy.ones(self.np, dtype=numpy.float64) * gamma_n
# Contact shear viscosity, Ns/m
- self.gamma_s = numpy.ones(self.np, dtype=numpy.float64) * gamma_s
+ self.gamma_t = numpy.ones(self.np, dtype=numpy.float64) * gamma_t
# Contact rolling viscosity, Ns/m?
self.gamma_r = numpy.ones(self.np, dtype=numpy.float64) * gamma_r
t@@ -692,7 +692,7 @@ class Spherebin:
# Wettability, 0=perfect
theta = 0.0;
if (capillaryCohesion == 1):
- self.kappa[0] = 2.0 * math.pi * gamma_s * numpy.cos(theta) # Prefactor
+ self.kappa[0] = 2.0 * math.pi * gamma_t * numpy.cos(theta) # Prefactor
self.V_b[0] = 1e-12 # Liquid volume at bond
else :
self.kappa[0] = 0.0; # Zero capillary force
(DIR) diff --git a/src/constants.cuh b/src/constants.cuh
t@@ -12,10 +12,10 @@ __constant__ unsigned int devC_np; // Number of particles
__constant__ char devC_nc; // Max. number of contacts a particle can have
__constant__ unsigned int devC_shearmodel; // Shear force model: 1: viscous, frictional, 2: elastic, frictional
__constant__ Float devC_k_n; // Material normal stiffness
-__constant__ Float devC_k_s; // Material shear stiffness
+__constant__ Float devC_k_t; // Material tangential stiffness
__constant__ Float devC_k_r; // Material rolling stiffness
__constant__ Float devC_gamma_n; // Material normal viscosity
-__constant__ Float devC_gamma_s; // Material shear viscosity
+__constant__ Float devC_gamma_t; // Material tangential viscosity
__constant__ Float devC_gamma_r; // Material rolling viscosity
__constant__ Float devC_gamma_wn; // Wall normal viscosity
__constant__ Float devC_gamma_ws; // Wall shear viscosity
(DIR) diff --git a/src/contactmodels.cuh b/src/contactmodels.cuh
t@@ -199,7 +199,7 @@ __device__ void contactLinearViscous(Float3* N, Float3* T, Float* es_dot, Float*
if (vel_t_ab_length > 0.f) {
// Tangential force by viscous model
- Float f_t_visc = devC_gamma_s * vel_t_ab_length;
+ Float f_t_visc = devC_gamma_t * vel_t_ab_length;
// Determine max. friction
Float f_t_limit;
t@@ -340,8 +340,8 @@ __device__ void contactLinear(Float3* N, Float3* T,
if (delta_t0_length > 0.f || vel_t_ab_length > 0.f) {
// Shear force: Visco-Elastic, limited by Coulomb friction
- Float3 f_t_elast = -1.0f * devC_k_s * delta_t0;
- Float3 f_t_visc = -1.0f * devC_gamma_s * vel_t_ab;
+ Float3 f_t_elast = -1.0f * devC_k_t * delta_t0;
+ Float3 f_t_visc = -1.0f * devC_gamma_t * vel_t_ab;
Float f_t_limit;
t@@ -369,12 +369,12 @@ __device__ void contactLinear(Float3* N, Float3* T,
// In a slip event, the tangential spring is adjusted to a
// length which is consistent with Coulomb's equation
// (Hinrichsen and Wolf, 2004)
- delta_t = -1.0f/devC_k_s * f_t + devC_gamma_s * vel_t_ab;
+ delta_t = -1.0f/devC_k_t * f_t + devC_gamma_t * vel_t_ab;
// Shear friction heat production rate:
// The energy lost from the tangential spring is dissipated as heat
//*es_dot += -dot(vel_t_ab, f_t);
- *es_dot += length(delta_t0 - delta_t) * devC_k_s / devC_dt; // Seen in EsyS-Particle
+ *es_dot += length(delta_t0 - delta_t) * devC_k_t / devC_dt; // Seen in EsyS-Particle
} else { // Static case
(DIR) diff --git a/src/datatypes.h b/src/datatypes.h
t@@ -68,10 +68,10 @@ const char NC = 32;
struct Particles {
Float *radius;
Float *k_n;
- Float *k_s;
+ Float *k_t;
Float *k_r;
Float *gamma_n;
- Float *gamma_s;
+ Float *gamma_t;
Float *gamma_r;
Float *mu_s;
Float *mu_d;
t@@ -111,10 +111,10 @@ struct Params {
unsigned int nw;
Float dt;
Float k_n;
- Float k_s;
+ Float k_t;
Float k_r;
Float gamma_n;
- Float gamma_s;
+ Float gamma_t;
Float gamma_r;
Float gamma_wn;
Float gamma_ws;
(DIR) diff --git a/src/device.cu b/src/device.cu
t@@ -110,20 +110,20 @@ __host__ void transferToConstantMemory(Particles* p,
// copy the values from the first particle into the designated constant memory.
//printf("(params.global == %d) ", params.global);
params->k_n = p->k_n[0];
- params->k_s = p->k_s[0];
+ params->k_t = p->k_t[0];
params->k_r = p->k_r[0];
params->gamma_n = p->gamma_n[0];
- params->gamma_s = p->gamma_s[0];
+ params->gamma_t = p->gamma_t[0];
params->gamma_r = p->gamma_r[0];
params->mu_s = p->mu_s[0];
params->mu_d = p->mu_d[0];
params->mu_r = p->mu_r[0];
params->rho = p->rho[0];
cudaMemcpyToSymbol("devC_k_n", ¶ms->k_n, sizeof(Float));
- cudaMemcpyToSymbol("devC_k_s", ¶ms->k_s, sizeof(Float));
+ cudaMemcpyToSymbol("devC_k_t", ¶ms->k_t, sizeof(Float));
cudaMemcpyToSymbol("devC_k_r", ¶ms->k_r, sizeof(Float));
cudaMemcpyToSymbol("devC_gamma_n", ¶ms->gamma_n, sizeof(Float));
- cudaMemcpyToSymbol("devC_gamma_s", ¶ms->gamma_s, sizeof(Float));
+ cudaMemcpyToSymbol("devC_gamma_t", ¶ms->gamma_t, sizeof(Float));
cudaMemcpyToSymbol("devC_gamma_r", ¶ms->gamma_r, sizeof(Float));
cudaMemcpyToSymbol("devC_gamma_wn", ¶ms->gamma_wn, sizeof(Float));
cudaMemcpyToSymbol("devC_gamma_ws", ¶ms->gamma_ws, sizeof(Float));
(DIR) diff --git a/src/file_io.cpp b/src/file_io.cpp
t@@ -96,10 +96,10 @@ int fwritebin(char *target,
fwrite(&p->radius[j], sizeof(p->radius[j]), 1, fp);
fwrite(&p->rho[j], sizeof(p->rho[j]), 1, fp);
fwrite(&p->k_n[j], sizeof(p->k_n[j]), 1, fp);
- fwrite(&p->k_s[j], sizeof(p->k_s[j]), 1, fp);
+ fwrite(&p->k_t[j], sizeof(p->k_t[j]), 1, fp);
fwrite(&p->k_r[j], sizeof(p->k_r[j]), 1, fp);
fwrite(&p->gamma_n[j], sizeof(p->gamma_n[j]), 1, fp);
- fwrite(&p->gamma_s[j], sizeof(p->gamma_s[j]), 1, fp);
+ fwrite(&p->gamma_t[j], sizeof(p->gamma_t[j]), 1, fp);
fwrite(&p->gamma_r[j], sizeof(p->gamma_r[j]), 1, fp);
fwrite(&p->mu_s[j], sizeof(p->mu_s[j]), 1, fp);
fwrite(&p->mu_d[j], sizeof(p->mu_d[j]), 1, fp);
t@@ -235,13 +235,13 @@ int fwritebin(char *target,
fwrite(&d, sizeof(d), 1, fp);
d = (double)p->k_n[j];
fwrite(&d, sizeof(d), 1, fp);
- d = (double)p->k_s[j];
+ d = (double)p->k_t[j];
fwrite(&d, sizeof(d), 1, fp);
d = (double)p->k_r[j];
fwrite(&d, sizeof(d), 1, fp);
d = (double)p->gamma_n[j];
fwrite(&d, sizeof(d), 1, fp);
- d = (double)p->gamma_s[j];
+ d = (double)p->gamma_t[j];
fwrite(&d, sizeof(d), 1, fp);
d = (double)p->gamma_r[j];
fwrite(&d, sizeof(d), 1, fp);
(DIR) diff --git a/src/main.cpp b/src/main.cpp
t@@ -167,10 +167,10 @@ int main(int argc, char *argv[])
p.m = new Float[p.np]; // Particle masses
p.I = new Float[p.np]; // Particle moment of inertia
p.k_n = new Float[p.np]; // Particle normal stiffnesses
- p.k_s = new Float[p.np]; // Particle shear stiffnesses
+ p.k_t = new Float[p.np]; // Particle shear stiffnesses
p.k_r = new Float[p.np]; // Particle rolling stiffnesses
p.gamma_n = new Float[p.np]; // Particle normal viscosity
- p.gamma_s = new Float[p.np]; // Particle shear viscosity
+ p.gamma_t = new Float[p.np]; // Particle shear viscosity
p.gamma_r = new Float[p.np]; // Particle rolling viscosity
p.mu_s = new Float[p.np]; // Inter-particle static shear contact friction coefficients
p.mu_d = new Float[p.np]; // Inter-particle dynamic shear contact friction coefficients
t@@ -254,13 +254,13 @@ int main(int argc, char *argv[])
exit(1);
if (fread(&p.k_n[j], sizeof(p.k_n[j]), 1, fp) != 1)
exit(1);
- if (fread(&p.k_s[j], sizeof(p.k_s[j]), 1, fp) != 1)
+ if (fread(&p.k_t[j], sizeof(p.k_t[j]), 1, fp) != 1)
exit(1);
if (fread(&p.k_r[j], sizeof(p.k_r[j]), 1, fp) != 1)
exit(1);
if (fread(&p.gamma_n[j], sizeof(p.gamma_n[j]), 1, fp) != 1)
exit(1);
- if (fread(&p.gamma_s[j], sizeof(p.gamma_s[j]), 1, fp) != 1)
+ if (fread(&p.gamma_t[j], sizeof(p.gamma_t[j]), 1, fp) != 1)
exit(1);
if (fread(&p.gamma_r[j], sizeof(p.gamma_r[j]), 1, fp) != 1)
exit(1);
t@@ -445,10 +445,10 @@ int main(int argc, char *argv[])
// Particle single-value parameters
delete[] p.radius;
delete[] p.k_n;
- delete[] p.k_s;
+ delete[] p.k_t;
delete[] p.k_r;
delete[] p.gamma_n;
- delete[] p.gamma_s;
+ delete[] p.gamma_t;
delete[] p.gamma_r;
delete[] p.mu_s;
delete[] p.mu_d;