tadd channel script - 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 e8e839893bdaf652fcc6ecec92ecefc7bd1e7885
(DIR) parent b2c51135074e3ea2410b9320d7a595ce2262ced4
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Sat, 5 Mar 2016 18:21:28 +0100
add channel script
Diffstat:
A python/channel.py | 155 +++++++++++++++++++++++++++++++
1 file changed, 155 insertions(+), 0 deletions(-)
---
(DIR) diff --git a/python/channel.py b/python/channel.py
t@@ -0,0 +1,155 @@
+#!/usr/bin/env python
+import sphere
+import numpy
+
+relaxation = False
+consolidation = True
+
+id_prefix = 'channel3'
+N = 40e3
+
+cube = sphere.sim('cube-init')
+cube.readlast()
+cube.adjustUpperWall(z_adjust=1.0)
+
+# Fill out grid with cubic packages
+grid = numpy.array((
+ [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+ [1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+ [1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+ [1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+ [1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+ [1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+ [1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+ [1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+ [1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+ [1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+ [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
+ [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
+ [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
+ [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
+ [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
+ [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
+ [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
+ [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
+ [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
+ [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
+ [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
+ [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
+ [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
+ [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
+ [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
+ [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
+ [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
+ [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
+ [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
+ [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
+))
+
+# World dimensions and cube grid
+nx = grid.shape[1] # horizontal cubes
+ny = 1 # horizontal (thickness) cubes
+nz = grid.shape[0] # vertical cubes
+dx = cube.L[0]
+dy = cube.L[1]
+dz = cube.L[2]
+Lx = dx*nx
+Ly = dy*ny
+Lz = dz*nz
+
+sim = sphere.sim(id_prefix + '-relaxation', nw=0)
+
+# insert particles into each cube
+for z in range(nz):
+ for y in range(ny):
+ for x in range(nx):
+
+ if (grid[z, x] == 0):
+ continue # skip to next iteration
+
+ for i in range(cube.np):
+ pos = [cube.x[i, 0] + x*dx,
+ cube.x[i, 1] + y*dy,
+ cube.x[i, 2] + (nz-z)*dz]
+
+ sim.addParticle(pos, radius=cube.radius[i], color=grid[z, x])
+
+# move to x=0
+min_x = numpy.min(sim.x[:, 0] - sim.radius[:])
+sim.x[:, 0] = sim.x[:, 0] - min_x
+
+# move to y=0
+min_y = numpy.min(sim.x[:, 1] - sim.radius[:])
+sim.x[:, 1] = sim.x[:, 1] - min_y
+
+# move to z=0
+min_z = numpy.min(sim.x[:, 2] - sim.radius[:])
+sim.x[:, 2] = sim.x[:, 2] - min_z
+
+sim.defineWorldBoundaries(L=[numpy.max(sim.x[:, 0] + sim.radius[:]),
+ numpy.max(sim.x[:, 1] + sim.radius[:]),
+ numpy.max(sim.x[:, 2] + sim.radius[:])*1.2])
+sim.k_t[0] = 2.0/3.0*sim.k_n[0]
+
+# sim.cleanup()
+sim.writeVTK()
+print("Number of particles: " + str(sim.np[0]))
+
+
+# Relaxation
+
+# Add gravitational acceleration
+sim.g[0] = 0.0
+sim.g[1] = 0.0
+sim.g[2] = -9.81
+
+sim.normalBoundariesXY()
+# sim.consolidate(normal_stress=0.0)
+
+# assign automatic colors, overwriting values from grid array
+sim.checkerboardColors(nx=grid.shape[1]/2, ny=2, nz=grid.shape[0])
+
+# Set duration of simulation, automatically determine timestep, etc.
+sim.initTemporal(total=3.0, file_dt=0.01, epsilon=0.07)
+sim.zeroKinematics()
+
+if relaxation:
+ sim.run(dry=True)
+ sim.run()
+ sim.writeVTKall()
+
+
+# Consolidation under constant normal stress
+sim.readlast()
+sim.id(id_prefix + '-' + str(int(N/100.)) + 'kPa')
+sim.cleanup()
+sim.initTemporal(current=0.0, total=10.0, file_dt=0.01, epsilon=0.07)
+
+# fix lowest plane of particles
+I = numpy.nonzero(sim.x[:, 2] < 1.5*numpy.mean(sim.radius))
+sim.fixvel[I] = -1
+sim.color[I] = 0
+
+sim.zeroKinematics()
+
+# Wall parameters
+sim.mu_ws[0] = 0.5
+sim.mu_wd[0] = 0.5
+sim.gamma_wn[0] = 1.0e2
+sim.gamma_wt[0] = 1.0e2
+# sim.gamma_wn[0] = 0.0
+# sim.gamma_wt[0] = 0.0
+
+# Particle parameters
+sim.mu_s[0] = 0.5
+sim.mu_d[0] = 0.5
+sim.gamma_n[0] = 0.0
+sim.gamma_t[0] = 0.0
+
+# apply effective normal stress from upper wall
+sim.consolidate(normal_stress=N)
+
+if consolidation:
+ sim.run(dry=True)
+ sim.run()
+ sim.writeVTKall()