tFix continued line indentation - 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 439019dc191d4e17968abd9d9ea70164a2489b84
(DIR) parent 84e33bd162733d251f944d84b1d799bc6951f4b7
(HTM) Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Mon, 2 Sep 2019 15:20:47 +0200
Fix continued line indentation
Diffstat:
M python/sphere.py | 241 ++++++++++++++++++-------------
1 file changed, 138 insertions(+), 103 deletions(-)
---
(DIR) diff --git a/python/sphere.py b/python/sphere.py
t@@ -313,19 +313,19 @@ class sim:
# Fluid velocities [m/s]
self.v_f = numpy.zeros((self.num[0], self.num[1], self.num[2], self.nd),
- dtype=numpy.float64)
+ dtype=numpy.float64)
# Fluid pressures [Pa]
self.p_f = numpy.zeros((self.num[0], self.num[1], self.num[2]),
- dtype=numpy.float64)
+ dtype=numpy.float64)
# Fluid cell porosities [-]
self.phi = numpy.zeros((self.num[0], self.num[1], self.num[2]),
- dtype=numpy.float64)
+ dtype=numpy.float64)
# Fluid cell porosity change [1/s]
self.dphi = numpy.zeros((self.num[0], self.num[1], self.num[2]),
- dtype=numpy.float64)
+ dtype=numpy.float64)
# Fluid density [kg/(m^3)]
self.rho_f = numpy.ones(1, dtype=numpy.float64) * 1.0e3
t@@ -365,8 +365,8 @@ class sim:
# Hold pressures constant in fluid cell (0: True, 1: False)
self.p_f_constant = numpy.zeros((self.num[0],
- self.num[1],
- self.num[2]), dtype=numpy.int32)
+ self.num[1],
+ self.num[2]), dtype=numpy.int32)
# Navier-Stokes
if self.cfd_solver[0] == 0:
t@@ -1113,15 +1113,15 @@ class sim:
self.bonds[i, 0] = numpy.fromfile(fh, dtype=numpy.uint32, count=1)
self.bonds[i, 1] = numpy.fromfile(fh, dtype=numpy.uint32, count=1)
self.bonds_delta_n = numpy.fromfile(fh, dtype=numpy.float64,
- count=self.nb0)
+ count=self.nb0)
self.bonds_delta_t = numpy.fromfile(fh, dtype=numpy.float64,
- count=self.nb0*self.nd)\
- .reshape(self.nb0, self.nd)
+ count=self.nb0*self.nd)\
+ .reshape(self.nb0, self.nd)
self.bonds_omega_n = numpy.fromfile(fh, dtype=numpy.float64,
- count=self.nb0)
+ count=self.nb0)
self.bonds_omega_t = numpy.fromfile(fh, dtype=numpy.float64,
- count=self.nb0*self.nd)\
- .reshape(self.nb0, self.nd)
+ count=self.nb0*self.nd)\
+ .reshape(self.nb0, self.nd)
else:
self.nb0 = 0
t@@ -1135,101 +1135,126 @@ class sim:
self.mu = numpy.fromfile(fh, dtype=numpy.float64, count=1)
self.v_f = numpy.empty((self.num[0],
- self.num[1],
- self.num[2],
- self.nd), dtype=numpy.float64)
+ self.num[1],
+ self.num[2],
+ self.nd), dtype=numpy.float64)
self.p_f = numpy.empty((self.num[0],
- self.num[1],
- self.num[2]), dtype=numpy.float64)
+ self.num[1],
+ self.num[2]), dtype=numpy.float64)
self.phi = numpy.empty((self.num[0],
- self.num[1],
- self.num[2]), dtype=numpy.float64)
+ self.num[1],
+ self.num[2]), dtype=numpy.float64)
self.dphi = numpy.empty((self.num[0],
- self.num[1],
- self.num[2]), dtype=numpy.float64)
+ self.num[1],
+ self.num[2]), dtype=numpy.float64)
for z in numpy.arange(self.num[2]):
for y in numpy.arange(self.num[1]):
for x in numpy.arange(self.num[0]):
self.v_f[x, y, z, 0] = numpy.fromfile(fh,
- dtype=numpy.float64,
- count=1)
+ dtype=numpy.float64,
+ count=1)
self.v_f[x, y, z, 1] = numpy.fromfile(fh,
+ dtype=numpy.float64,
+ count=1)
+ self.v_f[x, y, z, 2] = numpy.fromfile(fh,
+ dtype=numpy.float64,
+ count=1)
+ self.p_f[x, y, z] = numpy.fromfile(fh,
dtype=numpy.float64,
count=1)
- self.v_f[x, y, z, 2] = numpy.fromfile(fh,
+ self.phi[x, y, z] = numpy.fromfile(fh,
dtype=numpy.float64,
count=1)
- self.p_f[x, y, z] = numpy.fromfile(fh, dtype=numpy.float64,
- count=1)
- self.phi[x, y, z] = numpy.fromfile(fh, dtype=numpy.float64,
- count=1)
- self.dphi[x, y, z] = numpy.fromfile(fh, dtype=numpy.float64,
- count=1)\
+ self.dphi[x, y, z] = numpy.fromfile(fh,
+ dtype=numpy.float64,
+ count=1)\
/(self.time_dt*self.ndem)
if self.version >= 0.36:
- self.rho_f = numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.p_mod_A = numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.p_mod_f = numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.p_mod_phi = numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ self.rho_f = numpy.fromfile(fh, dtype=numpy.float64,
+ count=1)
+ self.p_mod_A = numpy.fromfile(fh, dtype=numpy.float64,
+ count=1)
+ self.p_mod_f = numpy.fromfile(fh, dtype=numpy.float64,
+ count=1)
+ self.p_mod_phi = numpy.fromfile(fh, dtype=numpy.float64,
+ count=1)
if self.version >= 2.12 and self.cfd_solver[0] == 1:
- self.bc_xn = numpy.fromfile(fh, dtype=numpy.int32, count=1)
- self.bc_xp = numpy.fromfile(fh, dtype=numpy.int32, count=1)
- self.bc_yn = numpy.fromfile(fh, dtype=numpy.int32, count=1)
- self.bc_yp = numpy.fromfile(fh, dtype=numpy.int32, count=1)
+ self.bc_xn = numpy.fromfile(fh, dtype=numpy.int32,
+ count=1)
+ self.bc_xp = numpy.fromfile(fh, dtype=numpy.int32,
+ count=1)
+ self.bc_yn = numpy.fromfile(fh, dtype=numpy.int32,
+ count=1)
+ self.bc_yp = numpy.fromfile(fh, dtype=numpy.int32,
+ count=1)
self.bc_bot = numpy.fromfile(fh, dtype=numpy.int32, count=1)
self.bc_top = numpy.fromfile(fh, dtype=numpy.int32, count=1)
- self.free_slip_bot = numpy.fromfile(fh, dtype=numpy.int32, count=1)
- self.free_slip_top = numpy.fromfile(fh, dtype=numpy.int32, count=1)
+ self.free_slip_bot = numpy.fromfile(fh, dtype=numpy.int32,
+ count=1)
+ self.free_slip_top = numpy.fromfile(fh, dtype=numpy.int32,
+ count=1)
if self.version >= 2.11:
- self.bc_bot_flux = numpy.fromfile(fh, dtype=numpy.float64,
- count=1)
- self.bc_top_flux = numpy.fromfile(fh, dtype=numpy.float64,
- count=1)
+ self.bc_bot_flux = numpy.fromfile(fh,
+ dtype=numpy.float64,
+ count=1)
+ self.bc_top_flux = numpy.fromfile(fh,
+ dtype=numpy.float64,
+ count=1)
else:
self.bc_bot_flux = numpy.zeros(1, dtype=numpy.float64)
self.bc_top_flux = numpy.zeros(1, dtype=numpy.float64)
if self.version >= 2.15:
self.p_f_constant = numpy.empty((self.num[0],
- self.num[1],
- self.num[2]),
- dtype=numpy.int32)
+ self.num[1],
+ self.num[2]),
+ dtype=numpy.int32)
for z in numpy.arange(self.num[2]):
for y in numpy.arange(self.num[1]):
for x in numpy.arange(self.num[0]):
self.p_f_constant[x, y, z] = \
- numpy.fromfile(fh, dtype=numpy.int32, count=1)
+ numpy.fromfile(fh, dtype=numpy.int32,
+ count=1)
else:
self.p_f_constant = numpy.zeros((self.num[0],
- self.num[1],
- self.num[2]),
- dtype=numpy.int32)
+ self.num[1],
+ self.num[2]),
+ dtype=numpy.int32)
if self.version >= 2.0 and self.cfd_solver == 0:
- self.gamma = numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.theta = numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.beta = numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.tolerance = numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.maxiter = numpy.fromfile(fh, dtype=numpy.uint32, count=1)
+ self.gamma = numpy.fromfile(fh, dtype=numpy.float64,
+ count=1)
+ self.theta = numpy.fromfile(fh, dtype=numpy.float64,
+ count=1)
+ self.beta = numpy.fromfile(fh, dtype=numpy.float64,
+ count=1)
+ self.tolerance = numpy.fromfile(fh, dtype=numpy.float64,
+ count=1)
+ self.maxiter = numpy.fromfile(fh, dtype=numpy.uint32,
+ count=1)
if self.version >= 1.01:
- self.ndem = numpy.fromfile(fh, dtype=numpy.uint32, count=1)
+ self.ndem = numpy.fromfile(fh, dtype=numpy.uint32,
+ count=1)
else:
self.ndem = 1
if self.version >= 1.04:
- self.c_phi = numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.c_v = numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ self.c_phi = numpy.fromfile(fh, dtype=numpy.float64,
+ count=1)
+ self.c_v = numpy.fromfile(fh, dtype=numpy.float64,
+ count=1)
if self.version == 1.06:
self.c_a = numpy.fromfile(fh, dtype=numpy.float64,
- count=1)
+ count=1)
elif self.version >= 1.07:
- self.dt_dem_fac = numpy.fromfile(fh, dtype=numpy.float64,
- count=1)
+ self.dt_dem_fac = numpy.fromfile(fh,
+ dtype=numpy.float64,
+ count=1)
else:
self.c_a = numpy.ones(1, dtype=numpy.float64)
else:
t@@ -1243,42 +1268,51 @@ class sim:
self.f_sum = numpy.empty_like(self.x)
for i in numpy.arange(self.np):
- self.f_d[i, :] = numpy.fromfile(fh, dtype=numpy.float64,
- count=self.nd)
+ self.f_d[i, :] = numpy.fromfile(fh,
+ dtype=numpy.float64,
+ count=self.nd)
for i in numpy.arange(self.np):
- self.f_p[i, :] = numpy.fromfile(fh, dtype=numpy.float64,
- count=self.nd)
+ self.f_p[i, :] = numpy.fromfile(fh,
+ dtype=numpy.float64,
+ count=self.nd)
for i in numpy.arange(self.np):
- self.f_v[i, :] = numpy.fromfile(fh, dtype=numpy.float64,
- count=self.nd)
+ self.f_v[i, :] = numpy.fromfile(fh,
+ dtype=numpy.float64,
+ count=self.nd)
for i in numpy.arange(self.np):
- self.f_sum[i, :] = numpy.fromfile(fh, dtype=numpy.float64,
- count=self.nd)
+ self.f_sum[i, :] = numpy.fromfile(fh,
+ dtype=numpy.float64,
+ count=self.nd)
else:
self.f_d = numpy.zeros((self.np, self.nd),
- dtype=numpy.float64)
+ dtype=numpy.float64)
self.f_p = numpy.zeros((self.np, self.nd),
- dtype=numpy.float64)
+ dtype=numpy.float64)
self.f_v = numpy.zeros((self.np, self.nd),
- dtype=numpy.float64)
+ dtype=numpy.float64)
self.f_sum = numpy.zeros((self.np, self.nd),
- dtype=numpy.float64)
+ dtype=numpy.float64)
elif self.version >= 2.0 and self.cfd_solver == 1:
- self.tolerance = numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.maxiter = numpy.fromfile(fh, dtype=numpy.uint32, count=1)
+ self.tolerance = numpy.fromfile(fh, dtype=numpy.float64,
+ count=1)
+ self.maxiter = numpy.fromfile(fh, dtype=numpy.uint32,
+ count=1)
self.ndem = numpy.fromfile(fh, dtype=numpy.uint32, count=1)
- self.c_phi = numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ self.c_phi = numpy.fromfile(fh, dtype=numpy.float64,
+ count=1)
self.f_p = numpy.empty_like(self.x)
for i in numpy.arange(self.np):
self.f_p[i, :] = numpy.fromfile(fh, dtype=numpy.float64,
- count=self.nd)
- self.beta_f = numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ count=self.nd)
+ self.beta_f = numpy.fromfile(fh, dtype=numpy.float64,
+ count=1)
self.k_c = numpy.fromfile(fh, dtype=numpy.float64, count=1)
if self.version >= 1.02:
- self.color = numpy.fromfile(fh, dtype=numpy.int32, count=self.np)
+ self.color = numpy.fromfile(fh, dtype=numpy.int32,
+ count=self.np)
else:
self.color = numpy.zeros(self.np, dtype=numpy.int32)
t@@ -2181,7 +2215,7 @@ class sim:
def red(ratio):
return numpy.fmin(1.0, 0.209*ratio**3. - 2.49*ratio**2. + 3.0*ratio
- + 0.0109)
+ + 0.0109)
def green(ratio):
return numpy.fmin(1.0, -2.44*ratio**2. + 2.15*ratio + 0.369)
def blue(ratio):
t@@ -2650,12 +2684,12 @@ class sim:
# Max. and min. coordinates of world
self.origo = numpy.array([numpy.amin(self.x[:, 0] - self.radius[:]),
- numpy.amin(self.x[:, 1] - self.radius[:]),
- numpy.amin(self.x[:, 2] - self.radius[:])]) \
+ numpy.amin(self.x[:, 1] - self.radius[:]),
+ numpy.amin(self.x[:, 2] - self.radius[:])]) \
- margin*r_max
self.L = numpy.array([numpy.amax(self.x[:, 0] + self.radius[:]),
- numpy.amax(self.x[:, 1] + self.radius[:]),
- numpy.amax(self.x[:, 2] + self.radius[:])]) \
+ numpy.amax(self.x[:, 1] + self.radius[:]),
+ numpy.amax(self.x[:, 2] + self.radius[:])]) \
+ margin*r_max
cellsize_min = 2.1 * r_max
t@@ -2864,7 +2898,7 @@ class sim:
'''
bondparticles = numpy.unique(numpy.random.random_integers(0, high=self.np-1,
- size=int(self.np*ratio)))
+ size=int(self.np*ratio)))
if bondparticles.size % 2 > 0:
bondparticles = bondparticles[:-1].copy()
bondparticles = bondparticles.reshape(int(bondparticles.size/2),
t@@ -3026,8 +3060,8 @@ class sim:
self.w_vel = numpy.array([1, 0, 0, 0, 0]) * wvel
self.w_sigma0 = numpy.array([0, 1, 1, 1, 1]) * normal_stress
self.w_n = numpy.array(([0, 0, -1], [-1, 0, 0],
- [1, 0, 0], [0, -1, 0], [0, 1, 0]),
- dtype=numpy.float64)
+ [1, 0, 0], [0, -1, 0], [0, 1, 0]),
+ dtype=numpy.float64)
self.w_x = numpy.zeros(5)
self.w_m = numpy.zeros(5)
self.w_force = numpy.zeros(5)
t@@ -3070,7 +3104,7 @@ class sim:
# Fix horizontal velocity to 0.0 of lowermost particles
d_max_below = numpy.max(self.radius[numpy.nonzero(self.x[:, 2] <
- (z_max-z_min)*0.3)])*2.0
+ (z_max-z_min)*0.3)])*2.0
I = numpy.nonzero(self.x[:, 2] < (z_min + d_max_below))
self.fixvel[I] = 1
self.angvel[I, 0] = 0.0
t@@ -3082,7 +3116,7 @@ class sim:
# Fix horizontal velocity to specific value of uppermost particles
d_max_top = numpy.max(self.radius[numpy.nonzero(self.x[:, 2] >
- (z_max-z_min)*0.7)])*2.0
+ (z_max-z_min)*0.7)])*2.0
I = numpy.nonzero(self.x[:, 2] > (z_max - d_max_top))
self.fixvel[I] = 1
self.angvel[I, 0] = 0.0
t@@ -3349,7 +3383,7 @@ class sim:
self.rho_f = numpy.ones(1, dtype=numpy.float64) * rho
self.p_f = numpy.ones((self.num[0], self.num[1], self.num[2]),
- dtype=numpy.float64) * p
+ dtype=numpy.float64) * p
if hydrostatic:
t@@ -3374,11 +3408,11 @@ class sim:
self.v_f = numpy.zeros((self.num[0], self.num[1], self.num[2], self.nd),
- dtype=numpy.float64)
+ dtype=numpy.float64)
self.phi = numpy.ones((self.num[0], self.num[1], self.num[2]),
- dtype=numpy.float64)
+ dtype=numpy.float64)
self.dphi = numpy.zeros((self.num[0], self.num[1], self.num[2]),
- dtype=numpy.float64)
+ dtype=numpy.float64)
self.p_mod_A = numpy.zeros(1, dtype=numpy.float64) # Amplitude [Pa]
self.p_mod_f = numpy.zeros(1, dtype=numpy.float64) # Frequency [Hz]
t@@ -3392,7 +3426,7 @@ class sim:
self.bc_top_flux = numpy.zeros(1, dtype=numpy.float64)
self.p_f_constant = numpy.zeros((self.num[0], self.num[1], self.num[2]),
- dtype=numpy.int32)
+ dtype=numpy.int32)
# Fluid solver type
# 0: Navier Stokes (fluid with inertia)
t@@ -4018,7 +4052,7 @@ class sim:
self.bonds_delta_t = numpy.array([[0.0, 0.0, 0.0]], dtype=numpy.uint32)
else:
self.bonds_delta_t = numpy.vstack((self.bonds_delta_t,
- [0.0, 0.0, 0.0]))
+ [0.0, 0.0, 0.0]))
if not hasattr(self, 'bonds_omega_n'):
self.bonds_omega_n = numpy.array([0.0], dtype=numpy.uint32)
t@@ -4027,10 +4061,11 @@ class sim:
self.bonds_omega_n = numpy.append(self.bonds_omega_n, [0.0])
if not hasattr(self, 'bonds_omega_t'):
- self.bonds_omega_t = numpy.array([[0.0, 0.0, 0.0]], dtype=numpy.uint32)
+ self.bonds_omega_t = numpy.array([[0.0, 0.0, 0.0]],
+ dtype=numpy.uint32)
else:
self.bonds_omega_t = numpy.vstack((self.bonds_omega_t,
- [0.0, 0.0, 0.0]))
+ [0.0, 0.0, 0.0]))
# Increment the number of bonds with one
self.nb0 += 1
t@@ -4715,7 +4750,7 @@ class sim:
shell=True)
contactdata = numpy.loadtxt('../output/' + self.sid + '.contacts.txt')
self.pairs = numpy.array((contactdata[:, 0], contactdata[:, 1]),
- dtype=numpy.int32)
+ dtype=numpy.int32)
self.overlaps = numpy.array(contactdata[:, 2])
def findCoordinationNumber(self):
t@@ -5520,11 +5555,11 @@ class sim:
# plot defined max compressive stress from tau/N ratio
ax.scatter(0., # prescribed stress
numpy.degrees(numpy.arctan(self.shearStress('defined')/
- self.currentNormalStress('defined'))),
+ self.currentNormalStress('defined'))),
marker='o', c='none', edgecolor='blue', s=300)
ax.scatter(0., # actual stress
numpy.degrees(numpy.arctan(self.shearStress('effective')/
- self.currentNormalStress('effective'))),
+ self.currentNormalStress('effective'))),
marker='+', color='blue', s=300)
ax.set_rmax(90)
t@@ -5561,7 +5596,7 @@ class sim:
# angle: 0 when vertical, 90 when horizontal
#hist, bins = numpy.histogram(datadata[:, 6], bins=10)
_, _, _ = plt.hist(90. - diplist, bins=range(0, 100, 10),
- alpha=0.75, facecolor='gray')
+ alpha=0.75, facecolor='gray')
theta_sigma1 = numpy.degrees(numpy.arctan(
self.currentNormalStress('defined')/\
self.shearStress('defined')))
t@@ -6747,7 +6782,7 @@ class sim:
raise ValueError
s = numpy.r_[2*tau[0]-tau[smoothing:1:-1], tau,
- 2*tau[-1]-tau[-1:-smoothing:-1]]
+ 2*tau[-1]-tau[-1:-smoothing:-1]]
if smoothing_window == 'flat': # moving average
w = numpy.ones(smoothing, 'd')