tsmall changes in scripts - 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 dcbfc0526787a023a06cc57200e81ff87a3eb709
(DIR) parent ef0e8ba6bae7051abf968c42029f18fa8352ba10
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Wed, 30 Oct 2013 09:12:01 +0100
small changes in scripts
Diffstat:
M python/darcy.py | 79 ++++++++++---------------------
M python/segregation.py | 64 +++++++++++++++++++++++++++++--
M python/uniaxial.py | 2 +-
3 files changed, 85 insertions(+), 60 deletions(-)
---
(DIR) diff --git a/python/darcy.py b/python/darcy.py
t@@ -14,7 +14,7 @@ plots = True
# Number of particles
-#np = 1e2
+#np = 2e2
np = 1e4
# Common simulation id
t@@ -37,6 +37,7 @@ init.defaultParams(mu_s = 0.4, mu_d = 0.4, nu = 8.9e-4)
# Initialize positions in random grid (also sets world size)
#init.initRandomGridPos(gridnum = numpy.array([9, 9, 1000]), periodic = 1, contactmodel = 2)
+#init.initRandomGridPos(gridnum = numpy.array([10, 10, 1000]), periodic = 1, contactmodel = 1)
#init.initRandomGridPos(gridnum = numpy.array([32, 32, 1000]), periodic = 1, contactmodel = 2)
init.initRandomGridPos(gridnum = numpy.array([32, 32, 1000]), periodic = 1, contactmodel = 1)
t@@ -44,14 +45,15 @@ init.initRandomGridPos(gridnum = numpy.array([32, 32, 1000]), periodic = 1, cont
#init.random2bonds(spacing=0.1)
# Set duration of simulation
-init.initTemporal(total = 7.0)
+init.initTemporal(total = 1.0)
+#init.initTemporal(total = 0.01)
init.time_file_dt[0] = 0.05
#init.time_file_dt[0] = init.time_dt[0]*0.99
#init.time_total[0] = init.time_dt[0]*2.0
#init.initTemporal(total = 0.5)
-#init.time_file_dt[0] = init.time_total[0]/5.0
+#init.time_dt[0] = 1.0e-5;
-#init.f_rho[2,2,4] = 5.0
+init.f_rho[2,2,4] = 1.1
#init.f_rho[6,6,10] = 1.1
#init.f_rho[:,:,-1] = 1.0001
t@@ -65,54 +67,21 @@ if (initialization == True):
init.run(darcyflow=True)
-### 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(darcyflow=True) # 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)
-
- project = cons.sid
- lastfile = status(cons.sid)
- sb = Spherebin()
- for i in range(lastfile+1):
- fn = "../output/{0}.output{1:0=5}.bin".format(project, i)
- sb.sid = project + ".output{:0=5}".format(i)
- sb.readbin(fn, verbose = False)
- for y in range(0,sb.num[1]):
- sb.plotFluidDensities(y = y)
- sb.plotFluidVelocities(y = y)
+ 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 = "pres", max_val = 2.0*devs, verbose = False)
+
+ project = init.sid
+ lastfile = status(init.sid)
+ sb = Spherebin()
+ for i in range(lastfile+1):
+ fn = "../output/{0}.output{1:0=5}.bin".format(project, i)
+ sb.sid = project + ".output{:0=5}".format(i)
+ sb.readbin(fn, verbose = False)
+ for y in range(0,sb.num[1]):
+ sb.plotFluidDensities(y = y)
+ sb.plotFluidVelocities(y = y)
(DIR) diff --git a/python/segregation.py b/python/segregation.py
t@@ -18,7 +18,7 @@ sim_id = "segregation"
# Deviatoric stress [Pa]
#devslist = [80e3, 10e3, 20e3, 40e3, 60e3, 120e3]
-devslist = [80e3]
+devslist = [120e3]
#devs = 0
### INITIALIZATION ###
t@@ -30,10 +30,10 @@ init = Spherebin(np = np, nd = 3, nw = 0, sid = sim_id + "-init")
init.generateRadii(radius_mean = 0.08)
# Use default params
-init.defaultParams(gamma_n = 0.0, mu_s = 0.3, mu_d = 0.3)
+init.defaultParams(gamma_n = 0.0, mu_s = 0.4, mu_d = 0.4)
# Initialize positions in random grid (also sets world size)
-hcells = np**(1.0/3.0)
+hcells = np**(1.0/3.0)*0.6
init.initRandomGridPos(gridnum = numpy.array([hcells, hcells, 1e9]), periodic = 1, contactmodel = 2)
# Decrease size of top particles
t@@ -115,7 +115,63 @@ for devs in devslist:
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)
+ #shear.shear(shear_strain_rate = 0.05, periodic = init.periodic)
+ shear_strain_rate = 0.05
+
+ ## Custom shear function
+ # Find lowest and heighest point
+ z_min = numpy.min(shear.x[:,2] - shear.radius)
+ z_max = numpy.max(shear.x[:,2] + shear.radius)
+
+ # the grid cell size is equal to the max. particle diameter
+ cellsize = shear.L[0] / shear.num[0]
+
+ # make grid one cell heigher to allow dilation
+ shear.num[2] += 1
+ shear.L[2] = shear.num[2] * cellsize
+
+ # zero kinematics
+ shear.zeroKinematics()
+
+ # set friction coefficients
+ shear.mu_s[0] = 0.5
+ shear.mu_d[0] = 0.5
+
+ # set the thickness of the horizons of fixed particles
+ #fixheight = 2*cellsize
+ #fixheight = cellsize
+
+ # Fix horizontal velocity to 0.0 of lowermost particles
+ d_max_below = numpy.max(shear.radius[numpy.nonzero(shear.x[:,2] <
+ (z_max-z_min)*0.3)])*2.0
+ #I = numpy.nonzero(shear.x[:,2] < (z_min + fixheight))
+ I = numpy.nonzero(shear.x[:,2] < (z_min + d_max_below))
+ shear.fixvel[I] = 1
+ shear.angvel[I,0] = 0.0
+ shear.angvel[I,1] = 0.0
+ shear.angvel[I,2] = 0.0
+ shear.vel[I,0] = 0.0 # x-dim
+ shear.vel[I,1] = 0.0 # y-dim
+
+ # Copy bottom fixed particles to top
+ z_offset = z_max-z_min
+ shearvel = (z_max-z_min)*shear_strain_rate
+ for i in I[0]:
+ x = shear.x[i,:] + numpy.array([0.0, 0.0, z_offset])
+ vel = numpy.array([shearvel, 0.0, 0.0])
+ shear.addParticle(x = x, radius = shear.radius[i], fixvel = 1, vel = vel)
+
+
+ # Set wall viscosities to zero
+ shear.gamma_wn[0] = 0.0
+ shear.gamma_wt[0] = 0.0
+
+ # Set wall friction coefficients to zero
+ shear.mu_ws[0] = 0.0
+ shear.mu_wd[0] = 0.0
+
+ # Readjust top wall
+ shear.adjustUpperWall()
# Set duration of simulation
#shear.initTemporal(total = 20.0)
(DIR) diff --git a/python/uniaxial.py b/python/uniaxial.py
t@@ -4,7 +4,7 @@
from sphere import *
### EXPERIMENT SETUP ###
-initialization = True
+initialization = False
consolidation = True
rendering = True
plots = True