-
Notifications
You must be signed in to change notification settings - Fork 1
/
particles.py
190 lines (156 loc) · 6.64 KB
/
particles.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
from numpy import *
from pylab import *
# PARTICLES MODULE:
# This set of classes was developed to model interacting spherical particles in 3
# dimensions. There are presently two types of interaction forces possible. Each
# is a class
#
# GranularMaterialsForce: This is a force that represents elastic forces between
# particles that are overlapping, and a damping, or dissipative force. It does
# not include the often used rotational degrees of freedom. Also note that the
# constraint forces are taken care of with an extra function call.
#
# LennardJonesForce: This force describes the weak attraction at large distances
# and strong repulsion experienced at close distances by mono-atomic gases.
#
# The equation of motion are integrated by the Verlet method, which is presently
# the only integration scheme supported.
#
# All particle data (position, velocity, acceleration, radii, and distances) are
# stored in the Particles class.
# 9/23/08 JVJ and Tim Bocek
class GranularMaterialForce:
def __init__(self, k = 1.5, gamma = .3, g = .25):
# parameters in force model
self.__k = k # Elastic 'bounce'
self.__gamma = gamma # Energy dissipation/loss
self.__g = g # Gravity
def __call__(self, p):
# Find position differences
d, dx, dy, dz = p.distanceMatrix(p.x,p.y,p.z)
# Compute overlap
dr = d - triu(p.sumOfRadii) - tril(p.sumOfRadii)
# No forces arising in no overlap cases
dr[dr>0]=0
# Compute elastic particle/particle forces
magnitude = -self.__k * dr
# Velocity differences
dv, dvx, dvy, dvz = p.distanceMatrix(p.vx,p.vy,p.vz)
# Damping terms
vijDotrij = dvx * dx + dvy * dy + dvz * dz
vijDotrij[dr==0] = 0.
# Damping is subtracted from force
magnitude = magnitude - self.__gamma * vijDotrij / d
magnitude[d==-1] = 0
cx, cy, cz = self.floorConstraint(p)
# Project onto components, sum all forces on each particle
p.ax = sum(magnitude * (-dx/d) * p.ratioOfRadii,axis=1) + cx
p.ay = sum(magnitude * (-dy/d)* p.ratioOfRadii,axis=1) - self.__g + cy
p.az = sum(magnitude * (-dz/d)* p.ratioOfRadii,axis=1) + cz
def floorConstraint(self,p):
" "" This is a highly specific function for a floor that responds (elasticity and damping) the same way a particle does. Presently, if constraints are to change, one would have to rewrite the function."" "
effectiveRadius = 3. # This is how 'hard' the floor is
floorDistance = p.y + p.L/2 - p.r
floorDistance[floorDistance > 0] = 0
lowerWallForce = -self.__k * floorDistance
lowerWallDamping = -self.__gamma * p.vy * floorDistance
lowerWallForce = lowerWallForce - lowerWallDamping
cx = 0
cy = lowerWallForce * effectiveRadius / p.r
cz = 0
return cx, cy, cz
class VerletIntegrator:
def __init__(self,dt=0.01):
# Time step
self.__dt = dt
def __call__(self,force,p):
# Position update
p.x = p.x + p.vx * self.__dt + .5 * p.ax * self.__dt ** 2
p.y = p.y + p.vy * self.__dt + .5 * p.ay * self.__dt ** 2
p.z = p.z + p.vz * self.__dt + .5 * p.az * self.__dt ** 2
# Update periodic BC
p.pbcUpdate()
# Store accelerations for averaging that is done
ax=p.ax
ay=p.ay
az=p.az
force(p) # Force update with new positions
# Velocity updates
p.vx = p.vx + 0.5 * (ax + p.ax) * self.__dt
p.vy = p.vy + 0.5 * (ay + p.ay) * self.__dt
p.vz = p.vz + 0.5 * (az + p.az) * self.__dt
class Particles:
def __init__(self,L,force,periodicX = 1,periodicY=1,periodicZ=1):
# Container size
self.L = L
# Total Number of particles
self.N = 0
# type
self.type = 'float32'
# Positions
self.x = array([],dtype=self.type)
self.y = array([],dtype=self.type)
self.z = array([],dtype=self.type)
# Velocities
self.vx = array([],dtype=self.type)
self.vy = array([],dtype=self.type)
self.vz = array([],dtype=self.type)
# Forces
self.ax = array([],dtype=self.type)
self.ay = array([],dtype=self.type)
self.az = array([],dtype=self.type)
# Radii
self.r = array([],dtype=self.type)
# Periodic on?
self.periodicX = periodicX
self.periodicY = periodicY
self.periodicZ = periodicZ
# Force function
self.f = force
def addParticle(self,x,y,z,vx,vy,vz,r):
self.x = hstack((self.x,x))
self.y = hstack((self.y,y))
self.z = hstack((self.z,z))
self.vx = hstack((self.vx,vx))
self.vy = hstack((self.vy,vy))
self.vz = hstack((self.vz,vz))
self.r = hstack((self.r,r))
self.N = self.N+1
temp = tile(self.r,(self.N,1))
self.sumOfRadii = temp + temp.T
self.ratioOfRadii = temp / temp.T
self.f(self)
def pbcUpdate(self):
" ""Moves paricles across periodic boundary"" "
if self.periodicX:
self.x[self.x > self.L/2] = self.x[self.x > self.L/2] - self.L
self.x[self.x < -self.L/2] = self.x[self.x < -self.L/2] + self.L
if self.periodicY:
self.y[self.y > self.L/2] = self.y[self.y > self.L/2] - self.L
self.y[self.y < -self.L/2] = self.y[self.y < -self.L/2] + self.L
if self.periodicZ:
self.z[self.z > self.L/2] = self.z[self.z > self.L/2] - self.L
self.z[self.z < -self.L/2] = self.z[self.z < -self.L/2] + self.L
def distanceMatrix(self,x,y,z):
" ""Computes distances between all particles and places the result in a matrix such that the ij th matrix entry corresponds to the distance between particle i and j"" "
xtemp = tile(x,(self.N,1))
dx = xtemp - xtemp.T
ytemp = tile(y,(self.N,1))
dy = ytemp - ytemp.T
ztemp = tile(z,(self.N,1))
dz = ztemp - ztemp.T
# Particles 'feel' each other across the periodic boundaries
if self.periodicX:
dx[dx>self.L/2]=dx[dx > self.L/2]-self.L
dx[dx<-self.L/2]=dx[dx < -self.L/2]+self.L
if self.periodicY:
dy[dy>self.L/2]=dy[dy>self.L/2]-self.L
dy[dy<-self.L/2]=dy[dy<-self.L/2]+self.L
if self.periodicZ:
dz[dz>self.L/2]=dz[dz>self.L/2]-self.L
dz[dz<-self.L/2]=dz[dz<-self.L/2]+self.L
# Total Distances
d = sqrt(dx**2+dy**2+dz**2)
# Mark zero entries with negative 1 to avoid divergences
d[d==0] = -1
return d, dx, dy, dz