tsupraglacial-master.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
---
tsupraglacial-master.py (3732B)
---
1 #!/usr/bin/env python
2
3 # sphere grain/fluid simulation: https://src.adamsgaard.dk/sphere
4 import sphere
5 import numpy
6
7 ### EXPERIMENT SETUP ###
8 initialization = False
9 creeping = True
10 plots = True
11
12 # Common simulation id
13 sim_id = "supraglacial"
14
15 # Fluid-pressure gradient [Pa/m]
16 dpdz = 0.0
17
18 # Grain density
19 rho_g = 3600.0
20
21 # Fluid density
22 rho_f = 1000.0
23
24 # Gravitational acceleration
25 g = 9.8
26
27 # Slope
28 slope_angle = 20.0
29
30 # Number of particles
31 np = 1e4
32
33 device = 0 # automatically choose best GPU
34
35
36 ### INITIALIZATION ###
37
38 # New class
39 init = sphere.sim(np = np, nd = 3, nw = 0, sid = sim_id + "-init")
40
41 # Uniform diameters from 0.3 cm to 0.7 cm
42 init.generateRadii(psd = 'uni', mean = 0.0025, variance = 0.001)
43
44 # Use default params
45 init.defaultParams(gamma_n = 100.0, mu_s = 0.6, mu_d = 0.6)
46 init.setYoungsModulus(1e8)
47
48 # Disable wall viscosities
49 init.gamma_wn[0] = 0.0
50 init.gamma_wt[0] = 0.0
51
52 # Add gravity
53 init.g[2] = -g
54
55 # Periodic x and y boundaries
56 init.periodicBoundariesXY()
57
58 # Initialize positions in random grid (also sets world size)
59 hcells = np**(1.0/3.0)
60 init.initRandomGridPos(gridnum = [hcells-2, hcells-2, 1e9])
61
62 # Set duration of simulation
63 init.initTemporal(total = 5.0)
64
65 if (initialization == True):
66
67 # Run sphere
68 init.run(dry = True)
69 init.run(device=device)
70
71 if (plots == True):
72 # Make a graph of energies
73 init.visualize('energy')
74
75 init.writeVTKall()
76
77
78 ### CREEP ###
79
80 # New class
81 creep = sphere.sim(np = init.np,
82 sid = sim_id + "-slope{}-dpdz{}".format(slope_angle, dpdz))
83
84 # Read last output file of initialization step
85 creep.readbin("../output/" + sim_id + "-init.output{:0=5}.bin"
86 .format(init.status()))
87
88 # Tilt gravity
89 creep.g[2] = -g*numpy.cos(numpy.deg2rad(slope_angle))
90 creep.g[0] = g*numpy.sin(numpy.deg2rad(slope_angle))
91
92 # Disable particle contact viscosities
93 creep.gamma_n[0] = 0.0
94 creep.gamma_t[0] = 0.0
95
96 # zero all velocities and accelerations
97 creep.zeroKinematics()
98
99 # Periodic x and y boundaries
100 creep.periodicBoundariesXY()
101
102 # Fit grid to grains
103 creep.adjustUpperWall(z_adjust=1.2)
104 creep.nw = 0 # no dynamic wall on top
105
106 # Fix bottom grains
107 z_min = numpy.min(creep.x[:,2] - creep.radius)
108 z_max = numpy.max(creep.x[:,2] + creep.radius)
109 d_max_below = numpy.max(creep.radius[numpy.nonzero(creep.x[:,2] <
110 (z_max-z_min)*0.3)])*2.0
111 I = numpy.nonzero(creep.x[:,2] < (z_min + d_max_below))
112 creep.fixvel[I] = 1
113
114 # set fluid and solver properties
115 creep.initFluid(mu=8.9e-4, p=0.0, rho=rho_f, cfd_solver=1) # water at 25 C
116 creep.setMaxIterations(2e5)
117 creep.setPermeabilityGrainSize()
118 creep.setFluidCompressibility(4.6e-10) # water at 25 C
119
120 # set fluid BCs
121 # creep.setFluidTopNoFlow()
122 # creep.setFluidBottomNoFlow()
123 # creep.setFluidXFixedPressure()
124 # creep.adaptiveGrid()
125 creep.setFluidTopFixedPressure()
126 creep.setFluidBottomFixedPressure()
127 creep.setFluidXPeriodic()
128 creep.setFluidYPeriodic()
129
130 # set fluid pressures at the boundaries and internally
131 dz = creep.L[2]/creep.num[2]
132 for iz in range(creep.num[2]):
133 z = iz + 0.5*dz
134 creep.p_f[:,:,iz] = numpy.abs(creep.L[2]*dpdz) + z*dpdz
135
136 # Remove fixvel constraint from uppermost grains
137 #creep.fixvel[numpy.nonzero(creep.x[:,2] > 0.5*creep.L[2])] = 0
138
139 # Produce regular coloring pattern
140 creep.checkerboardColors(creep.num[0], creep.num[1], creep.num[2])
141 creep.color[numpy.nonzero(creep.fixvel == 1)] == -1
142
143 # Adapt grid size during progressive deformation
144 #creep.adaptiveGrid()
145
146 # Set duration of simulation
147 creep.initTemporal(total=5.0, file_dt=0.01)
148
149 if (creeping == True):
150
151 # Run sphere
152 creep.run(dry = True)
153 creep.run(device=device)
154
155 if (plots == True):
156 # Make a graph of energies
157 creep.visualize('energy')
158
159 creep.writeVTKall()