tBegin implementing Python3 compability - 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 93256b5eb3dd7b1cee22f18d94892fe676eebe5d
(DIR) parent 116df3bcf276dec7fca4fc44a8b91c39db229f66
(HTM) Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Thu, 15 Aug 2019 15:49:32 +0200
Begin implementing Python3 compability
Diffstat:
M python/sphere.py | 358 ++++++++++++++++----------------
M src/CMakeLists.txt | 8 ++++----
2 files changed, 183 insertions(+), 183 deletions(-)
---
(DIR) diff --git a/python/sphere.py b/python/sphere.py
t@@ -1,4 +1,4 @@
-#!/usr/bin/env python2.7
+#!/usr/bin/env python
import math
import numpy
import matplotlib
t@@ -66,10 +66,10 @@ class sim:
self.version = numpy.ones(1, dtype=numpy.float64)*VERSION
# The number of spatial dimensions. Values other that 3 do not work
- self.nd = numpy.ones(1, dtype=numpy.int32)*int(nd)
+ self.nd = int(nd)
# The number of particles
- self.np = numpy.ones(1, dtype=numpy.uint32)*int(np)
+ self.np = int(np)
# The simulation id (text string)
self.sid = sid
t@@ -108,52 +108,52 @@ class sim:
## Particle data
# Particle position vectors [m]
- self.x = numpy.zeros((self.np[0], self.nd[0]),
+ self.x = numpy.zeros((self.np, self.nd),
dtype=numpy.float64)
# Particle radii [m]
- self.radius = numpy.ones(self.np[0], dtype=numpy.float64)
+ self.radius = numpy.ones(self.np, dtype=numpy.float64)
# The sums of x and y movement [m]
- self.xyzsum = numpy.zeros((self.np[0], 3), dtype=numpy.float64)
+ self.xyzsum = numpy.zeros((self.np, 3), dtype=numpy.float64)
# The linear velocities [m/s]
- self.vel = numpy.zeros((self.np[0], self.nd[0]),
+ self.vel = numpy.zeros((self.np, self.nd),
dtype=numpy.float64)
# Fix the particle horizontal velocities? 0: No, 1: Yes
- self.fixvel = numpy.zeros(self.np[0], dtype=numpy.float64)
+ self.fixvel = numpy.zeros(self.np, dtype=numpy.float64)
# The linear force vectors [N]
- self.force = numpy.zeros((self.np[0], self.nd[0]),
+ self.force = numpy.zeros((self.np, self.nd),
dtype=numpy.float64)
# The angular position vectors [rad]
- self.angpos = numpy.zeros((self.np[0], self.nd[0]),
+ self.angpos = numpy.zeros((self.np, self.nd),
dtype=numpy.float64)
# The angular velocity vectors [rad/s]
- self.angvel = numpy.zeros((self.np[0], self.nd[0]),
+ self.angvel = numpy.zeros((self.np, self.nd),
dtype=numpy.float64)
# The torque vectors [N*m]
- self.torque = numpy.zeros((self.np[0], self.nd[0]),
+ self.torque = numpy.zeros((self.np, self.nd),
dtype=numpy.float64)
# The shear friction energy dissipation rates [W]
- self.es_dot = numpy.zeros(self.np[0], dtype=numpy.float64)
+ self.es_dot = numpy.zeros(self.np, dtype=numpy.float64)
# The total shear energy dissipations [J]
- self.es = numpy.zeros(self.np[0], dtype=numpy.float64)
+ self.es = numpy.zeros(self.np, dtype=numpy.float64)
# The viscous energy dissipation rates [W]
- self.ev_dot = numpy.zeros(self.np[0], dtype=numpy.float64)
+ self.ev_dot = numpy.zeros(self.np, dtype=numpy.float64)
# The total viscois energy dissipation [J]
- self.ev = numpy.zeros(self.np[0], dtype=numpy.float64)
+ self.ev = numpy.zeros(self.np, dtype=numpy.float64)
# The total particle pressures [Pa]
- self.p = numpy.zeros(self.np[0], dtype=numpy.float64)
+ self.p = numpy.zeros(self.np, dtype=numpy.float64)
# The gravitational acceleration vector [N*m/s]
self.g = numpy.array([0.0, 0.0, 0.0], dtype=numpy.float64)
t@@ -227,43 +227,43 @@ class sim:
# nw = 1: Uniaxial (also used for shear experiments)
# nw = 2: Biaxial
# nw = 5: Triaxial
- self.nw = numpy.ones(1, dtype=numpy.uint32) * int(nw)
+ self.nw = int(nw)
# Wall modes
# 0: Fixed
# 1: Normal stress condition
# 2: Normal velocity condition
# 3: Normal stress and shear stress condition
- self.wmode = numpy.zeros(self.nw[0], dtype=numpy.int32)
+ self.wmode = numpy.zeros(self.nw, dtype=numpy.int32)
# Wall normals
- self.w_n = numpy.zeros((self.nw[0], self.nd[0]),
+ self.w_n = numpy.zeros((self.nw, self.nd),
dtype=numpy.float64)
- if self.nw[0] >= 1:
+ if self.nw >= 1:
self.w_n[0,2] = -1.0
- if self.nw[0] >= 2:
+ if self.nw >= 2:
self.w_n[1,0] = -1.0
- if self.nw[0] >= 3:
+ if self.nw >= 3:
self.w_n[2,0] = 1.0
- if self.nw[0] >= 4:
+ if self.nw >= 4:
self.w_n[3,1] = -1.0
- if self.nw[0] >= 5:
+ if self.nw >= 5:
self.w_n[4,1] = 1.0
# Wall positions on the axes that are parallel to the wall normal [m]
- self.w_x = numpy.ones(self.nw[0], dtype=numpy.float64)
+ self.w_x = numpy.ones(self.nw, dtype=numpy.float64)
# Wall masses [kg]
- self.w_m = numpy.zeros(self.nw[0], dtype=numpy.float64)
+ self.w_m = numpy.zeros(self.nw, dtype=numpy.float64)
# Wall velocities on the axes that are parallel to the wall normal [m/s]
- self.w_vel = numpy.zeros(self.nw[0], dtype=numpy.float64)
+ self.w_vel = numpy.zeros(self.nw, dtype=numpy.float64)
# Wall forces on the axes that are parallel to the wall normal [m/s]
- self.w_force = numpy.zeros(self.nw[0], dtype=numpy.float64)
+ self.w_force = numpy.zeros(self.nw, dtype=numpy.float64)
# Wall stress on the axes that are parallel to the wall normal [Pa]
- self.w_sigma0 = numpy.zeros(self.nw[0], dtype=numpy.float64)
+ self.w_sigma0 = numpy.zeros(self.nw, dtype=numpy.float64)
# Wall stress modulation amplitude [Pa]
self.w_sigma0_A = numpy.zeros(1, dtype=numpy.float64)
t@@ -279,7 +279,7 @@ class sim:
self.lambda_bar = numpy.ones(1, dtype=numpy.float64)
# Number of bonds
- self.nb0 = numpy.zeros(1, dtype=numpy.uint32)
+ self.nb0 = 0
# Bond tensile strength [Pa]
self.sigma_b = numpy.ones(1, dtype=numpy.float64) * numpy.infty
t@@ -288,20 +288,20 @@ class sim:
self.tau_b = numpy.ones(1, dtype=numpy.float64) * numpy.infty
# Bond pairs
- self.bonds = numpy.zeros((self.nb0[0], 2), dtype=numpy.uint32)
+ self.bonds = numpy.zeros((self.nb0, 2), dtype=numpy.uint32)
# Parallel bond movement
- self.bonds_delta_n = numpy.zeros(self.nb0[0], dtype=numpy.float64)
+ self.bonds_delta_n = numpy.zeros(self.nb0, dtype=numpy.float64)
# Shear bond movement
- self.bonds_delta_t = numpy.zeros((self.nb0[0], self.nd[0]),
+ self.bonds_delta_t = numpy.zeros((self.nb0, self.nd),
dtype=numpy.float64)
# Twisting bond movement
- self.bonds_omega_n = numpy.zeros(self.nb0[0], dtype=numpy.float64)
+ self.bonds_omega_n = numpy.zeros(self.nb0, dtype=numpy.float64)
# Bending bond movement
- self.bonds_omega_t = numpy.zeros((self.nb0[0], self.nd[0]),
+ self.bonds_omega_t = numpy.zeros((self.nb0, self.nd),
dtype=numpy.float64)
## Fluid parameters
t@@ -321,7 +321,7 @@ class sim:
# Fluid velocities [m/s]
self.v_f = numpy.zeros(
- (self.num[0], self.num[1], self.num[2], self.nd[0]),
+ (self.num[0], self.num[1], self.num[2], self.nd),
dtype=numpy.float64)
# Fluid pressures [Pa]
t@@ -411,13 +411,13 @@ class sim:
self.dt_dem_fac = numpy.ones(1, dtype=numpy.float64)
## Interaction forces
- self.f_d = numpy.zeros((self.np[0], self.nd[0]),
+ self.f_d = numpy.zeros((self.np, self.nd),
dtype=numpy.float64)
- self.f_p = numpy.zeros((self.np[0], self.nd[0]),
+ self.f_p = numpy.zeros((self.np, self.nd),
dtype=numpy.float64)
- self.f_v = numpy.zeros((self.np[0], self.nd[0]),
+ self.f_v = numpy.zeros((self.np, self.nd),
dtype=numpy.float64)
- self.f_sum = numpy.zeros((self.np[0], self.nd[0]),
+ self.f_sum = numpy.zeros((self.np, self.nd),
dtype=numpy.float64)
# Darcy
t@@ -436,7 +436,7 @@ class sim:
self.c_phi = numpy.ones(1, dtype=numpy.float64)
# Interaction forces
- self.f_p = numpy.zeros((self.np[0], self.nd[0]),
+ self.f_p = numpy.zeros((self.np, self.nd),
dtype=numpy.float64)
# Adiabatic fluid compressibility [1/Pa].
t@@ -451,7 +451,7 @@ class sim:
str(self.cfd_solver[0]) + ')')
# Particle color marker
- self.color = numpy.zeros(self.np[0], dtype=numpy.int32)
+ self.color = numpy.zeros(self.np, dtype=numpy.int32)
def __cmp__(self, other):
'''
t@@ -941,26 +941,26 @@ class sim:
'''
Deletes all particles in the simulation object.
'''
- self.np[0] = 0
- self.x = numpy.zeros((self.np[0], self.nd[0]), dtype=numpy.float64)
- self.radius = numpy.ones(self.np[0], dtype=numpy.float64)
- self.xyzsum = numpy.zeros((self.np[0], 3), dtype=numpy.float64)
- self.vel = numpy.zeros((self.np[0], self.nd[0]), dtype=numpy.float64)
- self.fixvel = numpy.zeros(self.np[0], dtype=numpy.float64)
- self.force = numpy.zeros((self.np[0], self.nd[0]), dtype=numpy.float64)
- self.angpos = numpy.zeros((self.np[0], self.nd[0]), dtype=numpy.float64)
- self.angvel = numpy.zeros((self.np[0], self.nd[0]), dtype=numpy.float64)
- self.torque = numpy.zeros((self.np[0], self.nd[0]), dtype=numpy.float64)
- self.es_dot = numpy.zeros(self.np[0], dtype=numpy.float64)
- self.es = numpy.zeros(self.np[0], dtype=numpy.float64)
- self.ev_dot = numpy.zeros(self.np[0], dtype=numpy.float64)
- self.ev = numpy.zeros(self.np[0], dtype=numpy.float64)
- self.p = numpy.zeros(self.np[0], dtype=numpy.float64)
- self.color = numpy.zeros(self.np[0], dtype=numpy.int32)
- self.f_d = numpy.zeros((self.np[0], self.nd[0]), dtype=numpy.float64)
- self.f_p = numpy.zeros((self.np[0], self.nd[0]), dtype=numpy.float64)
- self.f_v = numpy.zeros((self.np[0], self.nd[0]), dtype=numpy.float64)
- self.f_sum = numpy.zeros((self.np[0], self.nd[0]), dtype=numpy.float64)
+ self.np = 0
+ self.x = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
+ self.radius = numpy.ones(self.np, dtype=numpy.float64)
+ self.xyzsum = numpy.zeros((self.np, 3), dtype=numpy.float64)
+ self.vel = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
+ self.fixvel = numpy.zeros(self.np, dtype=numpy.float64)
+ self.force = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
+ self.angpos = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
+ self.angvel = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
+ self.torque = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
+ self.es_dot = numpy.zeros(self.np, dtype=numpy.float64)
+ self.es = numpy.zeros(self.np, dtype=numpy.float64)
+ self.ev_dot = numpy.zeros(self.np, dtype=numpy.float64)
+ self.ev = numpy.zeros(self.np, dtype=numpy.float64)
+ self.p = numpy.zeros(self.np, dtype=numpy.float64)
+ self.color = numpy.zeros(self.np, dtype=numpy.int32)
+ self.f_d = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
+ self.f_p = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
+ self.f_v = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
+ self.f_sum = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
def readbin(self, targetbin, verbose = True, bonds = True, sigma0mod = True,
esysparticle = False):
t@@ -997,8 +997,8 @@ class sim:
self.version = numpy.fromfile(fh, dtype=numpy.float64, count=1)
# Read the number of dimensions and particles
- self.nd = numpy.fromfile(fh, dtype=numpy.int32, count=1)
- self.np = numpy.fromfile(fh, dtype=numpy.uint32, count=1)
+ self.nd = int(numpy.fromfile(fh, dtype=numpy.int32, count=1))
+ self.np = int(numpy.fromfile(fh, dtype=numpy.uint32, count=1))
# Read the time variables
self.time_dt = \
t@@ -1013,21 +1013,21 @@ class sim:
numpy.fromfile(fh, dtype=numpy.uint32, count=1)
# Allocate array memory for particles
- self.x = numpy.empty((self.np[0], self.nd[0]), dtype=numpy.float64)
- self.radius = numpy.empty(self.np[0], dtype=numpy.float64)
- self.xyzsum = numpy.empty((self.np[0], 3), dtype=numpy.float64)
- self.vel = numpy.empty((self.np[0], self.nd[0]), dtype=numpy.float64)
- self.fixvel = numpy.empty(self.np[0], dtype=numpy.float64)
- self.es_dot = numpy.empty(self.np[0], dtype=numpy.float64)
- self.es = numpy.empty(self.np[0], dtype=numpy.float64)
- self.ev_dot = numpy.empty(self.np[0], dtype=numpy.float64)
- self.ev = numpy.empty(self.np[0], dtype=numpy.float64)
- self.p = numpy.empty(self.np[0], dtype=numpy.float64)
+ self.x = numpy.empty((self.np, self.nd), dtype=numpy.float64)
+ self.radius = numpy.empty(self.np, dtype=numpy.float64)
+ self.xyzsum = numpy.empty((self.np, 3), dtype=numpy.float64)
+ self.vel = numpy.empty((self.np, self.nd), dtype=numpy.float64)
+ self.fixvel = numpy.empty(self.np, dtype=numpy.float64)
+ self.es_dot = numpy.empty(self.np, dtype=numpy.float64)
+ self.es = numpy.empty(self.np, dtype=numpy.float64)
+ self.ev_dot = numpy.empty(self.np, dtype=numpy.float64)
+ self.ev = numpy.empty(self.np, dtype=numpy.float64)
+ self.p = numpy.empty(self.np, dtype=numpy.float64)
# Read remaining data from binary
- self.origo = numpy.fromfile(fh, dtype=numpy.float64, count=self.nd[0])
- self.L = numpy.fromfile(fh, dtype=numpy.float64, count=self.nd[0])
- self.num = numpy.fromfile(fh, dtype=numpy.uint32, count=self.nd[0])
+ self.origo = numpy.fromfile(fh, dtype=numpy.float64, count=self.nd)
+ self.L = numpy.fromfile(fh, dtype=numpy.float64, count=self.nd)
+ self.num = numpy.fromfile(fh, dtype=numpy.uint32, count=self.nd)
self.periodic = numpy.fromfile(fh, dtype=numpy.int32, count=1)
if self.version >= 2.14:
t@@ -1038,32 +1038,32 @@ class sim:
# Per-particle vectors
for i in numpy.arange(self.np):
self.x[i,:] =\
- numpy.fromfile(fh, dtype=numpy.float64, count=self.nd[0])
+ numpy.fromfile(fh, dtype=numpy.float64, count=self.nd)
self.radius[i] =\
numpy.fromfile(fh, dtype=numpy.float64, count=1)
if self.version >= 1.03:
self.xyzsum = numpy.fromfile(fh, dtype=numpy.float64,\
- count=self.np*3).reshape(self.np[0],3)
+ count=self.np*3).reshape(self.np,3)
else:
self.xyzsum = numpy.fromfile(fh, dtype=numpy.float64,\
- count=self.np*2).reshape(self.np[0],2)
+ count=self.np*2).reshape(self.np,2)
for i in numpy.arange(self.np):
self.vel[i,:] =\
- numpy.fromfile(fh, dtype=numpy.float64, count=self.nd[0])
+ numpy.fromfile(fh, dtype=numpy.float64, count=self.nd)
self.fixvel[i] =\
numpy.fromfile(fh, dtype=numpy.float64, count=1)
self.force = numpy.fromfile(fh, dtype=numpy.float64,\
- count=self.np*self.nd[0]).reshape(self.np[0], self.nd[0])
+ count=self.np*self.nd).reshape(self.np, self.nd)
self.angpos = numpy.fromfile(fh, dtype=numpy.float64,\
- count=self.np*self.nd[0]).reshape(self.np[0], self.nd[0])
+ count=self.np*self.nd).reshape(self.np, self.nd)
self.angvel = numpy.fromfile(fh, dtype=numpy.float64,\
- count=self.np*self.nd[0]).reshape(self.np[0], self.nd[0])
+ count=self.np*self.nd).reshape(self.np, self.nd)
self.torque = numpy.fromfile(fh, dtype=numpy.float64,\
- count=self.np*self.nd[0]).reshape(self.np[0], self.nd[0])
+ count=self.np*self.nd).reshape(self.np, self.nd)
if esysparticle:
return
t@@ -1076,7 +1076,7 @@ class sim:
self.p = numpy.fromfile(fh, dtype=numpy.float64, count=self.np)
# Constant, global physical parameters
- self.g = numpy.fromfile(fh, dtype=numpy.float64, count=self.nd[0])
+ self.g = numpy.fromfile(fh, dtype=numpy.float64, count=self.nd)
self.k_n = numpy.fromfile(fh, dtype=numpy.float64, count=1)
self.k_t = numpy.fromfile(fh, dtype=numpy.float64, count=1)
self.k_r = numpy.fromfile(fh, dtype=numpy.float64, count=1)
t@@ -1101,22 +1101,22 @@ class sim:
self.V_b = numpy.fromfile(fh, dtype=numpy.float64, count=1)
# Wall data
- self.nw = numpy.fromfile(fh, dtype=numpy.uint32, count=1)
- self.wmode = numpy.empty(self.nw[0], dtype=numpy.int32)
+ self.nw = int(numpy.fromfile(fh, dtype=numpy.uint32, count=1))
+ self.wmode = numpy.empty(self.nw, dtype=numpy.int32)
self.w_n = numpy.empty(self.nw*self.nd, dtype=numpy.float64)\
- .reshape(self.nw[0],self.nd[0])
- self.w_x = numpy.empty(self.nw[0], dtype=numpy.float64)
- self.w_m = numpy.empty(self.nw[0], dtype=numpy.float64)
- self.w_vel = numpy.empty(self.nw[0], dtype=numpy.float64)
- self.w_force = numpy.empty(self.nw[0], dtype=numpy.float64)
- self.w_sigma0 = numpy.empty(self.nw[0], dtype=numpy.float64)
-
- self.wmode = numpy.fromfile(fh, dtype=numpy.int32, count=self.nw[0])
- for i in numpy.arange(self.nw[0]):
+ .reshape(self.nw,self.nd)
+ self.w_x = numpy.empty(self.nw, dtype=numpy.float64)
+ self.w_m = numpy.empty(self.nw, dtype=numpy.float64)
+ self.w_vel = numpy.empty(self.nw, dtype=numpy.float64)
+ self.w_force = numpy.empty(self.nw, dtype=numpy.float64)
+ self.w_sigma0 = numpy.empty(self.nw, dtype=numpy.float64)
+
+ self.wmode = numpy.fromfile(fh, dtype=numpy.int32, count=self.nw)
+ for i in numpy.arange(self.nw):
self.w_n[i,:] =\
- numpy.fromfile(fh, dtype=numpy.float64, count=self.nd[0])
+ numpy.fromfile(fh, dtype=numpy.float64, count=self.nd)
self.w_x[i] = numpy.fromfile(fh, dtype=numpy.float64, count=1)
- for i in numpy.arange(self.nw[0]):
+ for i in numpy.arange(self.nw):
self.w_m[i] = numpy.fromfile(fh, dtype=numpy.float64, count=1)
self.w_vel[i] = numpy.fromfile(fh, dtype=numpy.float64, count=1)
self.w_force[i] =\
t@@ -1134,25 +1134,25 @@ class sim:
# Inter-particle bonds
self.lambda_bar =\
numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.nb0 = numpy.fromfile(fh, dtype=numpy.uint32, count=1)
+ self.nb0 = int(numpy.fromfile(fh, dtype=numpy.uint32, count=1))
self.sigma_b = numpy.fromfile(fh, dtype=numpy.float64, count=1)
self.tau_b = numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.bonds = numpy.empty((self.nb0[0], 2), dtype=numpy.uint32)
- for i in numpy.arange(self.nb0[0]):
+ self.bonds = numpy.empty((self.nb0, 2), dtype=numpy.uint32)
+ for i in numpy.arange(self.nb0):
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[0])
+ count=self.nb0)
self.bonds_delta_t = numpy.fromfile(fh, dtype=numpy.float64,
- count=self.nb0*self.nd[0]).reshape(self.nb0[0], self.nd[0])
+ count=self.nb0*self.nd).reshape(self.nb0, self.nd)
self.bonds_omega_n = numpy.fromfile(fh, dtype=numpy.float64,
- count=self.nb0[0])
+ count=self.nb0)
self.bonds_omega_t = numpy.fromfile(fh, dtype=numpy.float64,
- count=self.nb0*self.nd[0]).reshape(self.nb0[0], self.nd[0])
+ count=self.nb0*self.nd).reshape(self.nb0, self.nd)
else:
- self.nb0 = numpy.zeros(1, dtype=numpy.uint32)
+ self.nb0 = 0
if self.fluid:
t@@ -1165,7 +1165,7 @@ 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[0]),
+ (self.num[0], 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]),
t@@ -1295,30 +1295,30 @@ class sim:
self.f_v = numpy.empty_like(self.x)
self.f_sum = numpy.empty_like(self.x)
- for i in numpy.arange(self.np[0]):
+ for i in numpy.arange(self.np):
self.f_d[i,:] = \
numpy.fromfile(fh, dtype=numpy.float64,
- count=self.nd[0])
- for i in numpy.arange(self.np[0]):
+ count=self.nd)
+ for i in numpy.arange(self.np):
self.f_p[i,:] = \
numpy.fromfile(fh, dtype=numpy.float64,
- count=self.nd[0])
- for i in numpy.arange(self.np[0]):
+ count=self.nd)
+ for i in numpy.arange(self.np):
self.f_v[i,:] = \
numpy.fromfile(fh, dtype=numpy.float64,
- count=self.nd[0])
- for i in numpy.arange(self.np[0]):
+ count=self.nd)
+ for i in numpy.arange(self.np):
self.f_sum[i,:] = \
numpy.fromfile(fh, dtype=numpy.float64,
- count=self.nd[0])
+ count=self.nd)
else:
- self.f_d = numpy.zeros((self.np[0], self.nd[0]),
+ self.f_d = numpy.zeros((self.np, self.nd),
dtype=numpy.float64)
- self.f_p = numpy.zeros((self.np[0], self.nd[0]),
+ self.f_p = numpy.zeros((self.np, self.nd),
dtype=numpy.float64)
- self.f_v = numpy.zeros((self.np[0], self.nd[0]),
+ self.f_v = numpy.zeros((self.np, self.nd),
dtype=numpy.float64)
- self.f_sum = numpy.zeros((self.np[0], self.nd[0]),
+ self.f_sum = numpy.zeros((self.np, self.nd),
dtype=numpy.float64)
elif self.version >= 2.0 and self.cfd_solver == 1:
t@@ -1331,10 +1331,10 @@ class sim:
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[0]):
+ for i in numpy.arange(self.np):
self.f_p[i,:] = \
numpy.fromfile(fh, dtype=numpy.float64,
- count=self.nd[0])
+ count=self.nd)
self.beta_f = \
numpy.fromfile(fh, dtype=numpy.float64, count=1)
t@@ -1345,7 +1345,7 @@ class sim:
self.color =\
numpy.fromfile(fh, dtype=numpy.int32, count=self.np)
else:
- self.color = numpy.zeros(self.np[0], dtype=numpy.int32)
+ self.color = numpy.zeros(self.np, dtype=numpy.int32)
finally:
self.version[0] = VERSION
t@@ -1376,8 +1376,8 @@ class sim:
fh.write(self.version.astype(numpy.float64))
# Write the number of dimensions and particles
- fh.write(self.nd.astype(numpy.int32))
- fh.write(self.np.astype(numpy.uint32))
+ fh.write(numpy.array(self.nd).astype(numpy.int32))
+ fh.write(numpy.array(self.np).astype(numpy.uint32))
# Write the time variables
fh.write(self.time_dt.astype(numpy.float64))
t@@ -1398,14 +1398,14 @@ class sim:
fh.write(self.x[i,:].astype(numpy.float64))
fh.write(self.radius[i].astype(numpy.float64))
- if self.np[0] > 0:
+ if self.np > 0:
fh.write(self.xyzsum.astype(numpy.float64))
for i in numpy.arange(self.np):
fh.write(self.vel[i,:].astype(numpy.float64))
fh.write(self.fixvel[i].astype(numpy.float64))
- if self.np[0] > 0:
+ if self.np > 0:
fh.write(self.force.astype(numpy.float64))
fh.write(self.angpos.astype(numpy.float64))
t@@ -1440,14 +1440,14 @@ class sim:
fh.write(self.db.astype(numpy.float64))
fh.write(self.V_b.astype(numpy.float64))
- fh.write(self.nw.astype(numpy.uint32))
- for i in numpy.arange(self.nw[0]):
+ fh.write(numpy.array(self.nw).astype(numpy.uint32))
+ for i in numpy.arange(self.nw):
fh.write(self.wmode[i].astype(numpy.int32))
- for i in numpy.arange(self.nw[0]):
+ for i in numpy.arange(self.nw):
fh.write(self.w_n[i,:].astype(numpy.float64))
fh.write(self.w_x[i].astype(numpy.float64))
- for i in numpy.arange(self.nw[0]):
+ for i in numpy.arange(self.nw):
fh.write(self.w_m[i].astype(numpy.float64))
fh.write(self.w_vel[i].astype(numpy.float64))
fh.write(self.w_force[i].astype(numpy.float64))
t@@ -1457,10 +1457,10 @@ class sim:
fh.write(self.w_tau_x.astype(numpy.float64))
fh.write(self.lambda_bar.astype(numpy.float64))
- fh.write(self.nb0.astype(numpy.uint32))
+ fh.write(numpy.array(self.nb0).astype(numpy.uint32))
fh.write(self.sigma_b.astype(numpy.float64))
fh.write(self.tau_b.astype(numpy.float64))
- for i in numpy.arange(self.nb0[0]):
+ for i in numpy.arange(self.nb0):
fh.write(self.bonds[i,0].astype(numpy.uint32))
fh.write(self.bonds[i,1].astype(numpy.uint32))
fh.write(self.bonds_delta_n.astype(numpy.float64))
t@@ -1613,7 +1613,7 @@ class sim:
fn = "../output/{0}.output{1:0=5}.bin".format(self.sid, i)
sb.sid = self.sid + ".{:0=5}".format(i)
sb.readbin(fn, verbose = False)
- if sb.np[0] > 0:
+ if sb.np > 0:
if i == 0:
sb.writeVTK(verbose=verbose)
if forces:
t@@ -1692,7 +1692,7 @@ class sim:
+ 'byte_order="LittleEndian">\n') # VTK header
fh.write(' <UnstructuredGrid>\n')
fh.write(' <Piece NumberOfPoints="%d" NumberOfCells="0">\n' \
- % (self.np[0]))
+ % (self.np))
# Coordinates for each point (positions)
fh.write(' <Points>\n')
t@@ -2400,7 +2400,7 @@ class sim:
if histogram:
fig = plt.figure(figsize=(8,8))
figtitle = 'Particle size distribution, {0} particles'.format(\
- self.np[0])
+ self.np)
fig.text(0.5,0.95,figtitle,horizontalalignment='center',\
fontproperties=FontProperties(size=18))
bins = 20
t@@ -2438,7 +2438,7 @@ class sim:
V_small = V_sphere(r_small)
V_large = V_sphere(r_large)
- nlarge = int(V_small/V_large * ratio * self.np[0]) # ignore void volume
+ nlarge = int(V_small/V_large * ratio * self.np) # ignore void volume
self.radius[:] = r_small
self.radius[0:nlarge] = r_large
t@@ -2452,7 +2452,7 @@ class sim:
if verbose:
print("generateBimodalRadii created " + str(nlarge)
- + " large particles, and " + str(self.np[0] - nlarge)
+ + " large particles, and " + str(self.np - nlarge)
+ " small")
def checkerboardColors(self, nx = 6, ny = 6, nz = 6):
t@@ -2498,7 +2498,7 @@ class sim:
:returns: z cell index
:return type: int
'''
- if self.nw[0] > 0:
+ if self.nw > 0:
return int(self.w_x[0]/(self.L[2]/self.num[2]))
else:
raise Exception('No dynamic top wall present!')
t@@ -2585,7 +2585,7 @@ class sim:
overlaps = False
# Draw random position
- for d in range(self.nd[0]):
+ for d in range(self.nd):
self.x[i,d] = (self.L[d] - self.origo[d] - 2*r_max) \
* numpy.random.random_sample() \
+ self.origo[d] + r_max
t@@ -2598,7 +2598,7 @@ class sim:
if delta_len < 0.0:
overlaps = True
print("\rFinding non-overlapping particle positions, "
- + "{0} % complete".format(numpy.ceil(i/self.np[0]*100)))
+ + "{0} % complete".format(numpy.ceil(i/self.np*100)))
# Print newline
print()
t@@ -2625,7 +2625,7 @@ class sim:
if dx > 0.0:
cellsize_min = dx
else:
- if self.np[0] < 1:
+ if self.np < 1:
raise Exception('Error: You need to define dx in '
'defineWorldBoundaries if there are no particles in '
'the simulation.')
t@@ -2768,7 +2768,7 @@ class sim:
gridpos[2] = numpy.floor(i/((self.num[0])*(self.num[1]))) #\
#% ((self.num[0])*(self.num[1]))
- for d in range(self.nd[0]):
+ for d in range(self.nd):
self.x[i,d] = gridpos[d] * cellsize + 0.5*cellsize
# Allow pushing every 2.nd level out of lateral boundaries
t@@ -2829,7 +2829,7 @@ class sim:
# Place particles in grid structure, and randomly adjust the
# positions within the oversized cells (uniform distribution)
- for d in range(self.nd[0]):
+ for d in range(self.nd):
r = self.radius[i]*1.05
self.x[i,d] = gridpos[d] * cellsize \
+ ((cellsize-r) - r) * numpy.random.random_sample() + r
t@@ -2933,18 +2933,18 @@ class sim:
when output from one simulation is reused in another simulation.
'''
- self.force = numpy.zeros((self.np[0], self.nd[0]))
- self.torque = numpy.zeros((self.np[0], self.nd[0]))
+ self.force = numpy.zeros((self.np, self.nd))
+ self.torque = numpy.zeros((self.np, self.nd))
self.vel = numpy.zeros(self.np*self.nd, dtype=numpy.float64)\
- .reshape(self.np[0], self.nd[0])
+ .reshape(self.np, self.nd)
self.angvel = numpy.zeros(self.np*self.nd, dtype=numpy.float64)\
- .reshape(self.np[0], self.nd[0])
+ .reshape(self.np, self.nd)
self.angpos = numpy.zeros(self.np*self.nd, dtype=numpy.float64)\
- .reshape(self.np[0], self.nd[0])
- self.es = numpy.zeros(self.np[0], dtype=numpy.float64)
- self.ev = numpy.zeros(self.np[0], dtype=numpy.float64)
+ .reshape(self.np, self.nd)
+ self.es = numpy.zeros(self.np, dtype=numpy.float64)
+ self.ev = numpy.zeros(self.np, dtype=numpy.float64)
self.xyzsum = numpy.zeros(self.np*3, dtype=numpy.float64)\
- .reshape(self.np[0], 3)
+ .reshape(self.np, 3)
def adjustUpperWall(self, z_adjust = 1.1):
'''
t@@ -2959,7 +2959,7 @@ class sim:
self.nw = numpy.ones(1, dtype=numpy.int32)
self.wmode = numpy.zeros(1) # fixed BC
self.w_n = numpy.zeros(self.nw*self.nd, dtype=numpy.float64).reshape(\
- self.nw[0],self.nd[0])
+ self.nw,self.nd)
self.w_n[0,2] = -1.0
self.w_vel = numpy.zeros(1)
self.w_force = numpy.zeros(1)
t@@ -3025,7 +3025,7 @@ class sim:
:type normal_stress: float
'''
- self.nw[0] = 1
+ self.nw = 1
if normal_stress <= 0.0:
raise Exception('consolidate() error: The normal stress should be '
t@@ -3082,7 +3082,7 @@ class sim:
self.zeroKinematics()
# Initialize walls
- self.nw[0] = 5 # five dynamic walls
+ self.nw = 5 # five dynamic walls
self.wmode = numpy.array([2,1,1,1,1]) # BCs (vel, stress, stress, ...)
self.w_vel = numpy.array([1,0,0,0,0]) * wvel
self.w_sigma0 = numpy.array([0,1,1,1,1]) * normal_stress
t@@ -3109,7 +3109,7 @@ class sim:
:type shear_stress: float or bool
'''
- self.nw[0] = 1
+ self.nw = 1
# Find lowest and heighest point
z_min = numpy.min(self.x[:,2] - self.radius)
t@@ -3335,11 +3335,11 @@ class sim:
if dt > 0.0:
self.time_dt[0] = dt
- if self.np[0] > 0:
+ if self.np > 0:
print("Warning: Manually specifying the time step length when "
+ "simulating particles may produce instabilities.")
- elif self.np[0] > 0:
+ elif self.np > 0:
r_min = numpy.min(self.radius)
m_min = self.rho * 4.0/3.0*numpy.pi*r_min**3
t@@ -3444,7 +3444,7 @@ class sim:
self.p_f[:,:,iz] = p + (depth-dz) * rho * -self.g[2]
- self.v_f = numpy.zeros((self.num[0], self.num[1], self.num[2], self.nd[0]),
+ self.v_f = numpy.zeros((self.num[0], self.num[1], self.num[2], self.nd),
dtype=numpy.float64)
self.phi = numpy.ones((self.num[0], self.num[1], self.num[2]),
dtype=numpy.float64)
t@@ -3483,19 +3483,19 @@ class sim:
self.c_v = numpy.ones(1, dtype=numpy.float64)
self.dt_dem_fac = numpy.ones(1, dtype=numpy.float64)
- self.f_d = numpy.zeros((self.np[0], self.nd[0]), dtype=numpy.float64)
- self.f_p = numpy.zeros((self.np[0], self.nd[0]), dtype=numpy.float64)
- self.f_v = numpy.zeros((self.np[0], self.nd[0]), dtype=numpy.float64)
- self.f_sum = numpy.zeros((self.np[0], self.nd[0]), dtype=numpy.float64)
+ self.f_d = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
+ self.f_p = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
+ self.f_v = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
+ self.f_sum = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
elif self.cfd_solver[0] == 1:
self.tolerance = numpy.array(1.0e-3)
self.maxiter = numpy.array(1e4)
self.ndem = numpy.array(1)
self.c_phi = numpy.ones(1, dtype=numpy.float64)
- self.f_d = numpy.zeros((self.np[0], self.nd[0]), dtype=numpy.float64)
+ self.f_d = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
self.beta_f = numpy.ones(1, dtype=numpy.float64)*4.5e-10
- self.f_p = numpy.zeros((self.np[0], self.nd[0]), dtype=numpy.float64)
+ self.f_p = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
self.k_c = numpy.ones(1, dtype=numpy.float64)*4.6e-10
self.bc_xn = numpy.ones(1, dtype=numpy.int32)*2
t@@ -3822,7 +3822,7 @@ class sim:
# Contact normal viscosity. Critical damping: 2*sqrt(m*k_n).
# Normal force component elastic if nu = 0.0.
- #self.gamma_n = numpy.ones(self.np[0], dtype=numpy.float64) \
+ #self.gamma_n = numpy.ones(self.np, dtype=numpy.float64) \
# * nu_frac * 2.0 * math.sqrt(4.0/3.0 * math.pi \
# * numpy.amin(self.radius)**3 \
# * self.rho[0] * self.k_n[0])
t@@ -4888,7 +4888,7 @@ class sim:
particle). Requires a previous call to :func:`findOverlaps()`. Values
are stored in ``self.coordinationnumber``.
'''
- self.coordinationnumber = numpy.zeros(self.np[0], dtype=numpy.int)
+ self.coordinationnumber = numpy.zeros(self.np, dtype=numpy.int)
for i in numpy.arange(self.overlaps.size):
self.coordinationnumber[self.pairs[0,i]] += 1
self.coordinationnumber[self.pairs[1,i]] += 1
t@@ -5134,7 +5134,7 @@ class sim:
# contacts
strikelist = [] # strike direction of the normal vector, [0:360[
diplist = [] # dip of the normal vector, [0:90]
- for n in numpy.arange(self.nb0[0]):
+ for n in numpy.arange(self.nb0):
i = self.bonds[n,0]
j = self.bonds[n,1]
t@@ -5725,7 +5725,7 @@ class sim:
plt.axvline(90. - theta_sigma1, color='k', linestyle='dashed',
linewidth=1)
plt.xlim([0, 90.])
- plt.ylim([0, self.np[0]/10])
+ plt.ylim([0, self.np/10])
#plt.xlabel('$\\boldsymbol{f}_\text{n}$ [N]')
plt.xlabel('Contact angle [deg]')
plt.ylabel('Count $N$')
t@@ -6477,7 +6477,7 @@ class sim:
'''
lastfile = self.status()
- sb = sim(sid = self.sid, np = self.np[0], nw = self.nw[0], fluid = self.fluid)
+ sb = sim(sid = self.sid, np = self.np, nw = self.nw, fluid = self.fluid)
### Plotting
if outformat != 'txt':
t@@ -6618,14 +6618,14 @@ class sim:
# Allocate arrays on first run
if i == firststep:
- wforce = numpy.zeros((lastfile+1)*sb.nw[0],\
- dtype=numpy.float64).reshape((lastfile+1), sb.nw[0])
- wvel = numpy.zeros((lastfile+1)*sb.nw[0],\
- dtype=numpy.float64).reshape((lastfile+1), sb.nw[0])
- wpos = numpy.zeros((lastfile+1)*sb.nw[0],\
- dtype=numpy.float64).reshape((lastfile+1), sb.nw[0])
- wsigma0 = numpy.zeros((lastfile+1)*sb.nw[0],\
- dtype=numpy.float64).reshape((lastfile+1), sb.nw[0])
+ wforce = numpy.zeros((lastfile+1)*sb.nw,\
+ dtype=numpy.float64).reshape((lastfile+1), sb.nw)
+ wvel = numpy.zeros((lastfile+1)*sb.nw,\
+ dtype=numpy.float64).reshape((lastfile+1), sb.nw)
+ wpos = numpy.zeros((lastfile+1)*sb.nw,\
+ dtype=numpy.float64).reshape((lastfile+1), sb.nw)
+ wsigma0 = numpy.zeros((lastfile+1)*sb.nw,\
+ dtype=numpy.float64).reshape((lastfile+1), sb.nw)
maxpos = numpy.zeros((lastfile+1), dtype=numpy.float64)
logstress = numpy.zeros((lastfile+1), dtype=numpy.float64)
voidratio = numpy.zeros((lastfile+1), dtype=numpy.float64)
(DIR) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
t@@ -19,15 +19,15 @@ INCLUDE(FindCUDA)
# NOTE: Multiple arguments must be semi-colon selimited
IF (GPU_GENERATION EQUAL 1) # Kepler
SET(CUDA_NVCC_FLAGS
- #"--use_fast_math;-O3;-gencode=arch=compute_35,code=\"sm_35,compute_35\";--fmad=false -ccbin gcc-4.6")
+ "--use_fast_math;-O3;-gencode=arch=compute_35,code=\"sm_35,compute_35\";--fmad=false -ccbin gcc-8.3")
#"--use_fast_math;-O3;-gencode=arch=compute_35,code=\"sm_35,compute_35\";--fmad=false -ccbin gcc-4.6")
#"--use_fast_math;-O3;-gencode=arch=compute_35,code=\"sm_35,compute_35\";--fmad=false -ccbin gcc;-Xcompiler -fPIC")
- "--use_fast_math;-O3;-gencode=arch=compute_35,code=\"sm_35,compute_35\";--fmad=false;-ccbin clang-3.8")
+ #"--use_fast_math;-O3;-gencode=arch=compute_35,code=\"sm_35,compute_35\";--fmad=false;-ccbin clang-3.8")
ELSE() # Fermi
SET(CUDA_NVCC_FLAGS
- #"--use_fast_math;-O3;-gencode=arch=compute_20,code=\"sm_20,compute_20\";--fmad=false -ccbin gcc-4.6")
+ "--use_fast_math;-O3;-gencode=arch=compute_20,code=\"sm_20,compute_20\";--fmad=false -ccbin gcc-8.3")
#"--use_fast_math;-O3;-gencode=arch=compute_20,code=\"sm_20,compute_20\";--fmad=false -ccbin gcc;-Xcompiler -fPIC")
- "--use_fast_math;-O3;-gencode=arch=compute_20,code=\"sm_20,compute_20\";--fmad=false;-ccbin clang-3.8")
+ #"--use_fast_math;-O3;-gencode=arch=compute_20,code=\"sm_20,compute_20\";--fmad=false;-ccbin clang-3.8")
ENDIF (GPU_GENERATION EQUAL 1)
SET(CMAKE_CXX_FLAGS "-fPIC ${CMAKE_CXX_FLAGS}")