tRaytracer now working - 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 e956c0c522a59ce13be1db5c1fdfd3ff27279355
(DIR) parent 513ff2b8bafb8119bce86cf56b4d80d6b1e0519b
(HTM) Author: Anders Damsgaard <adc@geo.au.dk>
Date: Tue, 6 Nov 2012 14:26:53 +0100
Raytracer now working
Diffstat:
D img_out/about.txt | 2 --
D input/1e3-test-shear.bin | 0
D output/about.txt | 1 -
D python/1e4-largesize-uniaxial.py | 136 -------------------------------
D python/elastic_collision.py | 46 -------------------------------
D python/inelastic_collision.py | 49 -------------------------------
M python/sphere.py | 129 +++++++++----------------------
M python/tests.py | 6 +++++-
M src/device.cu | 21 ++++++---------------
M src/file_io.cpp | 7 +++++--
M src/sphere.cpp | 17 ++++++++++++++---
11 files changed, 68 insertions(+), 346 deletions(-)
---
(DIR) diff --git a/img_out/about.txt b/img_out/about.txt
t@@ -1,2 +0,0 @@
-This folder will contain images rendered from ../output binaries, and
-graphs plotted with ../python/sphere.py.
(DIR) diff --git a/input/1e3-test-shear.bin b/input/1e3-test-shear.bin
Binary files differ.
(DIR) diff --git a/output/about.txt b/output/about.txt
t@@ -1 +0,0 @@
-This folder will contain output binary files.
(DIR) diff --git a/python/1e4-largesize-uniaxial.py b/python/1e4-largesize-uniaxial.py
t@@ -1,136 +0,0 @@
-#!/usr/bin/env python
-
-# Import sphere functionality
-from sphere import *
-import subprocess
-
-# Determine system hostname
-import os;
-hostname = os.uname()[1]
-
-# Simulation ID
-sid = "1e4-largesize-uniaxial"
-
-# Delete previous output
-subprocess.call("rm ../{img_out,output}/" + sid + "*", shell=True);
-
-
-#### GRAVITATIONAL CONSOLIDATION
-
-# New class
-init = Spherebin(np = 1e4)
-
-# Generate random radii
-init.generateRadii(psd = 'logn', histogram = 1, radius_mean = 0.02)
-
-# Initialize particle parameters
-# Values from Yohannes et al. 2012
-init.defaultParams(ang_s = 28,
- ang_d = 28,
- ang_r = 0,
- rho = 2600,
- k_n = 1.16e9,
- k_t = 1.16e9,
- gamma_n = 0.0,
- gamma_t = 0.0,
- gamma_r = 0.0)
-
-# Place particles in grid-like arrangement
-# periodic: 0 = frictional x- and y-walls
-# periodic: 1 = periodic x- and y-boundaries
-# periodic: 2 = periodic x boundaries, frictional y-walls
-# shearmodel: 1 = viscous, frictional shear force model
-# shearmodel: 2 = elastic, frictional shear force model
-init.initRandomGridPos(gridnum = numpy.array([20,10,600]), periodic = 1, shearmodel = 2)
-
-# Initialize temporal parameters
-init.initTemporal(total = 5.0, file_dt = 0.10)
-
-# Write output binary for sphere
-init.writebin("../input/" + sid + "-initgrid.bin".format(sid))
-
-# Render start configuration
-render("../input/" + sid + "-initgrid.bin", out = sid + "-initgrid")
-
-# Run simulation
-subprocess.call("cd ..; ./sphere_*_X86_64 " + sid + "-initgrid", shell=True)
-
-# Plot energy
-visualize(sid + "-initgrid", "energy", savefig=True, outformat='png')
-
-
-#### CONSOLIDATION UNDER DEVIATORIC STRESS
-
-# New class
-cons = Spherebin(np = init.np, nd = init.nd)
-
-# Find out which output file was the last generated during gravitational consolidation
-lastfile = status(sid + "-initgrid", verbose = False)
-filename = "../output" + sid + "-initgrid.output{0}.bin".format(lastfile)
-
-
-## Uniaxial compression loop:
-# 1. Read last experiment binary
-# 2. Set new devs and zero current time
-# 3. Run sphere
-# 4. Visualize wall data and find void ratio (e)
-# 5. Save name of last binary written
-
-# Define deviatoric stresses to test
-stresses = numpy.array([10e3, 20e3, 40e3])
-voidratio = numpy.zeros(1, length(stresses))
-
-i = 0
-for devs in stresses:
-
- # Simulation ID for consolidation run
- cid = sid + "-uniaxial-{0}Pa".format(devs)
-
- # Read the last output file of from the gravitational consolidation experiment
- cons.readbin(filename)
-
- # Setup consolidation experiment
- cons.consolidate(deviatoric_stress = devs, periodic = init.periodic)
-
- # Zero time variables and set new total time
- cons.initTemporal(total = 3.0, file_dt = 0.05)
-
- # Write output binary for sphere
- cons.writebin("../input/" + cid + ".bin")
-
- # Run simulation
- subprocess.call("cd ..; ./sphere_*_X86_64 " + cid, shell=True)
-
- # Plot energy and wall data
- visualize(cid, "energy", savefig=True, outformat='png')
- visualize(cid, "walls", savefig=True, outformat='png')
-
- # Find void ratio at the end of the compressional period
- lastfile = status(cid, verbose = False)
- cons.readbin(cid)
- voidratio[i] = voidRatio()
- i = i+1
-
- # Name of last output file
- filename = "../output" + cid + ".output{0}.bin".format(lastfile)
-
-
-# Save plot of measured compression curve
-fig = plt.figure(figsize(15,10),dpi=300)
-figtitle = "{0}, measured compression curve".format(sid)
-fig.text(0.5,0.95,figtitle,horizontalalignment='center',fontproperties=FontProperties(size=18))
-ax1 = plt.subplot2grid((1,1),(0,0))
-ax1.set_xlabel('Stress [kPa]')
-ax1.set_ylabel('Void ratio')
-ax1.semilogx(stresses, voidratio, '+-')
-
-# Save values to text file
-fh = None
-try:
- fh = open(sid + "-uniaxial.txt", "w")
- for i in range(length(stresses)):
- fh.write("{0}\t{1}\n".format(stresses[i], voidratio[i]))
-finally:
- if fh is not None:
- fh.close()
-
(DIR) diff --git a/python/elastic_collision.py b/python/elastic_collision.py
t@@ -1,46 +0,0 @@
-#!/usr/bin/env python
-import subprocess
-import numpy
-import matplotlib.pyplot as plt
-
-# Import sphere functionality
-from sphere import *
-
-# New class
-init = Spherebin(np = 2, nd = 3)
-
-#init.generateRadii(radius_mean = 0.5, radius_variance = 1e-15, histogram = 0)
-init.radius[0] = numpy.ones(1, dtype=numpy.float64) * 0.5;
-init.radius[1] = numpy.ones(1, dtype=numpy.float64) * 0.52;
-
-init.defaultparams()
-
-# The world should never be less than 3 cells in ANY direction, due to contact searcg algorithm
-init.initsetup(gridnum = numpy.array([12,4,4]), periodic = 1, shearmodel = 2, g = numpy.array([0.0, 0.0, 0.0]))
-
-init.x[0] = numpy.array([1, 2, 2])
-init.x[1] = numpy.array([6, 2, 2])
-#init.x[2] = numpy.array([7, 2, 2])
-#init.x[3] = numpy.array([8, 2, 2])
-
-init.nu = numpy.zeros(init.np, dtype=numpy.float64)
-
-#for i in range(init.np):
-# init.x[i] = numpy.array([4+i*init.radius[i]*2, 2, 2])
-
-init.vel[0] = numpy.array([10.0, 0.0, 0.0])
-
-
-init.initTemporal(total = 1.0)
-
-init.writebin("../input/nc-test.bin")
-#render("../input/nc-test.bin", out = "~/Desktop/nc-test")
-
-subprocess.call("cd ..; rm output/*; ./sphere_darwin_X86_64 nc-test", shell=True)
-
-visualize("nc-test", "energy", savefig=True, outformat='png')
-#visualize("nc-test", "walls", savefig=True)
-
-subprocess.call("rm ../img_out/*; cd ../raytracer; ./render_all_outputs_GPU_clever.sh", shell=True)
-
-
(DIR) diff --git a/python/inelastic_collision.py b/python/inelastic_collision.py
t@@ -1,49 +0,0 @@
-#!/usr/bin/env python
-import subprocess
-import numpy
-import matplotlib.pyplot as plt
-
-# Import sphere functionality
-from sphere import *
-
-# New class
-init = Spherebin(np = 2, nd = 3)
-
-init.radius = numpy.ones(init.np, dtype=numpy.float64) * 0.5;
-
-init.defaultparams()
-
-# The world should never be less than 3 cells in ANY direction, due to contact searcg algorithm
-init.initsetup(gridnum = numpy.array([12,4,4]), periodic = 1, shearmodel = 2, g = numpy.array([0.0, 0.0, 0.0]))
-
-init.x[0] = numpy.array([1, 2, 2])
-init.x[1] = numpy.array([6, 2, 2])
-#init.x[2] = numpy.array([7, 2, 2])
-#init.x[3] = numpy.array([8, 2, 2])
-
-# Set fraction of critical damping (0 = elastic, 1 = completely inelastic)
-damping_fraction = 1.0
-init.nu = numpy.ones(init.np, dtype=numpy.float64) \
- * damping_fraction * 2.0 * math.sqrt(4.0/3.0 * math.pi * init.radius.min()**3 \
- * init.rho[0] * init.k_n[0])
-
-
-#for i in range(init.np):
-# init.x[i] = numpy.array([4+i*init.radius[i]*2, 2, 2])
-
-init.vel[0] = numpy.array([10.0, 0.0, 0.0])
-
-
-init.initTemporal(total = 1.0)
-
-init.writebin("../input/nc-test.bin")
-#render("../input/nc-test.bin", out = "~/Desktop/nc-test")
-
-subprocess.call("cd ..; rm output/*; ./sphere_darwin_X86_64 nc-test", shell=True)
-
-visualize("nc-test", "energy", savefig=True, outformat='png')
-#visualize("nc-test", "walls", savefig=True)
-
-subprocess.call("rm ../img_out/*; cd ../raytracer; ./render_all_outputs_GPU_clever.sh", shell=True)
-
-
(DIR) diff --git a/python/sphere.py b/python/sphere.py
t@@ -356,7 +356,7 @@ class Spherebin:
def generateRadii(self, psd = 'logn',
radius_mean = 440e-6,
radius_variance = 8.8e-9,
- histogram = 1):
+ histogram = True):
""" Draw random particle radii from the selected probability distribution.
Specify mean radius and variance. The variance should be kept at a
very low value.
t@@ -372,7 +372,7 @@ class Spherebin:
self.radius = numpy.random.uniform(radius_min, radius_max, self.np)
# Show radii as histogram
- if histogram == 1:
+ if (histogram == True):
fig = plt.figure(figsize=(15,10), dpi=300)
figtitle = 'Particle size distribution, {0} particles'.format(self.np[0])
fig.text(0.5,0.95,figtitle,horizontalalignment='center',fontproperties=FontProperties(size=18))
t@@ -432,17 +432,6 @@ class Spherebin:
print " "
self.contactmodel[0] = contactmodel
- # Initialize upper wall
- self.wmode[0] = 0 # 0: fixed, 1: devs, 2: vel
- self.w_n[0,2] = -1.0
- self.w_x[0] = self.L[2]
- self.w_m[0] = self.rho[0] * self.np * math.pi * r_max**3
- self.w_vel[0] = 0.0
- self.w_force[0] = 0.0
- self.w_devs[0] = 0.0
- #self.nw[0] = numpy.ones(1, dtype=numpy.uint32) * 1
- self.nw = numpy.ones(1, dtype=numpy.uint32) * 1
-
# Generate grid based on particle positions
def initGrid(self):
t@@ -566,16 +555,6 @@ class Spherebin:
self.num[0] += 1
self.num[1] += 1
- # Initialize upper wall
- self.wmode[0] = 0
- self.w_n[0,2] = -1.0
- self.w_x[0] = self.L[2]
- self.w_m[0] = self.rho[0] * self.np * math.pi * r_max**3
- self.w_vel[0] = 0.0
- self.w_force[0] = 0.0
- self.w_devs[0] = 0.0
- self.nw = numpy.ones(1, dtype=numpy.uint32) * 1
-
# Initialize particle positions to non-overlapping configuration
# in grid, with a certain element of randomness
t@@ -631,16 +610,6 @@ class Spherebin:
self.num[2] = numpy.ceil(z_max/cellsize)
self.L = self.num * cellsize
- # Initialize upper wall
- if (self.nw > 0):
- self.wmode[0] = 0
- self.w_n[0,2] = -1.0
- self.w_x[0] = self.L[2]
- self.w_m[0] = self.rho[0] * self.np * math.pi * r_max**3
- self.w_vel[0] = 0.0
- self.w_force[0] = 0.0
- self.w_devs[0] = 0.0
- self.nw = numpy.ones(1, dtype=numpy.uint32) * 1
# Adjust grid and upper wall for consolidation under deviatoric stress
def consolidate(self, deviatoric_stress = 10e3,
t@@ -663,14 +632,15 @@ class Spherebin:
self.angpos = numpy.zeros(self.np*self.nd, dtype=numpy.float64).reshape(self.np, self.nd)
# Initialize upper wall
- self.wmode[0] = 1 # devs
+ self.nw = numpy.array([1], dtype=numpy.uint32)
+ self.wmode = numpy.array([1]) # devs 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[0] = self.L[2]
- self.w_m[0] = self.rho[0] * self.np * math.pi * (cellsize/2.0)**3
- self.w_vel[0] = 0.0
- self.w_force[0] = 0.0
- self.w_devs[0] = deviatoric_stress
- self.nw = numpy.ones(1, dtype=numpy.uint32) * 1
+ self.w_x = numpy.array([self.L[2]])
+ self.w_m = numpy.array([self.rho[0] * self.np * math.pi * (cellsize/2.0)**3])
+ self.w_vel = numpy.array([0.0])
+ self.w_force = numpy.array([0.0])
+ self.w_devs = numpy.array([deviatoric_stress])
# Adjust grid and upper wall for consolidation under fixed upper wall velocity
def uniaxialStrainRate(self, wvel = -0.001,
t@@ -693,13 +663,15 @@ class Spherebin:
self.angpos = numpy.zeros(self.np*self.nd, dtype=numpy.float64).reshape(self.np, self.nd)
# Initialize upper wall
- self.wmode[0] = 2 # strain rate controlled
+ self.nw = numpy.array([1], dtype=numpy.uint32)
+ self.wmode = numpy.array([2]) # strain rate 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[0] = self.L[2]
- self.w_m[0] = self.rho[0] * self.np * math.pi * (cellsize/2.0)**3
- self.w_vel[0] = wvel
- self.w_force[0] = 0.0
- self.nw = numpy.ones(1, dtype=numpy.uint32) * 1
+ self.w_x = numpy.array([self.L[2]])
+ self.w_m = numpy.array([self.rho[0] * self.np * math.pi * (cellsize/2.0)**3])
+ self.w_vel = numpy.array([0.0])
+ self.w_force = numpy.array([0.0])
+ self.w_devs = numpy.array([0.0])
# Adjust grid and upper wall for shear, and fix boundary particle velocities
t@@ -978,60 +950,35 @@ class Spherebin:
return porosity_grid
+def convert(graphicsformat = "png",
+ folder = "../img_out"):
+ """ Converts all PPM images in img_out to graphicsformat,
+ using ImageMagick """
+ # Convert images
+ subprocess.call("for F in " + folder + "/*.ppm ; do BASE=`basename $F`; convert $F $F." + graphicsformat + " > /dev/null ; done", shell=True)
+
+ # Remove PPM files
+ subprocess.call("rm " + folder + "/*.ppm", shell=True)
+
def render(binary,
- out = '../img_out/out',
- graphicsformat = 'jpg',
- resolution = numpy.array([800, 800]),
- workhorse = 'GPU',
- method = 'pressure',
- max_val = 4e3,
- rt_bin = '../raytracer/rt',
+ graphicsformat = 'png',
verbose=True):
""" Render target binary using the sphere raytracer.
"""
- stdout = ""
+ quiet = ""
if (verbose == False):
- stdout = " > /dev/null"
-
- # Use raytracer to render the scene into a temporary PPM file
- if workhorse == 'GPU':
- subprocess.call(rt_bin + " GPU {0} {1} {2} {3}.ppm {4} {5}"\
- .format(binary, resolution[0], resolution[1], out, method, max_val) + stdout, shell=True)
- if workhorse == 'CPU':
- subprocess.call(rt_bin + " CPU {0} {1} {2} {3}.ppm"\
- .format(binary, resolution[0], resolution[1], out) + stdout, shell=True)
-
- # Use ImageMagick's convert command to change the graphics format
- subprocess.call("convert {0}.ppm {0}.{1}"\
- .format(out, graphicsformat), shell=True)
+ quiet = "-q"
- # Delete temporary PPM file
- subprocess.call("rm {0}.ppm".format(out), shell=True)
+ # Render images using sphere raytracer
+ subprocess.call("cd ..; ./sphere_* " + quiet + " -r " + binary, shell=True)
-def renderAll(project,
- method = "pressure",
- max_val = 10e3,
- out_folder = "../img_out",
- graphics_format = "png",
- workhorse = "GPU",
- resolution = numpy.array([800, 800]),
- rt_bin = "~/code/sphere/raytracer/rt",
- verbose = False):
- lastfile = status(project)
-
- for i in range(lastfile+1):
- # Input binary filename
- fn = "../output/{0}.output{1}.bin".format(project, i)
-
- # Output image name (without extension)
- out = out_folder + "/{0}.output{1}".format(project, i)
-
- # Call raytracer, also converts to format
- render(fn, out, graphics_format, resolution, workhorse, method, max_val, rt_bin, verbose)
+ # Convert images to compressed format
+ convert()
+
def video(project,
out_folder = "./",
video_format = "mp4",
t@@ -1286,10 +1233,10 @@ def run(project, verbose=True, hideinputfile=False):
quiet = ""
stdout = ""
if (verbose == False):
- quiet = "-q "
+ quiet = "-q"
if (hideinputfile == True):
stdout = " > /dev/null"
- subprocess.call("cd ..; ./sphere_*_X86_64 "+ quiet + project + stdout, shell=True)
+ subprocess.call("cd ..; ./sphere_* " + quiet + " input/" + project + ".bin" + stdout, shell=True)
def status(project):
""" Check the status.dat file for the target project,
(DIR) diff --git a/python/tests.py b/python/tests.py
t@@ -1,5 +1,6 @@
#!/usr/bin/env python
from sphere import *
+import subprocess
def compare(first, second, string):
if (first == second):
t@@ -13,7 +14,7 @@ print("### Input/output tests ###")
# Generate data in python
orig = Spherebin(np = 100, nw = 0)
-orig.generateRadii()
+orig.generateRadii(histogram = False)
orig.defaultParams()
orig.initRandomGridPos(g = numpy.zeros(orig.nd))
orig.initTemporal(current = 0.0, total = 0.0)
t@@ -39,3 +40,6 @@ cuda.time_current = orig.time_current
cuda.time_step_count = orig.time_step_count
compare(orig, cuda, "CUDA IO: ")
+# Remove temporary files
+subprocess.call("rm orig.bin; rm ../output/orig*.bin", shell=True)
+
(DIR) diff --git a/src/device.cu b/src/device.cu
t@@ -46,14 +46,15 @@ __host__ void DEM::initializeGPU(void)
cudaGetDeviceCount(&devicecount);
if(devicecount == 0) {
- std::cerr << "\nERROR: No CUDA-enabled devices availible. Bye.\n";
+ std::cerr << "\nERROR: No CUDA-enabled devices availible. Bye."
+ << std::endl;
exit(EXIT_FAILURE);
} else if (devicecount == 1) {
if (verbose == 1)
- cout << "\nSystem contains 1 CUDA compatible device.\n";
+ cout << " System contains 1 CUDA compatible device.\n";
} else {
if (verbose == 1)
- cout << "\nSystem contains " << devicecount << " CUDA compatible devices.\n";
+ cout << " System contains " << devicecount << " CUDA compatible devices.\n";
}
cudaGetDeviceProperties(&prop, cudadevice);
t@@ -61,13 +62,13 @@ __host__ void DEM::initializeGPU(void)
cudaRuntimeGetVersion(&cudaRuntimeVersion);
if (verbose == 1) {
- cout << "Using CUDA device ID: " << cudadevice << "\n";
+ cout << " Using CUDA device ID: " << cudadevice << "\n";
cout << " - Name: " << prop.name << ", compute capability: "
<< prop.major << "." << prop.minor << ".\n";
cout << " - CUDA Driver version: " << cudaDriverVersion/1000
<< "." << cudaDriverVersion%100
<< ", runtime version " << cudaRuntimeVersion/1000 << "."
- << cudaRuntimeVersion%100 << "\n\n";
+ << cudaRuntimeVersion%100 << std::endl;
}
// Comment following line when using a system only containing exclusive mode GPUs
t@@ -457,16 +458,6 @@ __host__ void DEM::startTime()
char file[200];
FILE *fp;
- // Copy data to constant global device memory
- transferToConstantDeviceMemory();
-
- // Allocate device memory for particle variables,
- // tied to previously declared pointers in structures
- allocateGlobalDeviceMemory();
-
- // Transfer data from host to gpu device memory
- transferToGlobalDeviceMemory();
-
// Synchronization point
cudaThreadSynchronize();
checkForCudaErrors("Start of startTime()");
(DIR) diff --git a/src/file_io.cpp b/src/file_io.cpp
t@@ -38,7 +38,7 @@ void DEM::readbin(const char *target)
std::ifstream ifs(target, std::ios_base::binary);
if (!ifs) {
cerr << "Could not read input binary file '"
- << target << endl;
+ << target << "'" << endl;
exit(1);
}
t@@ -351,11 +351,14 @@ void DEM::writePPM(const char *target)
// Open output file
std::ofstream ofs(target);
if (!ofs) {
- std::cerr << "could create output PPM file '"
+ std::cerr << "Could not create output PPM file '"
<< target << std::endl;
exit(1); // Return unsuccessful exit status
}
+ if (verbose == 1)
+ std::cout << " Saving image: " << target << std::endl;
+
// Write PPM header
ofs << "P6 " << width << " " << height << " 255\n";
(DIR) diff --git a/src/sphere.cpp b/src/sphere.cpp
t@@ -45,12 +45,23 @@ DEM::DEM(const std::string inputbin,
// Initialize CUDA
initializeGPU();
+ // Copy constant data to constant device memory
+ transferToConstantDeviceMemory();
+
+ // Allocate device memory for particle variables,
+ // tied to previously declared pointers in structures
+ allocateGlobalDeviceMemory();
+
+ // Transfer data from host to gpu device memory
+ transferToGlobalDeviceMemory();
+
// Render image using raytracer if requested
if (render_img == 1) {
- float3 eye = make_float3(0.6f * grid.L[0],
- -2.5f * grid.L[1],
- 0.52f * grid.L[2]);
+ float3 eye = make_float3( 2.5f * grid.L[0],
+ -5.0f * grid.L[1],
+ 0.5f * grid.L[2]);
//float focalLength = 0.8f*grid.L[0];
+ //float3 eye = make_float3(0.0f, 0.0f, 0.0f);
render(eye);
}