tMerge branch 'master' of github.com:anders-dc/sphere - sphere - GPU-based 3D d… | |
git clone git://src.adamsgaard.dk/sphere | |
Log | |
Files | |
Refs | |
LICENSE | |
--- | |
commit 3ed4555efc581a8084a47404fb7e6fb0a6d72258 | |
parent 69396d3dce34eea330a77c056986ea6f94f9dd20 | |
Author: Anders Damsgaard <[email protected]> | |
Date: Tue, 22 Apr 2014 09:36:04 +0200 | |
Merge branch 'master' of github.com:anders-dc/sphere | |
Diffstat: | |
M python/collapse.py | 36 +++++++++++++++++++++--------… | |
M python/collision.py | 4 ++-- | |
M python/fluidshear.py | 0 | |
D python/ns_cons_of_mass.py | 203 -----------------------------… | |
D python/ns_cons_of_mass_run.sh | 4 ---- | |
M python/segregation.py | 264 ++++++++++++++++-------------… | |
M python/shear-test.py | 3 +++ | |
D python/uniaxial.py | 85 -----------------------------… | |
8 files changed, 165 insertions(+), 434 deletions(-) | |
--- | |
diff --git a/python/collapse.py b/python/collapse.py | |
t@@ -13,23 +13,34 @@ plots = True | |
np = 1e4 | |
# Common simulation id | |
-sim_id = "collapse" | |
+sim_id = 'collapse' | |
### INITIALIZATION ### | |
# New class | |
-init = sphere.sim(np = np, nd = 3, nw = 0, sid = sim_id + "-init") | |
+init = sphere.sim(np = np, nd = 3, nw = 0, sid = sim_id + '-init') | |
-# Save radii | |
+# Set radii | |
init.generateRadii(radius_mean = 0.1) | |
# Use default params | |
-init.defaultParams(k_n = 1.0e6, k_t = 1.0e6, gamma_n = 100.0, gamma_t = 100.0,… | |
+init.defaultParams( | |
+ k_n = 1.0e6, k_t = 1.0e6, # normal and tangential stiffnesses | |
+ gamma_n = 100.0, gamma_t = 100.0, # normal and tangential viscosities | |
+ mu_s = 0.3, mu_d = 0.3) # static and dynamic frictions | |
# Initialize positions in random grid (also sets world size) | |
hcells = np**(1.0/3.0) | |
-init.initRandomGridPos(gridnum = numpy.array([hcells, hcells, 1e9]), periodic … | |
-#init.initRandomGridPos(gridnum = numpy.array([hcells, hcells, 1e9]), periodic… | |
+init.initRandomGridPos(gridnum = numpy.array([hcells, hcells, 1e9])) | |
+ | |
+# Choose the tangential contact model | |
+# 1) Visco-frictional (somewhat incorrect, fast computations) | |
+# 2) Elastic-viscous-frictional (more correct, slow computations in dense | |
+# packings) | |
+init.contactmodel[0] = 1 | |
+ | |
+# Add gravitational acceleration | |
+init.g[2] = -10.0 | |
# Set duration of simulation, automatically determine timestep, etc. | |
init.initTemporal(total = 10.0) | |
t@@ -42,7 +53,7 @@ if (initialization == True): | |
if (plots == True): | |
# Make a graph of energies | |
- visualize(init.sid, "energy", savefig=True, outformat='png') | |
+ init.visualize('energy', savefig=True, outformat='png') | |
### COLLAPSE ### | |
t@@ -50,8 +61,9 @@ if (initialization == True): | |
coll = sphere.sim(np = init.np, nw = init.nw, sid = sim_id) | |
# Read last output file of initialization step | |
-lastf = status(sim_id + "-init") | |
-coll.readbin("../output/" + sim_id + "-init.output{:0=5}.bin".format(lastf), v… | |
+lastf = status(sim_id + '-init') | |
+coll.readbin('../output/' + sim_id + '-init.output{:0=5}.bin'.format(lastf), | |
+ verbose = False) | |
# Setup collapse experiment by moving the +x boundary and resizing the grid | |
resizefactor = 3 | |
t@@ -73,10 +85,10 @@ if (collapse == True): | |
if (plots == True): | |
# Make a graph of the energies | |
- visualize(coll.sid, "energy", savefig=True, outformat='png') | |
+ init.visualize('energy', savefig=True, outformat='png') | |
if (rendering == True): | |
# Render images with raytracer with linear velocity as the color code | |
- print("Rendering images with raytracer") | |
- coll.render(method = "vel", max_val = 1.0, verbose = False) | |
+ print('Rendering images with raytracer') | |
+ coll.render(method = 'vel', max_val = 1.0, verbose = False) | |
coll.video() | |
diff --git a/python/collision.py b/python/collision.py | |
t@@ -1,8 +1,8 @@ | |
#!/usr/bin/env python | |
-""" | |
+''' | |
Example of two particles colliding. | |
Place script in sphere/python/ folder, and invoke with `python 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 | |
diff --git a/python/fluidshear.py b/python/fluidshear.py | |
diff --git a/python/ns_cons_of_mass.py b/python/ns_cons_of_mass.py | |
t@@ -1,203 +0,0 @@ | |
-#!/usr/bin/env python | |
- | |
-# Import sphere functionality | |
-import sphere | |
-from sphere import visualize, status | |
-import sys | |
-import numpy | |
-import pylab | |
- | |
-figformat = 'pdf' | |
- | |
-# Number of particles | |
-np = 1e4 | |
-#np = 50 | |
- | |
-# Common simulation id | |
-sim_id = "ns" | |
- | |
-# New class | |
-init = sphere.sim(np = np, nd = 3, nw = 0, sid = sim_id + "-init") | |
- | |
- | |
-# Save radii | |
-init.generateRadii(radius_mean = 0.05, histogram=False) | |
- | |
- | |
-# Use default params | |
-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([80, 80, 1000]), periodic = 1, c… | |
-init.initRandomGridPos(gridnum = numpy.array([40, 40, 1000]), periodic = 1, co… | |
-#init.initRandomGridPos(gridnum = numpy.array([6, 6, 1000]), periodic = 1, con… | |
- | |
-# Set duration of simulation | |
-#init.initTemporal(total = 2.5) | |
-#init.time_file_dt[0] = 0.01 | |
- | |
-#init.initTemporal(total = 0.01) | |
-init.initTemporal(total = 0.002) | |
-init.time_file_dt[0] = 0.001 | |
- | |
-#init.initTemporal(1) | |
-#init.time_file_dt[0] = init.time_dt[0]*0.9 | |
-#init.time_total[0] = init.time_file_dt[0]*10.5 | |
- | |
- | |
- | |
-# Small pertubation | |
-#init.p_f[init.num[0]/2,init.num[1]/2,init.num[2]/2] = 2.0 | |
-#init.p_f[:,:,init.num[2]-1] = 1.0 | |
-#init.p_f[:,:,0] = 2.0 | |
- | |
-#init.g[2] = -10.0 | |
-init.g[2] = 0.0 # uncomment to disable gravity | |
- | |
-# Run sphere | |
-init.run(dry=True) | |
-init.run(cfd=True) | |
- | |
-init.writeVTKall() | |
- | |
-#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 = sphere.sim() | |
-time = numpy.zeros(lastfile+1) | |
-sum_op_f = numpy.zeros(lastfile+1) | |
-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) | |
- | |
- time[i] = sb.time_current[0] | |
- sum_op_f[i] = sb.p_f.sum() | |
- | |
-#stack = numpy.vstack((time,sum_op_f)) | |
-#numpy.savetxt("sum_op_f", stack) | |
- | |
-# Set figure parameters | |
-fig_width_pt = 246.0 # Get this from LaTeX using \showthe\columnwidth | |
-#fig_width_pt = 50.0 # Get this from LaTeX using \showthe\columnwidth | |
-inches_per_pt = 1.0/72.27 # Convert pt to inch | |
-golden_mean = (pylab.sqrt(5)-1.0)/2.0 # Aesthetic ratio | |
-fig_width = fig_width_pt*inches_per_pt # width in inches | |
-fig_height = fig_width*golden_mean # height in inches | |
-fig_size = [fig_width,fig_height] | |
-params = {'backend': 'ps', | |
- 'axes.labelsize': 10, | |
- 'text.fontsize': 10, | |
- 'legend.fontsize': 10, | |
- 'xtick.labelsize': 8, | |
- 'ytick.labelsize': 8, | |
- 'text.usetex': | |
- True, | |
- 'figure.figsize': | |
- fig_size} | |
-pylab.rcParams.update(params) | |
- | |
-# generate data | |
-#x = pylab.arange(-2*pylab.pi,2*pylab.pi,0.01) | |
-#y1 = pylab.sin(x) | |
-#y2 = pylab.cos(x) | |
- | |
- | |
-y = init.num[1]/2 | |
- | |
-t = 0 | |
-fn = "../output/{0}.output{1:0=5}.bin".format(project, t) | |
-sb.sid = project + ".output{:0=5}".format(i) | |
-sb.readbin(fn, verbose = False) | |
-pylab.figure(1) | |
-pylab.clf() | |
-pylab.axes([0.125,0.2,0.95-0.125,0.95-0.2]) | |
-pylab.imshow(sb.p_f[:,y,:].T, origin='lower', interpolation='nearest', | |
- cmap=pylab.cm.Blues) | |
-#imgplt.set_interpolation('nearest') | |
-#pylab.set_interpolation('bicubic') | |
-#imgplt.set_cmap('hot') | |
-pylab.xlabel('$i_x$') | |
-pylab.ylabel('$i_z$') | |
-pylab.title('$t = {}$ s'.format(t*sb.time_file_dt[0])) | |
-#cb = pylab.colorbar(orientation = 'horizontal', shrink=0.8, pad=0.23) | |
-cb = pylab.colorbar(orientation = 'horizontal', shrink=0.8, pad=0.23, ticks=[s… | |
-cb.ax.set_xticklabels([1.0, sb.p_f[:,y,:].max()]) | |
-cb.set_label('H [m]') | |
- | |
-pylab.savefig('ns_cons_of_mass1.' + figformat) | |
-pylab.clf() | |
- | |
-t = 1 | |
-fn = "../output/{0}.output{1:0=5}.bin".format(project, t) | |
-sb.sid = project + ".output{:0=5}".format(i) | |
-sb.readbin(fn, verbose = False) | |
-pylab.figure(1) | |
-pylab.clf() | |
-pylab.axes([0.125,0.2,0.95-0.125,0.95-0.2]) | |
-pylab.imshow(sb.p_f[:,y,:].T, origin='lower', interpolation='nearest', | |
- cmap=pylab.cm.Blues) | |
-#pylab.set_interpolation('bicubic') | |
-pylab.xlabel('$i_x$') | |
-pylab.ylabel('$i_z$') | |
-pylab.title('$t = {}$ s'.format(t*sb.time_file_dt[0])) | |
-cb = pylab.colorbar(orientation = 'horizontal', shrink=0.8, pad=0.23, ticks=[s… | |
-#cb.ax.set_xticklabels([1.0, sb.p_f[:,y,:].max()]) | |
-cb.ax.set_xticklabels([sb.p_f[:,y,:].min(), sb.p_f[:,y,:].max()]) | |
-cb.set_label('H [m]') | |
-#pylab.tight_layout() | |
-pylab.savefig('ns_cons_of_mass2.' + figformat) | |
-pylab.clf() | |
- | |
-t = init.status() | |
-fn = "../output/{0}.output{1:0=5}.bin".format(project, t) | |
-sb.sid = project + ".output{:0=5}".format(i) | |
-sb.readbin(fn, verbose = False) | |
-pylab.figure(1) | |
-pylab.clf() | |
-pylab.axes([0.125,0.2,0.95-0.125,0.95-0.2]) | |
-pylab.imshow(sb.p_f[:,y,:].T, origin='lower', interpolation='nearest', | |
- cmap=pylab.cm.Blues) | |
-#pylab.set_interpolation('bicubic') | |
-pylab.xlabel('$i_x$') | |
-pylab.ylabel('$i_z$') | |
-pylab.title('$t = {}$ s'.format(t*sb.time_file_dt[0])) | |
-cb = pylab.colorbar(orientation = 'horizontal', shrink=0.8, pad=0.23, ticks=[s… | |
-#cb.ax.set_xticklabels([1.0, sb.p_f[:,y,:].max()]) | |
-cb.ax.set_xticklabels([sb.p_f[:,y,:].min(), sb.p_f[:,y,:].max()]) | |
-cb.set_label('H [m]') | |
-#pylab.tight_layout() | |
-pylab.savefig('ns_cons_of_mass3.' + figformat) | |
-pylab.clf() | |
- | |
-#pylab.axes([0.125,0.2,0.95-0.125,0.95-0.2]) | |
-pylab.axes([0.20,0.2,0.95-0.20,0.95-0.2]) | |
-pylab.plot(time, sum_op_f, '-k') | |
-pylab.xlabel('$t$ [s]') | |
-pylab.ylabel('$\sum H_i$ [m]') | |
-pylab.xlim([0,time.max()]) | |
-pylab.ylim([0,sum_op_f.max()*1.1]) | |
-#pylab.legend() | |
-#pylab.tight_layout() | |
-pylab.grid() | |
-pylab.savefig('ns_cons_of_mass4.' + figformat) | |
-pylab.clf() | |
- | |
-#pylab.tight_layout(h_pad=6.0) | |
-#pylab.savefig('ns_cons_of_mass.eps') | |
- | |
-#fig = matplotlib.pyplot.figure(figsize=(10,5),dpi=300) | |
-#ax = matplotlib.pyplot.subplot2grid((1, 1), (0, 0)) | |
-#ax.plot(time, sum_op_f) | |
-#fig.savefig("ns_cons_of_mass.png") | |
-#fig.clf() | |
diff --git a/python/ns_cons_of_mass_run.sh b/python/ns_cons_of_mass_run.sh | |
t@@ -1,4 +0,0 @@ | |
-#!/bin/sh | |
-rm ../output/*.vti; rm ../output/*.vtu; \ | |
- cd ~/code/sphere-cfd && cmake . && make -j2 &&\ | |
- cd ~/code/sphere-cfd/python && ipython -i ns_cons_of_mass.py | |
diff --git a/python/segregation.py b/python/segregation.py | |
t@@ -14,7 +14,7 @@ plots = True | |
np = 2e4 | |
# Common simulation id | |
-sim_id = "segregation" | |
+sim_id = 'segregation' | |
# Deviatoric stress [Pa] | |
#devslist = [80e3, 10e3, 20e3, 40e3, 60e3, 120e3] | |
t@@ -24,7 +24,7 @@ devslist = [120e3] | |
### INITIALIZATION ### | |
# New class | |
-init = sphere.sim(np = np, nd = 3, nw = 0, sid = sim_id + "-init") | |
+init = sphere.sim(np = np, nd = 3, nw = 0, sid = sim_id + '-init') | |
# Save radii | |
init.generateRadii(radius_mean = 0.08) | |
t@@ -34,7 +34,19 @@ init.defaultParams(gamma_n = 100.0, mu_s = 0.4, mu_d = 0.4) | |
# Initialize positions in random grid (also sets world size) | |
hcells = np**(1.0/3.0)*0.6 | |
-init.initRandomGridPos(gridnum = numpy.array([hcells, hcells, 1e9]), periodic … | |
+init.initRandomGridPos(gridnum = numpy.array([hcells, hcells, 1e9])) | |
+ | |
+# Choose the tangential contact model | |
+# 1) Visco-frictional (somewhat incorrect, fast computations) | |
+# 2) Elastic-viscous-frictional (more correct, slow computations in dense | |
+# packings) | |
+init.contactmodel[0] = 2 | |
+ | |
+# Set horizontal boundaries as periodic | |
+init.periodicBoundariesXY() | |
+ | |
+# Add gravity | |
+init.g[2] = -9.81 | |
# Decrease size of top particles | |
fraction = 0.80 | |
t@@ -48,143 +60,139 @@ init.initTemporal(total = 10.0) | |
if (initialization == True): | |
- # 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 = sphere.sim(np = init.np, nw = 1, sid = sim_id + "-cons-devs{}".format… | |
- | |
- # Read last output file of initialization step | |
- lastf = sphere.status(sim_id + "-init") | |
- cons.readbin("../output/" + sim_id + "-init.output{:0=5}.bin".format(lastf),… | |
- | |
- # 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): | |
- | |
# Run sphere | |
- cons.run(dry=True) # show values, don't run | |
- cons.run() # run | |
+ init.run(dry=True) | |
+ init.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') | |
+ # Make a graph of energies | |
+ init.visualize('energy', 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 = sphere.sim(np = cons.np, nw = cons.nw, sid = sim_id + "-shear-devs{}… | |
- | |
- # Read last output file of initialization step | |
- lastf = sphere.status(sim_id + "-cons-devs{}".format(devs)) | |
- shear.readbin("../output/" + sim_id + "-cons-devs{}.output{:0=5}.bin".format… | |
- | |
- # Setup shear experiment | |
- #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 | |
+ # Render images with raytracer | |
+ init.render(method = 'angvel', max_val = 0.3, verbose = False) | |
- # 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 | |
+### CONSOLIDATION ### | |
- # Readjust top wall | |
- shear.adjustUpperWall() | |
+for devs in devslist: | |
+ # New class | |
+ cons = sphere.sim(np = init.np, nw = 1, sid = sim_id | |
+ + '-cons-devs{}'.format(devs)) | |
+ | |
+ # Read last output file of initialization step | |
+ lastf = sphere.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) | |
+ | |
+ if (consolidation == True): | |
+ | |
+ # Run sphere | |
+ cons.run(dry=True) # show values, don't run | |
+ cons.run() # run | |
+ | |
+ if (plots == True): | |
+ # Make a graph of energies | |
+ cons.visualize('energy', savefig=True, outformat='png') | |
+ cons.visualize('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 = sphere.sim(np = cons.np, nw = cons.nw, | |
+ sid = sim_id + '-shear-devs{}'.format(devs)) | |
+ | |
+ # Read last output file of initialization step | |
+ lastf = sphere.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) | |
+ 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) | |
- shear.initTemporal(total = 400.0) | |
+ # Set duration of simulation | |
+ #shear.initTemporal(total = 20.0) | |
+ shear.initTemporal(total = 400.0) | |
- if (shearing == True): | |
+ if (shearing == True): | |
- # Run sphere | |
- shear.run(dry=True) | |
- shear.run() | |
+ # 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') | |
+ # Make a graph of energies | |
+ shear.visualize('energy', savefig=True, outformat='png') | |
+ shear.visualize('shear', savefig=True, outformat='png') | |
if (rendering == True): | |
- # Render images with raytracer | |
- shear.render(method = "pres", max_val = 2.0*devs, verbose = False) | |
- | |
+ # Render images with raytracer | |
+ shear.render(method = 'pres', max_val = 2.0*devs, verbose = False) | |
diff --git a/python/shear-test.py b/python/shear-test.py | |
t@@ -31,6 +31,9 @@ init.generateRadii(radius_mean = 0.02) | |
# Use default params | |
init.defaultParams(gamma_n = 100.0, mu_s = 0.6, mu_d = 0.6) | |
+# Add gravity | |
+init.g[2] = -9.81 | |
+ | |
# Initialize positions in random grid (also sets world size) | |
hcells = np**(1.0/3.0) | |
init.initRandomGridPos(gridnum = numpy.array([hcells, hcells, 1e9]), periodic … | |
diff --git a/python/uniaxial.py b/python/uniaxial.py | |
t@@ -1,85 +0,0 @@ | |
-#!/usr/bin/env python | |
- | |
-# Import sphere functionality | |
-import sphere | |
- | |
-### EXPERIMENT SETUP ### | |
-initialization = False | |
-consolidation = True | |
-rendering = True | |
-plots = True | |
- | |
-# Number of particles | |
-np = 2e3 | |
- | |
-# Common simulation id | |
-sim_id = "uniaxial-test" | |
- | |
-### INITIALIZATION ### | |
- | |
-# New class | |
-init = sphere.sim(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, co… | |
- | |
-# Set duration of simulation | |
-init.initTemporal(total = 5.0) | |
- | |
-if (initialization == True): | |
- | |
- # 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 = sphere.sim(np = np, nw = 1, sid = sim_id + "-cons") | |
- | |
-# Read last output file of initialization step | |
-lastf = sphere.status(sim_id + "-init") | |
-cons.readbin("../output/" + sim_id + "-init.output{:0=5}.bin".format(lastf), v… | |
- | |
-# Setup consolidation experiment | |
-cons.uniaxialStrainRate(wvel = -cons.L[2]*0.05) # five percent of height per s… | |
- | |
-# Set duration of simulation | |
-cons.initTemporal(total = 3.0) | |
-#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): | |
- | |
- # 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) | |
- |