tdeformationdepth.py - 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
---
tdeformationdepth.py (4949B)
---
1 #!/usr/bin/env python
2
3 # Import sphere functionality
4 import sphere
5 import numpy
6
7 ### EXPERIMENT SETUP ###
8 initialization = False
9 consolidation = True
10 shearing_norot = True
11 shearing = False
12 rendering = False
13 plots = True
14
15 # Number of particles
16 np = 1e4
17
18 # Common simulation id
19 sim_id = "defdepth"
20
21 # Deviatoric stress [Pa]
22 Nlist = [100e3, 10e3, 20e3, 50e3]
23
24 ### INITIALIZATION ###
25
26 # New class
27 init = sphere.sim(np = np, nd = 3, nw = 0, sid = sim_id + "-init")
28
29 # Save radii
30 init.generateRadii(mean = 0.02)
31
32 # Use default params
33 init.defaultParams(gamma_n = 100.0, mu_s = 0.6, mu_d = 0.6)
34
35 # GPU device (-1: auto select)
36 device = -1
37
38 # Add gravity
39 init.g[2] = -9.81
40
41 # Periodic x and y boundaries
42 init.periodicBoundariesXY()
43
44 # Initialize positions in random grid (also sets world size)
45 hcells = np**(1.0/3.0)
46 init.initRandomGridPos(gridnum = [hcells, hcells, 1e9])
47
48 # Set duration of simulation
49 init.initTemporal(total = 5.0)
50
51 if (initialization == True):
52
53 # Run sphere
54 init.run(dry = True)
55 init.run(device = device)
56
57 if (plots == True):
58 # Make a graph of energies
59 init.visualize('energy')
60
61 init.writeVTKall()
62
63 if (rendering == True):
64 # Render images with raytracer
65 init.render(method = "angvel", max_val = 0.3, verbose = False)
66
67
68 # For each normal stress, consolidate and subsequently shear the material
69 for N in Nlist:
70
71 ### CONSOLIDATION ###
72
73 # New class
74 cons = sphere.sim(np = init.np, nw = 1, sid = sim_id +
75 "-cons-N{}".format(N))
76
77 # Read last output file of initialization step
78 lastf = sphere.status(sim_id + "-init")
79 cons.readbin("../output/" + sim_id + "-init.output{:0=5}.bin".format(lastf), verbose=False)
80
81 # Periodic x and y boundaries
82 cons.periodicBoundariesXY()
83
84 # Setup consolidation experiment
85 cons.consolidate(normal_stress = N)
86 cons.adaptiveGrid()
87
88 # Set duration of simulation
89 cons.initTemporal(total = 1.5)
90
91 cons.w_m[0] *= 0.001
92 cons.mu_s[0] = 0.0
93 cons.mu_d[0] = 0.0
94 cons.gamma_wn[0] = 1e4
95 cons.gamma_wt[0] = 1e4
96 cons.contactmodel[0] = 1
97
98 if (consolidation == True):
99
100 # Run sphere
101 cons.run(dry = True) # show values, don't run
102 cons.run(device = device) # run
103
104 if (plots == True):
105 # Make a graph of energies
106 cons.visualize('energy')
107 cons.visualize('walls')
108
109 cons.writeVTKall()
110
111 if (rendering == True):
112 # Render images with raytracer
113 cons.render(method = "pres", max_val = 2.0*N, verbose = False)
114
115
116 ### SHEARING WITHOUT ROTATION ###
117
118 # New class
119 shear_nrt = sphere.sim(np = cons.np, nw = cons.nw, sid = sim_id +
120 "-shear_nrt-N{}".format(N))
121
122 # Read last output file of initialization step
123 lastf = sphere.status(sim_id + "-cons-N{}".format(N))
124 shear_nrt.readbin("../output/" + sim_id +
125 "-cons-N{}.output{:0=5}.bin".format(N, lastf),
126 verbose = False)
127
128 # Periodic x and y boundaries
129 shear_nrt.periodicBoundariesXY()
130
131 # Setup shear_nrt experiment
132 shear_nrt.shear(shear_strain_rate = 0.05)
133 shear_nrt.adaptiveGrid()
134
135 # Set duration of simulation
136 shear_nrt.initTemporal(total = 20.0)
137
138 # Disable rotation for regular grains
139 shear_nrt.fixvel[numpy.nonzero(shear_nrt.fixvel == 0.0)] = -10.0
140
141 if (shearing == True):
142
143 # Run sphere
144 shear_nrt.run(dry = True)
145 shear_nrt.run(device = device)
146
147 if (plots == True):
148 # Make a graph of energies
149 shear_nrt.visualize('energy')
150 shear_nrt.visualize('shear')
151
152 shear_nrt.writeVTKall()
153
154 if (rendering == True):
155 # Render images with raytracer
156 shear_nrt.render(method = "pres", max_val = 2.0*N, verbose = False)
157
158
159 ### SHEARING ###
160
161 # New class
162 shear = sphere.sim(np = cons.np, nw = cons.nw, sid = sim_id +
163 "-shear-N{}".format(N))
164
165 # Read last output file of initialization step
166 lastf = sphere.status(sim_id + "-cons-N{}".format(N))
167 shear.readbin("../output/" + sim_id +
168 "-cons-N{}.output{:0=5}.bin".format(N, lastf),
169 verbose = False)
170
171 # Periodic x and y boundaries
172 shear.periodicBoundariesXY()
173
174 # Setup shear experiment
175 shear.shear(shear_strain_rate = 0.05)
176 shear.adaptiveGrid()
177
178 # Set duration of simulation
179 shear.initTemporal(total = 20.0)
180
181 if (shearing == True):
182
183 # Run sphere
184 shear.run(dry = True)
185 shear.run(device = device)
186
187 if (plots == True):
188 # Make a graph of energies
189 shear.visualize('energy')
190 shear.visualize('shear')
191
192 shear.writeVTKall()
193
194 if (rendering == True):
195 # Render images with raytracer
196 shear.render(method = "pres", max_val = 2.0*N, verbose = False)