sphere

GPU-based 3D discrete element method algorithm with optional fluid coupling
git clone git://src.adamsgaard.dk/sphere
Log | Files | Refs | LICENSE

commit 7d8e9335c718e6e47e25c41b5cdd5d48e92b9e27
parent bd456ae010d2426342427ff84fb19d5f2609525e
Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date:   Sat, 28 Sep 2013 10:16:28 +0200

updated examples

Diffstat:
Apython/collision.py | 60++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Apython/segregation.py | 140+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Dpython/triaxial.py | 88-------------------------------------------------------------------------------
3 files changed, 200 insertions(+), 88 deletions(-)

diff --git a/python/collision.py b/python/collision.py @@ -0,0 +1,60 @@ +#!/usr/bin/env python +""" +Example of two particles colliding. +Place script in sphere/python folder, and invoke with `python 2.7 collision.py` +""" + +# Import the sphere module for setting up, running, and analyzing the +# experiment. We also need the numpy module when setting arrays in the sphere +# object. +import sphere +import numpy + + +### SIMULATION SETUP + +# Create a sphere object with two preallocated particles and a simulation ID +SB = sphere.Spherebin(np = 2, sid = 'collision') + +SB.radius[:] = 0.3 # set radii to 0.3 m + +# Define the positions of the two particles +SB.x[0, :] = numpy.array([0.0, 0.0, 0.0]) # particle 1 (idx 0) +SB.x[1, :] = numpy.array([1.0, 0.0, 0.0]) # particle 2 (idx 1) + +# The default velocity is [0,0,0]. Slam particle 1 into particle 2 by defining +# a positive x velocity for particle 1. +SB.vel[0, 0] = 1.0 + +# let's disable gravity in this simulation +GRAVITY = numpy.array([0.0, 0.0, 0.0]) + +# Set the world limits and the particle sorting grid. The particles need to stay +# within the world limits for the entire simulation, otherwise it will stop! +SB.initGridAndWorldsize(g = GRAVITY, margin = 10.0) + +# Define the temporal parameters, e.g. the total time (total) and the file +# output interval (file_dt), both in seconds +SB.initTemporal(total = 2.0, file_dt = 0.1) + +# Save the simulation as a input file for sphere +SB.writebin() + +# Using a 'dry' run, the sphere main program will display important parameters. +# sphere will end after displaying these values. +SB.run(dry = True) + + +### RUNNING THE SIMULATION + +# Start the simulation on the GPU from the sphere program +SB.run() + + +### ANALYSIS OF SIMULATION RESULTS + +# Plot the system energy through time, image saved as collision-energy.png +sphere.visualize(SB.sid, method = 'energy') + +# Render the particles using the built-in raytracer +SB.render() diff --git a/python/segregation.py b/python/segregation.py @@ -0,0 +1,140 @@ +#!/usr/bin/env python + +# Import sphere functionality +from sphere import * + +### EXPERIMENT SETUP ### +initialization = True +consolidation = True +shearing = True +rendering = True +plots = True + +# Number of particles +np = 2e4 + +# Common simulation id +sim_id = "segregation" + +# Deviatoric stress [Pa] +#devslist = [80e3, 10e3, 20e3, 40e3, 60e3, 120e3] +devslist = [80e3] +#devs = 0 + +### INITIALIZATION ### + +# New class +init = Spherebin(np = np, nd = 3, nw = 0, sid = sim_id + "-init") + +# Save radii +init.generateRadii(radius_mean = 0.08) + +# Use default params +init.defaultParams(gamma_n = 0.0, mu_s = 0.3, mu_d = 0.3) + +# Initialize positions in random grid (also sets world size) +hcells = np**(1.0/3.0) +init.initRandomGridPos(gridnum = numpy.array([hcells, hcells, 1e9]), periodic = 1, contactmodel = 2) + +# Decrease size of top particles +fraction = 0.80 +z_min = numpy.min(init.x[:,2] - init.radius) +z_max = numpy.max(init.x[:,2] + init.radius) +I = numpy.nonzero(init.x[:,2] > (z_max - (z_max-z_min)*fraction)) +init.radius[I] = init.radius[I] * 0.5 + +# Set duration of simulation +init.initTemporal(total = 10.0) + +if (initialization == True): + # Write input file for sphere + init.writebin() + + # Run sphere + init.run(dry=True) + init.run() + + if (plots == True): + # Make a graph of energies + visualize(init.sid, "energy", savefig=True, outformat='png') + + if (rendering == True): + # Render images with raytracer + init.render(method = "angvel", max_val = 0.3, verbose = False) + + +### CONSOLIDATION ### + +for devs in devslist: + # New class + cons = Spherebin(np = init.np, nw = 1, sid = sim_id + "-cons-devs{}".format(devs)) + + # Read last output file of initialization step + lastf = status(sim_id + "-init") + cons.readbin("../output/" + sim_id + "-init.output{:0=5}.bin".format(lastf), verbose=False) + + # Setup consolidation experiment + cons.consolidate(deviatoric_stress = devs, periodic = init.periodic) + + + # Set duration of simulation + cons.initTemporal(total = 1.5) + #cons.initTemporal(total = 0.0019, file_dt = 0.00009) + #cons.initTemporal(total = 0.0019, file_dt = 1e-6) + #cons.initTemporal(total = 0.19, file_dt = 0.019) + + cons.w_m[0] *= 0.001 + + + + if (consolidation == True): + # Write input file for sphere + cons.writebin() + + # Run sphere + cons.run(dry=True) # show values, don't run + cons.run() # run + + if (plots == True): + # Make a graph of energies + visualize(cons.sid, "energy", savefig=True, outformat='png') + visualize(cons.sid, "walls", savefig=True, outformat='png') + + if (rendering == True): + # Render images with raytracer + cons.render(method = "pres", max_val = 2.0*devs, verbose = False) + + +### SHEARING ### + + # New class + shear = Spherebin(np = cons.np, nw = cons.nw, sid = sim_id + "-shear-devs{}".format(devs)) + + # Read last output file of initialization step + lastf = status(sim_id + "-cons-devs{}".format(devs)) + shear.readbin("../output/" + sim_id + "-cons-devs{}.output{:0=5}.bin".format(devs, lastf), verbose = False) + + # Setup shear experiment + shear.shear(shear_strain_rate = 0.05, periodic = init.periodic) + + # Set duration of simulation + #shear.initTemporal(total = 20.0) + shear.initTemporal(total = 400.0) + + if (shearing == True): + # Write input file for sphere + shear.writebin() + + # Run sphere + shear.run(dry=True) + shear.run() + + if (plots == True): + # Make a graph of energies + visualize(shear.sid, "energy", savefig=True, outformat='png') + visualize(shear.sid, "shear", savefig=True, outformat='png') + + if (rendering == True): + # Render images with raytracer + shear.render(method = "pres", max_val = 2.0*devs, verbose = False) + diff --git a/python/triaxial.py b/python/triaxial.py @@ -1,88 +0,0 @@ -#!/usr/bin/env python - -# Import sphere functionality -from sphere import * - -### EXPERIMENT SETUP ### -initialization = False -consolidation = True -rendering = True -plots = True - -# Number of particles -np = 2e3 - -# Common simulation id -sim_id = "triaxial-test" - -# Normal stress (sigma_3) -devs = 10e3 - - -### INITIALIZATION ### - -# New class -init = Spherebin(np = np, nd = 3, nw = 0, sid = sim_id + "-init") - -# Save radii -init.generateRadii(radius_mean = 0.05) - -# Use default params -init.defaultParams(gamma_n = 0.0, mu_s = 0.4, mu_d = 0.4) - -# Initialize positions in random grid (also sets world size) -init.initRandomGridPos(gridnum = numpy.array([12, 12, 1000]), periodic = 0, contactmodel = 2) - -# Set duration of simulation -init.initTemporal(total = 5.0) - -if (initialization == True): - # Write input file for sphere - init.writebin() - - # Run sphere - init.run() - - if (plots == True): - # Make a graph of energies - visualize(init.sid, "energy", savefig=True, outformat='png') - - #if (rendering == True): - # Render images with raytracer - #init.render(method = "angvel", max_val = 0.3, verbose = False) - - -### CONSOLIDATION ### - -# New class -cons = Spherebin(np = np, nw = 1, sid = sim_id + "-cons") - -# Read last output file of initialization step -lastf = status(sim_id + "-init") -cons.readbin("../output/" + sim_id + "-init.output{:0=5}.bin".format(lastf), verbose=False) - -# Setup triaxial experiment -cons.triaxial(wvel = -cons.L[2]*0.01, deviatoric_stress = devs) # two percent of height per second - -# Set duration of simulation -cons.initTemporal(total = 5.0) - -cons.w_m[:] *= 0.001 - -if (consolidation == True): - # Write input file for sphere - cons.writebin() - - # Run sphere - cons.run(dry=True) # show values, don't run - cons.run() # run - - if (plots == True): - # Make a graph of energies - visualize(cons.sid, "energy", savefig=True, outformat='png') - visualize(cons.sid, "walls", savefig=True, outformat='png') - - if (rendering == True): - # Render images with raytracer - cons.render(method = "pres", max_val = 1e4, verbose = False) -