tMerge branch 'varN' Sinusoidal modulation of N, choose w_devs_A=0.0 to disable - 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 158d3182d4b897921eb8909106d2ab118a8c9a0f
(DIR) parent 5dbc0a8544a6dfa40505639271fb7ca1a1834be5
(HTM) Author: Anders Damsgaard <adc@geo.au.dk>
Date: Sun, 17 Mar 2013 21:25:27 +0100
Merge branch 'varN'
Sinusoidal modulation of N, choose w_devs_A=0.0 to disable
Diffstat:
M python/sphere.py | 19 +++++++++++++++++--
M src/datatypes.h | 2 ++
M src/device.cu | 3 ++-
M src/file_io.cpp | 9 ++++++++-
M src/integration.cuh | 11 +++++++----
M tests/pytestutils.py | 2 +-
D tests/sphere.py | 2 --
7 files changed, 37 insertions(+), 11 deletions(-)
---
(DIR) diff --git a/python/sphere.py b/python/sphere.py
t@@ -98,6 +98,8 @@ class Spherebin:
self.w_vel = numpy.zeros(self.nw, dtype=numpy.float64)
self.w_force = numpy.zeros(self.nw, dtype=numpy.float64)
self.w_devs = numpy.zeros(self.nw, dtype=numpy.float64)
+ self.w_devs_A = numpy.zeros(1, dtype=numpy.float64)
+ self.w_devs_f = numpy.zeros(1, dtype=numpy.float64)
self.lambda_bar = numpy.ones(1, dtype=numpy.float64)
self.nb0 = numpy.zeros(1, dtype=numpy.uint32)
t@@ -161,6 +163,8 @@ class Spherebin:
(self.w_vel == other.w_vel).all() and\
(self.w_force == other.w_force).all() and\
(self.w_devs == other.w_devs).all() and\
+ self.w_devs_A == other.w_devs_A and\
+ self.w_devs_f == other.w_devs_f and\
self.gamma_wn == other.gamma_wn and\
self.gamma_wt == other.gamma_wt and\
self.lambda_bar == other.lambda_bar and\
t@@ -179,7 +183,7 @@ class Spherebin:
- def readbin(self, targetbin, verbose = True, bonds = True):
+ def readbin(self, targetbin, verbose = True, bonds = True, devsmod = True):
'Reads a target SPHERE binary file'
fh = None
t@@ -285,6 +289,9 @@ class Spherebin:
self.w_vel[i] = numpy.fromfile(fh, dtype=numpy.float64, count=1)
self.w_force[i] = numpy.fromfile(fh, dtype=numpy.float64, count=1)
self.w_devs[i] = numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ if (devsmod == True):
+ self.w_devs_A = numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ self.w_devs_f = numpy.fromfile(fh, dtype=numpy.float64, count=1)
if (bonds == True):
# Inter-particle bonds
t@@ -391,6 +398,8 @@ class Spherebin:
fh.write(self.w_vel[i].astype(numpy.float64))
fh.write(self.w_force[i].astype(numpy.float64))
fh.write(self.w_devs[i].astype(numpy.float64))
+ fh.write(self.w_devs_A.astype(numpy.float64))
+ fh.write(self.w_devs_f.astype(numpy.float64))
fh.write(self.lambda_bar.astype(numpy.float64))
fh.write(self.nb0.astype(numpy.uint32))
t@@ -1043,6 +1052,10 @@ class Spherebin:
# Increment the number of bonds with one
self.nb0 += 1
+ def currentDevs(self):
+ ''' Return current magnitude of the deviatoric normal stress '''
+ return w_devs[0] + w_devs_A * numpy.sin(2.0 * numpy.pi * self.time_current)
+
def energy(self, method):
""" Calculate the sum of the energy components of all particles.
t@@ -1180,8 +1193,10 @@ class Spherebin:
if (cudamemcheck == True):
cudamemchk = "cuda-memcheck "
+ status = subprocess.call("cd ..; " + valgrindbin + cudamemchk + "./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)
+ if (status != 0):
+ raise Exception("Error, the sphere run returned with status " + str(status))
def torqueScript(self,
(DIR) diff --git a/src/datatypes.h b/src/datatypes.h
t@@ -94,6 +94,8 @@ struct Params {
unsigned int nb0; // Number of inter-particle bonds at t=0
Float sigma_b; // Bond tensile strength
Float tau_b; // Bond shear strength
+ Float devs_A; // Amplitude of modulations in deviatoric normal stress
+ Float devs_f; // Frequency of modulations in deviatoric normal stress
};
// Structure containing wall parameters
(DIR) diff --git a/src/device.cu b/src/device.cu
t@@ -801,7 +801,8 @@ __host__ void DEM::startTime()
dev_walls_wmode,
dev_walls_force_partial,
dev_walls_vel0,
- blocksPerGrid);
+ blocksPerGrid,
+ time.current);
}
// Synchronization point
(DIR) diff --git a/src/file_io.cpp b/src/file_io.cpp
t@@ -204,13 +204,18 @@ void DEM::readbin(const char *target)
ifs.read(as_bytes(walls.mvfd[i].z), sizeof(Float));
ifs.read(as_bytes(walls.mvfd[i].w), sizeof(Float));
}
+ ifs.read(as_bytes(params.devs_A), sizeof(params.devs_A));
+ ifs.read(as_bytes(params.devs_f), sizeof(params.devs_f));
+
// Read bond parameters
ifs.read(as_bytes(params.lambda_bar), sizeof(params.lambda_bar));
ifs.read(as_bytes(params.nb0), sizeof(params.nb0));
ifs.read(as_bytes(params.sigma_b), sizeof(params.sigma_b));
ifs.read(as_bytes(params.tau_b), sizeof(params.tau_b));
- k.bonds = new uint2[params.nb0];
+ std::cout << "\nparams.nb0 = " << params.nb0 << std::endl;
+ if (params.nb0 > 0)
+ 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) {
t@@ -363,6 +368,8 @@ void DEM::writebin(const char *target)
ofs.write(as_bytes(walls.mvfd[i].z), sizeof(Float));
ofs.write(as_bytes(walls.mvfd[i].w), sizeof(Float));
}
+ ofs.write(as_bytes(params.devs_A), sizeof(params.devs_A));
+ ofs.write(as_bytes(params.devs_f), sizeof(params.devs_f));
// Write bond parameters
ofs.write(as_bytes(params.lambda_bar), sizeof(params.lambda_bar));
(DIR) diff --git a/src/integration.cuh b/src/integration.cuh
t@@ -221,12 +221,14 @@ __global__ void summation(Float* in, Float *out)
}
// Update wall positions
-__global__ void integrateWalls(Float4* dev_walls_nx,
+__global__ void integrateWalls(
+ Float4* dev_walls_nx,
Float4* dev_walls_mvfd,
int* dev_walls_wmode,
Float* dev_walls_force_partial,
Float* dev_walls_vel0,
- unsigned int blocksPerGrid)
+ unsigned int blocksPerGrid,
+ Float t_current)
{
unsigned int idx = threadIdx.x + blockIdx.x * blockDim.x; // Thread id
t@@ -253,10 +255,11 @@ __global__ void integrateWalls(Float4* dev_walls_nx,
// Normal load = Deviatoric stress times wall surface area,
// directed downwards.
- Float N = -w_mvfd.w*devC_grid.L[0]*devC_grid.L[1];
+ Float sigma_0 = w_mvfd.w + devC_params.devs_A * sin(2.0 * 3.141596654 * devC_params.devs_f * t_current);
+ Float N = -sigma_0*devC_grid.L[0]*devC_grid.L[1];
// Calculate resulting acceleration of wall
- // (Wall mass is stored in w component of position Float4)
+ // (Wall mass is stored in x component of position Float4)
acc = (w_mvfd.z + N)/w_mvfd.x;
// If Wall BC is controlled by velocity, it should not change
(DIR) diff --git a/tests/pytestutils.py b/tests/pytestutils.py
t@@ -1,6 +1,6 @@
#!/usr/bin/env python
-from sphere import *
+from ../python/sphere import *
import subprocess
import sys
(DIR) diff --git a/tests/sphere.py b/tests/sphere.py
t@@ -1 +0,0 @@
-../python/sphere.py
-\ No newline at end of file