tbond_tests.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
---
tbond_tests.py (6294B)
---
1 #!/usr/bin/env python
2 from pytestutils import *
3 import sphere
4
5 def printKinematics(sb):
6 print('bonds_delta_n'); print(sb.bonds_delta_n)
7 print('bonds_delta_t'); print(sb.bonds_delta_t)
8 print('bonds_omega_n'); print(sb.bonds_omega_n)
9 print('bonds_omega_t'); print(sb.bonds_omega_t)
10 print('force'); print(sb.force)
11 print('torque'); print(sb.torque)
12 print('vel'); print(sb.vel)
13 print('angvel'); print(sb.angvel)
14
15 #### Bond tests ####
16 print("### Bond tests ###")
17
18 # Zero arrays
19 z2_1 = numpy.zeros((2,1))
20 z2_2 = numpy.zeros((2,2))
21 z1_3 = numpy.zeros((1,3))
22 z2_3 = numpy.zeros((2,3))
23
24 # Small value arrays
25 smallval = 1e-8
26 s2_1 = numpy.ones((2,1))*smallval
27
28 # Inter-particle distances to try (neg. for overlap)
29 #distances = [0.2, 0.0, -0.2]
30 #distances = [0.2, 0.0]
31 distances = []
32 #distances = [0.2]
33
34 for d in distances:
35
36 radii = 0.5
37 print("## Inter-particle distance: " + str(d/radii) + " radii")
38
39 sb = sphere.sim(np=2, sid='bondtest')
40 cleanup(sb)
41
42 # setup particles, bond, and simulation
43 sb.x[0,:] = numpy.array((10.0, 10.0, 10.0))
44 sb.x[1,:] = numpy.array((10.0+2.0*radii+d, 10.0, 10.0))
45 sb.radius = numpy.ones(sb.np)*radii
46 sb.initGridAndWorldsize(margin = 10, periodic = 1, contactmodel = 2, g = numpy.array([0.0, 0.0, 0.0]))
47 sb.bond(0, 1)
48 sb.defaultParams(gamma_n = 0.0, gamma_t = 0.0)
49 #sb.initTemporal(total=0.5, file_dt=0.01)
50 #sb.render(verbose=False)
51 #visualize(sb.sid, "energy")
52
53
54 print("# Stability test")
55 sb.x[0,:] = numpy.array((10.0, 10.0, 10.0))
56 sb.x[1,:] = numpy.array((10.0+2.0*radii+d, 10.0, 10.0))
57 sb.zeroKinematics()
58 sb.initTemporal(total=0.2, file_dt=0.01)
59 #sb.initTemporal(total=0.01, file_dt=0.0001)
60 sb.run(verbose=False)
61 #sb.run()
62 sb.readlast(verbose=False)
63 compareFloats(0.0, sb.energy("kin") + sb.energy("rot") + sb.energy("bondpot"), "Energy cons.")
64 compareNumpyArrays(sb.vel, z2_3, "vel\t")
65 compareNumpyArrays(sb.angvel, z2_3, "angvel\t")
66 #printKinematics(sb)
67 #visualize(sb.sid, "energy")
68 #sb.readbin('../output/' + sb.sid + '.output00001.bin')
69 #printKinematics(sb)
70
71 print("# Normal expansion")
72 sb.x[0,:] = numpy.array((10.0, 10.0, 10.0))
73 sb.x[1,:] = numpy.array((10.0+2.0*radii+d, 10.0, 10.0))
74 sb.zeroKinematics()
75 sb.initTemporal(total=0.2, file_dt=0.01)
76 sb.vel[1,0] = 1e-4
77 Ekinrot0 = sb.energy("kin") + sb.energy("rot")
78 sb.run(verbose=False)
79 sb.readlast(verbose=False)
80 compareFloats(Ekinrot0, sb.energy("kin") + sb.energy("rot") + sb.energy("bondpot"), "Energy cons.")
81 print("vel[:,0]"),
82 if ((sb.vel[0,0] <= 0.0) or (sb.vel[1,0] <= 0.0)):
83 print(failed())
84 else :
85 print(passed())
86 compareNumpyArrays(sb.vel[:,1:2], z2_2, "vel[:,1:2]")
87 compareNumpyArrays(sb.angvel, z2_3, "angvel\t")
88 printKinematics(sb)
89 #visualize(sb.sid, "energy")
90 continue
91
92 print("# Normal contraction")
93 sb.zeroKinematics()
94 sb.initTemporal(total=0.2, file_dt=0.01)
95 sb.vel[1,0] = -1e-4
96 Ekinrot0 = sb.energy("kin") + sb.energy("rot")
97 sb.run(verbose=False)
98 sb.readlast(verbose=False)
99 compareFloats(Ekinrot0, sb.energy("kin") + sb.energy("rot") + sb.energy("bondpot"), "Energy cons.")
100 print("vel[:,0]"),
101 if ((sb.vel[0,0] >= 0.0) or (sb.vel[1,0] >= 0.0)):
102 print(failed())
103 else :
104 print(passed())
105 compareNumpyArrays(sb.vel[:,1:2], z2_2, "vel[:,1:2]")
106 compareNumpyArrays(sb.angvel, z2_3, "angvel\t")
107 #printKinematics(sb)
108
109 print("# Shear")
110 sb.x[0,:] = numpy.array((10.0, 10.0, 10.0))
111 sb.x[1,:] = numpy.array((10.0+2.0*radii+d, 10.0, 10.0))
112 sb.zeroKinematics()
113 sb.initTemporal(total=0.2, file_dt=0.01)
114 sb.vel[1,2] = 1e-4
115 Ekinrot0 = sb.energy("kin") + sb.energy("rot")
116 sb.run(verbose=False)
117 sb.readlast(verbose=False)
118 compareFloats(Ekinrot0, sb.energy("kin") + sb.energy("rot") + sb.energy("bondpot"), "Energy cons.")
119 #print("vel[:,0]"),
120 #if ((sb.vel[0,0] <= 0.0) or (sb.vel[1,0] >= 0.0)):
121 # print(failed())
122 #else :
123 # print(passed())
124 compareNumpyArrays(sb.vel[:,1], z2_1, "vel[:,1]")
125 print("vel[:,2]"),
126 if ((sb.vel[0,2] <= 0.0) or (sb.vel[1,2] <= 0.0)):
127 print(failed())
128 else :
129 print(passed())
130 #compareNumpyArrays(sb.angvel[:,0:2:2], z2_2, "angvel[:,0:2:2]")
131 #print("angvel[:,1]"),
132 #if ((sb.angvel[0,1] >= 0.0) or (sb.angvel[1,1] >= 0.0)):
133 # print(failed())
134 #else :
135 # print(passed())
136 compareNumpyArrays(sb.angvel, z2_3, "angvel\t")
137 #printKinematics(sb)
138 #visualize(sb.sid, "energy")
139
140
141 #'''
142 print("# Twist")
143 sb.x[0,:] = numpy.array((10.0, 10.0, 10.0))
144 sb.x[1,:] = numpy.array((10.0+2.0*radii+d, 10.0, 10.0))
145 sb.zeroKinematics()
146 sb.initTemporal(total=0.2, file_dt=0.01)
147 #sb.initTemporal(total=0.001, file_dt=0.00001)
148 sb.angvel[1,0] = 1e-4
149 Ekinrot0 = sb.energy("kin") + sb.energy("rot")
150 sb.run(verbose=False)
151 sb.readlast(verbose=False)
152 compareFloats(Ekinrot0, sb.energy("kin") + sb.energy("rot") + sb.energy("bondpot"), "Energy cons.")
153 #compareNumpyArrays(sb.vel, z2_3, "vel\t")
154 print("angvel[:,0]"),
155 if ((sb.angvel[0,0] <= 0.0) or (sb.angvel[1,0] <= 0.0)):
156 raise Exception("Failed")
157 print(failed())
158 else :
159 print(passed())
160 compareNumpyArrays(sb.angvel[:,1:2], z2_2, "angvel[:,1:2]")
161 #printKinematics(sb)
162 #visualize(sb.sid, "energy")
163
164
165 #'''
166 print("# Bend")
167 sb.x[0,:] = numpy.array((10.0, 10.0, 10.0))
168 sb.x[1,:] = numpy.array((10.0+2.0*radii+d, 10.0, 10.0))
169 sb.zeroKinematics()
170 sb.initTemporal(total=0.2, file_dt=0.01)
171 sb.angvel[0,1] = -1e-4
172 sb.angvel[1,1] = 1e-4
173 Ekinrot0 = sb.energy("kin") + sb.energy("rot")
174 sb.run(verbose=False)
175 sb.readlast(verbose=False)
176 compareFloats(Ekinrot0, sb.energy("kin") + sb.energy("rot") + sb.energy("bondpot"), "Energy cons.")
177 #compareNumpyArrays(sb.vel, z2_3, "vel\t")
178 #compareNumpyArrays(sb.angvel[:,0:2:2], z2_2, "angvel[:,0:2:2]")
179 print("angvel[:,1]"),
180 if ((sb.angvel[0,1] == 0.0) or (sb.angvel[1,1] == 0.0)):
181 raise Exception("Failed")
182 print(failed())
183 else :
184 print(passed())
185 #printKinematics(sb)
186 #visualize(sb.sid, "energy")
187 #'''
188