tPerform many stylistic improvements based on pylint - 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 84e33bd162733d251f944d84b1d799bc6951f4b7
(DIR) parent 626a05132dc7966c24a0f1326342c2c738015d14
(HTM) Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Mon, 2 Sep 2019 14:45:51 +0200
Perform many stylistic improvements based on pylint
Diffstat:
M python/sphere.py | 3899 +++++++++++++++----------------
1 file changed, 1871 insertions(+), 2028 deletions(-)
---
(DIR) diff --git a/python/sphere.py b/python/sphere.py
t@@ -1,6 +1,8 @@
#!/usr/bin/env python
import math
import os
+import subprocess
+import pickle as pl
import numpy
try:
import matplotlib
t@@ -9,32 +11,30 @@ try:
import matplotlib.collections
matplotlib.rcParams.update({'font.size': 7, 'font.family': 'serif'})
matplotlib.rc('text', usetex=True)
- matplotlib.rcParams['text.latex.preamble']=[r"\usepackage{amsmath}"]
+ matplotlib.rcParams['text.latex.preamble'] = [r"\usepackage{amsmath}"]
from matplotlib.font_manager import FontProperties
- py_mpl=True
+ py_mpl = True
except ImportError:
print('Info: Could not find "matplotlib" python module. ' +
'Plotting functionality will be unavailable')
- py_mpl=False
-import subprocess
-import pickle as pl
+ py_mpl = False
try:
import vtk
- py_vtk=True
+ py_vtk = True
except ImportError:
print('Info: Could not find "vtk" python module. ' +
'Fluid VTK calls will be unavailable')
print('Consider installing with `pip install --user vtk`')
- py_vtk=False
+ py_vtk = False
numpy.seterr(all='warn', over='raise')
# Sphere version number. This field should correspond to the value in
# `../src/version.h`.
-VERSION=2.15
+VERSION = 2.15
# Transparency on plot legends
-legend_alpha=0.5
+legend_alpha = 0.5
class sim:
t@@ -61,295 +61,279 @@ class sim:
:type cfd_solver: int
'''
- def __init__(self,
- sid='unnamed',
- np=0,
- nd=3,
- nw=0,
- fluid=False,
- cfd_solver=0):
+ def __init__(self, sid='unnamed', np=0, nd=3, nw=0, fluid=False):
# Sphere version number
- self.version=numpy.ones(1, dtype=numpy.float64)*VERSION
+ self.version = numpy.ones(1, dtype=numpy.float64)*VERSION
# The number of spatial dimensions. Values other that 3 do not work
- self.nd=int(nd)
+ self.nd = int(nd)
# The number of particles
- self.np=int(np)
+ self.np = int(np)
# The simulation id (text string)
- self.sid=sid
+ self.sid = sid
## Time parameters
# Computational time step length [s]
- self.time_dt =numpy.zeros(1, dtype=numpy.float64)
+ self.time_dt = numpy.zeros(1, dtype=numpy.float64)
# Current time [s]
- self.time_current =numpy.zeros(1, dtype=numpy.float64)
+ self.time_current = numpy.zeros(1, dtype=numpy.float64)
# Total time [s]
- self.time_total =numpy.zeros(1, dtype=numpy.float64)
+ self.time_total = numpy.zeros(1, dtype=numpy.float64)
# File output interval [s]
- self.time_file_dt =numpy.zeros(1, dtype=numpy.float64)
+ self.time_file_dt = numpy.zeros(1, dtype=numpy.float64)
# The number of files written
- self.time_step_count=numpy.zeros(1, dtype=numpy.uint32)
+ self.time_step_count = numpy.zeros(1, dtype=numpy.uint32)
## World dimensions and grid data
# The Euclidean coordinate to the origo of the sorting grid
- self.origo =numpy.zeros(self.nd, dtype=numpy.float64)
+ self.origo = numpy.zeros(self.nd, dtype=numpy.float64)
- # The sorting grid size (x,y,z)
- self.L =numpy.zeros(self.nd, dtype=numpy.float64)
+ # The sorting grid size (x, y, z)
+ self.L = numpy.zeros(self.nd, dtype=numpy.float64)
# The number of sorting cells in each dimension
- self.num =numpy.zeros(self.nd, dtype=numpy.uint32)
+ self.num = numpy.zeros(self.nd, dtype=numpy.uint32)
# Whether to treat the lateral boundaries as periodic (1) or not (0)
- self.periodic=numpy.zeros(1, dtype=numpy.uint32)
+ self.periodic = numpy.zeros(1, dtype=numpy.uint32)
# Adaptively resize grid to assemblage height (0: no, 1: yes)
- self.adaptive=numpy.zeros(1, dtype=numpy.uint32)
+ self.adaptive = numpy.zeros(1, dtype=numpy.uint32)
## Particle data
# Particle position vectors [m]
- self.x =numpy.zeros((self.np, self.nd),
- dtype=numpy.float64)
+ self.x = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
# Particle radii [m]
- self.radius =numpy.ones(self.np, 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, 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, self.nd),
- dtype=numpy.float64)
+ 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, dtype=numpy.float64)
+ self.fixvel = numpy.zeros(self.np, dtype=numpy.float64)
# The linear force vectors [N]
- self.force =numpy.zeros((self.np, self.nd),
- dtype=numpy.float64)
+ self.force = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
# The angular position vectors [rad]
- self.angpos =numpy.zeros((self.np, self.nd),
- dtype=numpy.float64)
+ self.angpos = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
# The angular velocity vectors [rad/s]
- self.angvel =numpy.zeros((self.np, self.nd),
- dtype=numpy.float64)
+ self.angvel = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
# The torque vectors [N*m]
- self.torque =numpy.zeros((self.np, self.nd),
- dtype=numpy.float64)
+ 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, 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, 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, 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, dtype=numpy.float64)
+ self.ev = numpy.zeros(self.np, dtype=numpy.float64)
# The total particle pressures [Pa]
- self.p =numpy.zeros(self.np, 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)
+ self.g = numpy.array([0.0, 0.0, 0.0], dtype=numpy.float64)
# The Hookean coefficient for elastic stiffness normal to the contacts
# [N/m]
- self.k_n =numpy.ones(1, dtype=numpy.float64) * 1.16e9
+ self.k_n = numpy.ones(1, dtype=numpy.float64) * 1.16e9
# The Hookean coefficient for elastic stiffness tangential to the
# contacts [N/m]
- self.k_t =numpy.ones(1, dtype=numpy.float64) * 1.16e9
+ self.k_t = numpy.ones(1, dtype=numpy.float64) * 1.16e9
# The Hookean coefficient for elastic stiffness opposite of contact
# rotations. UNUSED
- self.k_r =numpy.zeros(1, dtype=numpy.float64)
+ self.k_r = numpy.zeros(1, dtype=numpy.float64)
# Young's modulus for contact stiffness [Pa]. This value is used
# instead of the Hookean stiffnesses (k_n, k_t) when self.E is larger
# than 0.0.
- self.E =numpy.zeros(1, dtype=numpy.float64)
+ self.E = numpy.zeros(1, dtype=numpy.float64)
# The viscosity normal to the contact [N/(m/s)]
- self.gamma_n =numpy.zeros(1, dtype=numpy.float64)
+ self.gamma_n = numpy.zeros(1, dtype=numpy.float64)
# The viscosity tangential to the contact [N/(m/s)]
- self.gamma_t =numpy.zeros(1, dtype=numpy.float64)
+ self.gamma_t = numpy.zeros(1, dtype=numpy.float64)
# The viscosity to contact rotation [N/(m/s)]
- self.gamma_r =numpy.zeros(1, dtype=numpy.float64)
+ self.gamma_r = numpy.zeros(1, dtype=numpy.float64)
# The coefficient of static friction on the contact [-]
- self.mu_s =numpy.ones(1, dtype=numpy.float64) * 0.5
+ self.mu_s = numpy.ones(1, dtype=numpy.float64) * 0.5
# The coefficient of dynamic friction on the contact [-]
- self.mu_d =numpy.ones(1, dtype=numpy.float64) * 0.5
+ self.mu_d = numpy.ones(1, dtype=numpy.float64) * 0.5
# The coefficient of rotational friction on the contact [-]
- self.mu_r =numpy.zeros(1, dtype=numpy.float64)
+ self.mu_r = numpy.zeros(1, dtype=numpy.float64)
# The viscosity normal to the walls [N/(m/s)]
- self.gamma_wn=numpy.zeros(1, dtype=numpy.float64)
+ self.gamma_wn = numpy.zeros(1, dtype=numpy.float64)
# The viscosity tangential to the walls [N/(m/s)]
- self.gamma_wt=numpy.zeros(1, dtype=numpy.float64)
+ self.gamma_wt = numpy.zeros(1, dtype=numpy.float64)
# The coeffient of static friction of the walls [-]
- self.mu_ws =numpy.ones(1, dtype=numpy.float64) * 0.5
+ self.mu_ws = numpy.ones(1, dtype=numpy.float64) * 0.5
# The coeffient of dynamic friction of the walls [-]
- self.mu_wd =numpy.ones(1, dtype=numpy.float64) * 0.5
+ self.mu_wd = numpy.ones(1, dtype=numpy.float64) * 0.5
# The particle density [kg/(m^3)]
- self.rho =numpy.ones(1, dtype=numpy.float64) * 2600.0
+ self.rho = numpy.ones(1, dtype=numpy.float64) * 2600.0
# The contact model to use
# 1: Normal: elasto-viscous, tangential: visco-frictional
# 2: Normal: elasto-viscous, tangential: elasto-visco-frictional
- self.contactmodel=numpy.ones(1, dtype=numpy.uint32) * 2 # lin-visc-el
+ self.contactmodel = numpy.ones(1, dtype=numpy.uint32) * 2 # lin-visc-el
# Capillary bond prefactor
- self.kappa =numpy.zeros(1, dtype=numpy.float64)
+ self.kappa = numpy.zeros(1, dtype=numpy.float64)
# Capillary bond debonding distance [m]
- self.db =numpy.zeros(1, dtype=numpy.float64)
+ self.db = numpy.zeros(1, dtype=numpy.float64)
# Capillary bond liquid volume [m^3]
- self.V_b =numpy.zeros(1, dtype=numpy.float64)
+ self.V_b = numpy.zeros(1, dtype=numpy.float64)
## Wall data
# Number of dynamic walls
# nw=1: Uniaxial (also used for shear experiments)
# nw=2: Biaxial
# nw=5: Triaxial
- self.nw =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, dtype=numpy.int32)
+ self.wmode = numpy.zeros(self.nw, dtype=numpy.int32)
# Wall normals
- self.w_n =numpy.zeros((self.nw, self.nd),
- dtype=numpy.float64)
+ self.w_n = numpy.zeros((self.nw, self.nd), dtype=numpy.float64)
if self.nw >= 1:
- self.w_n[0,2]=-1.0
+ self.w_n[0, 2] = -1.0
if self.nw >= 2:
- self.w_n[1,0]=-1.0
+ self.w_n[1, 0] = -1.0
if self.nw >= 3:
- self.w_n[2,0]= 1.0
+ self.w_n[2, 0] = 1.0
if self.nw >= 4:
- self.w_n[3,1]=-1.0
+ self.w_n[3, 1] = -1.0
if self.nw >= 5:
- self.w_n[4,1]= 1.0
+ 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, dtype=numpy.float64)
+ self.w_x = numpy.ones(self.nw, dtype=numpy.float64)
# Wall masses [kg]
- self.w_m =numpy.zeros(self.nw, 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, 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, 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, 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)
+ self.w_sigma0_A = numpy.zeros(1, dtype=numpy.float64)
# Wall stress modulation frequency [Hz]
- self.w_sigma0_f=numpy.zeros(1, dtype=numpy.float64)
+ self.w_sigma0_f = numpy.zeros(1, dtype=numpy.float64)
# Wall shear stress, enforced when wmode == 3
- self.w_tau_x=numpy.zeros(1, dtype=numpy.float64)
+ self.w_tau_x = numpy.zeros(1, dtype=numpy.float64)
## Bond parameters
# Radius multiplier to the parallel-bond radii
- self.lambda_bar=numpy.ones(1, dtype=numpy.float64)
+ self.lambda_bar = numpy.ones(1, dtype=numpy.float64)
# Number of bonds
- self.nb0=0
+ self.nb0 = 0
# Bond tensile strength [Pa]
- self.sigma_b=numpy.ones(1, dtype=numpy.float64) * numpy.infty
+ self.sigma_b = numpy.ones(1, dtype=numpy.float64) * numpy.infty
# Bond shear strength [Pa]
- self.tau_b=numpy.ones(1, dtype=numpy.float64) * numpy.infty
+ self.tau_b = numpy.ones(1, dtype=numpy.float64) * numpy.infty
# Bond pairs
- self.bonds=numpy.zeros((self.nb0, 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, 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, self.nd),
- dtype=numpy.float64)
+ self.bonds_delta_t = numpy.zeros((self.nb0, self.nd), dtype=numpy.float64)
# Twisting bond movement
- self.bonds_omega_n=numpy.zeros(self.nb0, 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, self.nd),
- dtype=numpy.float64)
+ self.bonds_omega_t = numpy.zeros((self.nb0, self.nd), dtype=numpy.float64)
## Fluid parameters
# Simulate fluid? True: Yes, False: no
- self.fluid=fluid
+ self.fluid = fluid
if self.fluid:
# Fluid solver type
# 0: Navier Stokes (fluid with inertia)
# 1: Stokes-Darcy (fluid without inertia)
- self.cfd_solver=numpy.zeros(1, dtype=numpy.int32)
+ self.cfd_solver = numpy.zeros(1, dtype=numpy.int32)
# Fluid dynamic viscosity [N/(m/s)]
- self.mu=numpy.zeros(1, dtype=numpy.float64)
+ self.mu = numpy.zeros(1, dtype=numpy.float64)
# Fluid velocities [m/s]
- self.v_f=numpy.zeros(
- (self.num[0], self.num[1], self.num[2], self.nd),
- dtype=numpy.float64)
+ self.v_f = numpy.zeros((self.num[0], self.num[1], self.num[2], self.nd),
+ dtype=numpy.float64)
# Fluid pressures [Pa]
- self.p_f=numpy.zeros((self.num[0], self.num[1], self.num[2]),
- dtype=numpy.float64)
+ self.p_f = numpy.zeros((self.num[0], self.num[1], self.num[2]),
+ dtype=numpy.float64)
# Fluid cell porosities [-]
- self.phi=numpy.zeros((self.num[0], self.num[1], self.num[2]),
- dtype=numpy.float64)
+ self.phi = numpy.zeros((self.num[0], self.num[1], self.num[2]),
+ 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)
+ self.dphi = numpy.zeros((self.num[0], self.num[1], self.num[2]),
+ dtype=numpy.float64)
# Fluid density [kg/(m^3)]
- self.rho_f=numpy.ones(1, dtype=numpy.float64) * 1.0e3
+ self.rho_f = numpy.ones(1, dtype=numpy.float64) * 1.0e3
# Pressure modulation at the top boundary
- self.p_mod_A=numpy.zeros(1, dtype=numpy.float64) # Amplitude [Pa]
- self.p_mod_f=numpy.zeros(1, dtype=numpy.float64) # Frequency [Hz]
- self.p_mod_phi=numpy.zeros(1, dtype=numpy.float64) # Shift [rad]
+ self.p_mod_A = numpy.zeros(1, dtype=numpy.float64) # Amplitude [Pa]
+ self.p_mod_f = numpy.zeros(1, dtype=numpy.float64) # Frequency [Hz]
+ self.p_mod_phi = numpy.zeros(1, dtype=numpy.float64) # Shift [rad]
## Fluid solver parameters
t@@ -358,10 +342,10 @@ class sim:
# 0: Dirichlet
# 1: Neumann
# 2: Periodic (default)
- self.bc_xn=numpy.ones(1, dtype=numpy.int32)*2 # Neg. x bc
- self.bc_xp=numpy.ones(1, dtype=numpy.int32)*2 # Pos. x bc
- self.bc_yn=numpy.ones(1, dtype=numpy.int32)*2 # Neg. y bc
- self.bc_yp=numpy.ones(1, dtype=numpy.int32)*2 # Pos. y bc
+ self.bc_xn = numpy.ones(1, dtype=numpy.int32)*2 # Neg. x bc
+ self.bc_xp = numpy.ones(1, dtype=numpy.int32)*2 # Pos. x bc
+ self.bc_yn = numpy.ones(1, dtype=numpy.int32)*2 # Neg. y bc
+ self.bc_yp = numpy.ones(1, dtype=numpy.int32)*2 # Pos. y bc
# Boundary conditions at the top and bottom of the fluid grid
# 0: Dirichlet (default)
t@@ -369,96 +353,91 @@ class sim:
# 2: Neumann no slip (Navier Stokes), Periodic (Darcy)
# 3: Periodic (Navier-Stokes solver only)
# 4: Constant flux (Darcy solver only)
- self.bc_bot=numpy.zeros(1, dtype=numpy.int32)
- self.bc_top=numpy.zeros(1, dtype=numpy.int32)
+ self.bc_bot = numpy.zeros(1, dtype=numpy.int32)
+ self.bc_top = numpy.zeros(1, dtype=numpy.int32)
# Free slip boundaries? 1: yes
- self.free_slip_bot=numpy.ones(1, dtype=numpy.int32)
- self.free_slip_top=numpy.ones(1, dtype=numpy.int32)
+ self.free_slip_bot = numpy.ones(1, dtype=numpy.int32)
+ self.free_slip_top = numpy.ones(1, dtype=numpy.int32)
# Boundary-normal flux (in case of bc_*=4)
- self.bc_bot_flux=numpy.zeros(1, dtype=numpy.float64)
- self.bc_top_flux=numpy.zeros(1, dtype=numpy.float64)
+ self.bc_bot_flux = numpy.zeros(1, dtype=numpy.float64)
+ self.bc_top_flux = numpy.zeros(1, dtype=numpy.float64)
# 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.p_f_constant = numpy.zeros((self.num[0],
+ self.num[1],
+ self.num[2]), dtype=numpy.int32)
# Navier-Stokes
if self.cfd_solver[0] == 0:
# Smoothing parameter, should be in the range [0.0;1.0[.
# 0.0=no smoothing.
- self.gamma=numpy.array(0.0)
+ self.gamma = numpy.array(0.0)
# Under-relaxation parameter, should be in the range ]0.0;1.0].
# 1.0=no under-relaxation
- self.theta=numpy.array(1.0)
+ self.theta = numpy.array(1.0)
# Velocity projection parameter, should be in the range
# [0.0;1.0]
- self.beta=numpy.array(0.0)
+ self.beta = numpy.array(0.0)
# Tolerance criteria for the normalized max. residual
- self.tolerance=numpy.array(1.0e-3)
+ self.tolerance = numpy.array(1.0e-3)
# The maximum number of iterations to perform per time step
- self.maxiter=numpy.array(1e4)
+ self.maxiter = numpy.array(1e4)
# The number of DEM time steps to perform between CFD updates
- self.ndem=numpy.array(1)
+ self.ndem = numpy.array(1)
# Porosity scaling factor
- self.c_phi=numpy.ones(1, dtype=numpy.float64)
+ self.c_phi = numpy.ones(1, dtype=numpy.float64)
# Fluid velocity scaling factor
- self.c_v=numpy.ones(1, dtype=numpy.float64)
+ self.c_v = numpy.ones(1, dtype=numpy.float64)
# DEM-CFD time scaling factor
- self.dt_dem_fac=numpy.ones(1, dtype=numpy.float64)
+ self.dt_dem_fac = numpy.ones(1, dtype=numpy.float64)
## Interaction forces
- 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)
+ 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)
# Darcy
elif self.cfd_solver[0] == 1:
# Tolerance criteria for the normalized max. residual
- self.tolerance=numpy.array(1.0e-3)
+ self.tolerance = numpy.array(1.0e-3)
# The maximum number of iterations to perform per time step
- self.maxiter=numpy.array(1e4)
+ self.maxiter = numpy.array(1e4)
# The number of DEM time steps to perform between CFD updates
- self.ndem=numpy.array(1)
+ self.ndem = numpy.array(1)
# Porosity scaling factor
- self.c_phi=numpy.ones(1, dtype=numpy.float64)
+ self.c_phi = numpy.ones(1, dtype=numpy.float64)
# Interaction forces
- self.f_p=numpy.zeros((self.np, self.nd),
- dtype=numpy.float64)
+ self.f_p = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
# Adiabatic fluid compressibility [1/Pa].
# Fluid bulk modulus=1/self.beta_f
- self.beta_f=numpy.ones(1, dtype=numpy.float64)*4.5e-10
+ self.beta_f = numpy.ones(1, dtype=numpy.float64)*4.5e-10
# Hydraulic permeability prefactor [m*m]
- self.k_c=numpy.ones(1, dtype=numpy.float64)*4.6e-10
+ self.k_c = numpy.ones(1, dtype=numpy.float64)*4.6e-10
else:
raise Exception('Value of cfd_solver not understood (' + \
- str(self.cfd_solver[0]) + ')')
+ str(self.cfd_solver[0]) + ')')
# Particle color marker
- self.color=numpy.zeros(self.np, dtype=numpy.int32)
+ self.color = numpy.zeros(self.np, dtype=numpy.int32)
def __eq__(self, other):
'''
t@@ -738,24 +717,24 @@ class sim:
return False
elif self.c_phi != other.c_phi:
print('c_phi')
- return(84)
+ return 84
elif self.c_v != other.c_v:
print('c_v')
elif self.dt_dem_fac != other.dt_dem_fac:
print('dt_dem_fac')
- return(85)
+ return 85
elif (self.f_d != other.f_d).any():
print('f_d')
- return(86)
+ return 86
elif (self.f_p != other.f_p).any():
print('f_p')
- return(87)
+ return 87
elif (self.f_v != other.f_v).any():
print('f_v')
- return(88)
+ return 88
elif (self.f_sum != other.f_sum).any():
print('f_sum')
- return(89)
+ return 89
if self.cfd_solver == 1:
if self.tolerance != other.tolerance:
t@@ -769,26 +748,26 @@ class sim:
return False
elif self.c_phi != other.c_phi:
print('c_phi')
- return(84)
+ return 84
elif (self.f_p != other.f_p).any():
print('f_p')
- return(86)
+ return 86
elif self.beta_f != other.beta_f:
print('beta_f')
- return(87)
+ return 87
elif self.k_c != other.k_c:
print('k_c')
- return(88)
- elif (self.bc_xn != other.bc_xn):
+ return 88
+ elif self.bc_xn != other.bc_xn:
print('bc_xn')
return False
- elif (self.bc_xp != other.bc_xp):
+ elif self.bc_xp != other.bc_xp:
print('bc_xp')
return False
- elif (self.bc_yn != other.bc_yn):
+ elif self.bc_yn != other.bc_yn:
print('bc_yn')
return False
- elif (self.bc_yp != other.bc_yp):
+ elif self.bc_yp != other.bc_yp:
print('bc_yp')
return False
t@@ -813,7 +792,7 @@ class sim:
if sid == '':
return self.sid
else:
- self.sid=sid
+ self.sid = sid
def idAppend(self, string):
'''
t@@ -825,22 +804,12 @@ class sim:
'''
self.sid += string
- def addParticle(self,
- x,
- radius,
- xyzsum=numpy.zeros(3),
- vel=numpy.zeros(3),
- fixvel=numpy.zeros(1),
- force=numpy.zeros(3),
- angpos=numpy.zeros(3),
- angvel=numpy.zeros(3),
- torque=numpy.zeros(3),
- es_dot=numpy.zeros(1),
- es=numpy.zeros(1),
- ev_dot=numpy.zeros(1),
- ev=numpy.zeros(1),
- p=numpy.zeros(1),
- color=0):
+ def addParticle(self, x, radius, xyzsum=numpy.zeros(3), vel=numpy.zeros(3),
+ fixvel=numpy.zeros(1), force=numpy.zeros(3),
+ angpos=numpy.zeros(3), angvel=numpy.zeros(3),
+ torque=numpy.zeros(3), es_dot=numpy.zeros(1),
+ es=numpy.zeros(1), ev_dot=numpy.zeros(1),
+ ev=numpy.zeros(1), p=numpy.zeros(1), color=0):
'''
Add a single particle to the simulation object. The only required
parameters are the position (x) and the radius (radius).
t@@ -849,17 +818,17 @@ class sim:
:type x: numpy.array
:param radius: The particle radius
:type radius: float
- :param vel: The particle linear velocity (default=[0,0,0])
+ :param vel: The particle linear velocity (default=[0, 0, 0])
:type vel: numpy.array
:param fixvel: 0: Do not fix particle velocity (default), 1: Fix
horizontal linear velocity, -1: Fix horizontal and vertical linear
velocity
:type fixvel: float
- :param angpos: The particle angular position (default=[0,0,0])
+ :param angpos: The particle angular position (default=[0, 0, 0])
:type angpos: numpy.array
- :param angvel: The particle angular velocity (default=[0,0,0])
+ :param angvel: The particle angular velocity (default=[0, 0, 0])
:type angvel: numpy.array
- :param torque: The particle torque (default=[0,0,0])
+ :param torque: The particle torque (default=[0, 0, 0])
:type torque: numpy.array
:param es_dot: The particle shear energy loss rate (default=0)
:type es_dot: float
t@@ -873,28 +842,28 @@ class sim:
:type p: float
'''
- self.np=self.np + 1
-
- self.x =numpy.append(self.x, [x], axis=0)
- self.radius=numpy.append(self.radius, radius)
- self.vel =numpy.append(self.vel, [vel], axis=0)
- self.xyzsum=numpy.append(self.xyzsum, [xyzsum], axis=0)
- self.fixvel=numpy.append(self.fixvel, fixvel)
- self.force =numpy.append(self.force, [force], axis=0)
- self.angpos=numpy.append(self.angpos, [angpos], axis=0)
- self.angvel=numpy.append(self.angvel, [angvel], axis=0)
- self.torque=numpy.append(self.torque, [torque], axis=0)
- self.es_dot=numpy.append(self.es_dot, es_dot)
- self.es =numpy.append(self.es, es)
- self.ev_dot=numpy.append(self.ev_dot, ev_dot)
- self.ev =numpy.append(self.ev, ev)
- self.p =numpy.append(self.p, p)
- self.color =numpy.append(self.color, color)
+ self.np += 1
+
+ self.x = numpy.append(self.x, [x], axis=0)
+ self.radius = numpy.append(self.radius, radius)
+ self.vel = numpy.append(self.vel, [vel], axis=0)
+ self.xyzsum = numpy.append(self.xyzsum, [xyzsum], axis=0)
+ self.fixvel = numpy.append(self.fixvel, fixvel)
+ self.force = numpy.append(self.force, [force], axis=0)
+ self.angpos = numpy.append(self.angpos, [angpos], axis=0)
+ self.angvel = numpy.append(self.angvel, [angvel], axis=0)
+ self.torque = numpy.append(self.torque, [torque], axis=0)
+ self.es_dot = numpy.append(self.es_dot, es_dot)
+ self.es = numpy.append(self.es, es)
+ self.ev_dot = numpy.append(self.ev_dot, ev_dot)
+ self.ev = numpy.append(self.ev, ev)
+ self.p = numpy.append(self.p, p)
+ self.color = numpy.append(self.color, color)
if self.fluid:
- self.f_d =numpy.append(self.f_d, [numpy.zeros(3)], axis=0)
- self.f_p =numpy.append(self.f_p, [numpy.zeros(3)], axis=0)
- self.f_v =numpy.append(self.f_v, [numpy.zeros(3)], axis=0)
- self.f_sum =numpy.append(self.f_sum, [numpy.zeros(3)], axis=0)
+ self.f_d = numpy.append(self.f_d, [numpy.zeros(3)], axis=0)
+ self.f_p = numpy.append(self.f_p, [numpy.zeros(3)], axis=0)
+ self.f_v = numpy.append(self.f_v, [numpy.zeros(3)], axis=0)
+ self.f_sum = numpy.append(self.f_sum, [numpy.zeros(3)], axis=0)
def deleteParticle(self, i):
'''
t@@ -918,59 +887,59 @@ class sim:
if type(i) == tuple:
raise Exception('Cannot parse tuples as index value. ' +
- 'Valid types are int, list and numpy.ndarray')
-
-
- self.x =numpy.delete(self.x, i, axis=0)
- self.radius=numpy.delete(self.radius, i)
- self.vel =numpy.delete(self.vel, i, axis=0)
- self.xyzsum=numpy.delete(self.xyzsum, i, axis=0)
- self.fixvel=numpy.delete(self.fixvel, i)
- self.force =numpy.delete(self.force, i, axis=0)
- self.angpos=numpy.delete(self.angpos, i, axis=0)
- self.angvel=numpy.delete(self.angvel, i, axis=0)
- self.torque=numpy.delete(self.torque, i, axis=0)
- self.es_dot=numpy.delete(self.es_dot, i)
- self.es =numpy.delete(self.es, i)
- self.ev_dot=numpy.delete(self.ev_dot, i)
- self.ev =numpy.delete(self.ev, i)
- self.p =numpy.delete(self.p, i)
- self.color =numpy.delete(self.color, i)
+ 'Valid types are int, list and numpy.ndarray')
+
+
+ self.x = numpy.delete(self.x, i, axis=0)
+ self.radius = numpy.delete(self.radius, i)
+ self.vel = numpy.delete(self.vel, i, axis=0)
+ self.xyzsum = numpy.delete(self.xyzsum, i, axis=0)
+ self.fixvel = numpy.delete(self.fixvel, i)
+ self.force = numpy.delete(self.force, i, axis=0)
+ self.angpos = numpy.delete(self.angpos, i, axis=0)
+ self.angvel = numpy.delete(self.angvel, i, axis=0)
+ self.torque = numpy.delete(self.torque, i, axis=0)
+ self.es_dot = numpy.delete(self.es_dot, i)
+ self.es = numpy.delete(self.es, i)
+ self.ev_dot = numpy.delete(self.ev_dot, i)
+ self.ev = numpy.delete(self.ev, i)
+ self.p = numpy.delete(self.p, i)
+ self.color = numpy.delete(self.color, i)
if self.fluid:
# Darcy and Navier-Stokes
- self.f_p =numpy.delete(self.f_p, i, axis=0)
+ self.f_p = numpy.delete(self.f_p, i, axis=0)
if self.cfd_solver[0] == 0: # Navier-Stokes
- self.f_d =numpy.delete(self.f_d, i, axis=0)
- self.f_v =numpy.delete(self.f_v, i, axis=0)
- self.f_sum =numpy.delete(self.f_sum, i, axis=0)
+ self.f_d = numpy.delete(self.f_d, i, axis=0)
+ self.f_v = numpy.delete(self.f_v, i, axis=0)
+ self.f_sum = numpy.delete(self.f_sum, i, axis=0)
def deleteAllParticles(self):
'''
Deletes all particles in the simulation object.
'''
- 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)
+ 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):
+ esysparticle=False):
'''
Reads a target ``sphere`` binary file.
t@@ -994,368 +963,327 @@ class sim:
:type esysparticle: bool
'''
- fh=None
+ fh = None
try:
if verbose:
print("Input file: {0}".format(targetbin))
- fh=open(targetbin, "rb")
+ fh = open(targetbin, "rb")
# Read the file version
- self.version=numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ self.version = numpy.fromfile(fh, dtype=numpy.float64, count=1)
# Read the number of dimensions and particles
- self.nd=int(numpy.fromfile(fh, dtype=numpy.int32, count=1))
- self.np=int(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=\
- numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.time_current =\
- numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.time_total =\
- numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.time_file_dt =\
- numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.time_step_count =\
- numpy.fromfile(fh, dtype=numpy.uint32, count=1)
+ self.time_dt = numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ self.time_current = numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ self.time_total = numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ self.time_file_dt = numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ self.time_step_count = numpy.fromfile(fh, dtype=numpy.uint32, count=1)
# Allocate array memory for particles
- 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)
+ 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)
- 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)
+ 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:
- self.adaptive=numpy.fromfile(fh, dtype=numpy.int32, count=1)
+ self.adaptive = numpy.fromfile(fh, dtype=numpy.int32, count=1)
else:
- self.adaptive=numpy.zeros(1, dtype=numpy.float64)
+ self.adaptive = numpy.zeros(1, dtype=numpy.float64)
# Per-particle vectors
for i in numpy.arange(self.np):
- self.x[i,:] =\
+ self.x[i, :] =\
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,3)
+ self.xyzsum = numpy.fromfile(fh, dtype=numpy.float64,\
+ count=self.np*3).reshape(self.np, 3)
else:
- self.xyzsum=numpy.fromfile(fh, dtype=numpy.float64,\
- count=self.np*2).reshape(self.np,2)
+ self.xyzsum = numpy.fromfile(fh, dtype=numpy.float64,\
+ 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)
- self.fixvel[i] =\
- numpy.fromfile(fh, dtype=numpy.float64, count=1)
-
- self.force=numpy.fromfile(fh, dtype=numpy.float64,\
- count=self.np*self.nd).reshape(self.np, self.nd)
-
- self.angpos=numpy.fromfile(fh, dtype=numpy.float64,\
- count=self.np*self.nd).reshape(self.np, self.nd)
- self.angvel=numpy.fromfile(fh, dtype=numpy.float64,\
- count=self.np*self.nd).reshape(self.np, self.nd)
- self.torque=numpy.fromfile(fh, dtype=numpy.float64,\
- count=self.np*self.nd).reshape(self.np, self.nd)
+ self.vel[i, :] = 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)\
+ .reshape(self.np, self.nd)
+
+ self.angpos = numpy.fromfile(fh, dtype=numpy.float64,\
+ count=self.np*self.nd)\
+ .reshape(self.np, self.nd)
+ self.angvel = numpy.fromfile(fh, dtype=numpy.float64,\
+ count=self.np*self.nd)\
+ .reshape(self.np, self.nd)
+ self.torque = numpy.fromfile(fh, dtype=numpy.float64,\
+ count=self.np*self.nd)\
+ .reshape(self.np, self.nd)
if esysparticle:
return
# Per-particle single-value parameters
- self.es_dot=numpy.fromfile(fh, dtype=numpy.float64, count=self.np)
- self.es =numpy.fromfile(fh, dtype=numpy.float64, count=self.np)
- self.ev_dot=numpy.fromfile(fh, dtype=numpy.float64, count=self.np)
- self.ev =numpy.fromfile(fh, dtype=numpy.float64, count=self.np)
- self.p =numpy.fromfile(fh, dtype=numpy.float64, count=self.np)
+ self.es_dot = numpy.fromfile(fh, dtype=numpy.float64, count=self.np)
+ self.es = numpy.fromfile(fh, dtype=numpy.float64, count=self.np)
+ self.ev_dot = numpy.fromfile(fh, dtype=numpy.float64, count=self.np)
+ self.ev = numpy.fromfile(fh, dtype=numpy.float64, count=self.np)
+ 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)
- 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)
+ 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)
if self.version >= 2.13:
- self.E=numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ self.E = numpy.fromfile(fh, dtype=numpy.float64, count=1)
else:
- self.E=numpy.zeros(1, dtype=numpy.float64)
- self.gamma_n =numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.gamma_t =numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.gamma_r =numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.mu_s =numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.mu_d =numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.mu_r =numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.gamma_wn =numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.gamma_wt =numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.mu_ws =numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.mu_wd =numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.rho =numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.contactmodel=numpy.fromfile(fh, dtype=numpy.uint32, count=1)
- self.kappa =numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.db =numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.V_b =numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ self.E = numpy.zeros(1, dtype=numpy.float64)
+ self.gamma_n = numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ self.gamma_t = numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ self.gamma_r = numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ self.mu_s = numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ self.mu_d = numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ self.mu_r = numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ self.gamma_wn = numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ self.gamma_wt = numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ self.mu_ws = numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ self.mu_wd = numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ self.rho = numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ self.contactmodel = numpy.fromfile(fh, dtype=numpy.uint32, count=1)
+ self.kappa = numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ self.db = numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ self.V_b = numpy.fromfile(fh, dtype=numpy.float64, count=1)
# Wall data
- 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,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)
+ 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, 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,:] =\
+ self.w_n[i, :] =\
numpy.fromfile(fh, dtype=numpy.float64, count=self.nd)
- self.w_x[i] =numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ self.w_x[i] = numpy.fromfile(fh, dtype=numpy.float64, count=1)
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] =\
- numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.w_sigma0[i]= numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ 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] = numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ self.w_sigma0[i] = numpy.fromfile(fh, dtype=numpy.float64, count=1)
if sigma0mod:
- self.w_sigma0_A=numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.w_sigma0_f=numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ self.w_sigma0_A = numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ self.w_sigma0_f = numpy.fromfile(fh, dtype=numpy.float64, count=1)
if self.version >= 2.1:
- self.w_tau_x=numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ self.w_tau_x = numpy.fromfile(fh, dtype=numpy.float64, count=1)
else:
- self.w_tau_x=numpy.zeros(1, dtype=numpy.float64)
+ self.w_tau_x = numpy.zeros(1, dtype=numpy.float64)
if bonds:
# Inter-particle bonds
- self.lambda_bar =\
- numpy.fromfile(fh, dtype=numpy.float64, 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, 2), dtype=numpy.uint32)
+ self.lambda_bar = numpy.fromfile(fh, dtype=numpy.float64, 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, 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)
- self.bonds_delta_t=numpy.fromfile(fh, dtype=numpy.float64,
- count=self.nb0*self.nd).reshape(self.nb0, self.nd)
- self.bonds_omega_n=numpy.fromfile(fh, dtype=numpy.float64,
- count=self.nb0)
- self.bonds_omega_t=numpy.fromfile(fh, dtype=numpy.float64,
- count=self.nb0*self.nd).reshape(self.nb0, self.nd)
+ 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)
+ self.bonds_delta_t = numpy.fromfile(fh, dtype=numpy.float64,
+ count=self.nb0*self.nd)\
+ .reshape(self.nb0, self.nd)
+ self.bonds_omega_n = numpy.fromfile(fh, dtype=numpy.float64,
+ count=self.nb0)
+ self.bonds_omega_t = numpy.fromfile(fh, dtype=numpy.float64,
+ count=self.nb0*self.nd)\
+ .reshape(self.nb0, self.nd)
else:
- self.nb0=0
+ self.nb0 = 0
if self.fluid:
if self.version >= 2.0:
- self.cfd_solver=numpy.fromfile(fh, dtype=numpy.int32,
- count=1)
+ self.cfd_solver = numpy.fromfile(fh, dtype=numpy.int32, count=1)
else:
- self.cfd_solver=numpy.zeros(1, dtype=numpy.int32)
-
- 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.p_f=\
- numpy.empty((self.num[0],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.dphi=\
- numpy.empty((self.num[0],self.num[1],self.num[2]),
- dtype=numpy.float64)
+ self.cfd_solver = numpy.zeros(1, dtype=numpy.int32)
+
+ 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.p_f = numpy.empty((self.num[0],
+ 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.dphi = numpy.empty((self.num[0],
+ 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)
- 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.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.time_dt*self.ndem)
+ self.v_f[x, y, z, 0] = numpy.fromfile(fh,
+ 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.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.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_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.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)
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)
+ 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.p_f_constant = numpy.empty((self.num[0],
+ 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)
+ self.p_f_constant[x, y, z] = \
+ 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.p_f_constant = numpy.zeros((self.num[0],
+ 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
+ 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)
+ self.c_a = numpy.fromfile(fh, dtype=numpy.float64,
+ 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)
+ self.c_a = numpy.ones(1, dtype=numpy.float64)
else:
- self.c_phi=numpy.ones(1, dtype=numpy.float64)
- self.c_v=numpy.ones(1, dtype=numpy.float64)
+ self.c_phi = numpy.ones(1, dtype=numpy.float64)
+ self.c_v = numpy.ones(1, dtype=numpy.float64)
if self.version >= 1.05:
- self.f_d=numpy.empty_like(self.x)
- self.f_p=numpy.empty_like(self.x)
- self.f_v=numpy.empty_like(self.x)
- self.f_sum=numpy.empty_like(self.x)
+ self.f_d = numpy.empty_like(self.x)
+ self.f_p = numpy.empty_like(self.x)
+ self.f_v = numpy.empty_like(self.x)
+ 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)
- 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)
+ 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.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.ndem=numpy.fromfile(fh, dtype=numpy.uint32, count=1)
- self.c_phi=\
- numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.f_p=numpy.empty_like(self.x)
+ 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.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)
-
- self.k_c=\
- numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ self.f_p[i, :] = numpy.fromfile(fh, dtype=numpy.float64,
+ 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)
+ self.color = numpy.zeros(self.np, dtype=numpy.int32)
finally:
- self.version[0]=VERSION
+ self.version[0] = VERSION
if fh is not None:
fh.close()
t@@ -1371,13 +1299,13 @@ class sim:
:param verbose: Show diagnostic information (default=True)
:type verbose: bool
'''
- fh=None
- try :
- targetbin=folder + "/" + self.sid + ".bin"
+ fh = None
+ try:
+ targetbin = folder + "/" + self.sid + ".bin"
if verbose:
print("Output file: {0}".format(targetbin))
- fh=open(targetbin, "wb")
+ fh = open(targetbin, "wb")
# Write the current version number
fh.write(self.version.astype(numpy.float64))
t@@ -1402,14 +1330,14 @@ class sim:
# Per-particle vectors
for i in numpy.arange(self.np):
- fh.write(self.x[i,:].astype(numpy.float64))
+ fh.write(self.x[i, :].astype(numpy.float64))
fh.write(self.radius[i].astype(numpy.float64))
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.vel[i, :].astype(numpy.float64))
fh.write(self.fixvel[i].astype(numpy.float64))
if self.np > 0:
t@@ -1451,7 +1379,7 @@ class sim:
for i in numpy.arange(self.nw):
fh.write(self.wmode[i].astype(numpy.int32))
for i in numpy.arange(self.nw):
- fh.write(self.w_n[i,:].astype(numpy.float64))
+ 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):
t@@ -1468,8 +1396,8 @@ class sim:
fh.write(self.sigma_b.astype(numpy.float64))
fh.write(self.tau_b.astype(numpy.float64))
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[i, 0].astype(numpy.uint32))
+ fh.write(self.bonds[i, 1].astype(numpy.uint32))
fh.write(self.bonds_delta_n.astype(numpy.float64))
fh.write(self.bonds_delta_t.astype(numpy.float64))
fh.write(self.bonds_omega_n.astype(numpy.float64))
t@@ -1482,13 +1410,13 @@ class sim:
for z in numpy.arange(self.num[2]):
for y in numpy.arange(self.num[1]):
for x in numpy.arange(self.num[0]):
- fh.write(self.v_f[x,y,z,0].astype(numpy.float64))
- fh.write(self.v_f[x,y,z,1].astype(numpy.float64))
- fh.write(self.v_f[x,y,z,2].astype(numpy.float64))
- fh.write(self.p_f[x,y,z].astype(numpy.float64))
- fh.write(self.phi[x,y,z].astype(numpy.float64))
- fh.write(self.dphi[x,y,z].astype(numpy.float64)*
- self.time_dt*self.ndem)
+ fh.write(self.v_f[x, y, z, 0].astype(numpy.float64))
+ fh.write(self.v_f[x, y, z, 1].astype(numpy.float64))
+ fh.write(self.v_f[x, y, z, 2].astype(numpy.float64))
+ fh.write(self.p_f[x, y, z].astype(numpy.float64))
+ fh.write(self.phi[x, y, z].astype(numpy.float64))
+ fh.write(self.dphi[x, y, z].astype(numpy.float64)*
+ self.time_dt*self.ndem)
fh.write(self.rho_f.astype(numpy.float64))
fh.write(self.p_mod_A.astype(numpy.float64))
t@@ -1511,7 +1439,7 @@ class sim:
for z in numpy.arange(self.num[2]):
for y in numpy.arange(self.num[1]):
for x in numpy.arange(self.num[0]):
- fh.write(self.p_f_constant[x,y,z].astype(
+ fh.write(self.p_f_constant[x, y, z].astype(
numpy.int32))
# Navier Stokes
t@@ -1528,13 +1456,13 @@ class sim:
fh.write(self.dt_dem_fac.astype(numpy.float64))
for i in numpy.arange(self.np):
- fh.write(self.f_d[i,:].astype(numpy.float64))
+ fh.write(self.f_d[i, :].astype(numpy.float64))
for i in numpy.arange(self.np):
- fh.write(self.f_p[i,:].astype(numpy.float64))
+ fh.write(self.f_p[i, :].astype(numpy.float64))
for i in numpy.arange(self.np):
- fh.write(self.f_v[i,:].astype(numpy.float64))
+ fh.write(self.f_v[i, :].astype(numpy.float64))
for i in numpy.arange(self.np):
- fh.write(self.f_sum[i,:].astype(numpy.float64))
+ fh.write(self.f_sum[i, :].astype(numpy.float64))
# Darcy
elif self.cfd_solver[0] == 1:
t@@ -1544,7 +1472,7 @@ class sim:
fh.write(self.ndem.astype(numpy.uint32))
fh.write(self.c_phi.astype(numpy.float64))
for i in numpy.arange(self.np):
- fh.write(self.f_p[i,:].astype(numpy.float64))
+ fh.write(self.f_p[i, :].astype(numpy.float64))
fh.write(self.beta_f.astype(numpy.float64))
fh.write(self.k_c.astype(numpy.float64))
t@@ -1627,8 +1555,8 @@ class sim:
print('skipping ' + fn_vtk +
': file exists and is newer than ' + fn)
if self.fluid:
- fn_vtk = "../output/fluid-{0}.{1:0=5}.vti".format(
- self.sid, i)
+ fn_vtk = "../output/fluid-{0}.{1:0=5}.vti" \
+ .format(self.sid, i)
if os.path.isfile(fn_vtk) and \
(os.path.getmtime(fn) < os.path.getmtime(fn_vtk)):
if verbose:
t@@ -1660,7 +1588,7 @@ class sim:
if verbose:
print("\tto")
sb.writeFluidVTK(verbose=verbose,
- cell_centered=cell_centered)
+ cell_centered=cell_centered)
else:
sb.writeFluidVTK(verbose=False, cell_centered=cell_centered)
t@@ -1697,20 +1625,20 @@ class sim:
:type verbose: bool
'''
- fh=None
- try :
- targetbin=folder + '/' + self.sid + '.vtu' # unstructured grid
+ fh = None
+ try:
+ targetbin = folder + '/' + self.sid + '.vtu' # unstructured grid
if verbose:
print('Output file: ' + targetbin)
- fh=open(targetbin, 'w')
+ fh = open(targetbin, 'w')
# the VTK data file format is documented in
# http://www.vtk.org/VTK/img/file-formats.pdf
fh.write('<?xml version="1.0"?>\n') # XML header
fh.write('<VTKFile type="UnstructuredGrid" version="0.1" '
- + 'byte_order="LittleEndian">\n') # VTK header
+ + 'byte_order="LittleEndian">\n') # VTK header
fh.write(' <UnstructuredGrid>\n')
fh.write(' <Piece NumberOfPoints="%d" NumberOfCells="0">\n' \
% (self.np))
t@@ -1718,10 +1646,10 @@ class sim:
# Coordinates for each point (positions)
fh.write(' <Points>\n')
fh.write(' <DataArray name="Position [m]" type="Float32" '
- + 'NumberOfComponents="3" format="ascii">\n')
+ + 'NumberOfComponents="3" format="ascii">\n')
fh.write(' ')
for i in range(self.np):
- fh.write('%f %f %f ' % (self.x[i,0], self.x[i,1], self.x[i,2]))
+ fh.write('%f %f %f ' % (self.x[i, 0], self.x[i, 1], self.x[i, 2]))
fh.write('\n')
fh.write(' </DataArray>\n')
fh.write(' </Points>\n')
t@@ -1731,7 +1659,7 @@ class sim:
# Radii
fh.write(' <DataArray type="Float32" Name="Diameter" '
- + 'format="ascii">\n')
+ + 'format="ascii">\n')
fh.write(' ')
for i in range(self.np):
fh.write('%f ' % (self.radius[i]*2.0))
t@@ -1740,21 +1668,21 @@ class sim:
# Displacements (xyzsum)
fh.write(' <DataArray type="Float32" Name="Displacement [m]" '
- + 'NumberOfComponents="3" format="ascii">\n')
+ + 'NumberOfComponents="3" format="ascii">\n')
fh.write(' ')
for i in range(self.np):
fh.write('%f %f %f ' % \
- (self.xyzsum[i,0], self.xyzsum[i,1], self.xyzsum[i,2]))
+ (self.xyzsum[i, 0], self.xyzsum[i, 1], self.xyzsum[i, 2]))
fh.write('\n')
fh.write(' </DataArray>\n')
# Velocity
fh.write(' <DataArray type="Float32" Name="Velocity [m/s]" '
- + 'NumberOfComponents="3" format="ascii">\n')
+ + 'NumberOfComponents="3" format="ascii">\n')
fh.write(' ')
for i in range(self.np):
fh.write('%f %f %f ' % \
- (self.vel[i,0], self.vel[i,1], self.vel[i,2]))
+ (self.vel[i, 0], self.vel[i, 1], self.vel[i, 2]))
fh.write('\n')
fh.write(' </DataArray>\n')
t@@ -1763,53 +1691,57 @@ class sim:
if self.cfd_solver == 0: # Navier Stokes
# Fluid interaction force
fh.write(' <DataArray type="Float32" '
- + 'Name="Fluid force total [N]" '
- + 'NumberOfComponents="3" format="ascii">\n')
+ + 'Name="Fluid force total [N]" '
+ + 'NumberOfComponents="3" format="ascii">\n')
fh.write(' ')
for i in range(self.np):
fh.write('%f %f %f ' % \
- (self.f_sum[i,0], self.f_sum[i,1], \
- self.f_sum[i,2]))
+ (self.f_sum[i, 0], self.f_sum[i, 1], \
+ self.f_sum[i, 2]))
fh.write('\n')
fh.write(' </DataArray>\n')
# Fluid drag force
fh.write(' <DataArray type="Float32" '
- + 'Name="Fluid drag force [N]" '
- + 'NumberOfComponents="3" format="ascii">\n')
+ + 'Name="Fluid drag force [N]" '
+ + 'NumberOfComponents="3" format="ascii">\n')
fh.write(' ')
for i in range(self.np):
fh.write('%f %f %f ' % \
- (self.f_d[i,0], self.f_d[i,1], self.f_d[i,2]))
+ (self.f_d[i, 0],
+ self.f_d[i, 1],
+ self.f_d[i, 2]))
fh.write('\n')
fh.write(' </DataArray>\n')
# Fluid pressure force
fh.write(' <DataArray type="Float32" '
- + 'Name="Fluid pressure force [N]" '
- + 'NumberOfComponents="3" format="ascii">\n')
+ + 'Name="Fluid pressure force [N]" '
+ + 'NumberOfComponents="3" format="ascii">\n')
fh.write(' ')
for i in range(self.np):
fh.write('%f %f %f ' % \
- (self.f_p[i,0], self.f_p[i,1], self.f_p[i,2]))
+ (self.f_p[i, 0], self.f_p[i, 1], self.f_p[i, 2]))
fh.write('\n')
fh.write(' </DataArray>\n')
if self.cfd_solver == 0: # Navier Stokes
# Fluid viscous force
fh.write(' <DataArray type="Float32" '
- + 'Name="Fluid viscous force [N]" '
- + 'NumberOfComponents="3" format="ascii">\n')
+ + 'Name="Fluid viscous force [N]" '
+ + 'NumberOfComponents="3" format="ascii">\n')
fh.write(' ')
for i in range(self.np):
fh.write('%f %f %f ' % \
- (self.f_v[i,0], self.f_v[i,1], self.f_v[i,2]))
+ (self.f_v[i, 0],
+ self.f_v[i, 1],
+ self.f_v[i, 2]))
fh.write('\n')
fh.write(' </DataArray>\n')
# fixvel
fh.write(' <DataArray type="Float32" Name="FixedVel" '
- + 'format="ascii">\n')
+ + 'format="ascii">\n')
fh.write(' ')
for i in range(self.np):
fh.write('%f ' % (self.fixvel[i]))
t@@ -1818,50 +1750,54 @@ class sim:
# Force
fh.write(' <DataArray type="Float32" Name="Force [N]" '
- + 'NumberOfComponents="3" format="ascii">\n')
+ + 'NumberOfComponents="3" format="ascii">\n')
fh.write(' ')
for i in range(self.np):
- fh.write('%f %f %f ' % \
- (self.force[i,0], self.force[i,1], self.force[i,2]))
+ fh.write('%f %f %f ' % (self.force[i, 0],
+ self.force[i, 1],
+ self.force[i, 2]))
fh.write('\n')
fh.write(' </DataArray>\n')
# Angular Position
fh.write(' <DataArray type="Float32" Name="Angular position'
- + '[rad]" '
- + 'NumberOfComponents="3" format="ascii">\n')
+ + '[rad]" '
+ + 'NumberOfComponents="3" format="ascii">\n')
fh.write(' ')
for i in range(self.np):
- fh.write('%f %f %f ' % \
- (self.angpos[i,0], self.angpos[i,1], self.angpos[i,2]))
+ fh.write('%f %f %f ' % (self.angpos[i, 0],
+ self.angpos[i, 1],
+ self.angpos[i, 2]))
fh.write('\n')
fh.write(' </DataArray>\n')
# Angular Velocity
fh.write(' <DataArray type="Float32" Name="Angular velocity'
- + ' [rad/s]" '
- + 'NumberOfComponents="3" format="ascii">\n')
+ + ' [rad/s]" '
+ + 'NumberOfComponents="3" format="ascii">\n')
fh.write(' ')
for i in range(self.np):
- fh.write('%f %f %f ' % \
- (self.angvel[i,0], self.angvel[i,1], self.angvel[i,2]))
+ fh.write('%f %f %f ' % (self.angvel[i, 0],
+ self.angvel[i, 1],
+ self.angvel[i, 2]))
fh.write('\n')
fh.write(' </DataArray>\n')
# Torque
fh.write(' <DataArray type="Float32" Name="Torque [Nm]" '
- + 'NumberOfComponents="3" format="ascii">\n')
+ + 'NumberOfComponents="3" format="ascii">\n')
fh.write(' ')
for i in range(self.np):
- fh.write('%f %f %f ' % \
- (self.torque[i,0], self.torque[i,1], self.torque[i,2]))
+ fh.write('%f %f %f ' % (self.torque[i, 0],
+ self.torque[i, 1],
+ self.torque[i, 2]))
fh.write('\n')
fh.write(' </DataArray>\n')
# Shear energy rate
fh.write(' <DataArray type="Float32" Name="Shear Energy '
- + 'Rate [J/s]" '
- + 'format="ascii">\n')
+ + 'Rate [J/s]" '
+ + 'format="ascii">\n')
fh.write(' ')
for i in range(self.np):
fh.write('%f ' % (self.es_dot[i]))
t@@ -1870,7 +1806,7 @@ class sim:
# Shear energy
fh.write(' <DataArray type="Float32" Name="Shear Energy [J]"'
- + ' format="ascii">\n')
+ + ' format="ascii">\n')
fh.write(' ')
for i in range(self.np):
fh.write('%f ' % (self.es[i]))
t@@ -1879,7 +1815,7 @@ class sim:
# Viscous energy rate
fh.write(' <DataArray type="Float32" '
- + 'Name="Viscous Energy Rate [J/s]" format="ascii">\n')
+ + 'Name="Viscous Energy Rate [J/s]" format="ascii">\n')
fh.write(' ')
for i in range(self.np):
fh.write('%f ' % (self.ev_dot[i]))
t@@ -1888,8 +1824,8 @@ class sim:
# Shear energy
fh.write(' <DataArray type="Float32" '
- + 'Name="Viscous Energy [J]" '
- + 'format="ascii">\n')
+ + 'Name="Viscous Energy [J]" '
+ + 'format="ascii">\n')
fh.write(' ')
for i in range(self.np):
fh.write('%f ' % (self.ev[i]))
t@@ -1898,7 +1834,7 @@ class sim:
# Pressure
fh.write(' <DataArray type="Float32" Name="Pressure [Pa]" '
- + 'format="ascii">\n')
+ + 'format="ascii">\n')
fh.write(' ')
for i in range(self.np):
fh.write('%f ' % (self.p[i]))
t@@ -1907,7 +1843,7 @@ class sim:
# Color
fh.write(' <DataArray type="Int32" Name="Type color" '
- + 'format="ascii">\n')
+ + 'format="ascii">\n')
fh.write(' ')
for i in range(self.np):
fh.write('%d ' % (self.color[i]))
t@@ -1918,13 +1854,13 @@ class sim:
fh.write(' </PointData>\n')
fh.write(' <Cells>\n')
fh.write(' <DataArray type="Int32" Name="connectivity" '
- + 'format="ascii">\n')
+ + 'format="ascii">\n')
fh.write(' </DataArray>\n')
fh.write(' <DataArray type="Int32" Name="offsets" '
- + 'format="ascii">\n')
+ + 'format="ascii">\n')
fh.write(' </DataArray>\n')
fh.write(' <DataArray type="UInt8" Name="types" '
- + 'format="ascii">\n')
+ + 'format="ascii">\n')
fh.write(' </DataArray>\n')
fh.write(' </Cells>\n')
fh.write(' </Piece>\n')
t@@ -1955,45 +1891,45 @@ class sim:
print('Error: vtk module not found, cannot writeVTKforces.')
return
- filename=folder + '/forces-' + self.sid + '.vtp' # Polygon data
+ filename = folder + '/forces-' + self.sid + '.vtp' # Polygon data
# points mark the particle centers
- points=vtk.vtkPoints()
+ points = vtk.vtkPoints()
# lines mark the particle connectivity
- lines=vtk.vtkCellArray()
+ lines = vtk.vtkCellArray()
# colors
- #colors=vtk.vtkUnsignedCharArray()
+ #colors = vtk.vtkUnsignedCharArray()
#colors.SetNumberOfComponents(3)
#colors.SetName('Colors')
#colors.SetNumberOfTuples(self.overlaps.size)
# scalars
- forces=vtk.vtkDoubleArray()
+ forces = vtk.vtkDoubleArray()
forces.SetName("Force [N]")
forces.SetNumberOfComponents(1)
#forces.SetNumberOfTuples(self.overlaps.size)
forces.SetNumberOfValues(self.overlaps.size)
- stresses=vtk.vtkDoubleArray()
+ stresses = vtk.vtkDoubleArray()
stresses.SetName("Stress [Pa]")
stresses.SetNumberOfComponents(1)
stresses.SetNumberOfValues(self.overlaps.size)
for i in numpy.arange(self.overlaps.size):
- points.InsertNextPoint(self.x[self.pairs[0,i],:])
- points.InsertNextPoint(self.x[self.pairs[1,i],:])
- line=vtk.vtkLine()
+ points.InsertNextPoint(self.x[self.pairs[0, i], :])
+ points.InsertNextPoint(self.x[self.pairs[1, i], :])
+ line = vtk.vtkLine()
line.GetPointIds().SetId(0, 2*i) # index of particle 1
line.GetPointIds().SetId(1, 2*i + 1) # index of particle 2
lines.InsertNextCell(line)
- #colors.SetTupleValue(i, [100,100,100])
+ #colors.SetTupleValue(i, [100, 100, 100])
forces.SetValue(i, self.f_n_magn[i])
stresses.SetValue(i, self.sigma_contacts[i])
# initalize VTK data structure
- polydata=vtk.vtkPolyData()
+ polydata = vtk.vtkPolyData()
polydata.SetPoints(points)
polydata.SetLines(lines)
t@@ -2006,7 +1942,7 @@ class sim:
#polydata.GetPointData().SetScalars(stresses) # default scalar
# write VTK XML image data file
- writer=vtk.vtkXMLPolyDataWriter()
+ writer = vtk.vtkXMLPolyDataWriter()
writer.SetFileName(filename)
if vtk.VTK_MAJOR_VERSION <= 5:
writer.SetInput(polydata)
t@@ -2019,7 +1955,7 @@ class sim:
def writeFluidVTK(self, folder='../output/', cell_centered=True,
- verbose=True):
+ verbose=True):
'''
Writes a VTK file for the fluid grid to the ``../output/`` folder by
default. The file name will be in the format ``fluid-<self.sid>.vti``.
t@@ -2075,11 +2011,11 @@ class sim:
print('Error: vtk module not found, cannot writeFluidVTK.')
return
- filename=folder + '/fluid-' + self.sid + '.vti' # image grid
+ filename = folder + '/fluid-' + self.sid + '.vti' # image grid
# initalize VTK data structure
- grid=vtk.vtkImageData()
- dx=(self.L-self.origo)/self.num # cell center spacing
+ grid = vtk.vtkImageData()
+ dx = (self.L-self.origo)/self.num # cell center spacing
if cell_centered:
grid.SetOrigin(self.origo)
else:
t@@ -2091,7 +2027,7 @@ class sim:
grid.SetDimensions(self.num) # no. of points in each direction
# array of scalars: hydraulic pressures
- pres=vtk.vtkDoubleArray()
+ pres = vtk.vtkDoubleArray()
pres.SetName("Pressure [Pa]")
pres.SetNumberOfComponents(1)
if cell_centered:
t@@ -2100,7 +2036,7 @@ class sim:
pres.SetNumberOfTuples(grid.GetNumberOfPoints())
# array of vectors: hydraulic velocities
- vel=vtk.vtkDoubleArray()
+ vel = vtk.vtkDoubleArray()
vel.SetName("Velocity [m/s]")
vel.SetNumberOfComponents(3)
if cell_centered:
t@@ -2109,7 +2045,7 @@ class sim:
vel.SetNumberOfTuples(grid.GetNumberOfPoints())
# array of scalars: porosities
- poros=vtk.vtkDoubleArray()
+ poros = vtk.vtkDoubleArray()
poros.SetName("Porosity [-]")
poros.SetNumberOfComponents(1)
if cell_centered:
t@@ -2118,7 +2054,7 @@ class sim:
poros.SetNumberOfTuples(grid.GetNumberOfPoints())
# array of scalars: porosity change
- dporos=vtk.vtkDoubleArray()
+ dporos = vtk.vtkDoubleArray()
dporos.SetName("Porosity change [1/s]")
dporos.SetNumberOfComponents(1)
if cell_centered:
t@@ -2128,7 +2064,7 @@ class sim:
# array of scalars: Reynold's number
self.ReynoldsNumber()
- Re=vtk.vtkDoubleArray()
+ Re = vtk.vtkDoubleArray()
Re.SetName("Reynolds number [-]")
Re.SetNumberOfComponents(1)
if cell_centered:
t@@ -2139,7 +2075,7 @@ class sim:
# Find permeabilities if the Darcy solver is used
if self.cfd_solver[0] == 1:
self.findPermeabilities()
- k=vtk.vtkDoubleArray()
+ k = vtk.vtkDoubleArray()
k.SetName("Permeability [m*m]")
k.SetNumberOfComponents(1)
if cell_centered:
t@@ -2148,7 +2084,7 @@ class sim:
k.SetNumberOfTuples(grid.GetNumberOfPoints())
self.findHydraulicConductivities()
- K=vtk.vtkDoubleArray()
+ K = vtk.vtkDoubleArray()
K.SetName("Conductivity [m/s]")
K.SetNumberOfComponents(1)
if cell_centered:
t@@ -2156,7 +2092,7 @@ class sim:
else:
K.SetNumberOfTuples(grid.GetNumberOfPoints())
- p_f_constant=vtk.vtkDoubleArray()
+ p_f_constant = vtk.vtkDoubleArray()
p_f_constant.SetName("Constant pressure [-]")
p_f_constant.SetNumberOfComponents(1)
if cell_centered:
t@@ -2168,16 +2104,16 @@ class sim:
for z in range(self.num[2]):
for y in range(self.num[1]):
for x in range(self.num[0]):
- idx=x + self.num[0]*y + self.num[0]*self.num[1]*z;
- pres.SetValue(idx, self.p_f[x,y,z])
- vel.SetTuple(idx, self.v_f[x,y,z,:])
- poros.SetValue(idx, self.phi[x,y,z])
- dporos.SetValue(idx, self.dphi[x,y,z])
- Re.SetValue(idx, self.Re[x,y,z])
+ idx = x + self.num[0]*y + self.num[0]*self.num[1]*z
+ pres.SetValue(idx, self.p_f[x, y, z])
+ vel.SetTuple(idx, self.v_f[x, y, z, :])
+ poros.SetValue(idx, self.phi[x, y, z])
+ dporos.SetValue(idx, self.dphi[x, y, z])
+ Re.SetValue(idx, self.Re[x, y, z])
if self.cfd_solver[0] == 1:
- k.SetValue(idx, self.k[x,y,z])
- K.SetValue(idx, self.K[x,y,z])
- p_f_constant.SetValue(idx, self.p_f_constant[x,y,z])
+ k.SetValue(idx, self.k[x, y, z])
+ K.SetValue(idx, self.K[x, y, z])
+ p_f_constant.SetValue(idx, self.p_f_constant[x, y, z])
# add pres array to grid
if cell_centered:
t@@ -2202,7 +2138,7 @@ class sim:
grid.GetPointData().AddArray(p_f_constant)
# write VTK XML image data file
- writer=vtk.vtkXMLImageDataWriter()
+ writer = vtk.vtkXMLImageDataWriter()
writer.SetFileName(filename)
#writer.SetInput(grid) # deprecated from VTK 6
writer.SetInputData(grid)
t@@ -2227,25 +2163,25 @@ class sim:
return
# create a rendering window and renderer
- ren=vtk.vtkRenderer()
- renWin=vtk.vtkRenderWindow()
+ ren = vtk.vtkRenderer()
+ renWin = vtk.vtkRenderWindow()
renWin.AddRenderer(ren)
# create a renderwindowinteractor
- iren=vtk.vtkRenderWindowInteractor()
+ iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(renWin)
if coloring.any():
- #min_value=numpy.min(coloring)
- max_value=numpy.max(coloring)
- #min_rgb=numpy.array([50, 50, 50])
- #max_rgb=numpy.array([255, 255, 255])
+ #min_value = numpy.min(coloring)
+ max_value = numpy.max(coloring)
+ #min_rgb = numpy.array([50, 50, 50])
+ #max_rgb = numpy.array([255, 255, 255])
#def color(value):
#return (max_rgb - min_rgb) * (value - min_value)
def red(ratio):
- return numpy.fmin(1.0, 0.209*ratio**3. - 2.49*ratio**2.
- + 3.0*ratio + 0.0109)
+ return numpy.fmin(1.0, 0.209*ratio**3. - 2.49*ratio**2. + 3.0*ratio
+ + 0.0109)
def green(ratio):
return numpy.fmin(1.0, -2.44*ratio**2. + 2.15*ratio + 0.369)
def blue(ratio):
t@@ -2254,28 +2190,28 @@ class sim:
for i in numpy.arange(self.np):
# create source
- source=vtk.vtkSphereSource()
- source.SetCenter(self.x[i,:])
+ source = vtk.vtkSphereSource()
+ source.SetCenter(self.x[i, :])
source.SetRadius(self.radius[i])
source.SetThetaResolution(resolution)
source.SetPhiResolution(resolution)
# mapper
- mapper=vtk.vtkPolyDataMapper()
+ mapper = vtk.vtkPolyDataMapper()
if vtk.VTK_MAJOR_VERSION <= 5:
mapper.SetInput(source.GetOutput())
else:
mapper.SetInputConnection(source.GetOutputPort())
# actor
- actor=vtk.vtkActor()
+ actor = vtk.vtkActor()
actor.SetMapper(mapper)
# color
if coloring.any():
- ratio=coloring[i]/max_value
- r,g,b=red(ratio), green(ratio), blue(ratio)
- actor.GetProperty().SetColor(r,g,b)
+ ratio = coloring[i]/max_value
+ r, g, b = red(ratio), green(ratio), blue(ratio)
+ actor.GetProperty().SetColor(r, g, b)
# assign actor to the renderer
ren.AddActor(actor)
t@@ -2299,7 +2235,7 @@ class sim:
:func:`readstep`.
'''
- fn='../output/' + self.sid + '.output00000.bin'
+ fn = '../output/' + self.sid + '.output00000.bin'
self.readbin(fn, verbose)
def readsecond(self, verbose=True):
t@@ -2313,7 +2249,7 @@ class sim:
See also :func:`readbin()`, :func:`readfirst()`, :func:`readlast()`,
and :func:`readstep`.
'''
- fn='../output/' + self.sid + '.output00001.bin'
+ fn = '../output/' + self.sid + '.output00001.bin'
self.readbin(fn, verbose)
def readstep(self, step, verbose=True):
t@@ -2329,7 +2265,7 @@ class sim:
See also :func:`readbin()`, :func:`readfirst()`, :func:`readlast()`,
and :func:`readsecond`.
'''
- fn="../output/{0}.output{1:0=5}.bin".format(self.sid, step)
+ fn = "../output/{0}.output{1:0=5}.bin".format(self.sid, step)
self.readbin(fn, verbose)
def readlast(self, verbose=True):
t@@ -2343,8 +2279,8 @@ class sim:
See also :func:`readbin()`, :func:`readfirst()`, :func:`readsecond`, and
:func:`readstep`.
'''
- lastfile=status(self.sid)
- fn="../output/{0}.output{1:0=5}.bin".format(self.sid, lastfile)
+ lastfile = status(self.sid)
+ fn = "../output/{0}.output{1:0=5}.bin".format(self.sid, lastfile)
self.readbin(fn, verbose)
def readTime(self, time, verbose=True):
t@@ -2360,20 +2296,20 @@ class sim:
'''
self.readfirst(verbose=False)
- t_first=self.currentTime()
- n_first=self.time_step_count[0]
+ t_first = self.currentTime()
+ n_first = self.time_step_count[0]
self.readlast(verbose=False)
- t_last=self.currentTime()
- n_last=self.time_step_count[0]
+ t_last = self.currentTime()
+ n_last = self.time_step_count[0]
if time < t_first or time > t_last:
raise Exception('Error: The specified time {} s is outside the ' +
- 'range of output files [{}; {}] s.'.format(time, \
- t_first, t_last))
+ 'range of output files [{}; {}] s.'
+ .format(time, t_first, t_last))
- dt_dn=(t_last - t_first)/(n_last - n_first)
- step=int((time - t_first)/dt_dn) + n_first + 1
+ dt_dn = (t_last - t_first)/(n_last - n_first)
+ step = int((time - t_first)/dt_dn) + n_first + 1
self.readstep(step, verbose=verbose)
def generateRadii(self, psd='logn', mean=440e-6, variance=8.8e-9,
t@@ -2403,26 +2339,26 @@ class sim:
'''
if psd == 'logn': # Log-normal probability distribution
- mu=math.log(\
- (mean**2)/math.sqrt(variance+mean**2))
- sigma=math.sqrt(math.log(variance/(mean**2)+1))
- self.radius=numpy.random.lognormal(mu, sigma, self.np)
+ mu = math.log((mean**2)/math.sqrt(variance+mean**2))
+ sigma = math.sqrt(math.log(variance/(mean**2)+1))
+ self.radius = numpy.random.lognormal(mu, sigma, self.np)
elif psd == 'uni': # Uniform distribution
- radius_min=mean - variance
- radius_max=mean + variance
- self.radius=numpy.random.uniform(radius_min, radius_max, self.np)
+ radius_min = mean - variance
+ radius_max = mean + variance
+ self.radius = numpy.random.uniform(radius_min, radius_max, self.np)
else:
raise Exception('Particle size distribution type not understood ('
- + str(psd) + '). Valid values are \'uni\' or \'logn\'')
+ + str(psd) + '). '
+ + 'Valid values are \'uni\' or \'logn\'')
# Show radii as histogram
if histogram and py_mpl:
- fig=plt.figure(figsize=(8,8))
- figtitle='Particle size distribution, {0} particles'.format(\
- self.np)
- fig.text(0.5,0.95,figtitle,horizontalalignment='center',\
- fontproperties=FontProperties(size=18))
- bins=20
+ fig = plt.figure(figsize=(8, 8))
+ figtitle = 'Particle size distribution, {0} particles'\
+ .format(self.np)
+ fig.text(0.5, 0.95, figtitle, horizontalalignment='center',
+ fontproperties=FontProperties(size=18))
+ bins = 20
# Create histogram
plt.hist(self.radius, bins)
t@@ -2452,24 +2388,24 @@ class sim:
if r_small >= r_large:
raise Exception("r_large should be larger than r_small")
- V_small=V_sphere(r_small)
- V_large=V_sphere(r_large)
- nlarge=int(V_small/V_large * ratio * self.np) # ignore void volume
+ V_small = V_sphere(r_small)
+ V_large = V_sphere(r_large)
+ nlarge = int(V_small/V_large * ratio * self.np) # ignore void volume
- self.radius[:]=r_small
- self.radius[0:nlarge]=r_large
+ self.radius[:] = r_small
+ self.radius[0:nlarge] = r_large
numpy.random.shuffle(self.radius)
# Test volumetric ratio
- V_small_total=V_small * (self.np - nlarge)
- V_large_total=V_large * nlarge
+ V_small_total = V_small * (self.np - nlarge)
+ V_large_total = V_large * nlarge
if abs(V_large_total/V_small_total - ratio) > 1.0e5:
raise Exception("Volumetric ratio seems wrong")
if verbose:
print("generateBimodalRadii created " + str(nlarge)
- + " large particles, and " + str(self.np - nlarge)
- + " small")
+ + " large particles, and " + str(self.np - nlarge)
+ + " small")
def checkerboardColors(self, nx=6, ny=6, nz=6):
'''
t@@ -2482,17 +2418,17 @@ class sim:
:param nz: Number of color values along the z azis
:type nz: int
'''
- x_min=numpy.min(self.x[:,0])
- x_max=numpy.max(self.x[:,0])
- y_min=numpy.min(self.x[:,1])
- y_max=numpy.max(self.x[:,1])
- z_min=numpy.min(self.x[:,2])
- z_max=numpy.max(self.x[:,2])
+ x_min = numpy.min(self.x[:, 0])
+ x_max = numpy.max(self.x[:, 0])
+ y_min = numpy.min(self.x[:, 1])
+ y_max = numpy.max(self.x[:, 1])
+ z_min = numpy.min(self.x[:, 2])
+ z_max = numpy.max(self.x[:, 2])
for i in numpy.arange(self.np):
- ix=numpy.floor((self.x[i,0] - x_min)/(x_max/nx))
- iy=numpy.floor((self.x[i,1] - y_min)/(y_max/ny))
- iz=numpy.floor((self.x[i,2] - z_min)/(z_max/nz))
- self.color[i]=(-1)**ix + (-1)**iy + (-1)**iz
+ ix = numpy.floor((self.x[i, 0] - x_min)/(x_max/nx))
+ iy = numpy.floor((self.x[i, 1] - y_min)/(y_max/ny))
+ iz = numpy.floor((self.x[i, 2] - z_min)/(z_max/nz))
+ self.color[i] = (-1)**ix + (-1)**iy + (-1)**iz
def contactModel(self, contactmodel):
'''
t@@ -2505,7 +2441,7 @@ class sim:
(visco-frictional=1, elasto-visco-frictional=2)
:type contactmodel: int
'''
- self.contactmodel[0]=contactmodel
+ self.contactmodel[0] = contactmodel
def wall0iz(self):
'''
t@@ -2526,7 +2462,7 @@ class sim:
See also :func:`periodicBoundariesXY()` and
:func:`periodicBoundariesX()`
'''
- self.periodic[0]=0
+ self.periodic[0] = 0
def periodicBoundariesXY(self):
'''
t@@ -2535,7 +2471,7 @@ class sim:
See also :func:`normalBoundariesXY()` and
:func:`periodicBoundariesX()`
'''
- self.periodic[0]=1
+ self.periodic[0] = 1
def periodicBoundariesX(self):
'''
t@@ -2544,7 +2480,7 @@ class sim:
See also :func:`normalBoundariesXY()` and
:func:`periodicBoundariesXY()`
'''
- self.periodic[0]=2
+ self.periodic[0] = 2
def adaptiveGrid(self):
'''
t@@ -2555,7 +2491,7 @@ class sim:
See also :func:`staticGrid()`
'''
- self.adaptive[0]=1
+ self.adaptive[0] = 1
def staticGrid(self):
'''
t@@ -2563,9 +2499,9 @@ class sim:
See also :func:`adaptiveGrid()`
'''
- self.adaptive[0]=0
+ self.adaptive[0] = 0
- def initRandomPos(self, gridnum=numpy.array([12, 12, 36])):
+ def initRandomPos(self, gridnum=numpy.array([12, 12, 36]), dx=-1.0):
'''
Initialize particle positions in completely random configuration. Radii
*must* be set beforehand. If the x and y boundaries are set as periodic,
t@@ -2582,39 +2518,39 @@ class sim:
'''
# Calculate cells in grid
- self.num=gridnum
+ self.num = gridnum
+ r_max = numpy.max(self.radius)
# Cell configuration
if dx > 0.0:
- cellsize=dx
+ cellsize = dx
else:
- cellsize=2.1 * numpy.amax(self.radius)
+ cellsize = 2.1 * numpy.amax(self.radius)
# World size
- self.L=self.num * cellsize
-
+ self.L = self.num * cellsize
# Particle positions randomly distributed without overlap
for i in range(self.np):
- overlaps=True
+ overlaps = True
while overlaps:
- overlaps=False
+ overlaps = False
# Draw random position
for d in range(self.nd):
- self.x[i,d]=(self.L[d] - self.origo[d] - 2*r_max) \
+ self.x[i, d] = (self.L[d] - self.origo[d] - 2*r_max) \
* numpy.random.random_sample() \
+ self.origo[d] + r_max
# Check other particles for overlaps
for j in range(i-1):
- delta=self.x[i] - self.x[j]
- delta_len=math.sqrt(numpy.dot(delta,delta)) \
- - (self.radius[i] + self.radius[j])
+ delta = self.x[i] - self.x[j]
+ delta_len = math.sqrt(numpy.dot(delta, delta)) \
+ - (self.radius[i] + self.radius[j])
if delta_len < 0.0:
- overlaps=True
+ overlaps = True
print("\rFinding non-overlapping particle positions, "
- + "{0} % complete".format(numpy.ceil(i/self.np*100)))
+ + "{0} % complete".format(numpy.ceil(i/self.np*100)))
# Print newline
print()
t@@ -2639,34 +2575,33 @@ class sim:
# Cell configuration
if dx > 0.0:
- cellsize_min=dx
+ cellsize_min = dx
else:
if self.np < 1:
- raise Exception('Error: You need to define dx in '
- 'defineWorldBoundaries if there are no particles in '
- 'the simulation.')
- cellsize_min=2.1 * numpy.amax(self.radius)
+ raise Exception('Error: You need to define dx in ' +
+ 'defineWorldBoundaries if there are no ' +
+ 'particles in the simulation.')
+ cellsize_min = 2.1 * numpy.amax(self.radius)
# Lower boundary of the sorting grid
- self.origo[:]=origo[:]
+ self.origo[:] = origo[:]
# Upper boundary of the sorting grid
- self.L[:]=L[:]
+ self.L[:] = L[:]
# Adjust the number of sorting cells along each axis to fit the largest
# particle size and the world size
- self.num[0]=numpy.ceil((self.L[0]-self.origo[0])/cellsize_min)
- self.num[1]=numpy.ceil((self.L[1]-self.origo[1])/cellsize_min)
- self.num[2]=numpy.ceil((self.L[2]-self.origo[2])/cellsize_min)
+ self.num[0] = numpy.ceil((self.L[0]-self.origo[0])/cellsize_min)
+ self.num[1] = numpy.ceil((self.L[1]-self.origo[1])/cellsize_min)
+ self.num[2] = numpy.ceil((self.L[2]-self.origo[2])/cellsize_min)
#if (self.num.any() < 4):
#if (self.num[0] < 4 or self.num[1] < 4 or self.num[2] < 4):
if self.num[0] < 3 or self.num[1] < 3 or self.num[2] < 3:
raise Exception("Error: The grid must be at least 3 cells in each "
- + "direction\nGrid: x={}, y={}, z={}\n".format(\
- self.num[0], self.num[1], self.num[2])
- + "Please increase the world size.")
-
+ + "direction\nGrid: x={}, y={}, z={}\n"
+ .format(self.num[0], self.num[1], self.num[2])
+ + "Please increase the world size.")
def initGrid(self, dx=-1):
'''
t@@ -2683,22 +2618,21 @@ class sim:
# Cell configuration
if dx > 0.0:
- cellsize_min=dx
+ cellsize_min = dx
else:
- cellsize_min=2.1 * numpy.amax(self.radius)
- self.num[0]=numpy.ceil((self.L[0]-self.origo[0])/cellsize_min)
- self.num[1]=numpy.ceil((self.L[1]-self.origo[1])/cellsize_min)
- self.num[2]=numpy.ceil((self.L[2]-self.origo[2])/cellsize_min)
+ cellsize_min = 2.1 * numpy.amax(self.radius)
+ self.num[0] = numpy.ceil((self.L[0]-self.origo[0])/cellsize_min)
+ self.num[1] = numpy.ceil((self.L[1]-self.origo[1])/cellsize_min)
+ self.num[2] = numpy.ceil((self.L[2]-self.origo[2])/cellsize_min)
if self.num[0] < 4 or self.num[1] < 4 or self.num[2] < 4:
raise Exception("Error: The grid must be at least 3 cells in each "
- + "direction\nGrid: x={}, y={}, z={}".format(\
- self.num[0], self.num[1], self.num[2]))
+ + "direction\nGrid: x={}, y={}, z={}"
+ .format(self.num[0], self.num[1], self.num[2]))
# Put upper wall at top boundary
if self.nw > 0:
- self.w_x[0]=self.L[0]
-
+ self.w_x[0] = self.L[0]
def initGridAndWorldsize(self, margin=2.0):
'''
t@@ -2712,31 +2646,30 @@ class sim:
'''
# Cell configuration
- r_max=numpy.amax(self.radius)
+ r_max = numpy.amax(self.radius)
# 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[:])]) \
- - 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[:])]) \
- + margin*r_max
-
- cellsize_min=2.1 * r_max
- self.num[0]=numpy.ceil((self.L[0]-self.origo[0])/cellsize_min)
- self.num[1]=numpy.ceil((self.L[1]-self.origo[1])/cellsize_min)
- self.num[2]=numpy.ceil((self.L[2]-self.origo[2])/cellsize_min)
+ 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[:])]) \
+ - 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[:])]) \
+ + margin*r_max
+
+ cellsize_min = 2.1 * r_max
+ self.num[0] = numpy.ceil((self.L[0]-self.origo[0])/cellsize_min)
+ self.num[1] = numpy.ceil((self.L[1]-self.origo[1])/cellsize_min)
+ self.num[2] = numpy.ceil((self.L[2]-self.origo[2])/cellsize_min)
if self.num[0] < 4 or self.num[1] < 4 or self.num[2] < 4:
raise Exception("Error: The grid must be at least 3 cells in each "
- + "direction, num=" + str(self.num))
+ + "direction, num=" + str(self.num))
# Put upper wall at top boundary
if self.nw > 0:
- self.w_x[0]=self.L[0]
-
+ self.w_x[0] = self.L[0]
def initGridPos(self, gridnum=numpy.array([12, 12, 36])):
'''
t@@ -2750,19 +2683,19 @@ class sim:
'''
# Calculate cells in grid
- self.num=numpy.asarray(gridnum)
+ self.num = numpy.asarray(gridnum)
# World size
- r_max=numpy.amax(self.radius)
- cellsize=2.1 * r_max
- self.L=self.num * cellsize
+ r_max = numpy.amax(self.radius)
+ cellsize = 2.1 * r_max
+ self.L = self.num * cellsize
# Check whether there are enough grid cells
if (self.num[0]*self.num[1]*self.num[2]-(2**3)) < self.np:
print("Error! The grid is not sufficiently large.")
raise NameError('Error! The grid is not sufficiently large.')
- gridpos=numpy.zeros(self.nd, dtype=numpy.uint32)
+ gridpos = numpy.zeros(self.nd, dtype=numpy.uint32)
# Make sure grid is sufficiently large if every second level is moved
if self.periodic[0] == 1:
t@@ -2774,32 +2707,30 @@ class sim:
print("Error! The grid is not sufficiently large.")
raise NameError('Error! The grid is not sufficiently large.')
-
# Particle positions randomly distributed without overlap
for i in range(self.np):
# Find position in 3d mesh from linear index
- gridpos[0]=(i % (self.num[0]))
- gridpos[1]=numpy.floor(i/(self.num[0])) % (self.num[0])
- gridpos[2]=numpy.floor(i/((self.num[0])*(self.num[1]))) #\
+ gridpos[0] = (i % (self.num[0]))
+ gridpos[1] = numpy.floor(i/(self.num[0])) % (self.num[0])
+ gridpos[2] = numpy.floor(i/((self.num[0])*(self.num[1]))) #\
#% ((self.num[0])*(self.num[1]))
for d in range(self.nd):
- self.x[i,d]=gridpos[d] * cellsize + 0.5*cellsize
+ self.x[i, d] = gridpos[d] * cellsize + 0.5*cellsize
# Allow pushing every 2.nd level out of lateral boundaries
if self.periodic[0] == 1:
# Offset every second level
if gridpos[2] % 2:
- self.x[i,0] += 0.5*cellsize
- self.x[i,1] += 0.5*cellsize
+ self.x[i, 0] += 0.5*cellsize
+ self.x[i, 1] += 0.5*cellsize
# Readjust grid to correct size
if self.periodic[0] == 1:
self.num[0] += 1
self.num[1] += 1
-
def initRandomGridPos(self, gridnum=numpy.array([12, 12, 32]),
padding=2.1):
'''
t@@ -2819,46 +2750,47 @@ class sim:
'''
# Calculate cells in grid
- coarsegrid=numpy.floor(numpy.asarray(gridnum)/2)
+ coarsegrid = numpy.floor(numpy.asarray(gridnum)/2)
# World size
- r_max=numpy.amax(self.radius)
+ r_max = numpy.amax(self.radius)
# Cells in grid 2*size to make space for random offset
- cellsize=padding * r_max * 2
+ cellsize = padding * r_max * 2
# Check whether there are enough grid cells
if ((coarsegrid[0]-1)*(coarsegrid[1]-1)*(coarsegrid[2]-1)) < self.np:
print("Error! The grid is not sufficiently large.")
raise NameError('Error! The grid is not sufficiently large.')
- gridpos=numpy.zeros(self.nd, dtype=numpy.uint32)
+ gridpos = numpy.zeros(self.nd, dtype=numpy.uint32)
# Particle positions randomly distributed without overlap
for i in range(self.np):
# Find position in 3d mesh from linear index
- gridpos[0]=(i % (coarsegrid[0]))
- gridpos[1]=numpy.floor(i/(coarsegrid[0]))%(coarsegrid[1]) # Thanks Horacio!
- gridpos[2]=numpy.floor(i/((coarsegrid[0])*(coarsegrid[1])))
+ gridpos[0] = (i % (coarsegrid[0]))
+ gridpos[1] = numpy.floor(i/(coarsegrid[0]))%(coarsegrid[1]) # Thanks Horacio!
+ gridpos[2] = numpy.floor(i/((coarsegrid[0])*(coarsegrid[1])))
# Place particles in grid structure, and randomly adjust the
# positions within the oversized cells (uniform distribution)
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
+ r = self.radius[i]*1.05
+ self.x[i, d] = gridpos[d] * cellsize \
+ + ((cellsize-r) - r) \
+ * numpy.random.random_sample() + r
# Calculate new grid with cell size equal to max. particle diameter
- x_max=numpy.max(self.x[:,0] + self.radius)
- y_max=numpy.max(self.x[:,1] + self.radius)
- z_max=numpy.max(self.x[:,2] + self.radius)
+ x_max = numpy.max(self.x[:, 0] + self.radius)
+ y_max = numpy.max(self.x[:, 1] + self.radius)
+ z_max = numpy.max(self.x[:, 2] + self.radius)
# Adjust size of world
- self.num[0]=numpy.ceil(x_max/cellsize)
- self.num[1]=numpy.ceil(y_max/cellsize)
- self.num[2]=numpy.ceil(z_max/cellsize)
- self.L=self.num * cellsize
+ self.num[0] = numpy.ceil(x_max/cellsize)
+ self.num[1] = numpy.ceil(y_max/cellsize)
+ self.num[2] = numpy.ceil(z_max/cellsize)
+ self.L = self.num * cellsize
def createBondPair(self, i, j, spacing=-0.1):
'''
t@@ -2875,41 +2807,41 @@ class sim:
:type spacing: float
'''
- x_i=self.x[i]
- r_i=self.radius[i]
- r_j=self.radius[j]
- dist_ij=(r_i + r_j)*(1.0 + spacing)
+ x_i = self.x[i]
+ r_i = self.radius[i]
+ r_j = self.radius[j]
+ dist_ij = (r_i + r_j)*(1.0 + spacing)
- dazi=numpy.random.rand(1) * 360.0 # azimuth
- azi=numpy.radians(dazi)
- dang=numpy.random.rand(1) * 180.0 - 90.0 # angle
- ang=numpy.radians(dang)
+ dazi = numpy.random.rand(1) * 360.0 # azimuth
+ azi = numpy.radians(dazi)
+ dang = numpy.random.rand(1) * 180.0 - 90.0 # angle
+ ang = numpy.radians(dang)
- x_j=numpy.copy(x_i)
- x_j[0]=x_j[0] + dist_ij * numpy.cos(azi) * numpy.cos(ang)
- x_j[1]=x_j[1] + dist_ij * numpy.sin(azi) * numpy.cos(ang)
- x_j[2]=x_j[2] + dist_ij * numpy.sin(ang) * numpy.cos(azi)
- self.x[j]=x_j
+ x_j = numpy.copy(x_i)
+ x_j[0] = x_j[0] + dist_ij * numpy.cos(azi) * numpy.cos(ang)
+ x_j[1] = x_j[1] + dist_ij * numpy.sin(azi) * numpy.cos(ang)
+ x_j[2] = x_j[2] + dist_ij * numpy.sin(ang) * numpy.cos(azi)
+ self.x[j] = x_j
- if self.x[j,0] < self.origo[0]:
- self.x[j,0] += x_i[0] - x_j[0]
- if self.x[j,1] < self.origo[1]:
- self.x[j,1] += x_i[1] - x_j[1]
- if self.x[j,2] < self.origo[2]:
- self.x[j,2] += x_i[2] - x_j[2]
+ if self.x[j, 0] < self.origo[0]:
+ self.x[j, 0] += x_i[0] - x_j[0]
+ if self.x[j, 1] < self.origo[1]:
+ self.x[j, 1] += x_i[1] - x_j[1]
+ if self.x[j, 2] < self.origo[2]:
+ self.x[j, 2] += x_i[2] - x_j[2]
- if self.x[j,0] > self.L[0]:
- self.x[j,0] -= abs(x_j[0] - x_i[0])
- if self.x[j,1] > self.L[1]:
- self.x[j,1] -= abs(x_j[1] - x_i[1])
- if self.x[j,2] > self.L[2]:
- self.x[j,2] -= abs(x_j[2] - x_i[2])
+ if self.x[j, 0] > self.L[0]:
+ self.x[j, 0] -= abs(x_j[0] - x_i[0])
+ if self.x[j, 1] > self.L[1]:
+ self.x[j, 1] -= abs(x_j[1] - x_i[1])
+ if self.x[j, 2] > self.L[2]:
+ self.x[j, 2] -= abs(x_j[2] - x_i[2])
- self.bond(i,j) # register bond
+ self.bond(i, j) # register bond
# Check that the spacing is correct
- x_ij=self.x[i] - self.x[j]
- x_ij_length=numpy.sqrt(x_ij.dot(x_ij))
+ x_ij = self.x[i] - self.x[j]
+ x_ij_length = numpy.sqrt(x_ij.dot(x_ij))
if (x_ij_length - dist_ij) > dist_ij*0.01:
print(x_i); print(r_i)
print(x_j); print(r_j)
t@@ -2931,16 +2863,16 @@ class sim:
:type spacing: float
'''
- bondparticles=numpy.unique(\
- numpy.random.random_integers(0, high=self.np-1,\
- size=int(self.np*ratio)))
+ bondparticles = numpy.unique(numpy.random.random_integers(0, high=self.np-1,
+ size=int(self.np*ratio)))
if bondparticles.size % 2 > 0:
- bondparticles=bondparticles[:-1].copy()
- bondparticles =\
- bondparticles.reshape(int(bondparticles.size/2), 2).copy()
+ bondparticles = bondparticles[:-1].copy()
+ bondparticles = bondparticles.reshape(int(bondparticles.size/2),
+ 2).copy()
for n in numpy.arange(bondparticles.shape[0]):
- self.createBondPair(bondparticles[n,0], bondparticles[n,1], spacing)
+ self.createBondPair(bondparticles[n, 0], bondparticles[n, 1],
+ spacing)
def zeroKinematics(self):
'''
t@@ -2948,18 +2880,17 @@ class sim:
when output from one simulation is reused in another simulation.
'''
- 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, self.nd)
- self.angvel=numpy.zeros(self.np*self.nd, dtype=numpy.float64)\
- .reshape(self.np, self.nd)
- self.angpos=numpy.zeros(self.np*self.nd, 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, 3)
+ 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, self.nd)
+ self.angvel = numpy.zeros(self.np*self.nd, dtype=numpy.float64)\
+ .reshape(self.np, self.nd)
+ self.angpos = numpy.zeros(self.np*self.nd, 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, 3)
def adjustUpperWall(self, z_adjust=1.1):
'''
t@@ -2971,17 +2902,17 @@ class sim:
'''
# Initialize upper wall
- self.nw=1
- self.wmode=numpy.zeros(1) # fixed BC
- self.w_n=numpy.zeros(self.nw*self.nd, dtype=numpy.float64).reshape(\
- self.nw,self.nd)
- self.w_n[0,2]=-1.0
- self.w_vel=numpy.zeros(1)
- self.w_force=numpy.zeros(1)
- self.w_sigma0=numpy.zeros(1)
-
- self.w_x=numpy.zeros(1)
- self.w_m=numpy.zeros(1)
+ self.nw = 1
+ self.wmode = numpy.zeros(1) # fixed BC
+ self.w_n = numpy.zeros(self.nw*self.nd, dtype=numpy.float64)\
+ .reshape(self.nw, self.nd)
+ self.w_n[0, 2] = -1.0
+ self.w_vel = numpy.zeros(1)
+ self.w_force = numpy.zeros(1)
+ self.w_sigma0 = numpy.zeros(1)
+
+ self.w_x = numpy.zeros(1)
+ self.w_m = numpy.zeros(1)
self.adjustWall(idx=0, adjust=z_adjust)
def adjustWall(self, idx, adjust=1.1):
t@@ -3001,32 +2932,27 @@ class sim:
'''
if idx == 0:
- dim=2
+ dim = 2
elif idx == 1 or idx == 2:
- dim=0
+ dim = 0
elif idx == 3 or idx == 4:
- dim=1
+ dim = 1
else:
print("adjustWall: idx value not understood")
- xmin=numpy.min(self.x[:,dim] - self.radius)
- xmax=numpy.max(self.x[:,dim] + self.radius)
-
- cellsize=self.L[0] / self.num[0]
- d_max=numpy.max(self.radius)*2.0
+ xmin = numpy.min(self.x[:, dim] - self.radius)
+ xmax = numpy.max(self.x[:, dim] + self.radius)
- self.num[dim]=numpy.ceil(((xmax-xmin)*adjust + xmin)/cellsize)
- self.L[dim]=(xmax-xmin)*adjust + xmin
+ cellsize = self.L[0] / self.num[0]
+ self.num[dim] = numpy.ceil(((xmax-xmin)*adjust + xmin)/cellsize)
+ self.L[dim] = (xmax-xmin)*adjust + xmin
# Initialize upper wall
if idx == 0 or idx == 1 or idx == 3:
- self.w_x[idx]=numpy.array([xmax])
+ self.w_x[idx] = numpy.array([xmax])
else:
- self.w_x[idx]=numpy.array([xmin])
- #self.w_m[idx]=numpy.array([self.rho[0]*self.np*math.pi \
- # *(cellsize/2.0)**3])
- #self.w_m[idx]=numpy.array([self.rho*self.L[0]*self.L[1]*d_max])
- self.w_m[idx]=self.totalMass()
+ self.w_x[idx] = numpy.array([xmin])
+ self.w_m[idx] = self.totalMass()
def consolidate(self, normal_stress=10e3):
'''
t@@ -3037,11 +2963,12 @@ class sim:
:type normal_stress: float
'''
- self.nw=1
+ self.nw = 1
if normal_stress <= 0.0:
raise Exception('consolidate() error: The normal stress should be '
- 'a positive value, but is ' + str(normal_stress) + ' Pa')
+ 'a positive value, but is ' + str(normal_stress) +
+ ' Pa')
# Zero the kinematics of all particles
self.zeroKinematics()
t@@ -3050,14 +2977,14 @@ class sim:
self.adjustUpperWall()
# Set the top wall BC to a value of normal stress
- self.wmode=numpy.array([1])
- self.w_sigma0=numpy.ones(1) * normal_stress
+ self.wmode = numpy.array([1])
+ self.w_sigma0 = numpy.ones(1) * normal_stress
# Set top wall to a certain mass corresponding to the selected normal
# stress
- #self.w_sigma0=numpy.zeros(1)
- #self.w_m[0]=numpy.abs(normal_stress*self.L[0]*self.L[1]/self.g[2])
- self.w_m[0]=self.totalMass()
+ #self.w_sigma0 = numpy.zeros(1)
+ #self.w_m[0] = numpy.abs(normal_stress*self.L[0]*self.L[1]/self.g[2])
+ self.w_m[0] = self.totalMass()
def uniaxialStrainRate(self, wvel=-0.001):
'''
t@@ -3074,8 +3001,8 @@ class sim:
# Initialize upper wall
self.adjustUpperWall()
- self.wmode=numpy.array([2]) # strain rate BC
- self.w_vel=numpy.array([wvel])
+ self.wmode = numpy.array([2]) # strain rate BC
+ self.w_vel = numpy.array([wvel])
def triaxial(self, wvel=-0.001, normal_stress=10.0e3):
'''
t@@ -3094,15 +3021,16 @@ class sim:
self.zeroKinematics()
# Initialize 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
- self.w_n=numpy.array(([0,0,-1], [-1,0,0], [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)
+ 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
+ self.w_n = numpy.array(([0, 0, -1], [-1, 0, 0],
+ [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)
for i in range(5):
self.adjustWall(idx=i)
t@@ -3121,60 +3049,60 @@ class sim:
:type shear_stress: float or bool
'''
- self.nw=1
+ self.nw = 1
# Find lowest and heighest point
- z_min=numpy.min(self.x[:,2] - self.radius)
- z_max=numpy.max(self.x[:,2] + self.radius)
+ z_min = numpy.min(self.x[:, 2] - self.radius)
+ z_max = numpy.max(self.x[:, 2] + self.radius)
# the grid cell size is equal to the max. particle diameter
- cellsize=self.L[0] / self.num[0]
+ cellsize = self.L[0] / self.num[0]
# make grid one cell heigher to allow dilation
self.num[2] += 1
- self.L[2]=self.num[2] * cellsize
+ self.L[2] = self.num[2] * cellsize
# zero kinematics
self.zeroKinematics()
# Adjust grid and placement of upper wall
- self.wmode=numpy.array([1])
+ self.wmode = numpy.array([1])
# 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
- I=numpy.nonzero(self.x[:,2] < (z_min + d_max_below))
- self.fixvel[I]=1
- self.angvel[I,0]=0.0
- self.angvel[I,1]=0.0
- self.angvel[I,2]=0.0
- self.vel[I,0]=0.0 # x-dim
- self.vel[I,1]=0.0 # y-dim
- self.color[I]=-1
+ d_max_below = numpy.max(self.radius[numpy.nonzero(self.x[:, 2] <
+ (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
+ self.angvel[I, 1] = 0.0
+ self.angvel[I, 2] = 0.0
+ self.vel[I, 0] = 0.0 # x-dim
+ self.vel[I, 1] = 0.0 # y-dim
+ self.color[I] = -1
# 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
- I=numpy.nonzero(self.x[:,2] > (z_max - d_max_top))
- self.fixvel[I]=1
- self.angvel[I,0]=0.0
- self.angvel[I,1]=0.0
- self.angvel[I,2]=0.0
+ d_max_top = numpy.max(self.radius[numpy.nonzero(self.x[:, 2] >
+ (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
+ self.angvel[I, 1] = 0.0
+ self.angvel[I, 2] = 0.0
if not shear_stress:
- self.vel[I,0]=(z_max-z_min)*shear_strain_rate
+ self.vel[I, 0] = (z_max-z_min)*shear_strain_rate
else:
- self.vel[I,0]=0.0
- self.wmode[0]=3
- self.w_tau_x[0]=float(shear_stress)
- self.vel[I,1]=0.0 # y-dim
- self.color[I]=-1
+ self.vel[I, 0] = 0.0
+ self.wmode[0] = 3
+ self.w_tau_x[0] = float(shear_stress)
+ self.vel[I, 1] = 0.0 # y-dim
+ self.color[I] = -1
# Set wall tangential viscosity to zero
- self.gamma_wt[0]=0.0
+ self.gamma_wt[0] = 0.0
# Set wall friction coefficients to zero
- self.mu_ws[0]=0.0
- self.mu_wd[0]=0.0
+ self.mu_ws[0] = 0.0
+ self.mu_wd[0] = 0.0
def largestFluidTimeStep(self, safety=0.5, v_max=-1.0):
'''
t@@ -3187,8 +3115,8 @@ class sim:
(`self.mu`), and increases with fluid cell size (`self.L/self.num`)
and fluid velocities (`self.v_f`).
- NOTE: The fluid time step with the Darcy solver is an arbitrarily
- large value. In practice, this is not a problem since the short
+ NOTE: The fluid time step with the Darcy solver is an arbitrarily
+ large value. In practice, this is not a problem since the short
DEM time step is stable for fluid computations.
:param safety: Safety factor which is multiplied to the largest time
t@@ -3204,45 +3132,43 @@ class sim:
if self.fluid:
# Normalized velocities
- v_norm=numpy.empty(self.num[0]*self.num[1]*self.num[2])
- idx=0
+ v_norm = numpy.empty(self.num[0]*self.num[1]*self.num[2])
+ idx = 0
for x in numpy.arange(self.num[0]):
for y in numpy.arange(self.num[1]):
for z in numpy.arange(self.num[2]):
- v_norm[idx]=numpy.sqrt(self.v_f[x,y,z,:].dot(\
- self.v_f[x,y,z,:]))
+ v_norm[idx] = numpy.sqrt(self.v_f[x, y, z, :]\
+ .dot(self.v_f[x, y, z, :]))
idx += 1
- v_max_obs=numpy.amax(v_norm)
+ v_max_obs = numpy.amax(v_norm)
if v_max_obs == 0:
- v_max_obs=1.0e-7
+ v_max_obs = 1.0e-7
if v_max < 0.0:
- v_max=v_max_obs
+ v_max = v_max_obs
- dx_min=numpy.min(self.L/self.num)
- dt_min_cfl=dx_min/v_max
+ dx_min = numpy.min(self.L/self.num)
+ dt_min_cfl = dx_min/v_max
# Navier-Stokes
if self.cfd_solver[0] == 0:
- dt_min_von_neumann=0.5*dx_min**2/(self.mu[0] + 1.0e-16)
-
+ dt_min_von_neumann = 0.5*dx_min**2/(self.mu[0] + 1.0e-16)
return numpy.min([dt_min_von_neumann, dt_min_cfl])*safety
# Darcy
elif self.cfd_solver[0] == 1:
- self.cellSize()
return dt_min_cfl
'''
# Determine on the base of the diffusivity coefficient
# components
#self.hydraulicPermeability()
- #alpha_max=numpy.max(self.k/(self.beta_f*0.9*self.mu))
- k_max=2.7e-10 # hardcoded in darcy.cuh
- phi_min=0.1 # hardcoded in darcy.cuh
- alpha_max=k_max/(self.beta_f*phi_min*self.mu)
+ #alpha_max = numpy.max(self.k/(self.beta_f*0.9*self.mu))
+ k_max = 2.7e-10 # hardcoded in darcy.cuh
+ phi_min = 0.1 # hardcoded in darcy.cuh
+ alpha_max = k_max/(self.beta_f*phi_min*self.mu)
print(alpha_max)
return safety * 1.0/(2.0*alpha_max)*1.0/(
1.0/(self.dx[0]**2) + \
t@@ -3252,13 +3178,12 @@ class sim:
'''
# Determine value on the base of the hydraulic conductivity
- g=numpy.max(numpy.abs(self.g))
+ g = numpy.max(numpy.abs(self.g))
# Bulk modulus of fluid
- K=1.0/self.beta_f[0]
+ K = 1.0/self.beta_f[0]
self.hydraulicDiffusivity()
- self.cellSize()
return safety * 1.0/(2.0*self.D)*1.0/( \
1.0/(self.dx[0]**2) + \
t@@ -3266,13 +3191,6 @@ class sim:
1.0/(self.dx[2]**2))
'''
- def cellSize(self):
- '''
- Calculate the particle sorting (and fluid) cell dimensions.
- These values are stored in `self.dx` and are NOT returned.
- '''
- self.dx=self.L/self.num
-
def hydraulicConductivity(self, phi=0.35):
'''
Determine the hydraulic conductivity (K) [m/s] from the permeability
t@@ -3285,13 +3203,13 @@ class sim:
:return type: float
'''
if self.cfd_solver[0] == 1:
- k=self.k_c * phi**3/(1.0 - phi**2)
- self.K_c=k*self.rho_f*numpy.abs(self.g[2])/self.mu
+ k = self.k_c * phi**3/(1.0 - phi**2)
+ self.K_c = k*self.rho_f*numpy.abs(self.g[2])/self.mu
return self.K_c[0]
else:
raise Exception('This function only works for the Darcy solver')
- def hydraulicPermeability(self, phi_min=0.3, phi_max=0.99):
+ def hydraulicPermeability(self):
'''
Determine the hydraulic permeability (k) [m*m] from the Kozeny-Carman
relationship, using the permeability prefactor (`self.k_c`), and the
t@@ -3313,9 +3231,10 @@ class sim:
(`self.cfd_solver[0] == 1`)
'''
if self.cfd_solver[0] == 1:
- self.hydraulicConductivity()
- phi_bar=numpy.mean(self.phi)
- self.D=self.K_c/(self.rho_f*g*(self.k_n[0] + phi_bar*K))
+ self.hydraulicConductivity()
+ phi_bar = numpy.mean(self.phi)
+ self.D = self.K_c/(self.rho_f*self.g[2]
+ *(self.k_n[0] + phi_bar*self.K))
else:
raise Exception('This function only works for the Darcy solver')
t@@ -3345,29 +3264,29 @@ class sim:
'''
if dt > 0.0:
- self.time_dt[0]=dt
+ self.time_dt[0] = dt
if self.np > 0:
print("Warning: Manually specifying the time step length when "
- + "simulating particles may produce instabilities.")
+ "simulating particles may produce instabilities.")
elif self.np > 0:
- r_min=numpy.min(self.radius)
- m_min=self.rho * 4.0/3.0*numpy.pi*r_min**3
+ r_min = numpy.min(self.radius)
+ m_min = self.rho * 4.0/3.0*numpy.pi*r_min**3
if self.E > 0.001:
- k_max=numpy.max(numpy.pi/2.0*self.E*self.radius)
+ k_max = numpy.max(numpy.pi/2.0*self.E*self.radius)
else:
- k_max=numpy.max([self.k_n[:], self.k_t[:]])
+ k_max = numpy.max([self.k_n[:], self.k_t[:]])
# Radjaii et al 2011
- self.time_dt[0]=epsilon/(numpy.sqrt(k_max/m_min))
+ self.time_dt[0] = epsilon/(numpy.sqrt(k_max/m_min))
# Zhang and Campbell, 1992
- #self.time_dt[0]=0.075*math.sqrt(m_min/k_max)
+ #self.time_dt[0] = 0.075*math.sqrt(m_min/k_max)
# Computational time step (O'Sullivan et al, 2003)
- #self.time_dt[0]=0.17*math.sqrt(m_min/k_max)
+ #self.time_dt[0] = 0.17*math.sqrt(m_min/k_max)
elif not self.fluid:
raise Exception('Error: Could not automatically set a time step.')
t@@ -3376,15 +3295,14 @@ class sim:
# by von Neumann stability analysis of the diffusion and advection
# terms
if self.fluid:
- fluid_time_dt=self.largestFluidTimeStep()
- self.time_dt[0]=numpy.min([fluid_time_dt, self.time_dt[0]])
-
+ fluid_time_dt = self.largestFluidTimeStep()
+ self.time_dt[0] = numpy.min([fluid_time_dt, self.time_dt[0]])
# Time at start
- self.time_current[0]=current
- self.time_total[0]=total
- self.time_file_dt[0]=file_dt
- self.time_step_count[0]=0
+ self.time_current[0] = current
+ self.time_total[0] = total
+ self.time_file_dt[0] = file_dt
+ self.time_step_count[0] = step_count
def dry(self):
'''
t@@ -3392,7 +3310,7 @@ class sim:
See also :func:`wet()`
'''
- self.fluid=False
+ self.fluid = False
def wet(self):
'''
t@@ -3400,7 +3318,7 @@ class sim:
See also :func:`dry()`
'''
- self.fluid=True
+ self.fluid = True
self.initFluid()
def initFluid(self, mu=8.9e-4, rho=1.0e3, p=0.0, hydrostatic=False,
t@@ -3426,93 +3344,92 @@ class sim:
Accepted values: 0 (Navier Stokes, default) and 1 (Darcy).
:type cfd_solver: int
'''
- self.fluid=True
- self.mu=numpy.ones(1, dtype=numpy.float64) * mu
- self.rho_f=numpy.ones(1, dtype=numpy.float64) * rho
+ self.fluid = True
+ self.mu = numpy.ones(1, dtype=numpy.float64) * mu
+ 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
+ self.p_f = numpy.ones((self.num[0], self.num[1], self.num[2]),
+ dtype=numpy.float64) * p
if hydrostatic:
- dz=self.L[2]/self.num[2]
+ dz = self.L[2]/self.num[2]
# Zero pressure gradient from grid top to top wall, linear pressure
# distribution from top wall to grid bottom
if self.nw == 1:
- wall0_iz=int(self.w_x[0]/(self.L[2]/self.num[2]))
- self.p_f[:,:,wall0_iz:]=p
+ wall0_iz = int(self.w_x[0]/(self.L[2]/self.num[2]))
+ self.p_f[:, :, wall0_iz:] = p
for iz in numpy.arange(wall0_iz - 1):
- z=dz*iz + 0.5*dz
- depth=self.w_x[0] - z
- self.p_f[:,:,iz]=p + (depth-dz) * rho * -self.g[2]
+ z = dz*iz + 0.5*dz
+ depth = self.w_x[0] - z
+ self.p_f[:, :, iz] = p + (depth-dz) * rho * -self.g[2]
# Linear pressure distribution from grid top to grid bottom
else:
for iz in numpy.arange(self.num[2] - 1):
- z=dz*iz + 0.5*dz
- depth=self.L[2] - z
- self.p_f[:,:,iz]=p + (depth-dz) * rho * -self.g[2]
+ z = dz*iz + 0.5*dz
+ depth = self.L[2] - z
+ 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),
- dtype=numpy.float64)
- self.phi=numpy.ones((self.num[0], self.num[1], self.num[2]),
- dtype=numpy.float64)
- self.dphi=numpy.zeros((self.num[0], self.num[1], self.num[2]),
- dtype=numpy.float64)
+ 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)
+ self.dphi = numpy.zeros((self.num[0], self.num[1], self.num[2]),
+ 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]
- self.p_mod_phi=numpy.zeros(1, dtype=numpy.float64) # Shift [rad]
+ self.p_mod_A = numpy.zeros(1, dtype=numpy.float64) # Amplitude [Pa]
+ self.p_mod_f = numpy.zeros(1, dtype=numpy.float64) # Frequency [Hz]
+ self.p_mod_phi = numpy.zeros(1, dtype=numpy.float64) # Shift [rad]
- self.bc_bot=numpy.zeros(1, dtype=numpy.int32)
- self.bc_top=numpy.zeros(1, dtype=numpy.int32)
- self.free_slip_bot=numpy.ones(1, dtype=numpy.int32)
- self.free_slip_top=numpy.ones(1, dtype=numpy.int32)
- self.bc_bot_flux=numpy.zeros(1, dtype=numpy.float64)
- self.bc_top_flux=numpy.zeros(1, dtype=numpy.float64)
+ self.bc_bot = numpy.zeros(1, dtype=numpy.int32)
+ self.bc_top = numpy.zeros(1, dtype=numpy.int32)
+ self.free_slip_bot = numpy.ones(1, dtype=numpy.int32)
+ self.free_slip_top = numpy.ones(1, dtype=numpy.int32)
+ self.bc_bot_flux = numpy.zeros(1, dtype=numpy.float64)
+ 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)
+ self.p_f_constant = numpy.zeros((self.num[0], self.num[1], self.num[2]),
+ dtype=numpy.int32)
# Fluid solver type
# 0: Navier Stokes (fluid with inertia)
# 1: Stokes-Darcy (fluid without inertia)
- self.cfd_solver=numpy.ones(1)*cfd_solver
+ self.cfd_solver = numpy.ones(1)*cfd_solver
if self.cfd_solver[0] == 0:
- self.gamma=numpy.array(0.0)
- self.theta=numpy.array(1.0)
- self.beta=numpy.array(0.0)
- 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.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, 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)
+ self.gamma = numpy.array(0.0)
+ self.theta = numpy.array(1.0)
+ self.beta = numpy.array(0.0)
+ 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.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, 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, self.nd), dtype=numpy.float64)
- self.beta_f=numpy.ones(1, dtype=numpy.float64)*4.5e-10
- 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
- self.bc_xp=numpy.ones(1, dtype=numpy.int32)*2
- self.bc_yn=numpy.ones(1, dtype=numpy.int32)*2
- self.bc_yp=numpy.ones(1, dtype=numpy.int32)*2
+ 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, self.nd), dtype=numpy.float64)
+ self.beta_f = numpy.ones(1, dtype=numpy.float64)*4.5e-10
+ 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
+ self.bc_xp = numpy.ones(1, dtype=numpy.int32)*2
+ self.bc_yn = numpy.ones(1, dtype=numpy.int32)*2
+ self.bc_yp = numpy.ones(1, dtype=numpy.int32)*2
else:
raise Exception('Value of cfd_solver not understood (' + \
t@@ -3531,7 +3448,7 @@ class sim:
:return type: float
'''
if value != -1:
- self.time_current[0]=value
+ self.time_current[0] = value
else:
return self.time_current[0]
t@@ -3543,7 +3460,7 @@ class sim:
The default behavior for the boundary is fixed value (Dirichlet), see
:func:`setFluidBottomFixedPressure()`.
'''
- self.bc_bot[0]=1
+ self.bc_bot[0] = 1
def setFluidBottomNoFlowNoSlip(self):
'''
t@@ -3553,7 +3470,7 @@ class sim:
The default behavior for the boundary is fixed value (Dirichlet), see
:func:`setFluidBottomFixedPressure()`.
'''
- self.bc_bot[0]=2
+ self.bc_bot[0] = 2
def setFluidBottomFixedPressure(self):
'''
t@@ -3563,7 +3480,7 @@ class sim:
This is the default behavior for the boundary. See also
:func:`setFluidBottomNoFlow()`
'''
- self.bc_bot[0]=0
+ self.bc_bot[0] = 0
def setFluidBottomFixedFlux(self, specific_flux):
'''
t@@ -3575,8 +3492,8 @@ class sim:
:param specific_flux: Specific flux values across boundary (positive
values upwards), [m/s]
'''
- self.bc_bot[0]=4
- self.bc_bot_flux[0]=specific_flux
+ self.bc_bot[0] = 4
+ self.bc_bot_flux[0] = specific_flux
def setFluidTopNoFlow(self):
'''
t@@ -3586,7 +3503,7 @@ class sim:
The default behavior for the boundary is fixed value (Dirichlet), see
:func:`setFluidTopFixedPressure()`.
'''
- self.bc_top[0]=1
+ self.bc_top[0] = 1
def setFluidTopNoFlowNoSlip(self):
'''
t@@ -3596,7 +3513,7 @@ class sim:
The default behavior for the boundary is fixed value (Dirichlet), see
:func:`setFluidTopFixedPressure()`.
'''
- self.bc_top[0]=2
+ self.bc_top[0] = 2
def setFluidTopFixedPressure(self):
'''
t@@ -3606,7 +3523,7 @@ class sim:
This is the default behavior for the boundary. See also
:func:`setFluidTopNoFlow()`
'''
- self.bc_top[0]=0
+ self.bc_top[0] = 0
def setFluidTopFixedFlux(self, specific_flux):
'''
t@@ -3618,8 +3535,8 @@ class sim:
:param specific_flux: Specific flux values across boundary (positive
values upwards), [m/s]
'''
- self.bc_top[0]=4
- self.bc_top_flux[0]=specific_flux
+ self.bc_top[0] = 4
+ self.bc_top_flux[0] = specific_flux
def setFluidXFixedPressure(self):
'''
t@@ -3631,8 +3548,8 @@ class sim:
:func:`setFluidXNoFlow()`, and
:func:`setFluidXPeriodic()` (default)
'''
- self.bc_xn[0]=0
- self.bc_xp[0]=0
+ self.bc_xn[0] = 0
+ self.bc_xp[0] = 0
def setFluidXNoFlow(self):
'''
t@@ -3644,8 +3561,8 @@ class sim:
:func:`setFluidXNoFlow()`, and
:func:`setFluidXPeriodic()` (default)
'''
- self.bc_xn[0]=1
- self.bc_xp[0]=1
+ self.bc_xn[0] = 1
+ self.bc_xp[0] = 1
def setFluidXPeriodic(self):
'''
t@@ -3656,8 +3573,8 @@ class sim:
:func:`setFluidXFixedPressure()` and
:func:`setFluidXNoFlow()`
'''
- self.bc_xn[0]=2
- self.bc_xp[0]=2
+ self.bc_xn[0] = 2
+ self.bc_xp[0] = 2
def setFluidYFixedPressure(self):
'''
t@@ -3668,8 +3585,8 @@ class sim:
:func:`setFluidYNoFlow()` and
:func:`setFluidYPeriodic()` (default)
'''
- self.bc_yn[0]=0
- self.bc_yp[0]=0
+ self.bc_yn[0] = 0
+ self.bc_yp[0] = 0
def setFluidYNoFlow(self):
'''
t@@ -3680,8 +3597,8 @@ class sim:
:func:`setFluidYFixedPressure()` and
:func:`setFluidYPeriodic()` (default)
'''
- self.bc_yn[0]=1
- self.bc_yp[0]=1
+ self.bc_yn[0] = 1
+ self.bc_yp[0] = 1
def setFluidYPeriodic(self):
'''
t@@ -3692,8 +3609,8 @@ class sim:
:func:`setFluidYFixedPressure()` and
:func:`setFluidYNoFlow()`
'''
- self.bc_yn[0]=2
- self.bc_yp[0]=2
+ self.bc_yn[0] = 2
+ self.bc_yp[0] = 2
def setPermeabilityGrainSize(self, verbose=True):
'''
t@@ -3705,9 +3622,8 @@ class sim:
prefactor.
:type verbose: bool
'''
- self.setPermeabilityPrefactor(
- k_c=numpy.mean(self.radius*2.0)**2.0/180.0,
- verbose=verbose)
+ self.setPermeabilityPrefactor(k_c=numpy.mean(self.radius*2.0)**2.0/180.0,
+ verbose=verbose)
def setPermeabilityPrefactor(self, k_c, verbose=True):
'''
t@@ -3723,11 +3639,11 @@ class sim:
:type verbose: bool
'''
if self.cfd_solver[0] == 1:
- self.k_c[0]=k_c
+ self.k_c[0] = k_c
if verbose:
- phi=numpy.array([0.1, 0.35, 0.9])
- k=self.k_c * phi**3/(1.0 - phi**2)
- K=k * self.rho_f*numpy.abs(self.g[2])/self.mu
+ phi = numpy.array([0.1, 0.35, 0.9])
+ k = self.k_c * phi**3/(1.0 - phi**2)
+ K = k * self.rho_f*numpy.abs(self.g[2])/self.mu
print('Hydraulic permeability limits for porosity phi=' + \
str(phi) + ':')
print('\tk=' + str(k) + ' m*m')
t@@ -3736,7 +3652,7 @@ class sim:
print('\tK=' + str(K) + ' m/s')
else:
raise Exception('setPermeabilityPrefactor() only relevant for the '
- + 'Darcy solver (cfd_solver=1)')
+ 'Darcy solver (cfd_solver=1)')
def findPermeabilities(self):
'''
t@@ -3748,11 +3664,11 @@ class sim:
`self.readbin`.
'''
if self.cfd_solver[0] == 1:
- phi=numpy.clip(self.phi, 0.1, 0.9)
- self.k=self.k_c * phi**3/(1.0 - phi**2)
+ phi = numpy.clip(self.phi, 0.1, 0.9)
+ self.k = self.k_c * phi**3/(1.0 - phi**2)
else:
raise Exception('findPermeabilities() only relevant for the '
- + 'Darcy solver (cfd_solver=1)')
+ 'Darcy solver (cfd_solver=1)')
def findHydraulicConductivities(self):
'''
t@@ -3765,10 +3681,10 @@ class sim:
'''
if self.cfd_solver[0] == 1:
self.findPermeabilities()
- self.K=self.k*self.rho_f*numpy.abs(self.g[2])/self.mu
+ self.K = self.k*self.rho_f*numpy.abs(self.g[2])/self.mu
else:
raise Exception('findPermeabilities() only relevant for the '
- + 'Darcy solver (cfd_solver=1)')
+ 'Darcy solver (cfd_solver=1)')
def defaultParams(self, mu_s=0.5, mu_d=0.5, mu_r=0.0, rho=2600, k_n=1.16e9,
k_t=1.16e9, k_r=0, gamma_n=0.0, gamma_t=0.0, gamma_r=0.0,
t@@ -3806,19 +3722,19 @@ class sim:
'''
# Particle material density, kg/m^3
- self.rho=numpy.ones(1, dtype=numpy.float64) * rho
+ self.rho = numpy.ones(1, dtype=numpy.float64) * rho
### Dry granular material parameters
# Contact normal elastic stiffness, N/m
- self.k_n=numpy.ones(1, dtype=numpy.float64) * k_n
+ self.k_n = numpy.ones(1, dtype=numpy.float64) * k_n
# Contact shear elastic stiffness (for contactmodel=2), N/m
- self.k_t=numpy.ones(1, dtype=numpy.float64) * k_t
+ self.k_t = numpy.ones(1, dtype=numpy.float64) * k_t
# Contact rolling elastic stiffness (for contactmodel=2), N/m
- self.k_r=numpy.ones(1, dtype=numpy.float64) * k_r
+ self.k_r = numpy.ones(1, dtype=numpy.float64) * k_r
# Contact normal viscosity. Critical damping: 2*sqrt(m*k_n).
# Normal force component elastic if nu=0.0.
t@@ -3826,51 +3742,51 @@ class sim:
# * nu_frac * 2.0 * math.sqrt(4.0/3.0 * math.pi \
# * numpy.amin(self.radius)**3 \
# * self.rho[0] * self.k_n[0])
- self.gamma_n=numpy.ones(1, dtype=numpy.float64) * gamma_n
+ self.gamma_n = numpy.ones(1, dtype=numpy.float64) * gamma_n
# Contact shear viscosity, Ns/m
- self.gamma_t=numpy.ones(1, dtype=numpy.float64) * gamma_t
+ self.gamma_t = numpy.ones(1, dtype=numpy.float64) * gamma_t
# Contact rolling viscosity, Ns/m?
- self.gamma_r=numpy.ones(1, dtype=numpy.float64) * gamma_r
+ self.gamma_r = numpy.ones(1, dtype=numpy.float64) * gamma_r
# Contact static shear friction coefficient
- #self.mu_s=numpy.ones(1, dtype=numpy.float64) * \
+ #self.mu_s = numpy.ones(1, dtype=numpy.float64) * \
#numpy.tan(numpy.radians(ang_s))
- self.mu_s=numpy.ones(1, dtype=numpy.float64) * mu_s
+ self.mu_s = numpy.ones(1, dtype=numpy.float64) * mu_s
# Contact dynamic shear friction coefficient
- #self.mu_d=numpy.ones(1, dtype=numpy.float64) * \
+ #self.mu_d = numpy.ones(1, dtype=numpy.float64) * \
#numpy.tan(numpy.radians(ang_d))
- self.mu_d=numpy.ones(1, dtype=numpy.float64) * mu_d
+ self.mu_d = numpy.ones(1, dtype=numpy.float64) * mu_d
# Contact rolling friction coefficient
- #self.mu_r=numpy.ones(1, dtype=numpy.float64) * \
+ #self.mu_r = numpy.ones(1, dtype=numpy.float64) * \
#numpy.tan(numpy.radians(ang_r))
- self.mu_r=numpy.ones(1, dtype=numpy.float64) * mu_r
+ self.mu_r = numpy.ones(1, dtype=numpy.float64) * mu_r
# Wall viscosities
- self.gamma_wn[0]=gamma_wn # normal
- self.gamma_wt[0]=gamma_wt # sliding
+ self.gamma_wn[0] = gamma_wn # normal
+ self.gamma_wt[0] = gamma_wt # sliding
# Wall friction coefficients
- self.mu_ws=self.mu_s # static
- self.mu_wd=self.mu_d # dynamic
+ self.mu_ws = self.mu_s # static
+ self.mu_wd = self.mu_d # dynamic
### Parameters related to capillary bonds
# Wettability, 0=perfect
- theta=0.0;
+ theta = 0.0
if capillaryCohesion == 1:
# Prefactor
- self.kappa[0]=2.0 * math.pi * gamma_t * numpy.cos(theta)
- self.V_b[0]=1e-12 # Liquid volume at bond
- else :
- self.kappa[0]=0.0; # Zero capillary force
- self.V_b[0]=0.0; # Zero liquid volume at bond
+ self.kappa[0] = 2.0 * math.pi * gamma_t * numpy.cos(theta)
+ self.V_b[0] = 1e-12 # Liquid volume at bond
+ else:
+ self.kappa[0] = 0.0 # Zero capillary force
+ self.V_b[0] = 0.0 # Zero liquid volume at bond
# Debonding distance
- self.db[0]=(1.0 + theta/2.0) * self.V_b**(1.0/3.0)
+ self.db[0] = (1.0 + theta/2.0) * self.V_b**(1.0/3.0)
def setStiffnessNormal(self, k_n):
'''
t@@ -3880,7 +3796,7 @@ class sim:
:param k_n: The elastic stiffness coefficient [N/m]
:type k_n: float
'''
- self.k_n[0]=k_n
+ self.k_n[0] = k_n
def setStiffnessTangential(self, k_t):
'''
t@@ -3890,7 +3806,7 @@ class sim:
:param k_t: The elastic stiffness coefficient [N/m]
:type k_t: float
'''
- self.k_t[0]=k_t
+ self.k_t[0] = k_t
def setYoungsModulus(self, E):
'''
t@@ -3905,7 +3821,7 @@ class sim:
:param E: The elastic modulus [Pa]
:type E: float
'''
- self.E[0]=E
+ self.E[0] = E
def setDampingNormal(self, gamma, over_damping=False):
'''
t@@ -3920,9 +3836,9 @@ class sim:
See also: :func:`setDampingTangential(gamma)`
'''
- self.gamma_n[0]=gamma
- critical_gamma=2.0*numpy.sqrt(self.smallestMass()*self.k_n[0])
- damping_ratio=gamma/critical_gamma
+ self.gamma_n[0] = gamma
+ critical_gamma = 2.0*numpy.sqrt(self.smallestMass()*self.k_n[0])
+ damping_ratio = gamma/critical_gamma
if damping_ratio < 1.0:
print('Info: The system is under-dampened (ratio='
+ str(damping_ratio)
t@@ -3931,18 +3847,20 @@ class sim:
elif damping_ratio > 1.0:
if over_damping:
print('Warning: The system is over-dampened (ratio='
- + str(damping_ratio) + ') in the normal component. '
- + '\nCritical damping=' + str(critical_gamma) + '.')
+ + str(damping_ratio) + ') in the normal component. '
+ '\nCritical damping=' + str(critical_gamma) + '.')
else:
raise Exception('Warning: The system is over-dampened (ratio='
- + str(damping_ratio) + ') in the normal component.\n'
- + 'Call this function once more with `over_damping=True` '
- + 'if this is what you want. \nCritical damping='
- + str(critical_gamma) + '.')
+ + str(damping_ratio) + ') in the normal '
+ 'component.\n'
+ 'Call this function once more with '
+ '`over_damping=True` if this is what you want.'
+ '\nCritical damping=' + str(critical_gamma) +
+ '.')
else:
- print('Warning: The system is critically dampened (ratio='
- + str(damping_ratio) + ') in the normal component. '
- + '\nCritical damping=' + str(critical_gamma) + '.')
+ print('Warning: The system is critically dampened (ratio=' +
+ str(damping_ratio) + ') in the normal component. '
+ '\nCritical damping=' + str(critical_gamma) + '.')
def setDampingTangential(self, gamma, over_damping=False):
'''
t@@ -3957,8 +3875,8 @@ class sim:
See also: :func:`setDampingNormal(gamma)`
'''
- self.gamma_t[0]=gamma
- damping_ratio=gamma/(2.0*numpy.sqrt(self.smallestMass()*self.k_t[0]))
+ self.gamma_t[0] = gamma
+ damping_ratio = gamma/(2.0*numpy.sqrt(self.smallestMass()*self.k_t[0]))
if damping_ratio < 1.0:
print('Info: The system is under-dampened (ratio='
+ str(damping_ratio)
t@@ -3966,12 +3884,13 @@ class sim:
elif damping_ratio > 1.0:
if over_damping:
print('Warning: The system is over-dampened (ratio='
- + str(damping_ratio) + ') in the tangential component.')
+ + str(damping_ratio) + ') in the tangential component.')
else:
raise Exception('Warning: The system is over-dampened (ratio='
- + str(damping_ratio) + ') in the tangential component.\n'
- + 'Call this function once more with `over_damping=True` '
- + 'if this is what you want.')
+ + str(damping_ratio) + ') in the tangential '
+ 'component.\n'
+ 'Call this function once more with '
+ '`over_damping=True` if this is what you want.')
else:
print('Warning: The system is critically dampened (ratio='
+ str(damping_ratio) + ') in the tangential component.')
t@@ -3988,7 +3907,7 @@ class sim:
See also: :func:`setDynamicFriction(mu_d)`
'''
- self.mu_s[0]=mu_s
+ self.mu_s[0] = mu_s
def setDynamicFriction(self, mu_d):
'''
t@@ -4005,7 +3924,7 @@ class sim:
See also: :func:`setStaticFriction(mu_s)`
'''
- self.mu_d[0]=mu_d
+ self.mu_d[0] = mu_d
def setFluidCompressibility(self, beta_f):
'''
t@@ -4020,10 +3939,10 @@ class sim:
See also: :func:`setFluidDensity()` and :func:`setFluidViscosity()`
'''
if self.cfd_solver[0] == 1:
- self.beta_f[0]=beta_f
+ self.beta_f[0] = beta_f
else:
raise Exception('setFluidCompressibility() only relevant for the '
- + 'Darcy solver (cfd_solver=1)')
+ 'Darcy solver (cfd_solver=1)')
def setFluidViscosity(self, mu):
'''
t@@ -4037,7 +3956,7 @@ class sim:
See also: :func:`setFluidDensity()` and
:func:`setFluidCompressibility()`
'''
- self.mu[0]=mu
+ self.mu[0] = mu
def setFluidDensity(self, rho_f):
'''
t@@ -4050,7 +3969,7 @@ class sim:
See also: :func:`setFluidViscosity()` and
:func:`setFluidCompressibility()`
'''
- self.rho_f[0]=rho_f
+ self.rho_f[0] = rho_f
def scaleSize(self, factor):
'''
t@@ -4082,38 +4001,36 @@ class sim:
:type j: int
'''
- self.lambda_bar[0]=1.0 # Radius multiplier to parallel-bond radii
+ self.lambda_bar[0] = 1.0 # Radius multiplier to parallel-bond radii
if not hasattr(self, 'bonds'):
- self.bonds=numpy.array([[i,j]], dtype=numpy.uint32)
- else :
- self.bonds=numpy.vstack((self.bonds, [i,j]))
+ self.bonds = numpy.array([[i, j]], dtype=numpy.uint32)
+ else:
+ self.bonds = numpy.vstack((self.bonds, [i, j]))
if not hasattr(self, 'bonds_delta_n'):
- self.bonds_delta_n=numpy.array([0.0], dtype=numpy.uint32)
- else :
- #self.bonds_delta_n=numpy.vstack((self.bonds_delta_n, [0.0]))
- self.bonds_delta_n=numpy.append(self.bonds_delta_n, [0.0])
+ self.bonds_delta_n = numpy.array([0.0], dtype=numpy.uint32)
+ else:
+ #self.bonds_delta_n = numpy.vstack((self.bonds_delta_n, [0.0]))
+ self.bonds_delta_n = numpy.append(self.bonds_delta_n, [0.0])
if not hasattr(self, 'bonds_delta_t'):
- 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]))
+ 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]))
if not hasattr(self, 'bonds_omega_n'):
- self.bonds_omega_n=numpy.array([0.0], dtype=numpy.uint32)
- else :
- #self.bonds_omega_n=numpy.vstack((self.bonds_omega_n, [0.0]))
- self.bonds_omega_n=numpy.append(self.bonds_omega_n, [0.0])
+ self.bonds_omega_n = numpy.array([0.0], dtype=numpy.uint32)
+ else:
+ #self.bonds_omega_n = numpy.vstack((self.bonds_omega_n, [0.0]))
+ 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)
- else :
- self.bonds_omega_t=numpy.vstack((self.bonds_omega_t,\
- [0.0, 0.0, 0.0]))
+ 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]))
# Increment the number of bonds with one
self.nb0 += 1
t@@ -4178,7 +4095,7 @@ class sim:
:returns: The total mass in [kg]
'''
- m=0.0
+ m = 0.0
for i in range(self.np):
m += self.mass(i)
return m
t@@ -4226,7 +4143,7 @@ class sim:
:return type: float
'''
return 0.5*self.mass(idx) \
- *numpy.sqrt(numpy.dot(self.vel[idx,:], self.vel[idx,:]))**2
+ *numpy.sqrt(numpy.dot(self.vel[idx, :], self.vel[idx, :]))**2
def totalKineticEnergy(self):
'''
t@@ -4234,7 +4151,7 @@ class sim:
:returns: The kinetic energy of all particles [J]
'''
- esum=0.0
+ esum = 0.0
for i in range(self.np):
esum += self.kineticEnergy(i)
return esum
t@@ -4249,7 +4166,7 @@ class sim:
:return type: float
'''
return 0.5*self.momentOfInertia(idx) \
- *numpy.sqrt(numpy.dot(self.angvel[idx,:], self.angvel[idx,:]))**2
+ *numpy.sqrt(numpy.dot(self.angvel[idx, :], self.angvel[idx, :]))**2
def totalRotationalEnergy(self):
'''
t@@ -4257,7 +4174,7 @@ class sim:
:returns: The rotational energy of all particles [J]
'''
- esum=0.0
+ esum = 0.0
for i in range(self.np):
esum += self.rotationalEnergy(i)
return esum
t@@ -4280,7 +4197,7 @@ class sim:
:returns: The normal viscous energy lost by all particles [J]
:return type: float
'''
- esum=0.0
+ esum = 0.0
for i in range(self.np):
esum += self.viscousEnergy(i)
return esum
t@@ -4303,7 +4220,7 @@ class sim:
:returns: The total frictional energy lost of all particles [J]
:return type: float
'''
- esum=0.0
+ esum = 0.0
for i in range(self.np):
esum += self.frictionalEnergy(i)
return esum
t@@ -4325,24 +4242,24 @@ class sim:
'''
if method == 'pot':
- m=numpy.ones(self.np)*4.0/3.0*math.pi*self.radius**3*self.rho
- return numpy.sum(m*math.sqrt(numpy.dot(self.g,self.g))*self.x[:,2])
+ m = numpy.ones(self.np)*4.0/3.0*math.pi*self.radius**3*self.rho
+ return numpy.sum(m*math.sqrt(numpy.dot(self.g, self.g))*self.x[:, 2])
elif method == 'kin':
- m=numpy.ones(self.np)*4.0/3.0*math.pi*self.radius**3*self.rho
- esum=0.0
+ m = numpy.ones(self.np)*4.0/3.0*math.pi*self.radius**3*self.rho
+ esum = 0.0
for i in range(self.np):
esum += 0.5*m[i]*math.sqrt(\
- numpy.dot(self.vel[i,:],self.vel[i,:]))**2
+ numpy.dot(self.vel[i, :], self.vel[i, :]))**2
return esum
elif method == 'rot':
- m=numpy.ones(self.np)*4.0/3.0*math.pi*self.radius**3*self.rho
- esum=0.0
+ m = numpy.ones(self.np)*4.0/3.0*math.pi*self.radius**3*self.rho
+ esum = 0.0
for i in range(self.np):
esum += 0.5*2.0/5.0*m[i]*self.radius[i]**2 \
*math.sqrt(\
- numpy.dot(self.angvel[i,:],self.angvel[i,:]))**2
+ numpy.dot(self.angvel[i, :], self.angvel[i, :]))**2
return esum
elif method == 'shear':
t@@ -4359,22 +4276,22 @@ class sim:
elif method == 'bondpot':
if self.nb0 > 0:
- R_bar=self.lambda_bar*numpy.minimum(\
- self.radius[self.bonds[:,0]],\
- self.radius[self.bonds[:,1]])
- A=numpy.pi*R_bar**2
- I=0.25*numpy.pi*R_bar**4
- J=I*2.0
- bondpot_fn=numpy.sum(\
+ R_bar = self.lambda_bar*numpy.minimum(\
+ self.radius[self.bonds[:, 0]],\
+ self.radius[self.bonds[:, 1]])
+ A = numpy.pi*R_bar**2
+ I = 0.25*numpy.pi*R_bar**4
+ J = I*2.0
+ bondpot_fn = numpy.sum(\
0.5*A*self.k_n*numpy.abs(self.bonds_delta_n)**2)
- bondpot_ft=numpy.sum(\
+ bondpot_ft = numpy.sum(\
0.5*A*self.k_t*numpy.linalg.norm(self.bonds_delta_t)**2)
- bondpot_tn=numpy.sum(\
+ bondpot_tn = numpy.sum(\
0.5*J*self.k_t*numpy.abs(self.bonds_omega_n)**2)
- bondpot_tt=numpy.sum(\
+ bondpot_tt = numpy.sum(\
0.5*I*self.k_n*numpy.linalg.norm(self.bonds_omega_t)**2)
return bondpot_fn + bondpot_ft + bondpot_tn + bondpot_tt
- else :
+ else:
return 0.0
else:
raise Exception('Unknownw energy() method "' + method + '"')
t@@ -4388,15 +4305,15 @@ class sim:
'''
# Find the bulk volume
- V_t=(self.L[0] - self.origo[0]) \
+ V_t = (self.L[0] - self.origo[0]) \
*(self.L[1] - self.origo[1]) \
*(self.w_x[0] - self.origo[2])
# Find the volume of solids
- V_s=numpy.sum(4.0/3.0 * math.pi * self.radius**3)
+ V_s = numpy.sum(4.0/3.0 * math.pi * self.radius**3)
# Return the void ratio
- e=(V_t - V_s)/V_s
+ e = (V_t - V_s)/V_s
return e
def bulkPorosity(self, trim=True):
t@@ -4411,27 +4328,26 @@ class sim:
:return type: float
'''
- V_total=0.0
+ V_total = 0.0
if trim:
- min_x=numpy.min(self.x[:,0] - self.radius)
- min_y=numpy.min(self.x[:,1] - self.radius)
- min_z=numpy.min(self.x[:,2] - self.radius)
- max_x=numpy.max(self.x[:,0] + self.radius)
- max_y=numpy.max(self.x[:,1] + self.radius)
- max_z=numpy.max(self.x[:,2] + self.radius)
- V_total=(max_x - min_x)*(max_y - min_y)*(max_z - min_z)
+ min_x = numpy.min(self.x[:, 0] - self.radius)
+ min_y = numpy.min(self.x[:, 1] - self.radius)
+ min_z = numpy.min(self.x[:, 2] - self.radius)
+ max_x = numpy.max(self.x[:, 0] + self.radius)
+ max_y = numpy.max(self.x[:, 1] + self.radius)
+ max_z = numpy.max(self.x[:, 2] + self.radius)
+ V_total = (max_x - min_x)*(max_y - min_y)*(max_z - min_z)
else:
if self.nw == 0:
- V_total=self.L[0] * self.L[1] * self.L[2]
+ V_total = self.L[0] * self.L[1] * self.L[2]
elif self.nw == 1:
- V_total=self.L[0] * self.L[1] * self.w_x[0]
+ V_total = self.L[0] * self.L[1] * self.w_x[0]
if V_total <= 0.0:
raise Exception("Could not determine total volume")
-
# Find the volume of solids
- V_solid=numpy.sum(V_sphere(self.radius))
+ V_solid = numpy.sum(V_sphere(self.radius))
return (V_total - V_solid) / V_total
def porosity(self, slices=10, verbose=False):
t@@ -4453,24 +4369,23 @@ class sim:
self.writebin(verbose=False)
# Run porosity program on binary
- pipe=subprocess.Popen(\
- ["../porosity",\
- "-s","{}".format(slices),\
- "../input/" + self.sid + ".bin"],\
- stdout=subprocess.PIPE)
- output, err=pipe.communicate()
+ pipe = subprocess.Popen(["../porosity",\
+ "-s", "{}".format(slices),
+ "../input/" + self.sid + ".bin"],
+ stdout=subprocess.PIPE)
+ output, err = pipe.communicate()
if err:
print(err)
raise Exception("Could not run external 'porosity' program")
# read one line of output at a time
- s2=output.split(b'\n')
- depth=[]
- porosity=[]
+ s2 = output.split(b'\n')
+ depth = []
+ porosity = []
for row in s2:
if row != '\n' or row != '' or row != ' ': # skip blank lines
- s3=row.split(b'\t')
+ s3 = row.split(b'\t')
if s3 != '' and len(s3) == 2: # make sure line has two vals
depth.append(float(s3[0]))
porosity.append(float(s3[1]))
t@@ -4505,34 +4420,34 @@ class sim:
self.writebin(verbose=False)
- quiet=""
- stdout=""
- dryarg=""
- fluidarg=""
- devicearg=""
- valgrindbin=""
- cudamemchk=""
- binary="sphere"
+ quiet = ""
+ stdout = ""
+ dryarg = ""
+ fluidarg = ""
+ devicearg = ""
+ valgrindbin = ""
+ cudamemchk = ""
+ binary = "sphere"
if not verbose:
- quiet="-q "
+ quiet = "-q "
if hideinputfile:
- stdout=" > /dev/null"
+ stdout = " > /dev/null"
if dry:
- dryarg="--dry "
+ dryarg = "--dry "
if valgrind:
- valgrindbin="valgrind -q --track-origins=yes "
+ valgrindbin = "valgrind -q --track-origins=yes "
if cudamemcheck:
- cudamemchk="cuda-memcheck --leak-check full "
+ cudamemchk = "cuda-memcheck --leak-check full "
if self.fluid:
- fluidarg="--fluid "
+ fluidarg = "--fluid "
if device != -1:
- devicearg="-d " + str(device) + " "
+ devicearg = "-d " + str(device) + " "
- cmd="cd ..; " + valgrindbin + cudamemchk + "./" + binary + " " \
+ cmd = "cd ..; " + valgrindbin + cudamemchk + "./" + binary + " " \
+ quiet + dryarg + fluidarg + devicearg + \
"input/" + self.sid + ".bin " + stdout
#print(cmd)
- status=subprocess.call(cmd, shell=True)
+ status = subprocess.call(cmd, shell=True)
if status != 0:
print("Warning: the sphere run returned with status " + str(status))
t@@ -4582,10 +4497,10 @@ class sim:
'''
- filename=self.sid + ".sh"
- fh=None
- try :
- fh=open(filename, "w")
+ filename = self.sid + ".sh"
+ fh = None
+ try:
+ fh = open(filename, "w")
fh.write('#!/bin/sh\n')
fh.write('#PBS -N ' + self.sid + '\n')
t@@ -4597,7 +4512,7 @@ class sim:
fh.write('CUDAPATH=' + cudapath + '\n')
fh.write('export PATH=$CUDAPATH/bin:$PATH\n')
fh.write('export LD_LIBRARY_PATH=$CUDAPATH/lib64'
- + ':$CUDAPATH/lib:$LD_LIBRARY_PATH\n')
+ + ':$CUDAPATH/lib:$LD_LIBRARY_PATH\n')
fh.write('echo "`whoami`@`hostname`"\n')
fh.write('echo "Start at `date`"\n')
fh.write('ORIGDIR=' + spheredir + '\n')
t@@ -4614,7 +4529,7 @@ class sim:
fh.write('cp $WORKDIR/output/* $ORIGDIR/output/\n')
fh.write('echo "End at `date`"\n')
- finally :
+ finally:
if fh is not None:
fh.close()
t@@ -4645,21 +4560,21 @@ class sim:
print("Rendering {} images with the raytracer".format(self.sid))
- quiet=""
+ quiet = ""
if not verbose:
- quiet="-q"
+ quiet = "-q"
# Render images using sphere raytracer
if method == "normal":
subprocess.call("cd ..; for F in `ls output/" + self.sid
- + "*.bin`; do ./sphere " + quiet
- + " --render $F; done", shell=True)
- else :
+ + "*.bin`; do ./sphere " + quiet
+ + " --render $F; done", shell=True)
+ else:
subprocess.call("cd ..; for F in `ls output/" + self.sid
- + "*.bin`; do ./sphere " + quiet
- + " --method " + method + " {}".format(max_val)
- + " -l {}".format(lower_cutoff)
- + " --render $F; done", shell=True)
+ + "*.bin`; do ./sphere " + quiet
+ + " --method " + method + " {}".format(max_val)
+ + " -l {}".format(lower_cutoff)
+ + " --render $F; done", shell=True)
# Convert images to compressed format
if verbose:
t@@ -4668,7 +4583,7 @@ class sim:
def video(self, out_folder="./", video_format="mp4",
graphics_folder="../img_out/", graphics_format="png", fps=25,
- qscale=1, bitrate=1800, verbose=False):
+ verbose=False):
'''
Uses ffmpeg to combine images to animation. All images should be
rendered beforehand using :func:`render()`.
t@@ -4707,9 +4622,9 @@ class sim:
'''
# Displacement of the upper, fixed particles in the shear direction
- #xdisp=self.time_current[0] * self.shearVel()
- fixvel=numpy.nonzero(self.fixvel > 0.0)
- return numpy.max(self.xyzsum[fixvel,0])
+ #xdisp = self.time_current[0] * self.shearVel()
+ fixvel = numpy.nonzero(self.fixvel > 0.0)
+ return numpy.max(self.xyzsum[fixvel, 0])
def shearVelocity(self):
'''
t@@ -4723,9 +4638,9 @@ class sim:
See also: :func:`shearStrainRate()` and :func:`shearDisplacement()`
'''
# Displacement of the upper, fixed particles in the shear direction
- #xdisp=self.time_current[0] * self.shearVel()
- fixvel=numpy.nonzero(self.fixvel > 0.0)
- return numpy.max(self.vel[fixvel,0])
+ #xdisp = self.time_current[0] * self.shearVel()
+ fixvel = numpy.nonzero(self.fixvel > 0.0)
+ return numpy.max(self.vel[fixvel, 0])
def shearVel(self):
'''
t@@ -4746,10 +4661,10 @@ class sim:
'''
# Current height
- w_x0=self.w_x[0]
+ w_x0 = self.w_x[0]
# Displacement of the upper, fixed particles in the shear direction
- xdisp=self.shearDisplacement()
+ xdisp = self.shearDisplacement()
# Return shear strain
return xdisp/w_x0
t@@ -4766,8 +4681,8 @@ class sim:
#return self.shearStrain()/self.time_current[1]
# Current height
- w_x0=self.w_x[0]
- v=self.shearVelocity()
+ w_x0 = self.w_x[0]
+ v = self.shearVelocity()
# Return shear strain rate
return v/w_x0
t@@ -4796,11 +4711,12 @@ class sim:
'''
self.writebin(verbose=False)
subprocess.call('cd .. && ./sphere --contacts input/' + self.sid
- + '.bin > output/' + self.sid + '.contacts.txt', shell=True)
- contactdata=numpy.loadtxt('../output/' + self.sid + '.contacts.txt')
- self.pairs=numpy.array((contactdata[:,0], contactdata[:,1]),
- dtype=numpy.int32)
- self.overlaps=numpy.array(contactdata[:,2])
+ + '.bin > output/' + self.sid + '.contacts.txt',
+ shell=True)
+ contactdata = numpy.loadtxt('../output/' + self.sid + '.contacts.txt')
+ self.pairs = numpy.array((contactdata[:, 0], contactdata[:, 1]),
+ dtype=numpy.int32)
+ self.overlaps = numpy.array(contactdata[:, 2])
def findCoordinationNumber(self):
'''
t@@ -4808,10 +4724,10 @@ class sim:
particle). Requires a previous call to :func:`findOverlaps()`. Values
are stored in ``self.coordinationnumber``.
'''
- self.coordinationnumber=numpy.zeros(self.np, 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
+ self.coordinationnumber[self.pairs[0, i]] += 1
+ self.coordinationnumber[self.pairs[1, i]] += 1
def findMeanCoordinationNumber(self):
'''
t@@ -4834,7 +4750,7 @@ class sim:
See also: :func:`findOverlaps()` and :func:`findContactStresses()`
'''
self.findOverlaps()
- self.f_n_magn=self.k_n * numpy.abs(self.overlaps)
+ self.f_n_magn = self.k_n * numpy.abs(self.overlaps)
def contactSurfaceArea(self, i, j, overlap):
'''
t@@ -4849,13 +4765,12 @@ class sim:
:returns: Contact area [m*m]
:return type: float or array of floats
'''
- r_i=self.radius[i]
- r_j=self.radius[j]
- d=r_i + r_j + overlap
- contact_radius=1./(2.*d)*(
- (-d + r_i - r_j)*(-d - r_i + r_j)*
- (-d + r_i + r_j)*( d + r_i + r_j)
- )**0.5
+ r_i = self.radius[i]
+ r_j = self.radius[j]
+ d = r_i + r_j + overlap
+ contact_radius = 1./(2.*d)*((-d + r_i - r_j)*(-d - r_i + r_j)*
+ (-d + r_i + r_j)*(d + r_i + r_j)
+ )**0.5
return numpy.pi*contact_radius**2.
def contactParticleArea(self, i, j):
t@@ -4871,7 +4786,7 @@ class sim:
:returns: Contact area [m*m]
:return type: float or array of floats
'''
- r_bar=(self.radius[i] + self.radius[j])*0.5
+ r_bar = (self.radius[i] + self.radius[j])*0.5
return numpy.pi*r_bar**2.
def findAllContactSurfaceAreas(self):
t@@ -4883,8 +4798,8 @@ class sim:
:returns: Array of contact surface areas
:return type: array of floats
'''
- return self.contactSurfaceArea(self.pairs[0,:], self.pairs[1,:],
- self.overlaps)
+ return self.contactSurfaceArea(self.pairs[0, :], self.pairs[1, :],
+ self.overlaps)
def findAllAverageParticlePairAreas(self):
'''
t@@ -4895,7 +4810,7 @@ class sim:
:returns: Array of contact surface areas
:return type: array of floats
'''
- return self.contactParticleArea(self.pairs[0,:], self.pairs[1,:])
+ return self.contactParticleArea(self.pairs[0, :], self.pairs[1, :])
def findContactStresses(self, area='average'):
'''
t@@ -4913,13 +4828,13 @@ class sim:
'''
self.findNormalForces()
if area == 'average':
- areas=self.findAllAverageParticlePairAreas()
+ areas = self.findAllAverageParticlePairAreas()
elif area == 'contact':
- areas=self.findAllContactSurfaceAreas()
+ areas = self.findAllContactSurfaceAreas()
else:
raise Exception('Contact area type "' + area + '" not understood')
- self.sigma_contacts=self.f_n_magn/areas
+ self.sigma_contacts = self.f_n_magn/areas
def findLoadedContacts(self, threshold):
'''
t@@ -4956,13 +4871,14 @@ class sim:
self.writebin(verbose=False)
- nd=''
+ nd = ''
if disp == '2d':
- nd='-2d '
+ nd = '-2d '
- subprocess.call("cd .. && ./forcechains " + nd + "-f " + outformat \
- + " -lc " + str(lc) + " -uc " + str(uc) + " input/" + self.sid \
- + ".bin > python/tmp.gp", shell=True)
+ subprocess.call("cd .. && ./forcechains " + nd + "-f " + outformat
+ + " -lc " + str(lc) + " -uc " + str(uc)
+ + " input/" + self.sid + ".bin > python/tmp.gp",
+ shell=True)
subprocess.call("gnuplot tmp.gp && rm tmp.gp", shell=True)
t@@ -4984,42 +4900,41 @@ class sim:
+ ".bin > python/fc-tmp.txt", shell=True)
# data will have the shape (numcontacts, 7)
- data=numpy.loadtxt("fc-tmp.txt", skiprows=1)
+ data = numpy.loadtxt("fc-tmp.txt", skiprows=1)
# find the max. value of the normal force
- f_n_max=numpy.amax(data[:,6])
+ f_n_max = numpy.amax(data[:, 6])
# specify the lower limit of force chains to do statistics on
- f_n_lim=lower_limit * f_n_max * 0.6
+ f_n_lim = lower_limit * f_n_max * 0.6
# find the indexes of these contacts
- I=numpy.nonzero(data[:,6] > f_n_lim)
+ I = numpy.nonzero(data[:, 6] > f_n_lim)
# loop through these contacts and find the strike and dip of the
# contacts
- strikelist=[] # strike direction of the normal vector, [0:360[
- diplist=[] # dip of the normal vector, [0:90]
+ strikelist = [] # strike direction of the normal vector, [0:360[
+ diplist = [] # dip of the normal vector, [0:90]
for i in I[0]:
- x1=data[i,0]
- y1=data[i,1]
- z1=data[i,2]
- x2=data[i,3]
- y2=data[i,4]
- z2=data[i,5]
+ x1 = data[i, 0]
+ y1 = data[i, 1]
+ z1 = data[i, 2]
+ x2 = data[i, 3]
+ y2 = data[i, 4]
+ z2 = data[i, 5]
if z1 < z2:
- xlower=x1; ylower=y1; zlower=z1
- xupper=x2; yupper=y2; zupper=z2
- else :
- xlower=x2; ylower=y2; zlower=z2
- xupper=x1; yupper=y1; zupper=z1
+ xlower = x1; ylower = y1; zlower = z1
+ xupper = x2; yupper = y2; zupper = z2
+ else:
+ xlower = x2; ylower = y2; zlower = z2
+ xupper = x1; yupper = y1; zupper = z1
# Vector pointing downwards
- dx=xlower - xupper
- dy=ylower - yupper
- dz=zlower - zupper
- dhoriz=numpy.sqrt(dx**2 + dy**2)
+ dx = xlower - xupper
+ dy = ylower - yupper
+ dhoriz = numpy.sqrt(dx**2 + dy**2)
# Find dip angle
diplist.append(math.degrees(math.atan((zupper - zlower)/dhoriz)))
t@@ -5027,17 +4942,17 @@ class sim:
# Find strike angle
if ylower >= yupper: # in first two quadrants
strikelist.append(math.acos(dx/dhoriz))
- else :
+ else:
strikelist.append(2.0*numpy.pi - math.acos(dx/dhoriz))
- plt.figure(figsize=[4,4])
- ax=plt.subplot(111, polar=True)
+ plt.figure(figsize=[4, 4])
+ ax = plt.subplot(111, polar=True)
ax.scatter(strikelist, diplist, c='k', marker='+')
ax.set_rmax(90)
ax.set_rticks([])
plt.savefig('fc-' + self.sid + '-rose.' + graphics_format,\
- transparent=True)
+ transparent=True)
subprocess.call('rm fc-tmp.txt', shell=True)
t@@ -5055,32 +4970,31 @@ class sim:
return
# loop through these contacts and find the strike and dip of the
# contacts
- strikelist=[] # strike direction of the normal vector, [0:360[
- diplist=[] # dip of the normal vector, [0:90]
+ strikelist = [] # strike direction of the normal vector, [0:360[
+ diplist = [] # dip of the normal vector, [0:90]
for n in numpy.arange(self.nb0):
- i=self.bonds[n,0]
- j=self.bonds[n,1]
+ i = self.bonds[n, 0]
+ j = self.bonds[n, 1]
- x1=self.x[i,0]
- y1=self.x[i,1]
- z1=self.x[i,2]
- x2=self.x[j,0]
- y2=self.x[j,1]
- z2=self.x[j,2]
+ x1 = self.x[i, 0]
+ y1 = self.x[i, 1]
+ z1 = self.x[i, 2]
+ x2 = self.x[j, 0]
+ y2 = self.x[j, 1]
+ z2 = self.x[j, 2]
if z1 < z2:
- xlower=x1; ylower=y1; zlower=z1
- xupper=x2; yupper=y2; zupper=z2
- else :
- xlower=x2; ylower=y2; zlower=z2
- xupper=x1; yupper=y1; zupper=z1
+ xlower = x1; ylower = y1; zlower = z1
+ xupper = x2; yupper = y2; zupper = z2
+ else:
+ xlower = x2; ylower = y2; zlower = z2
+ xupper = x1; yupper = y1; zupper = z1
# Vector pointing downwards
- dx=xlower - xupper
- dy=ylower - yupper
- dz=zlower - zupper
- dhoriz=numpy.sqrt(dx**2 + dy**2)
+ dx = xlower - xupper
+ dy = ylower - yupper
+ dhoriz = numpy.sqrt(dx**2 + dy**2)
# Find dip angle
diplist.append(math.degrees(math.atan((zupper - zlower)/dhoriz)))
t@@ -5088,16 +5002,16 @@ class sim:
# Find strike angle
if ylower >= yupper: # in first two quadrants
strikelist.append(math.acos(dx/dhoriz))
- else :
+ else:
strikelist.append(2.0*numpy.pi - math.acos(dx/dhoriz))
- plt.figure(figsize=[4,4])
- ax=plt.subplot(111, polar=True)
+ plt.figure(figsize=[4, 4])
+ ax = plt.subplot(111, polar=True)
ax.scatter(strikelist, diplist, c='k', marker='+')
ax.set_rmax(90)
ax.set_rticks([])
plt.savefig('bonds-' + self.sid + '-rose.' + graphics_format,\
- transparent=True)
+ transparent=True)
def status(self):
'''
t@@ -5118,7 +5032,7 @@ class sim:
:returns: The particle momentum [N*s]
:return type: numpy.array
'''
- return self.rho*V_sphere(self.radius[idx])*self.vel[idx,:]
+ return self.rho*V_sphere(self.radius[idx])*self.vel[idx, :]
def totalMomentum(self):
'''
t@@ -5127,7 +5041,7 @@ class sim:
:returns: The sum of particle momentums (m*v) [N*s]
:return type: numpy.array
'''
- m_sum=numpy.zeros(3)
+ m_sum = numpy.zeros(3)
for i in range(self.np):
m_sum += self.momentum(i)
return m_sum
t@@ -5146,40 +5060,40 @@ class sim:
return
# Bin data and error bars for alternative visualization
- h_total=numpy.max(self.x[:,2]) - numpy.min(self.x[:,2])
- h_slice=h_total / zslices
+ h_total = numpy.max(self.x[:, 2]) - numpy.min(self.x[:, 2])
+ h_slice = h_total / zslices
- zpos=numpy.zeros(zslices)
- xdisp=numpy.zeros(zslices)
- err=numpy.zeros(zslices)
+ zpos = numpy.zeros(zslices)
+ xdisp = numpy.zeros(zslices)
+ err = numpy.zeros(zslices)
for iz in range(zslices):
# Find upper and lower boundaries of bin
- zlower=iz * h_slice
- zupper=zlower + h_slice
+ zlower = iz * h_slice
+ zupper = zlower + h_slice
# Save depth
- zpos[iz]=zlower + 0.5*h_slice
+ zpos[iz] = zlower + 0.5*h_slice
# Find particle indexes within that slice
- I=numpy.nonzero((self.x[:,2] > zlower) & (self.x[:,2] < zupper))
+ I = numpy.nonzero((self.x[:, 2] > zlower) & (self.x[:, 2] < zupper))
# Save mean x displacement
- xdisp[iz]=numpy.mean(self.xyzsum[I,0])
+ xdisp[iz] = numpy.mean(self.xyzsum[I, 0])
# Save x displacement standard deviation
- err[iz]=numpy.std(self.xyzsum[I,0])
+ err[iz] = numpy.std(self.xyzsum[I, 0])
plt.figure(figsize=[4, 4])
- ax=plt.subplot(111)
- ax.scatter(self.xyzsum[:,0], self.x[:,2], c='gray', marker='+')
+ ax = plt.subplot(111)
+ ax.scatter(self.xyzsum[:, 0], self.x[:, 2], c='gray', marker='+')
ax.errorbar(xdisp, zpos, xerr=err,
c='black', linestyle='-', linewidth=1.4)
ax.set_xlabel("Horizontal particle displacement, [m]")
ax.set_ylabel("Vertical position, [m]")
plt.savefig(self.sid + '-sheardisp.' + graphics_format,
- transparent=True)
+ transparent=True)
def porosities(self, graphics_format='pdf', zslices=16):
'''
t@@ -5196,16 +5110,15 @@ class sim:
print('Error: matplotlib module not found, cannot sheardisp.')
return
- porosity, depth=self.porosity(zslices)
+ porosity, depth = self.porosity(zslices)
plt.figure(figsize=[4, 4])
- ax=plt.subplot(111)
- ax.plot(porosity, depth,
- c='black', linestyle='-', linewidth=1.4)
+ ax = plt.subplot(111)
+ ax.plot(porosity, depth, c='black', linestyle='-', linewidth=1.4)
ax.set_xlabel('Horizontally averaged porosity, [-]')
ax.set_ylabel('Vertical position, [m]')
plt.savefig(self.sid + '-porositiy.' + graphics_format,
- transparent=True)
+ transparent=True)
def thinsection_x1x3(self, x2='center', graphics_format='png', cbmax=None,
arrowscale=0.01, velarrowscale=1.0, slipscale=1.0,
t@@ -5246,32 +5159,32 @@ class sim:
print('Error: matplotlib module not found (thinsection_x1x3).')
return
- if x2 == 'center' :
- x2=(self.L[1] - self.origo[1]) / 2.0
+ if x2 == 'center':
+ x2 = (self.L[1] - self.origo[1]) / 2.0
# Initialize plot circle positionsr, radii and pressures
- ilist=[]
- xlist=[]
- ylist=[]
- rlist=[]
- plist=[]
- pmax=0.0
- rmax=0.0
- axlist=[]
- aylist=[]
- daxlist=[]
- daylist=[]
- dvxlist=[]
- dvylist=[]
- # Black circle at periphery of particles with angvel[:,1] > 0.0
- cxlist=[]
- cylist=[]
- crlist=[]
+ ilist = []
+ xlist = []
+ ylist = []
+ rlist = []
+ plist = []
+ pmax = 0.0
+ rmax = 0.0
+ axlist = []
+ aylist = []
+ daxlist = []
+ daylist = []
+ dvxlist = []
+ dvylist = []
+ # Black circle at periphery of particles with angvel[:, 1] > 0.0
+ cxlist = []
+ cylist = []
+ crlist = []
# Loop over all particles, find intersections
for i in range(self.np):
- delta=abs(self.x[i,1] - x2) # distance between centre and plane
+ delta = abs(self.x[i, 1] - x2) # distance between centre and plane
if delta < self.radius[i]: # if the sphere intersects the plane
t@@ -5279,60 +5192,60 @@ class sim:
ilist.append(i)
# Store position on plane
- xlist.append(self.x[i,0])
- ylist.append(self.x[i,2])
+ xlist.append(self.x[i, 0])
+ ylist.append(self.x[i, 2])
# Store radius of intersection
- r_circ=math.sqrt(self.radius[i]**2 - delta**2)
+ r_circ = math.sqrt(self.radius[i]**2 - delta**2)
if r_circ > rmax:
- rmax=r_circ
+ rmax = r_circ
rlist.append(r_circ)
# Store pos. and radius if it is spinning around pos. y
- if self.angvel[i,1] > 0.0:
- cxlist.append(self.x[i,0])
- cylist.append(self.x[i,2])
+ if self.angvel[i, 1] > 0.0:
+ cxlist.append(self.x[i, 0])
+ cylist.append(self.x[i, 2])
crlist.append(r_circ)
# Store pressure
- pval=self.p[i]
+ pval = self.p[i]
if cbmax != None:
if pval > cbmax:
- pval=cbmax
+ pval = cbmax
plist.append(pval)
# Store rotational velocity data for arrows
# Save two arrows per particle
- axlist.append(self.x[i,0]) # x starting point of arrow
- axlist.append(self.x[i,0]) # x starting point of arrow
+ axlist.append(self.x[i, 0]) # x starting point of arrow
+ axlist.append(self.x[i, 0]) # x starting point of arrow
# y starting point of arrow
- aylist.append(self.x[i,2] + r_circ*0.5)
+ aylist.append(self.x[i, 2] + r_circ*0.5)
# y starting point of arrow
- aylist.append(self.x[i,2] - r_circ*0.5)
+ aylist.append(self.x[i, 2] - r_circ*0.5)
# delta x for arrow end point
- daxlist.append(self.angvel[i,1]*arrowscale)
+ daxlist.append(self.angvel[i, 1]*arrowscale)
# delta x for arrow end point
- daxlist.append(-self.angvel[i,1]*arrowscale)
+ daxlist.append(-self.angvel[i, 1]*arrowscale)
daylist.append(0.0) # delta y for arrow end point
daylist.append(0.0) # delta y for arrow end point
# Store linear velocity data
# delta x for arrow end point
- dvxlist.append(self.vel[i,0]*velarrowscale)
+ dvxlist.append(self.vel[i, 0]*velarrowscale)
# delta y for arrow end point
- dvylist.append(self.vel[i,2]*velarrowscale)
+ dvylist.append(self.vel[i, 2]*velarrowscale)
if r_circ > self.radius[i]:
raise Exception("Error, circle radius is larger than the "
- + "particle radius")
+ "particle radius")
if self.p[i] > pmax:
- pmax=self.p[i]
+ pmax = self.p[i]
if verbose:
print("Max. pressure of intersecting spheres: " + str(pmax) + " Pa")
t@@ -5340,28 +5253,28 @@ class sim:
print("Value limited to: " + str(cbmax) + " Pa")
# Save circle data
- filename='../gnuplot/data/' + self.sid + '-ts-x1x3.txt'
- fh=None
- try :
- fh=open(filename, 'w')
+ filename = '../gnuplot/data/' + self.sid + '-ts-x1x3.txt'
+ fh = None
+ try:
+ fh = open(filename, 'w')
for (x, y, r, p) in zip(xlist, ylist, rlist, plist):
fh.write("{}\t{}\t{}\t{}\n".format(x, y, r, p))
- finally :
+ finally:
if fh is not None:
fh.close()
# Save circle data for articles spinning with pos. y
- filename='../gnuplot/data/' + self.sid + '-ts-x1x3-circ.txt'
- fh=None
- try :
- fh=open(filename, 'w')
+ filename = '../gnuplot/data/' + self.sid + '-ts-x1x3-circ.txt'
+ fh = None
+ try:
+ fh = open(filename, 'w')
for (x, y, r) in zip(cxlist, cylist, crlist):
fh.write("{}\t{}\t{}\n".format(x, y, r))
- finally :
+ finally:
if fh is not None:
fh.close()
t@@ -5369,41 +5282,41 @@ class sim:
# radius
# Output format: x, y, deltax, deltay
# gnuplot> plot '-' using 1:2:3:4 with vectors head filled lt 2
- filename='../gnuplot/data/' + self.sid + '-ts-x1x3-arrows.txt'
- fh=None
- try :
- fh=open(filename, 'w')
+ filename = '../gnuplot/data/' + self.sid + '-ts-x1x3-arrows.txt'
+ fh = None
+ try:
+ fh = open(filename, 'w')
for (ax, ay, dax, day) in zip(axlist, aylist, daxlist, daylist):
fh.write("{}\t{}\t{}\t{}\n".format(ax, ay, dax, day))
- finally :
+ finally:
if fh is not None:
fh.close()
# Save linear velocity data
# Output format: x, y, deltax, deltay
# gnuplot> plot '-' using 1:2:3:4 with vectors head filled lt 2
- filename='../gnuplot/data/' + self.sid + '-ts-x1x3-velarrows.txt'
- fh=None
- try :
- fh=open(filename, 'w')
+ filename = '../gnuplot/data/' + self.sid + '-ts-x1x3-velarrows.txt'
+ fh = None
+ try:
+ fh = open(filename, 'w')
for (x, y, dvx, dvy) in zip(xlist, ylist, dvxlist, dvylist):
fh.write("{}\t{}\t{}\t{}\n".format(x, y, dvx, dvy))
- finally :
+ finally:
if fh is not None:
fh.close()
# Check whether there are slips between the particles intersecting the
# plane
- sxlist=[]
- sylist=[]
- dsxlist=[]
- dsylist=[]
- anglelist=[] # angle of the slip vector
- slipvellist=[] # velocity of the slip
+ sxlist = []
+ sylist = []
+ dsxlist = []
+ dsylist = []
+ anglelist = [] # angle of the slip vector
+ slipvellist = [] # velocity of the slip
for i in ilist:
# Loop through other particles, and check whether they are in
t@@ -5413,43 +5326,43 @@ class sim:
if i != j:
# positions
- x_i=self.x[i,:]
- x_j=self.x[j,:]
+ x_i = self.x[i, :]
+ x_j = self.x[j, :]
# radii
- r_i=self.radius[i]
- r_j=self.radius[j]
+ r_i = self.radius[i]
+ r_j = self.radius[j]
# Inter-particle vector
- x_ij=x_i - x_j
- x_ij_length=numpy.sqrt(x_ij.dot(x_ij))
+ x_ij = x_i - x_j
+ x_ij_length = numpy.sqrt(x_ij.dot(x_ij))
# Check for overlap
if x_ij_length - (r_i + r_j) < 0.0:
# contact plane normal vector
- n_ij=x_ij / x_ij_length
+ n_ij = x_ij / x_ij_length
- vel_i=self.vel[i,:]
- vel_j=self.vel[j,:]
- angvel_i=self.angvel[i,:]
- angvel_j=self.angvel[j,:]
+ vel_i = self.vel[i, :]
+ vel_j = self.vel[j, :]
+ angvel_i = self.angvel[i, :]
+ angvel_j = self.angvel[j, :]
# Determine the tangential contact surface velocity in
# the x,z plane
- dot_delta=(vel_i - vel_j) \
+ dot_delta = (vel_i - vel_j) \
+ r_i * numpy.cross(n_ij, angvel_i) \
+ r_j * numpy.cross(n_ij, angvel_j)
# Subtract normal component to get tangential velocity
- dot_delta_n=n_ij * numpy.dot(dot_delta, n_ij)
- dot_delta_t=dot_delta - dot_delta_n
+ dot_delta_n = n_ij * numpy.dot(dot_delta, n_ij)
+ dot_delta_t = dot_delta - dot_delta_n
# Save slip velocity data for gnuplot
if dot_delta_t[0] != 0.0 or dot_delta_t[2] != 0.0:
# Center position of the contact
- cpos=x_i - x_ij * 0.5
+ cpos = x_i - x_ij * 0.5
sxlist.append(cpos[0])
sylist.append(cpos[2])
t@@ -5464,51 +5377,51 @@ class sim:
# Write slip lines to text file
- filename='../gnuplot/data/' + self.sid + '-ts-x1x3-slips.txt'
- fh=None
- try :
- fh=open(filename, 'w')
+ filename = '../gnuplot/data/' + self.sid + '-ts-x1x3-slips.txt'
+ fh = None
+ try:
+ fh = open(filename, 'w')
for (sx, sy, dsx, dsy) in zip(sxlist, sylist, dsxlist, dsylist):
fh.write("{}\t{}\t{}\t{}\n".format(sx, sy, dsx, dsy))
- finally :
+ finally:
if fh is not None:
fh.close()
# Plot thinsection with gnuplot script
- gamma=self.shearstrain()
+ gamma = self.shearStrain()
subprocess.call('''cd ../gnuplot/scripts && gnuplot -e "sid='{}'; ''' \
+ '''gamma='{:.4}'; xmin='{}'; xmax='{}'; ymin='{}'; ''' \
+ '''ymax='{}'" plotts.gp'''.format(\
- self.sid, self.shearstrain(), self.origo[0], self.L[0], \
+ self.sid, self.shearStrain(), self.origo[0], self.L[0], \
self.origo[2], self.L[2]), shell=True)
# Find all particles who have a slip velocity higher than slipvel
- slipvellimit=0.01
- slipvels=numpy.nonzero(numpy.array(slipvellist) > slipvellimit)
+ slipvellimit = 0.01
+ slipvels = numpy.nonzero(numpy.array(slipvellist) > slipvellimit)
# Bin slip angle data for histogram
- binno=36/2
- hist_ang, bins_ang=numpy.histogram(numpy.array(anglelist)[slipvels],\
+ binno = 36/2
+ hist_ang, bins_ang = numpy.histogram(numpy.array(anglelist)[slipvels],\
bins=binno, density=False)
- center_ang=(bins_ang[:-1] + bins_ang[1:]) / 2.0
+ center_ang = (bins_ang[:-1] + bins_ang[1:]) / 2.0
- center_ang_mirr=numpy.concatenate((center_ang, center_ang + math.pi))
- hist_ang_mirr=numpy.tile(hist_ang, 2)
+ center_ang_mirr = numpy.concatenate((center_ang, center_ang + math.pi))
+ hist_ang_mirr = numpy.tile(hist_ang, 2)
# Write slip angles to text file
#numpy.savetxt(self.sid + '-ts-x1x3-slipangles.txt', zip(center_ang,\
#hist_ang), fmt="%f\t%f")
- fig=plt.figure()
- ax=fig.add_subplot(111, polar=True)
+ fig = plt.figure()
+ ax = fig.add_subplot(111, polar=True)
ax.bar(center_ang_mirr, hist_ang_mirr, width=30.0/180.0)
fig.savefig('../img_out/' + self.sid + '-ts-x1x3-slipangles.' +
- graphics_format)
+ graphics_format)
fig.clf()
- def plotContacts(self, graphics_format='png', figsize=[4,4], title=None,
+ def plotContacts(self, graphics_format='png', figsize=[4, 4], title=None,
lower_limit=0.0, upper_limit=1.0, alpha=1.0,
return_data=False, outfolder='.',
f_min=None, f_max=None, histogram=True,
t@@ -5536,86 +5449,83 @@ class sim:
+ ".bin > python/contacts-tmp.txt", shell=True)
# data will have the shape (numcontacts, 7)
- data=numpy.loadtxt('contacts-tmp.txt', skiprows=1)
+ data = numpy.loadtxt('contacts-tmp.txt', skiprows=1)
# find the max. value of the normal force
- f_n_max=numpy.amax(data[:,6])
+ f_n_max = numpy.amax(data[:, 6])
# specify the lower limit of force chains to do statistics on
- f_n_lim=lower_limit * f_n_max
+ f_n_lim = lower_limit * f_n_max
if f_min:
- f_n_lim=f_min
+ f_n_lim = f_min
if f_max:
- f_n_max=f_max
+ f_n_max = f_max
# find the indexes of these contacts
- I=numpy.nonzero(data[:,6] >= f_n_lim)
+ I = numpy.nonzero(data[:, 6] >= f_n_lim)
# loop through these contacts and find the strike and dip of the
# contacts
# strike direction of the normal vector, [0:360[
- strikelist=numpy.empty(len(I[0]))
- diplist=numpy.empty(len(I[0])) # dip of the normal vector, [0:90]
- forcemagnitude=data[I,6]
- j=0
+ strikelist = numpy.empty(len(I[0]))
+ diplist = numpy.empty(len(I[0])) # dip of the normal vector, [0:90]
+ forcemagnitude = data[I, 6]
+ j = 0
for i in I[0]:
- x1=data[i,0]
- y1=data[i,1]
- z1=data[i,2]
- x2=data[i,3]
- y2=data[i,4]
- z2=data[i,5]
+ x1 = data[i, 0]
+ y1 = data[i, 1]
+ z1 = data[i, 2]
+ x2 = data[i, 3]
+ y2 = data[i, 4]
+ z2 = data[i, 5]
if z1 < z2:
- xlower=x1; ylower=y1; zlower=z1
- xupper=x2; yupper=y2; zupper=z2
- else :
- xlower=x2; ylower=y2; zlower=z2
- xupper=x1; yupper=y1; zupper=z1
+ xlower = x1; ylower = y1; zlower = z1
+ xupper = x2; yupper = y2; zupper = z2
+ else:
+ xlower = x2; ylower = y2; zlower = z2
+ xupper = x1; yupper = y1; zupper = z1
# Vector pointing downwards
- dx=xlower - xupper
- dy=ylower - yupper
- dz=zlower - zupper
- dhoriz=numpy.sqrt(dx**2 + dy**2)
+ dx = xlower - xupper
+ dy = ylower - yupper
+ dhoriz = numpy.sqrt(dx**2 + dy**2)
# Find dip angle
- diplist[j]=numpy.degrees(numpy.arctan((zupper - zlower)/dhoriz))
+ diplist[j] = numpy.degrees(numpy.arctan((zupper - zlower)/dhoriz))
# Find strike angle
if ylower >= yupper: # in first two quadrants
- strikelist[j]=numpy.arccos(dx/dhoriz)
- else :
- strikelist[j]=2.0*numpy.pi - numpy.arccos(dx/dhoriz)
+ strikelist[j] = numpy.arccos(dx/dhoriz)
+ else:
+ strikelist[j] = 2.0*numpy.pi - numpy.arccos(dx/dhoriz)
j += 1
- fig=plt.figure(figsize=figsize)
- ax=plt.subplot(111, polar=True)
- cs=ax.scatter(strikelist, 90. - diplist, marker='o',
- c=forcemagnitude,
- s=forcemagnitude/f_n_max*40.,
- alpha=alpha,
- edgecolors='none',
- vmin=f_n_max*lower_limit,
- vmax=f_n_max*upper_limit,
- cmap=matplotlib.cm.get_cmap('afmhot_r'))
+ fig = plt.figure(figsize=figsize)
+ ax = plt.subplot(111, polar=True)
+ cs = ax.scatter(strikelist, 90. - diplist, marker='o',
+ c=forcemagnitude,
+ s=forcemagnitude/f_n_max*40.,
+ alpha=alpha,
+ edgecolors='none',
+ vmin=f_n_max*lower_limit,
+ vmax=f_n_max*upper_limit,
+ cmap=matplotlib.cm.get_cmap('afmhot_r'))
plt.colorbar(cs, extend='max')
# plot defined max compressive stress from tau/N ratio
ax.scatter(0., # prescribed stress
- numpy.degrees(numpy.arctan(
- self.shearStress('defined')/
- self.currentNormalStress('defined'))),
- marker='o', c='none', edgecolor='blue', s=300)
+ numpy.degrees(numpy.arctan(self.shearStress('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'))),
- marker='+', color='blue', s=300)
+ numpy.degrees(numpy.arctan(self.shearStress('effective')/
+ self.currentNormalStress('effective'))),
+ marker='+', color='blue', s=300)
ax.set_rmax(90)
ax.set_rticks([])
t@@ -5635,8 +5545,8 @@ class sim:
fig.clf()
if histogram:
- #hist, bins=numpy.histogram(datadata[:,6], bins=10)
- n, bins, patches=plt.hist(data[:,6], alpha=0.75, facecolor='gray')
+ #hist, bins = numpy.histogram(datadata[:, 6], bins=10)
+ _, _, _ = plt.hist(data[:, 6], alpha=0.75, facecolor='gray')
#plt.xlabel('$\\boldsymbol{f}_\text{n}$ [N]')
plt.yscale('log', nonposy='clip')
plt.xlabel('Contact load [N]')
t@@ -5649,10 +5559,10 @@ class sim:
plt.clf()
# angle: 0 when vertical, 90 when horizontal
- #hist, bins=numpy.histogram(datadata[:,6], bins=10)
- n, bins, patches=plt.hist(90. - diplist, bins=range(0, 100, 10),
+ #hist, bins = numpy.histogram(datadata[:, 6], bins=10)
+ _, _, _ = plt.hist(90. - diplist, bins=range(0, 100, 10),
alpha=0.75, facecolor='gray')
- theta_sigma1=numpy.degrees(numpy.arctan(
+ theta_sigma1 = numpy.degrees(numpy.arctan(
self.currentNormalStress('defined')/\
self.shearStress('defined')))
plt.axvline(90. - theta_sigma1, color='k', linestyle='dashed',
t@@ -5671,27 +5581,27 @@ class sim:
if forcechains:
- #color=matplotlib.cm.spectral(data[:,6]/f_n_max)
+ #color = matplotlib.cm.spectral(data[:, 6]/f_n_max)
for i in I[0]:
- x1=data[i,0]
- #y1=data[i,1]
- z1=data[i,2]
- x2=data[i,3]
- #y2=data[i,4]
- z2=data[i,5]
- f_n=data[i,6]
+ x1 = data[i, 0]
+ #y1 = data[i, 1]
+ z1 = data[i, 2]
+ x2 = data[i, 3]
+ #y2 = data[i, 4]
+ z2 = data[i, 5]
+ f_n = data[i, 6]
- lw_max=1.0
+ lw_max = 1.0
if f_n >= f_n_max:
- lw=lw_max
+ lw = lw_max
else:
- lw=(f_n - f_n_lim)/(f_n_max - f_n_lim)*lw_max
+ lw = (f_n - f_n_lim)/(f_n_max - f_n_lim)*lw_max
#print lw
- plt.plot([x1,x2], [z1,z2], '-k', linewidth=lw)
+ plt.plot([x1, x2], [z1, z2], '-k', linewidth=lw)
- axfc1=plt.gca()
+ axfc1 = plt.gca()
axfc1.spines['right'].set_visible(False)
axfc1.spines['left'].set_visible(False)
# Only show ticks on the left and bottom spines
t@@ -5738,18 +5648,18 @@ class sim:
return
if y == -1:
- y=self.num[1]/2
+ y = self.num[1]/2
- plt.figure(figsize=[8,8])
+ plt.figure(figsize=[8, 8])
plt.title('Fluid pressures')
- imgplt=plt.imshow(self.p_f[:,y,:].T, origin='lower')
+ imgplt = plt.imshow(self.p_f[:, y, :].T, origin='lower')
imgplt.set_interpolation('nearest')
#imgplt.set_interpolation('bicubic')
#imgplt.set_cmap('hot')
plt.xlabel('$x_1$')
plt.ylabel('$x_3$')
plt.colorbar()
- filename='p_f-' + self.sid + '-y' + str(y) + '.' + graphics_format
+ filename = 'p_f-' + self.sid + '-y' + str(y) + '.' + graphics_format
plt.savefig(filename, transparent=False)
if verbose:
print('saved to ' + filename)
t@@ -5778,18 +5688,18 @@ class sim:
return
if z == -1:
- z=self.num[2]/2
+ z = self.num[2]/2
- plt.figure(figsize=[8,8])
+ plt.figure(figsize=[8, 8])
plt.title('Fluid pressures')
- imgplt=plt.imshow(self.p_f[:,:,z].T, origin='lower')
+ imgplt = plt.imshow(self.p_f[:, :, z].T, origin='lower')
imgplt.set_interpolation('nearest')
#imgplt.set_interpolation('bicubic')
#imgplt.set_cmap('hot')
plt.xlabel('$x_1$')
plt.ylabel('$x_2$')
plt.colorbar()
- filename='p_f-' + self.sid + '-z' + str(z) + '.' + graphics_format
+ filename = 'p_f-' + self.sid + '-z' + str(z) + '.' + graphics_format
plt.savefig(filename, transparent=False)
if verbose:
print('saved to ' + filename)
t@@ -5818,13 +5728,13 @@ class sim:
return
if y == -1:
- y=self.num[1]/2
+ y = self.num[1]/2
plt.title('Fluid velocities')
- plt.figure(figsize=[8,8])
+ plt.figure(figsize=[8, 8])
plt.subplot(131)
- imgplt=plt.imshow(self.v_f[:,y,:,0].T, origin='lower')
+ imgplt = plt.imshow(self.v_f[:, y, :, 0].T, origin='lower')
imgplt.set_interpolation('nearest')
#imgplt.set_interpolation('bicubic')
#imgplt.set_cmap('hot')
t@@ -5834,7 +5744,7 @@ class sim:
plt.colorbar(orientation='horizontal')
plt.subplot(132)
- imgplt=plt.imshow(self.v_f[:,y,:,1].T, origin='lower')
+ imgplt = plt.imshow(self.v_f[:, y, :, 1].T, origin='lower')
imgplt.set_interpolation('nearest')
#imgplt.set_interpolation('bicubic')
#imgplt.set_cmap('hot')
t@@ -5844,7 +5754,7 @@ class sim:
plt.colorbar(orientation='horizontal')
plt.subplot(133)
- imgplt=plt.imshow(self.v_f[:,y,:,2].T, origin='lower')
+ imgplt = plt.imshow(self.v_f[:, y, :, 2].T, origin='lower')
imgplt.set_interpolation('nearest')
#imgplt.set_interpolation('bicubic')
#imgplt.set_cmap('hot')
t@@ -5853,7 +5763,7 @@ class sim:
plt.ylabel('$x_3$')
plt.colorbar(orientation='horizontal')
- filename='v_f-' + self.sid + '-y' + str(y) + '.' + graphics_format
+ filename = 'v_f-' + self.sid + '-y' + str(y) + '.' + graphics_format
plt.savefig(filename, transparent=False)
if verbose:
print('saved to ' + filename)
t@@ -5881,13 +5791,13 @@ class sim:
return
if z == -1:
- z=self.num[2]/2
+ z = self.num[2]/2
plt.title("Fluid velocities")
- plt.figure(figsize=[8,8])
+ plt.figure(figsize=[8, 8])
plt.subplot(131)
- imgplt=plt.imshow(self.v_f[:,:,z,0].T, origin='lower')
+ imgplt = plt.imshow(self.v_f[:, :, z, 0].T, origin='lower')
imgplt.set_interpolation('nearest')
#imgplt.set_interpolation('bicubic')
#imgplt.set_cmap('hot')
t@@ -5897,7 +5807,7 @@ class sim:
plt.colorbar(orientation='horizontal')
plt.subplot(132)
- imgplt=plt.imshow(self.v_f[:,:,z,1].T, origin='lower')
+ imgplt = plt.imshow(self.v_f[:, :, z, 1].T, origin='lower')
imgplt.set_interpolation('nearest')
#imgplt.set_interpolation('bicubic')
#imgplt.set_cmap('hot')
t@@ -5907,7 +5817,7 @@ class sim:
plt.colorbar(orientation='horizontal')
plt.subplot(133)
- imgplt=plt.imshow(self.v_f[:,:,z,2].T, origin='lower')
+ imgplt = plt.imshow(self.v_f[:, :, z, 2].T, origin='lower')
imgplt.set_interpolation('nearest')
#imgplt.set_interpolation('bicubic')
#imgplt.set_cmap('hot')
t@@ -5916,7 +5826,7 @@ class sim:
plt.ylabel('$x_2$')
plt.colorbar(orientation='horizontal')
- filename='v_f-' + self.sid + '-z' + str(z) + '.' + graphics_format
+ filename = 'v_f-' + self.sid + '-z' + str(z) + '.' + graphics_format
plt.savefig(filename, transparent=False)
if verbose:
print('saved to ' + filename)
t@@ -5941,56 +5851,58 @@ class sim:
print('Error: matplotlib module not found (plotFluidDiffAdvPresZ).')
return
- # The v_z values are read from self.v_f[0,0,:,2]
- dz=self.L[2]/self.num[2]
- rho=self.rho_f
+ # The v_z values are read from self.v_f[0, 0, :, 2]
+ dz = self.L[2]/self.num[2]
+ rho = self.rho_f
# Central difference gradients
- dvz_dz=(self.v_f[0,0,1:,2] - self.v_f[0,0,:-1,2])/(2.0*dz)
- dvzvz_dz=(self.v_f[0,0,1:,2]**2 - self.v_f[0,0,:-1,2]**2)/(2.0*dz)
+ dvz_dz = (self.v_f[0, 0, 1:, 2] - self.v_f[0, 0, :-1, 2])/(2.0*dz)
+ dvzvz_dz = (self.v_f[0, 0, 1:, 2]**2 - self.v_f[0, 0, :-1, 2]**2)\
+ /(2.0*dz)
# Diffusive contribution to velocity change
- dvz_diff=2.0*self.mu/rho*dvz_dz*self.time_dt
+ dvz_diff = 2.0*self.mu/rho*dvz_dz*self.time_dt
# Advective contribution to velocity change
- dvz_adv=dvzvz_dz*self.time_dt
+ dvz_adv = dvzvz_dz*self.time_dt
# Pressure gradient
- dp_dz=(self.p_f[0,0,1:] - self.p_f[0,0,:-1])/(2.0*dz)
+ dp_dz = (self.p_f[0, 0, 1:] - self.p_f[0, 0, :-1])/(2.0*dz)
- cellno=numpy.arange(1, self.num[2])
+ cellno = numpy.arange(1, self.num[2])
- fig=plt.figure()
- titlesize=12
+ fig = plt.figure()
+ titlesize = 12
- plt.subplot(1,3,1)
+ plt.subplot(1, 3, 1)
plt.title('Pressure', fontsize=titlesize)
plt.ylabel('$i_z$')
plt.xlabel('$p_z$')
- plt.plot(self.p_f[0,0,:], numpy.arange(self.num[2]))
+ plt.plot(self.p_f[0, 0, :], numpy.arange(self.num[2]))
plt.grid()
- plt.subplot(1,3,2)
+ plt.subplot(1, 3, 2)
plt.title('Pressure gradient', fontsize=titlesize)
plt.ylabel('$i_z$')
plt.xlabel('$\Delta p_z$')
plt.plot(dp_dz, cellno)
plt.grid()
- plt.subplot(1,3,3)
+ plt.subplot(1, 3, 3)
plt.title('Velocity prediction terms', fontsize=titlesize)
plt.ylabel('$i_z$')
plt.xlabel('$\Delta v_z$')
plt.plot(dvz_diff, cellno, label='Diffusion')
plt.plot(dvz_adv, cellno, label='Advection')
plt.plot(dvz_diff+dvz_adv, cellno, '--', label='Sum')
- leg=plt.legend(loc='best', prop={'size':8})
+ leg = plt.legend(loc='best', prop={'size':8})
leg.get_frame().set_alpha(0.5)
plt.grid()
plt.tight_layout()
- filename='../output/{}-diff_adv-t={:.2e}s-mu={:.2e}Pa-s.{}'.format(
- self.sid, self.time_current[0], self.mu[0], graphics_format)
+ filename = '../output/{}-diff_adv-t={:.2e}s-mu={:.2e}Pa-s.{}'\
+ .format(self.sid, self.time_current[0], self.mu[0],
+ graphics_format)
plt.savefig(filename)
if verbose:
print('saved to ' + filename)
t@@ -6007,16 +5919,16 @@ class sim:
'''
# find magnitude of fluid velocity vectors
- self.v_f_magn=numpy.empty_like(self.p_f)
+ self.v_f_magn = numpy.empty_like(self.p_f)
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_magn[x,y,z]=\
- self.v_f[x,y,z,:].dot(self.v_f[x,y,z,:])
+ self.v_f_magn[x, y, z] = \
+ self.v_f[x, y, z, :].dot(self.v_f[x, y, z, :])
- self.Re=self.rho_f*self.v_f_magn*self.L[0]/self.num[0]/(self.mu + \
+ Re = self.rho_f*self.v_f_magn*self.L[0]/self.num[0]/(self.mu + \
1.0e-16)
- return self.Re
+ return Re
def plotLoadCurve(self, graphics_format='png', verbose=True):
'''
t@@ -6037,57 +5949,56 @@ class sim:
print('Error: matplotlib module not found (plotLoadCurve).')
return
- t=numpy.empty(self.status())
- H=numpy.empty_like(t)
- sb=sim(self.sid, fluid=self.fluid)
+ t = numpy.empty(self.status())
+ H = numpy.empty_like(t)
+ sb = sim(self.sid, fluid=self.fluid)
sb.readfirst(verbose=False)
- load=0.0
for i in numpy.arange(1, self.status()+1):
sb.readstep(i, verbose=False)
if i == 0:
- load=sb.w_sigma0[0]
- t[i-1] =sb.time_current[0]
- H[i-1]=sb.w_x[0]
+ load = sb.w_sigma0[0]
+ t[i-1] = sb.time_current[0]
+ H[i-1] = sb.w_x[0]
# find consolidation parameters
- self.H0=H[0]
- self.H100=H[-1]
- self.H50=(self.H0 + self.H100)/2.0
- T50=0.197 # case I
+ H0 = H[0]
+ H100 = H[-1]
+ H50 = (H0 + H100)/2.0
+ T50 = 0.197 # case I
# find the time where 50% of the consolidation (H50) has happened by
# linear interpolation. The values in H are expected to be
# monotonically decreasing. See Numerical Recipies p. 115
- i_lower=0
- i_upper=self.status()-1
- while (i_upper - i_lower > 1):
- i_midpoint=int((i_upper + i_lower)/2)
- if self.H50 < H[i_midpoint]:
- i_lower=i_midpoint
+ i_lower = 0
+ i_upper = self.status()-1
+ while i_upper - i_lower > 1:
+ i_midpoint = int((i_upper + i_lower)/2)
+ if H50 < H[i_midpoint]:
+ i_lower = i_midpoint
else:
- i_upper=i_midpoint
- self.t50=t[i_lower] + (t[i_upper] - t[i_lower]) * \
- (self.H50 - H[i_lower])/(H[i_upper] - H[i_lower])
+ i_upper = i_midpoint
+ t50 = t[i_lower] + (t[i_upper] - t[i_lower]) * \
+ (H50 - H[i_lower])/(H[i_upper] - H[i_lower])
- self.c_coeff=T50*self.H50**2.0/(self.t50)
+ c_coeff = T50*H50**2.0/(t50)
if self.fluid:
- e=numpy.mean(sb.phi[:,:,3:-8]) # ignore boundaries
+ e = numpy.mean(sb.phi[:, :, 3:-8]) # ignore boundaries
else:
- e=sb.voidRatio()
+ e = sb.voidRatio()
- self.phi_bar=e
- fig=plt.figure()
+ phi_bar = e
+ fig = plt.figure()
plt.xlabel('Time [s]')
plt.ylabel('Height [m]')
plt.title('$c_v$=%.2e m$^2$ s$^{-1}$ at %.1f kPa and $e$=%.2f' \
- % (self.c_coeff, sb.w_sigma0[0]/1000.0, e))
+ % (c_coeff, sb.w_sigma0[0]/1000.0, e))
plt.semilogx(t, H, '+-')
- plt.axhline(y=self.H0, color='gray')
- plt.axhline(y=self.H50, color='gray')
- plt.axhline(y=self.H100, color='gray')
- plt.axvline(x=self.t50, color='red')
+ plt.axhline(y=H0, color='gray')
+ plt.axhline(y=H50, color='gray')
+ plt.axhline(y=H100, color='gray')
+ plt.axvline(x=t50, color='red')
plt.grid()
- filename=self.sid + '-loadcurve.' + graphics_format
+ filename = self.sid + '-loadcurve.' + graphics_format
plt.savefig(filename)
if verbose:
print('saved to ' + filename)
t@@ -6102,8 +6013,7 @@ class sim:
See also: :func:`plotConvergence()`
'''
- self.conv=numpy.loadtxt('../output/' + self.sid + '-conv.log',
- dtype=numpy.int32)
+ return numpy.loadtxt('../output/' + self.sid + '-conv.log', dtype=numpy.int32)
def plotConvergence(self, graphics_format='png', verbose=True):
'''
t@@ -6122,15 +6032,15 @@ class sim:
print('Error: matplotlib module not found (plotConvergence).')
return
- fig=plt.figure()
- self.convergence()
+ fig = plt.figure()
+ conv = self.convergence()
plt.title('Convergence evolution in CFD solver in "' + self.sid + '"')
plt.xlabel('Time step')
plt.ylabel('Jacobi iterations')
- plt.plot(self.conv[:,0], self.conv[:,1])
+ plt.plot(conv[:, 0], conv[:, 1])
plt.grid()
- filename=self.sid + '-conv.' + graphics_format
+ filename = self.sid + '-conv.' + graphics_format
plt.savefig(filename)
if verbose:
print('saved to ' + filename)
t@@ -6167,17 +6077,17 @@ class sim:
print('Error: matplotlib module not found (plotSinFunction).')
return
- fig=plt.figure(figsize=[8,6])
- steps_left=(self.time_total[0] - self.time_current[0]) \
+ fig = plt.figure(figsize=[8, 6])
+ steps_left = (self.time_total[0] - self.time_current[0]) \
/self.time_file_dt[0]
- t=numpy.linspace(self.time_current[0], self.time_total[0], steps_left)
- f=baseval + A*numpy.sin(2.0*numpy.pi*f*t + phi)
+ t = numpy.linspace(self.time_current[0], self.time_total[0], steps_left)
+ f = baseval + A*numpy.sin(2.0*numpy.pi*f*t + phi)
plt.plot(t, f, plotstyle)
plt.grid()
plt.xlabel(xlabel)
plt.ylabel(ylabel)
plt.tight_layout()
- filename=self.sid + '-sin.' + outformat
+ filename = self.sid + '-sin.' + outformat
plt.savefig(filename)
if verbose:
print(filename)
t@@ -6199,12 +6109,12 @@ class sim:
See also: :func:`setFluidPressureModulation()` and
:func:`disableTopWallNormalStressModulation()`
'''
- self.w_sigma0_A[0]=A
- self.w_sigma0_f[0]=f
+ self.w_sigma0_A[0] = A
+ self.w_sigma0_f[0] = f
if plot and py_mpl:
self.plotSinFunction(self.w_sigma0[0], A, f, phi=0.0,
- xlabel='$t$ [s]', ylabel='$\\sigma_0$ [Pa]')
+ xlabel='$t$ [s]', ylabel='$\\sigma_0$ [Pa]')
def disableTopWallNormalStressModulation(self):
'''
t@@ -6232,13 +6142,13 @@ class sim:
See also: :func:`setTopWallNormalStressModulation()` and
:func:`disableFluidPressureModulation()`
'''
- self.p_mod_A[0]=A
- self.p_mod_f[0]=f
- self.p_mod_phi[0]=phi
+ self.p_mod_A[0] = A
+ self.p_mod_f[0] = f
+ self.p_mod_phi[0] = phi
if plot:
- self.plotSinFunction(self.p_f[0,0,-1], A, f, phi=0.0,
- xlabel='$t$ [s]', ylabel='$p_f$ [kPa]')
+ self.plotSinFunction(self.p_f[0, 0, -1], A, f, phi=0.0,
+ xlabel='$t$ [s]', ylabel='$p_f$ [kPa]')
def disableFluidPressureModulation(self):
'''
t@@ -6262,20 +6172,17 @@ class sim:
'(plotPrescribedFluidPressures).')
return
- fig=plt.figure()
- conv=numpy.loadtxt('../output/' + self.sid + '-conv.log')
+ fig = plt.figure()
plt.title('Prescribed fluid pressures at the top in "' + self.sid + '"')
plt.xlabel('Time [s]')
plt.ylabel('Pressure [Pa]')
- t=numpy.linspace(0,self.time_total,\
- self.time_total/self.time_file_dt)
- p=self.p_f[0,0,-1] + \
- self.p_mod_A * \
- numpy.sin(2.0*numpy.pi*self.p_mod_f*t + self.p_mod_phi)
+ t = numpy.linspace(0, self.time_total, self.time_total/self.time_file_dt)
+ p = self.p_f[0, 0, -1] + self.p_mod_A * \
+ numpy.sin(2.0*numpy.pi*self.p_mod_f*t + self.p_mod_phi)
plt.plot(t, p, '.-')
plt.grid()
- filename='../output/' + self.sid + '-pres.' + graphics_format
+ filename = '../output/' + self.sid + '-pres.' + graphics_format
plt.savefig(filename)
if verbose:
print('saved to ' + filename)
t@@ -6294,8 +6201,8 @@ class sim:
:return type: numpy.array
'''
if idx == -1:
- idx=range(self.np)
- return self.force[idx,:]/(V_sphere(self.radius[idx])*self.rho[0]) + \
+ idx = range(self.np)
+ return self.force[idx, :]/(V_sphere(self.radius[idx])*self.rho[0]) + \
self.g
def setGamma(self, gamma):
t@@ -6315,7 +6222,7 @@ class sim:
:func:`setBeta()`, :func:`setTolerance()`,
:func:`setDEMstepsPerCFDstep()` and :func:`setMaxIterations()`
'''
- self.gamma=numpy.asarray(gamma)
+ self.gamma = numpy.asarray(gamma)
def setTheta(self, theta):
'''
t@@ -6335,7 +6242,7 @@ class sim:
:func:`setBeta()`, :func:`setTolerance()`,
:func:`setDEMstepsPerCFDstep()` and :func:`setMaxIterations()`
'''
- self.theta=numpy.asarray(theta)
+ self.theta = numpy.asarray(theta)
def setBeta(self, beta):
t@@ -6355,7 +6262,7 @@ class sim:
:func:`setDEMstepsPerCFDstep()` and
:func:`setMaxIterations()`
'''
- self.beta=numpy.asarray(beta)
+ self.beta = numpy.asarray(beta)
def setTolerance(self, tolerance):
'''
t@@ -6373,7 +6280,7 @@ class sim:
:func:`setTheta()`, :func:`setBeta()`, :func:`setDEMstepsPerCFDstep()` and
:func:`setMaxIterations()`
'''
- self.tolerance=numpy.asarray(tolerance)
+ self.tolerance = numpy.asarray(tolerance)
def setMaxIterations(self, maxiter):
'''
t@@ -6393,7 +6300,7 @@ class sim:
:func:`setTheta()`, :func:`setBeta()`, :func:`setDEMstepsPerCFDstep()`
and :func:`setTolerance()`
'''
- self.maxiter=numpy.asarray(maxiter)
+ self.maxiter = numpy.asarray(maxiter)
def setDEMstepsPerCFDstep(self, ndem):
'''
t@@ -6409,7 +6316,7 @@ class sim:
:func:`setTheta()`, :func:`setBeta()`, :func:`setTolerance()` and
:func:`setMaxIterations()`.
'''
- self.ndem=numpy.asarray(ndem)
+ self.ndem = numpy.asarray(ndem)
def shearStress(self, type='effective'):
'''
t@@ -6428,13 +6335,13 @@ class sim:
elif type == 'effective':
- fixvel=numpy.nonzero(self.fixvel > 0.0)
- force=numpy.zeros(3)
+ fixvel = numpy.nonzero(self.fixvel > 0.0)
+ force = numpy.zeros(3)
# Summation of shear stress contributions
for i in fixvel[0]:
- if self.vel[i,0] > 0.0:
- force += -self.force[i,:]
+ if self.vel[i, 0] > 0.0:
+ force += -self.force[i, :]
return force[0]/(self.L[0]*self.L[1])
t@@ -6484,8 +6391,8 @@ class sim:
:type smoothing_window: str
'''
- lastfile=self.status()
- sb=sim(sid=self.sid, np=self.np, nw=self.nw, fluid=self.fluid)
+ lastfile = self.status()
+ sb = sim(sid=self.sid, np=self.np, nw=self.nw, fluid=self.fluid)
if not py_mpl:
print('Error: matplotlib module not found (visualize).')
t@@ -6494,118 +6401,118 @@ class sim:
### Plotting
if outformat != 'txt':
if figsize:
- fig=plt.figure(figsize=figsize)
+ fig = plt.figure(figsize=figsize)
else:
- fig=plt.figure(figsize=(8,8))
+ fig = plt.figure(figsize=(8, 8))
if method == 'energy':
if figsize:
- fig=plt.figure(figsize=figsize)
+ fig = plt.figure(figsize=figsize)
else:
- fig=plt.figure(figsize=(20,8))
+ fig = plt.figure(figsize=(20, 8))
# Allocate arrays
- t=numpy.zeros(lastfile-firststep + 1)
- Epot=numpy.zeros_like(t)
- Ekin=numpy.zeros_like(t)
- Erot=numpy.zeros_like(t)
- Es =numpy.zeros_like(t)
- Ev =numpy.zeros_like(t)
- Es_dot=numpy.zeros_like(t)
- Ev_dot=numpy.zeros_like(t)
- Ebondpot=numpy.zeros_like(t)
- Esum=numpy.zeros_like(t)
+ t = numpy.zeros(lastfile-firststep + 1)
+ Epot = numpy.zeros_like(t)
+ Ekin = numpy.zeros_like(t)
+ Erot = numpy.zeros_like(t)
+ Es = numpy.zeros_like(t)
+ Ev = numpy.zeros_like(t)
+ Es_dot = numpy.zeros_like(t)
+ Ev_dot = numpy.zeros_like(t)
+ Ebondpot = numpy.zeros_like(t)
+ Esum = numpy.zeros_like(t)
# Read energy values from simulation binaries
for i in numpy.arange(firststep, lastfile+1):
sb.readstep(i, verbose=False)
- Epot[i]=sb.energy("pot")
- Ekin[i]=sb.energy("kin")
- Erot[i]=sb.energy("rot")
- Es[i] =sb.energy("shear")
- Ev[i] =sb.energy("visc_n")
- Es_dot[i]=sb.energy("shearrate")
- Ev_dot[i]=sb.energy("visc_n_rate")
- Ebondpot[i]=sb.energy("bondpot")
- Esum[i]=Epot[i] + Ekin[i] + Erot[i] + Es[i] + Ev[i] +\
+ Epot[i] = sb.energy("pot")
+ Ekin[i] = sb.energy("kin")
+ Erot[i] = sb.energy("rot")
+ Es[i] = sb.energy("shear")
+ Ev[i] = sb.energy("visc_n")
+ Es_dot[i] = sb.energy("shearrate")
+ Ev_dot[i] = sb.energy("visc_n_rate")
+ Ebondpot[i] = sb.energy("bondpot")
+ Esum[i] = Epot[i] + Ekin[i] + Erot[i] + Es[i] + Ev[i] +\
Ebondpot[i]
- t[i]=sb.currentTime()
+ t[i] = sb.currentTime()
if outformat != 'txt':
# Potential energy
- ax1=plt.subplot2grid((2,5),(0,0))
+ ax1 = plt.subplot2grid((2, 5), (0, 0))
ax1.set_xlabel('Time [s]')
ax1.set_ylabel('Total potential energy [J]')
ax1.plot(t, Epot, '+-')
ax1.grid()
# Kinetic energy
- ax2=plt.subplot2grid((2,5),(0,1))
+ ax2 = plt.subplot2grid((2, 5), (0, 1))
ax2.set_xlabel('Time [s]')
ax2.set_ylabel('Total kinetic energy [J]')
ax2.plot(t, Ekin, '+-')
ax2.grid()
# Rotational energy
- ax3=plt.subplot2grid((2,5),(0,2))
+ ax3 = plt.subplot2grid((2, 5), (0, 2))
ax3.set_xlabel('Time [s]')
ax3.set_ylabel('Total rotational energy [J]')
ax3.plot(t, Erot, '+-')
ax3.grid()
# Bond energy
- ax4=plt.subplot2grid((2,5),(0,3))
+ ax4 = plt.subplot2grid((2, 5), (0, 3))
ax4.set_xlabel('Time [s]')
ax4.set_ylabel('Bond energy [J]')
ax4.plot(t, Ebondpot, '+-')
ax4.grid()
# Total energy
- ax5=plt.subplot2grid((2,5),(0,4))
+ ax5 = plt.subplot2grid((2, 5), (0, 4))
ax5.set_xlabel('Time [s]')
ax5.set_ylabel('Total energy [J]')
ax5.plot(t, Esum, '+-')
ax5.grid()
# Shear energy rate
- ax6=plt.subplot2grid((2,5),(1,0))
+ ax6 = plt.subplot2grid((2, 5), (1, 0))
ax6.set_xlabel('Time [s]')
ax6.set_ylabel('Frictional dissipation rate [W]')
ax6.plot(t, Es_dot, '+-')
ax6.grid()
# Shear energy
- ax7=plt.subplot2grid((2,5),(1,1))
+ ax7 = plt.subplot2grid((2, 5), (1, 1))
ax7.set_xlabel('Time [s]')
ax7.set_ylabel('Total frictional dissipation [J]')
ax7.plot(t, Es, '+-')
ax7.grid()
# Visc_n energy rate
- ax8=plt.subplot2grid((2,5),(1,2))
+ ax8 = plt.subplot2grid((2, 5), (1, 2))
ax8.set_xlabel('Time [s]')
ax8.set_ylabel('Viscous dissipation rate [W]')
ax8.plot(t, Ev_dot, '+-')
ax8.grid()
# Visc_n energy
- ax9=plt.subplot2grid((2,5),(1,3))
+ ax9 = plt.subplot2grid((2, 5), (1, 3))
ax9.set_xlabel('Time [s]')
ax9.set_ylabel('Total viscous dissipation [J]')
ax9.plot(t, Ev, '+-')
ax9.grid()
# Combined view
- ax10=plt.subplot2grid((2,5),(1,4))
+ ax10 = plt.subplot2grid((2, 5), (1, 4))
ax10.set_xlabel('Time [s]')
ax10.set_ylabel('Energy [J]')
ax10.plot(t, Epot, '+-g')
ax10.plot(t, Ekin, '+-b')
ax10.plot(t, Erot, '+-r')
- ax10.legend(('$\sum E_{pot}$','$\sum E_{kin}$',\
- '$\sum E_{rot}$'), 'upper right', shadow=True)
+ ax10.legend(('$\sum E_{pot}$', '$\sum E_{kin}$',
+ '$\sum E_{rot}$'), 'upper right', shadow=True)
ax10.grid()
if xlim:
t@@ -6630,33 +6537,32 @@ class sim:
# Allocate arrays on first run
if i == firststep:
- wforce=numpy.zeros((lastfile+1)*sb.nw,\
+ wforce = numpy.zeros((lastfile+1)*sb.nw,\
dtype=numpy.float64).reshape((lastfile+1), sb.nw)
- wvel =numpy.zeros((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,\
+ wpos = numpy.zeros((lastfile+1)*sb.nw,\
dtype=numpy.float64).reshape((lastfile+1), sb.nw)
- wsigma0 =numpy.zeros((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)
+ maxpos = numpy.zeros((lastfile+1), dtype=numpy.float64)
+ logstress = numpy.zeros((lastfile+1), dtype=numpy.float64)
+ voidratio = numpy.zeros((lastfile+1), dtype=numpy.float64)
- wforce[i]=sb.w_force[0]
- wvel[i] =sb.w_vel[0]
- wpos[i] =sb.w_x[0]
- wsigma0[i] =sb.w_sigma0[0]
- maxpos[i]=numpy.max(sb.x[:,2]+sb.radius)
- logstress[i] =\
- numpy.log((sb.w_force[0]/(sb.L[0]*sb.L[1]))/1000.0)
- voidratio[i]=sb.voidRatio()
+ wforce[i] = sb.w_force[0]
+ wvel[i] = sb.w_vel[0]
+ wpos[i] = sb.w_x[0]
+ wsigma0[i] = sb.w_sigma0[0]
+ maxpos[i] = numpy.max(sb.x[:, 2]+sb.radius)
+ logstress[i] = numpy.log((sb.w_force[0]/(sb.L[0]*sb.L[1]))/1000.0)
+ voidratio[i] = sb.voidRatio()
- t=numpy.linspace(0.0, sb.time_current, lastfile+1)
+ t = numpy.linspace(0.0, sb.time_current, lastfile+1)
# Plotting
if outformat != 'txt':
# linear plot of time vs. wall position
- ax1=plt.subplot2grid((2,2),(0,0))
+ ax1 = plt.subplot2grid((2, 2), (0, 0))
ax1.set_xlabel('Time [s]')
ax1.set_ylabel('Position [m]')
ax1.plot(t, wpos, '+-', label="upper wall")
t@@ -6664,27 +6570,27 @@ class sim:
ax1.legend()
ax1.grid()
- #ax2=plt.subplot2grid((2,2),(1,0))
+ #ax2 = plt.subplot2grid((2, 2), (1, 0))
#ax2.set_xlabel('Time [s]')
#ax2.set_ylabel('Force [N]')
#ax2.plot(t, wforce, '+-')
# semilog plot of log stress vs. void ratio
- ax2=plt.subplot2grid((2,2),(1,0))
+ ax2 = plt.subplot2grid((2, 2), (1, 0))
ax2.set_xlabel('log deviatoric stress [kPa]')
ax2.set_ylabel('Void ratio [-]')
ax2.plot(logstress, voidratio, '+-')
ax2.grid()
# linear plot of time vs. wall velocity
- ax3=plt.subplot2grid((2,2),(0,1))
+ ax3 = plt.subplot2grid((2, 2), (0, 1))
ax3.set_xlabel('Time [s]')
ax3.set_ylabel('Velocity [m/s]')
ax3.plot(t, wvel, '+-')
ax3.grid()
# linear plot of time vs. deviatoric stress
- ax4=plt.subplot2grid((2,2),(1,1))
+ ax4 = plt.subplot2grid((2, 2), (1, 1))
ax4.set_xlabel('Time [s]')
ax4.set_ylabel('Deviatoric stress [Pa]')
ax4.plot(t, wsigma0, '+-', label="$\sigma_0$")
t@@ -6704,26 +6610,26 @@ class sim:
for i in numpy.arange(firststep, lastfile+1):
sb.readstep(i, verbose=False)
- vol=(sb.w_x[0]-sb.origo[2]) * (sb.w_x[1]-sb.w_x[2]) \
+ vol = (sb.w_x[0]-sb.origo[2]) * (sb.w_x[1]-sb.w_x[2]) \
* (sb.w_x[3] - sb.w_x[4])
# Allocate arrays on first run
if i == firststep:
- axial_strain=numpy.zeros(lastfile+1, dtype=numpy.float64)
+ axial_strain = numpy.zeros(lastfile+1, dtype=numpy.float64)
deviatoric_stress =\
numpy.zeros(lastfile+1, dtype=numpy.float64)
volumetric_strain =\
numpy.zeros(lastfile+1, dtype=numpy.float64)
- w0pos0=sb.w_x[0]
- vol0=vol
+ w0pos0 = sb.w_x[0]
+ vol0 = vol
- sigma1=sb.w_force[0]/\
+ sigma1 = sb.w_force[0]/\
((sb.w_x[1]-sb.w_x[2])*(sb.w_x[3]-sb.w_x[4]))
- axial_strain[i]=(w0pos0 - sb.w_x[0])/w0pos0
- volumetric_strain[i]=(vol0-vol)/vol0
- deviatoric_stress[i]=sigma1 / sb.w_sigma0[1]
+ axial_strain[i] = (w0pos0 - sb.w_x[0])/w0pos0
+ volumetric_strain[i] = (vol0-vol)/vol0
+ deviatoric_stress[i] = sigma1 / sb.w_sigma0[1]
#print(lastfile)
#print(axial_strain)
t@@ -6734,20 +6640,20 @@ class sim:
if outformat != 'txt':
# linear plot of deviatoric stress
- ax1=plt.subplot2grid((2,1),(0,0))
+ ax1 = plt.subplot2grid((2, 1), (0, 0))
ax1.set_xlabel('Axial strain, $\gamma_1$, [-]')
ax1.set_ylabel('Deviatoric stress, $\sigma_1 - \sigma_3$, [Pa]')
ax1.plot(axial_strain, deviatoric_stress, '+-')
#ax1.legend()
ax1.grid()
- #ax2=plt.subplot2grid((2,2),(1,0))
+ #ax2 = plt.subplot2grid((2, 2), (1, 0))
#ax2.set_xlabel('Time [s]')
#ax2.set_ylabel('Force [N]')
#ax2.plot(t, wforce, '+-')
# semilog plot of log stress vs. void ratio
- ax2=plt.subplot2grid((2,1),(1,0))
+ ax2 = plt.subplot2grid((2, 1), (1, 0))
ax2.set_xlabel('Axial strain, $\gamma_1$ [-]')
ax2.set_ylabel('Volumetric strain, $\gamma_v$, [-]')
ax2.plot(axial_strain, volumetric_strain, '+-')
t@@ -6766,130 +6672,119 @@ class sim:
# First iteration: Allocate arrays and find constant values
if i == firststep:
# Shear displacement
- self.xdisp =numpy.zeros(lastfile+1, dtype=numpy.float64)
+ xdisp = numpy.zeros(lastfile+1, dtype=numpy.float64)
# Normal stress
- self.sigma_eff= numpy.zeros(lastfile+1, dtype=numpy.float64)
+ sigma_eff = numpy.zeros(lastfile+1, dtype=numpy.float64)
# Normal stress
- self.sigma_def= numpy.zeros(lastfile+1, dtype=numpy.float64)
+ sigma_def = numpy.zeros(lastfile+1, dtype=numpy.float64)
# Shear stress
- self.tau =numpy.zeros(lastfile+1, dtype=numpy.float64)
+ tau = numpy.zeros(lastfile+1, dtype=numpy.float64)
# Upper wall position
- self.dilation=numpy.zeros(lastfile+1, dtype=numpy.float64)
+ dilation = numpy.zeros(lastfile+1, dtype=numpy.float64)
+
+ # Peak shear stress
+ tau_p = 0.0
- # Upper wall position
- self.tau_p=0.0 # Peak shear stress
# Shear strain value of peak sh. stress
- self.tau_p_shearstrain=0.0
+ tau_p_shearstrain = 0.0
- fixvel=numpy.nonzero(sb.fixvel > 0.0)
- #fixvel_upper=numpy.nonzero(sb.vel[fixvel,0] > 0.0)
- shearvel=sb.vel[fixvel,0].max()
- w_x0=sb.w_x[0] # Original height
- A=sb.L[0] * sb.L[1] # Upper surface area
+ fixvel = numpy.nonzero(sb.fixvel > 0.0)
+ #fixvel_upper = numpy.nonzero(sb.vel[fixvel, 0] > 0.0)
+ shearvel = sb.vel[fixvel, 0].max()
+ w_x0 = sb.w_x[0] # Original height
+ A = sb.L[0] * sb.L[1] # Upper surface area
if i == firststep+1:
- w_x0=sb.w_x[0] # Original height
+ w_x0 = sb.w_x[0] # Original height
# Summation of shear stress contributions
for j in fixvel[0]:
- if sb.vel[j,0] > 0.0:
- self.tau[i] += -sb.force[j,0]/A
+ if sb.vel[j, 0] > 0.0:
+ tau[i] += -sb.force[j, 0]/A
if i > 0:
- self.xdisp[i]=self.xdisp[i-1] +sb.time_file_dt[0]*shearvel
- self.sigma_eff[i]=sb.w_force[0] / A
- self.sigma_def[i]=sb.w_sigma0[0]
+ xdisp[i] = self.xdisp[i-1] + \
+ sb.time_file_dt[0]*shearvel
+ sigma_eff[i] = sb.w_force[0] / A
+ sigma_def[i] = sb.w_sigma0[0]
# dilation in meters
- #dilation[i]=sb.w_x[0] - w_x0
+ #dilation[i] = sb.w_x[0] - w_x0
# dilation in percent
- #dilation[i]=(sb.w_x[0] - w_x0)/w_x0 * 100.0 # dilation in percent
+ #dilation[i] = (sb.w_x[0] - w_x0)/w_x0 * 100.0 # dilation in percent
# dilation in number of mean particle diameters
- d_bar=numpy.mean(self.radius)*2.0
+ d_bar = numpy.mean(self.radius)*2.0
if numpy.isnan(d_bar):
- #raise Exception("Error, d_bar is NaN. Please check that the"
- # + " radii are initialized.")
print('No radii in self.radius, attempting to read first '
- + 'file')
+ + 'file')
self.readfirst()
- d_bar=numpy.mean(self.radius)*2.0
- self.dilation[i]=(sb.w_x[0] - w_x0)/d_bar
+ d_bar = numpy.mean(self.radius)*2.0
+ dilation[i] = (sb.w_x[0] - w_x0)/d_bar
# Test if this was the max. shear stress
- if self.tau[i] > self.tau_p:
- self.tau_p=self.tau[i]
- self.tau_p_shearstrain=self.xdisp[i]/w_x0
+ if tau[i] > tau_p:
+ tau_p = tau[i]
+ tau_p_shearstrain = xdisp[i]/w_x0
- self.shear_strain=self.xdisp/w_x0
+ shear_strain = xdisp/w_x0
# Copy values so they can be modified during smoothing
- shear_strain=self.shear_strain
- tau=self.tau
- sigma_def=self.sigma_def
+ shear_strain_smooth = self.shear_strain
+ tau_smooth = tau
+ sigma_def_smooth = sigma_def
# Optionally smooth the shear stress
if smoothing > 2:
- if not smoothing_window in ['flat', 'hanning', 'hamming',
- 'bartlett', 'blackman']:
+ if smoothing_window not in ['flat', 'hanning', 'hamming',
+ 'bartlett', 'blackman']:
raise ValueError
- s=numpy.r_[2*tau[0]-tau[smoothing:1:-1], tau,
- 2*tau[-1]-tau[-1:-smoothing:-1]]
+ s = numpy.r_[2*tau[0]-tau[smoothing:1:-1], tau,
+ 2*tau[-1]-tau[-1:-smoothing:-1]]
if smoothing_window == 'flat': # moving average
- w=numpy.ones(smoothing, 'd')
+ w = numpy.ones(smoothing, 'd')
else:
- w=getattr(numpy, smoothing_window)(smoothing)
- y=numpy.convolve(w/w.sum(), s, mode='same')
- tau=y[smoothing-1:-smoothing+1]
+ w = getattr(np, smoothing_window)(smoothing)
+ y = numpy.convolve(w/w.sum(), s, mode='same')
+ tau_smooth = y[smoothing-1:-smoothing+1]
# Plot stresses
if outformat != 'txt':
- shearinfo="$\\tau_p$={:.3} Pa at $\gamma$={:.3}".format(\
+ shearinfo = "$\\tau_p$={:.3} Pa at $\gamma$={:.3}".format(\
self.tau_p, self.tau_p_shearstrain)
fig.text(0.01, 0.01, shearinfo, horizontalalignment='left',
fontproperties=FontProperties(size=14))
- ax1=plt.subplot2grid((2,1), (0,0))
+ ax1 = plt.subplot2grid((2, 1), (0, 0))
ax1.set_xlabel('Shear strain [-]')
- #ax1.set_ylabel('Stress [Pa]')
- #ax1.set_ylabel('Shear friction $\\tau/\\sigma\'$ [-]')
ax1.set_ylabel('Shear friction $\\tau/\\sigma_0$ [-]')
- #ax1.plot(xdisp / w_x0, sigma_eff, '+-g', label="$\sigma'$")
- #ax1.plot(xdisp / w_x0, sigma_def, '+-b', label="$\sigma_0$")
- #ax1.plot(xdisp / w_x0, tau, '+-r', label="$\\tau$")
if smoothing > 2:
- ax1.plot(shear_strain[1:-(smoothing+1)/2],
- tau[1:-(smoothing+1)/2] /
- sigma_def[1:-(smoothing+1)/2],
- '-', label="$\\tau/\\sigma_0$")
+ ax1.plot(shear_strain_smooth[1:-(smoothing+1)/2],
+ tau_smooth[1:-(smoothing+1)/2] /
+ sigma_def_smooth[1:-(smoothing+1)/2],
+ '-', label="$\\tau/\\sigma_0$")
else:
ax1.plot(shear_strain[1:],\
- tau[1:]/sigma_def[1:],\
- #self.tau[1:]/self.sigma_eff[1:],\
- '-', label="$\\tau/\\sigma_0$")
- #'.-', label="$\\tau$")
- #ax1.legend(loc=4)
+ tau[1:]/sigma_def[1:],\
+ '-', label="$\\tau/\\sigma_0$")
ax1.grid()
# Plot dilation
- ax2=plt.subplot2grid((2,1),(1,0))
+ ax2 = plt.subplot2grid((2, 1), (1, 0))
ax2.set_xlabel('Shear strain [-]')
- #ax2.set_ylabel('Dilation [m]')
- #ax2.set_ylabel('Dilation [%]')
ax2.set_ylabel('Dilation, $\Delta h/(2\\bar{r})$ [m]')
- #ax2.plot(self.shear_strain, self.dilation, '.-')
if smoothing > 2:
- ax2.plot(self.shear_strain[1:-(smoothing+1)/2],
- self.dilation[1:-(smoothing+1)/2], '-')
+ ax2.plot(shear_strain_smooth[1:-(smoothing+1)/2],
+ dilation_smooth[1:-(smoothing+1)/2], '-')
else:
- ax2.plot(self.shear_strain, self.dilation, '-')
+ ax2.plot(shear_strain, self.dilation, '-')
ax2.grid()
if xlim:
t@@ -6898,28 +6793,26 @@ class sim:
fig.tight_layout()
- else :
+ else:
# Write values to textfile
- filename="shear-stresses-{0}.txt".format(self.sid)
+ filename = "shear-stresses-{0}.txt".format(self.sid)
#print("Writing stress data to " + filename)
- fh=None
- try :
- fh=open(filename, "w")
- L=sb.L[2] - sb.origo[2] # Initial height
+ fh = None
+ try:
+ fh = open(filename, "w")
for i in numpy.arange(firststep, lastfile+1):
# format: shear distance [mm], sigma [kPa], tau [kPa],
# Dilation [%]
- fh.write("{0}\t{1}\t{2}\t{3}\n".format(xdisp[i],
- sigma_eff[i]/1000.0,
- tau[i]/1000.0,
- dilation[i]))
- finally :
+ fh.write("{0}\t{1}\t{2}\t{3}\n"
+ .format(xdisp[i], sigma_eff[i]/1000.0,
+ tau[i]/1000.0, dilation[i]))
+ finally:
if fh is not None:
fh.close()
elif method == 'shear-displacement':
- time=numpy.zeros(lastfile+1, dtype=numpy.float64)
+ time = numpy.zeros(lastfile+1, dtype=numpy.float64)
# Read stress values from simulation binaries
for i in numpy.arange(firststep, lastfile+1):
sb.readstep(i, verbose=False)
t@@ -6928,176 +6821,162 @@ class sim:
if i == firststep:
# Shear displacement
- self.xdisp =numpy.zeros(lastfile+1, dtype=numpy.float64)
+ xdisp = numpy.zeros(lastfile+1, dtype=numpy.float64)
# Normal stress
- self.sigma_eff= numpy.zeros(lastfile+1, dtype=numpy.float64)
+ sigma_eff = numpy.zeros(lastfile+1, dtype=numpy.float64)
# Normal stress
- self.sigma_def= numpy.zeros(lastfile+1, dtype=numpy.float64)
+ sigma_def = numpy.zeros(lastfile+1, dtype=numpy.float64)
# Shear stress
- self.tau_eff =numpy.zeros(lastfile+1, dtype=numpy.float64)
+ tau_eff = numpy.zeros(lastfile+1, dtype=numpy.float64)
# Upper wall position
- self.dilation=numpy.zeros(lastfile+1, dtype=numpy.float64)
+ dilation = numpy.zeros(lastfile+1, dtype=numpy.float64)
# Mean porosity
- self.phi_bar=numpy.zeros(lastfile+1, dtype=numpy.float64)
+ phi_bar = numpy.zeros(lastfile+1, dtype=numpy.float64)
# Mean fluid pressure
- self.p_f_bar=numpy.zeros(lastfile+1, dtype=numpy.float64)
- self.p_f_top=numpy.zeros(lastfile+1, dtype=numpy.float64)
+ p_f_bar = numpy.zeros(lastfile+1, dtype=numpy.float64)
+ p_f_top = numpy.zeros(lastfile+1, dtype=numpy.float64)
# Upper wall position
- self.tau_p=0.0 # Peak shear stress
+ tau_p = 0.0 # Peak shear stress
# Shear strain value of peak sh. stress
- self.tau_p_shearstrain=0.0
+ tau_p_shearstrain = 0.0
- fixvel=numpy.nonzero(sb.fixvel > 0.0)
- #fixvel_upper=numpy.nonzero(sb.vel[fixvel,0] > 0.0)
- w_x0=sb.w_x[0] # Original height
- A=sb.L[0]*sb.L[1] # Upper surface area
+ fixvel = numpy.nonzero(sb.fixvel > 0.0)
+ #fixvel_upper=numpy.nonzero(sb.vel[fixvel, 0] > 0.0)
+ w_x0 = sb.w_x[0] # Original height
+ A = sb.L[0]*sb.L[1] # Upper surface area
- d_bar=numpy.mean(sb.radius)*2.0
+ d_bar = numpy.mean(sb.radius)*2.0
# Shear velocity
- self.v=numpy.zeros(lastfile+1, dtype=numpy.float64)
+ v = numpy.zeros(lastfile+1, dtype=numpy.float64)
- time[i]=sb.time_current[0]
+ time[i] = sb.time_current[0]
if i == firststep+1:
- w_x0=sb.w_x[0] # Original height
+ w_x0 = sb.w_x[0] # Original height
# Summation of shear stress contributions
for j in fixvel[0]:
- if sb.vel[j,0] > 0.0:
- self.tau_eff[i] += -sb.force[j,0]/A
+ if sb.vel[j, 0] > 0.0:
+ tau_eff[i] += -sb.force[j, 0]/A
if i > 0:
- self.xdisp[i]=sb.xyzsum[fixvel,0].max()
- self.v[i] =sb.vel[fixvel,0].max()
+ xdisp[i] = sb.xyzsum[fixvel, 0].max()
+ v[i] = sb.vel[fixvel, 0].max()
- self.sigma_eff[i]=sb.w_force[0]/A
- self.sigma_def[i]=sb.currentNormalStress()
+ sigma_eff[i] = sb.w_force[0]/A
+ sigma_def[i] = sb.currentNormalStress()
# dilation in number of mean particle diameters
- self.dilation[i]=(sb.w_x[0] - w_x0)/d_bar
+ dilation[i] = (sb.w_x[0] - w_x0)/d_bar
- wall0_iz=int(sb.w_x[0]/(sb.L[2]/sb.num[2]))
+ wall0_iz = int(sb.w_x[0]/(sb.L[2]/sb.num[2]))
if self.fluid:
if i > 0:
- self.phi_bar[i]=numpy.mean(sb.phi[:,:,0:wall0_iz])
+ phi_bar[i] = numpy.mean(sb.phi[:, :, 0:wall0_iz])
if i == firststep+1:
- self.phi_bar[0]=self.phi_bar[1]
- self.p_f_bar[i]=numpy.mean(sb.p_f[:,:,0:wall0_iz])
- self.p_f_top[i]=sb.p_f[0,0,-1]
+ phi_bar[0] = phi_bar[1]
+ p_f_bar[i] = numpy.mean(sb.p_f[:, :, 0:wall0_iz])
+ p_f_top[i] = sb.p_f[0, 0, -1]
# Test if this was the max. shear stress
- if self.tau_eff[i] > self.tau_p:
- self.tau_p=self.tau_eff[i]
- self.tau_p_shearstrain=self.xdisp[i]/w_x0
+ if tau_eff[i] > tau_p:
+ tau_p = tau_eff[i]
+ tau_p_shearstrain = xdisp[i]/w_x0
- self.shear_strain=self.xdisp/w_x0
+ shear_strain = self.xdisp/w_x0
# Plot stresses
if outformat != 'txt':
if figsize:
- fig=plt.figure(figsize=figsize)
+ fig = plt.figure(figsize=figsize)
else:
- fig=plt.figure(figsize=(8,12))
-
- #shearinfo="$\\tau_p$={:.3} Pa at $\gamma$={:.3}".format(\
- #self.tau_p, self.tau_p_shearstrain)
- #fig.text(0.01, 0.01, shearinfo, horizontalalignment='left',
- #fontproperties=FontProperties(size=14))
+ fig = plt.figure(figsize=(8, 12))
# Upper plot
- ax1=plt.subplot(3,1,1)
- ax1.plot(time, self.xdisp, 'k', label='Displacement')
+ ax1 = plt.subplot(3, 1, 1)
+ ax1.plot(time, xdisp, 'k', label='Displacement')
ax1.set_ylabel('Horizontal displacement [m]')
- ax2=ax1.twinx()
+ ax2 = ax1.twinx()
- #ax2color='#666666'
- ax2color='blue'
+ #ax2color = '#666666'
+ ax2color = 'blue'
if self.fluid:
- ax2.plot(time, self.phi_bar, color=ax2color,
- label='Porosity')
+ ax2.plot(time, phi_bar, color=ax2color, label='Porosity')
ax2.set_ylabel('Mean porosity $\\bar{\\phi}$ [-]')
else:
- ax2.plot(time, self.dilation, color=ax2color,
- label='Dilation')
+ ax2.plot(time, dilation, color=ax2color, label='Dilation')
ax2.set_ylabel('Dilation, $\Delta h/(2\\bar{r})$ [-]')
for tl in ax2.get_yticklabels():
tl.set_color(ax2color)
# Middle plot
- ax5=plt.subplot(3, 1, 2, sharex=ax1)
- ax5.semilogy(time[1:], self.v[1:], label='Shear velocity')
+ ax5 = plt.subplot(3, 1, 2, sharex=ax1)
+ ax5.semilogy(time[1:], v[1:], label='Shear velocity')
ax5.set_ylabel('Shear velocity [ms$^{-1}$]')
# shade stick periods
- collection=\
+ collection = \
matplotlib.collections.BrokenBarHCollection.span_where(
- time, ymin=1.0e-7, ymax=1.0,
- where=numpy.isclose(self.v, 0.0),
- facecolor='black', alpha=0.2,
- linewidth=0)
+ time, ymin=1.0e-7, ymax=1.0,
+ where=numpy.isclose(v, 0.0),
+ facecolor='black', alpha=0.2,
+ linewidth=0)
ax5.add_collection(collection)
# Lower plot
- ax3=plt.subplot(3, 1, 3, sharex=ax1)
+ ax3 = plt.subplot(3, 1, 3, sharex=ax1)
if sb.w_sigma0_A > 1.0e-3:
- lns0=ax3.plot(time, self.sigma_def/1000.0,
- '-k', label="$\\sigma_0$")
- lns1=ax3.plot(time, self.sigma_eff/1000.0,
- '--k', label="$\\sigma'$")
- lns2=ax3.plot(time,
- numpy.ones_like(time)*sb.w_tau_x/1000.0,
- '-r', label="$\\tau$")
- lns3=ax3.plot(time, self.tau_eff/1000.0,
- '--r', label="$\\tau'$")
+ lns0 = ax3.plot(time, sigma_def/1000.0,
+ '-k', label="$\\sigma_0$")
+ lns1 = ax3.plot(time, sigma_eff/1000.0,
+ '--k', label="$\\sigma'$")
+ lns2 = ax3.plot(time, numpy.ones_like(time)*sb.w_tau_x/1000.0,
+ '-r', label="$\\tau$")
+ lns3 = ax3.plot(time, tau_eff/1000.0,
+ '--r', label="$\\tau'$")
ax3.set_ylabel('Stress [kPa]')
- #ax3.legend(loc='upper left')
else:
- ax3.plot(time, self.tau_eff/sb.w_sigma0[0],
- '-k', label="$Shear friction$")
+ ax3.plot(time, tau_eff/sb.w_sigma0[0],
+ '-k', label="$Shear friction$")
ax3.plot([0, time[-1]],
- [sb.w_tau_x/self.sigma_def,
- sb.w_tau_x/self.sigma_def],
- '--k', label="$Applied shear friction$")
+ [sb.w_tau_x/sigma_def, sb.w_tau_x/sigma_def],
+ '--k', label="$Applied shear friction$")
ax3.set_ylabel('Shear friction $\\tau\'/\\sigma_0$ [-]')
# axis limits
- ax3.set_ylim([sb.w_tau_x/self.sigma_def[0]*0.5,
- sb.w_tau_x/self.sigma_def[0]*1.5])
-
+ ax3.set_ylim([sb.w_tau_x/sigma_def[0]*0.5,
+ sb.w_tau_x/sigma_def[0]*1.5])
if self.fluid:
- ax4=ax3.twinx()
- #ax4color='#666666'
- ax4color=ax2color
- lns4=ax4.plot(time, self.p_f_top/1000.0, '-',
- color=ax4color,
- label='$p_\\text{f}^\\text{forcing}$')
- lns5=ax4.plot(time, self.p_f_bar/1000.0, '--',
- color=ax4color,
- label='$\\bar{p}_\\text{f}$')
+ ax4 = ax3.twinx()
+ #ax4color = '#666666'
+ ax4color = ax2color
+ lns4 = ax4.plot(time, p_f_top/1000.0, '-', color=ax4color,
+ label='$p_\\text{f}^\\text{forcing}$')
+ lns5 = ax4.plot(time, p_f_bar/1000.0, '--', color=ax4color,
+ label='$\\bar{p}_\\text{f}$')
ax4.set_ylabel('Mean fluid pressure '
- + '$\\bar{p_\\text{f}}$ [kPa]')
+ + '$\\bar{p_\\text{f}}$ [kPa]')
for tl in ax4.get_yticklabels():
tl.set_color(ax4color)
if sb.w_sigma0_A > 1.0e-3:
#ax4.legend(loc='upper right')
- lns=lns0+lns1+lns2+lns3+lns4+lns5
- labs=[l.get_label() for l in lns]
+ lns = lns0+lns1+lns2+lns3+lns4+lns5
+ labs = [l.get_label() for l in lns]
ax4.legend(lns, labs, loc='upper right',
- fancybox=True, framealpha=legend_alpha)
+ fancybox=True, framealpha=legend_alpha)
if xlim:
ax4.set_xlim(xlim)
-
# aesthetics
ax3.set_xlabel('Time [s]')
t@@ -7119,83 +6998,68 @@ class sim:
elif method == 'rate-dependence':
if figsize:
- fig=plt.figure(figsize=figsize)
+ fig = plt.figure(figsize=figsize)
else:
- fig=plt.figure(figsize=(8,6))
+ fig = plt.figure(figsize=(8, 6))
- tau=numpy.empty(sb.status())
- N=numpy.empty(sb.status())
- #v=numpy.empty(sb.status())
- shearstrainrate=numpy.empty(sb.status())
- shearstrain=numpy.empty(sb.status())
+ tau = numpy.empty(sb.status())
+ N = numpy.empty(sb.status())
+ #v = numpy.empty(sb.status())
+ shearstrainrate = numpy.empty(sb.status())
+ shearstrain = numpy.empty(sb.status())
for i in numpy.arange(firststep, sb.status()):
sb.readstep(i+1, verbose=False)
- #tau=sb.shearStress()
- tau[i]=sb.w_tau_x # defined shear stress
- #tau[i]=sb.shearStress()[0] # measured shear stress along x
- N[i]=sb.currentNormalStress() # defined normal stress
- #v[i]=sb.shearVel()
- shearstrainrate[i]=sb.shearStrainRate()
- shearstrain[i]=sb.shearStrain()
+ #tau = sb.shearStress()
+ tau[i] = sb.w_tau_x # defined shear stress
+ N[i] = sb.currentNormalStress() # defined normal stress
+ shearstrainrate[i] = sb.shearStrainRate()
+ shearstrain[i] = sb.shearStrain()
# remove nonzero sliding velocities and their associated values
- #idx=numpy.nonzero(v)
- idx=numpy.nonzero(shearstrainrate)
- #v_nonzero=v[idx]
- shearstrainrate_nonzero=shearstrainrate[idx]
- tau_nonzero=tau[idx]
- N_nonzero=N[idx]
- shearstrain_nonzero=shearstrain[idx]
-
- ax1=plt.subplot(111)
+ idx = numpy.nonzero(shearstrainrate)
+ shearstrainrate_nonzero = shearstrainrate[idx]
+ tau_nonzero = tau[idx]
+ N_nonzero = N[idx]
+ shearstrain_nonzero = shearstrain[idx]
+
+ ax1 = plt.subplot(111)
#ax1.semilogy(N/1000., v)
#ax1.semilogy(tau_nonzero/N_nonzero, v_nonzero, '+k')
#ax1.plot(tau/N, v, '.')
- friction=tau_nonzero/N_nonzero
- #CS=ax1.scatter(friction, v_nonzero, c=shearstrain_nonzero,
+ friction = tau_nonzero/N_nonzero
+ #CS = ax1.scatter(friction, v_nonzero, c=shearstrain_nonzero,
#linewidth=0)
if cmap:
- CS=ax1.scatter(friction, shearstrainrate_nonzero,
- c=shearstrain_nonzero, linewidth=0.1,
- cmap=cmap)
+ CS = ax1.scatter(friction, shearstrainrate_nonzero,
+ c=shearstrain_nonzero, linewidth=0.1,
+ cmap=cmap)
else:
- CS=ax1.scatter(friction, shearstrainrate_nonzero,
- c=shearstrain_nonzero, linewidth=0.1,
- cmap=matplotlib.cm.get_cmap('afmhot'))
+ CS = ax1.scatter(friction, shearstrainrate_nonzero,
+ c=shearstrain_nonzero, linewidth=0.1,
+ cmap=matplotlib.cm.get_cmap('afmhot'))
ax1.set_yscale('log')
- x_min=numpy.floor(numpy.min(friction))
- x_max=numpy.ceil(numpy.max(friction))
+ x_min = numpy.floor(numpy.min(friction))
+ x_max = numpy.ceil(numpy.max(friction))
ax1.set_xlim([x_min, x_max])
- #y_min=numpy.min(v_nonzero)*0.5
- #y_max=numpy.max(v_nonzero)*2.0
- y_min=numpy.min(shearstrainrate_nonzero)*0.5
- y_max=numpy.max(shearstrainrate_nonzero)*2.0
+ y_min = numpy.min(shearstrainrate_nonzero)*0.5
+ y_max = numpy.max(shearstrainrate_nonzero)*2.0
ax1.set_ylim([y_min, y_max])
- cb=plt.colorbar(CS)
+ cb = plt.colorbar(CS)
cb.set_label('Shear strain $\\gamma$ [-]')
- #ax1.set_xlabel('Effective normal stress [kPa]')
ax1.set_xlabel('Friction $\\tau/N$ [-]')
- #ax1.set_ylabel('Shear velocity [m/s]')
ax1.set_ylabel('Shear strain rate $\\dot{\\gamma}$ [s$^{-1}$]')
- '''
- ax2=plt.subplot(212)
- ax2.plot(tau/N, v, '.')
- ax1.set_xlabel('Friction $\\tau/N$ [-]')
- ax1.set_ylabel('Shear velocity [m/s]')
- '''
-
elif method == 'inertia':
- t=numpy.zeros(sb.status())
- I=numpy.zeros(sb.status())
+ t = numpy.zeros(sb.status())
+ I = numpy.zeros(sb.status())
for i in numpy.arange(firststep, sb.status()):
sb.readstep(i, verbose=False)
- t[i]=sb.currentTime()
- I[i]=sb.inertiaParameterPlanarShear()
+ t[i] = sb.currentTime()
+ I[i] = sb.inertiaParameterPlanarShear()
# Plotting
if outformat != 'txt':
t@@ -7204,7 +7068,7 @@ class sim:
ax1.set_xlim(xlim)
# linear plot of deviatoric stress
- ax1=plt.subplot2grid((1,1),(0,0))
+ ax1 = plt.subplot2grid((1, 1), (0, 0))
ax1.set_xlabel('Time $t$ [s]')
ax1.set_ylabel('Inertia parameter $I$ [-]')
ax1.semilogy(t, I)
t@@ -7219,11 +7083,11 @@ class sim:
# Allocate arrays on first run
if i == firststep:
- p_mean=numpy.zeros(lastfile+1, dtype=numpy.float64)
+ p_mean = numpy.zeros(lastfile+1, dtype=numpy.float64)
- p_mean[i]=numpy.mean(sb.p_f)
+ p_mean[i] = numpy.mean(sb.p_f)
- t=numpy.linspace(0.0, sb.time_current, lastfile+1)
+ t = numpy.linspace(0.0, sb.time_current, lastfile+1)
# Plotting
if outformat != 'txt':
t@@ -7232,7 +7096,7 @@ class sim:
ax1.set_xlim(xlim)
# linear plot of deviatoric stress
- ax1=plt.subplot2grid((1,1),(0,0))
+ ax1 = plt.subplot2grid((1, 1), (0, 0))
ax1.set_xlabel('Time $t$, [s]')
ax1.set_ylabel('Mean fluid pressure, $\\bar{p}_f$, [kPa]')
ax1.plot(t, p_mean/1000.0, '+-')
t@@ -7242,70 +7106,63 @@ class sim:
elif method == 'fluid-pressure':
if figsize:
- fig=plt.figure(figsize=figsize)
+ fig = plt.figure(figsize=figsize)
else:
- fig=plt.figure(figsize=(8,6))
+ fig = plt.figure(figsize=(8, 6))
sb.readfirst(verbose=False)
# cell midpoint cell positions
- zpos_c=numpy.zeros(sb.num[2])
- dz=sb.L[2]/sb.num[2]
+ zpos_c = numpy.zeros(sb.num[2])
+ dz = sb.L[2]/sb.num[2]
for i in numpy.arange(sb.num[2]):
- zpos_c[i]=i*dz + 0.5*dz
+ zpos_c[i] = i*dz + 0.5*dz
- shear_strain=numpy.zeros(sb.status())
- pres=numpy.zeros((sb.num[2], sb.status()))
+ shear_strain = numpy.zeros(sb.status())
+ pres = numpy.zeros((sb.num[2], sb.status()))
# Read pressure values from simulation binaries
for i in numpy.arange(firststep, sb.status()):
sb.readstep(i, verbose=False)
- pres[:,i]=numpy.average(numpy.average(sb.p_f, axis=0), axis=0)
- shear_strain[i]=sb.shearStrain()
- t=numpy.linspace(0.0, sb.time_current, lastfile+1)
+ pres[:, i] = numpy.average(numpy.average(sb.p_f, axis=0), axis=0)
+ shear_strain[i] = sb.shearStrain()
+ t = numpy.linspace(0.0, sb.time_current, lastfile+1)
# Plotting
if outformat != 'txt':
- ax=plt.subplot(1,1,1)
+ ax = plt.subplot(1, 1, 1)
pres /= 1000.0 # Pa to kPa
if xlim:
- sb.readstep(10,verbose=False)
- gamma_per_i=sb.shearStrain()/10.0
- i_min=int(xlim[0]/gamma_per_i)
- i_max=int(xlim[1]/gamma_per_i)
- pres=pres[:,i_min:i_max]
+ sb.readstep(10, verbose=False)
+ gamma_per_i = sb.shearStrain()/10.0
+ i_min = int(xlim[0]/gamma_per_i)
+ i_max = int(xlim[1]/gamma_per_i)
+ pres = pres[:, i_min:i_max]
else:
- i_min=0
- i_max=sb.status()
+ i_min = 0
+ i_max = sb.status()
# use largest difference in p from 0 as +/- limit on colormap
#print i_min, i_max
- #p_ext=numpy.max(numpy.abs(pres))
- p_ext=numpy.max(numpy.abs(pres[0:9,:])) # for article2
+ p_ext = numpy.max(numpy.abs(pres))
if sb.wmode[0] == 3:
- x=t
+ x = t
else:
- x=shear_strain
+ x = shear_strain
if xlim:
- x=x[i_min:i_max]
+ x = x[i_min:i_max]
if cmap:
- im1=ax.pcolormesh(
- x, zpos_c, pres,
- #cmap=matplotlib.cm.get_cmap('bwr'),
- cmap=cmap,
- vmin=-p_ext, vmax=p_ext,
- rasterized=True)
+ im1 = ax.pcolormesh(x, zpos_c, pres, cmap=cmap,
+ vmin=-p_ext, vmax=p_ext,
+ rasterized=True)
else:
- im1=ax.pcolormesh(
- x, zpos_c, pres,
- #cmap=matplotlib.cm.get_cmap('bwr'),
- cmap=matplotlib.cm.get_cmap('RdBu_r'),
- #cmap=matplotlib.cm.get_cmap('coolwarm'),
- vmin=-p_ext, vmax=p_ext,
- rasterized=True)
+ im1 = ax.pcolormesh(x, zpos_c, pres,
+ cmap=matplotlib.cm.get_cmap('RdBu_r'),
+ vmin=-p_ext, vmax=p_ext,
+ rasterized=True)
ax.set_xlim([0, numpy.max(x)])
if sb.w_x[0] < sb.L[2]:
ax.set_ylim([zpos_c[0], sb.w_x[0]])
t@@ -7315,18 +7172,15 @@ class sim:
ax.set_xlabel('Time $t$ [s]')
else:
ax.set_xlabel('Shear strain $\\gamma$ [-]')
- #ax.set_xlabel('Time $t$ [s]')
ax.set_ylabel('Vertical position $z$ [m]')
- #ax.set_title(sb.id())
-
if xlim:
ax.set_xlim([x[0], x[-1]])
# for article2
ax.set_ylim([zpos_c[0], zpos_c[9]])
- cb=plt.colorbar(im1)
+ cb = plt.colorbar(im1)
cb.set_label('$p_\\text{f}$ [kPa]')
cb.solids.set_rasterized(True)
plt.tight_layout()
t@@ -7336,53 +7190,48 @@ class sim:
sb.readfirst(verbose=False)
if not sb.fluid:
raise Exception('Porosities can only be visualized in wet ' +
- 'simulations')
+ 'simulations')
- wall0_iz=int(sb.w_x[0]/(sb.L[2]/sb.num[2]))
+ wall0_iz = int(sb.w_x[0]/(sb.L[2]/sb.num[2]))
# cell midpoint cell positions
- zpos_c=numpy.zeros(sb.num[2])
- dz=sb.L[2]/sb.num[2]
+ zpos_c = numpy.zeros(sb.num[2])
+ dz = sb.L[2]/sb.num[2]
for i in numpy.arange(firststep, sb.num[2]):
- zpos_c[i]=i*dz + 0.5*dz
+ zpos_c[i] = i*dz + 0.5*dz
- shear_strain=numpy.zeros(sb.status())
- poros=numpy.zeros((sb.num[2], sb.status()))
+ shear_strain = numpy.zeros(sb.status())
+ poros = numpy.zeros((sb.num[2], sb.status()))
# Read pressure values from simulation binaries
for i in numpy.arange(firststep, sb.status()):
sb.readstep(i, verbose=False)
- poros[:,i]=numpy.average(numpy.average(sb.phi, axis=0),axis=0)
- shear_strain[i]=sb.shearStrain()
- t=numpy.linspace(0.0, sb.time_current, lastfile+1)
+ poros[:, i] = numpy.average(numpy.average(sb.phi, axis=0), axis=0)
+ shear_strain[i] = sb.shearStrain()
+ t = numpy.linspace(0.0, sb.time_current, lastfile+1)
# Plotting
if outformat != 'txt':
- ax=plt.subplot(1,1,1)
+ ax = plt.subplot(1, 1, 1)
- poros_max=numpy.max(poros[0:wall0_iz-1,1:])
- poros_min=numpy.min(poros)
+ poros_max = numpy.max(poros[0:wall0_iz-1, 1:])
+ poros_min = numpy.min(poros)
if sb.wmode[0] == 3:
- x=t
+ x = t
else:
- x=shear_strain
+ x = shear_strain
if cmap:
- im1=ax.pcolormesh(
- x, zpos_c, poros,
- cmap=cmap,
- vmin=poros_min, vmax=poros_max,
- rasterized=True)
+ im1 = ax.pcolormesh(x, zpos_c, poros,
+ cmap=cmap,
+ vmin=poros_min, vmax=poros_max,
+ rasterized=True)
else:
- im1=ax.pcolormesh(
- x, zpos_c, poros,
- cmap=matplotlib.cm.get_cmap('Blues_r'),
- #cmap=matplotlib.cm.get_cmap('bwr'),
- #cmap=matplotlib.cm.get_cmap('coolwarm'),
- #vmin=-p_ext, vmax=p_ext,
- vmin=poros_min, vmax=poros_max,
- rasterized=True)
+ im1 = ax.pcolormesh(x, zpos_c, poros,
+ cmap=matplotlib.cm.get_cmap('Blues_r'),
+ vmin=poros_min, vmax=poros_max,
+ rasterized=True)
ax.set_xlim([0, numpy.max(x)])
if sb.w_x[0] < sb.L[2]:
ax.set_ylim([zpos_c[0], sb.w_x[0]])
t@@ -7392,15 +7241,12 @@ class sim:
ax.set_xlabel('Time $t$ [s]')
else:
ax.set_xlabel('Shear strain $\\gamma$ [-]')
- #ax.set_xlabel('Time $t$ [s]')
ax.set_ylabel('Vertical position $z$ [m]')
if xlim:
ax.set_xlim(xlim)
- #ax.set_title(sb.id())
-
- cb=plt.colorbar(im1)
+ cb = plt.colorbar(im1)
cb.set_label('Mean horizontal porosity $\\bar{\phi}$ [-]')
cb.solids.set_rasterized(True)
plt.tight_layout()
t@@ -7409,31 +7255,30 @@ class sim:
elif method == 'contacts':
for i in numpy.arange(sb.status()+1):
- fn="../output/{0}.output{1:0=5}.bin".format(self.sid, i)
- sb.sid=self.sid + ".{:0=5}".format(i)
+ fn = "../output/{0}.output{1:0=5}.bin".format(self.sid, i)
+ sb.sid = self.sid + ".{:0=5}".format(i)
sb.readbin(fn, verbose=True)
if f_min and f_max:
sb.plotContacts(lower_limit=0.25, upper_limit=0.75,
- outfolder='../img_out/',
- f_min=f_min, f_max=f_max,
- title="t={:.2f} s, $N$={:.0f} kPa".format(
- sb.currentTime(),
- sb.currentNormalStress('defined')/1000.)
- )
+ outfolder='../img_out/',
+ f_min=f_min, f_max=f_max,
+ title="t={:.2f} s, $N$={:.0f} kPa"
+ .format(sb.currentTime(),
+ sb.currentNormalStress('defined')
+ /1000.))
else:
sb.plotContacts(lower_limit=0.25, upper_limit=0.75,
- title="t={:.2f} s, $N$={:.0f} kPa".format(
- sb.currentTime(),
- sb.currentNormalStress('defined')/1000.),
- outfolder='../img_out/')
+ title="t={:.2f} s, $N$={:.0f} kPa"
+ .format(sb.currentTime(),
+ sb.currentNormalStress('defined')
+ /1000.), outfolder='../img_out/')
# render images to movie
subprocess.call('cd ../img_out/ && ' +
- 'ffmpeg -sameq -i {}.%05d-contacts.png '.format(self.sid) +
- '{}-contacts.mp4'.format(self.sid),
- #'convert -quality 100 {}.*.png {}-contacts.avi'.format(
- #self.sid, self.sid),
- shell=True)
+ 'ffmpeg -sameq -i {}.%05d-contacts.png '
+ .format(self.sid) +
+ '{}-contacts.mp4'.format(self.sid),
+ shell=True)
else:
print("Visualization type '" + method + "' not understood")
t@@ -7441,10 +7286,10 @@ class sim:
# Optional save of figure content
if xlim:
- filename='{0}-{1}-{3}.{2}'.format(self.sid, method, outformat,
- xlim[-1])
+ filename = '{0}-{1}-{3}.{2}'.format(self.sid, method, outformat,
+ xlim[-1])
else:
- filename='{0}-{1}.{2}'.format(self.sid, method, outformat)
+ filename = '{0}-{1}.{2}'.format(self.sid, method, outformat)
if pickle:
pl.dump(fig, file(filename + '.pickle', 'w'))
t@@ -7472,8 +7317,8 @@ def convert(graphics_format='png', folder='../img_out', remove_ppm=False):
:type remove_ppm: bool
'''
- #quiet=' > /dev/null'
- quiet=''
+ #quiet = ' > /dev/null'
+ quiet = ''
# Convert images
subprocess.call('for F in ' + folder \
+ '/*.ppm ; do BASE=`basename $F .ppm`; convert $F ' \
t@@ -7507,15 +7352,15 @@ def render(binary, method='pres', max_val=1e3, lower_cutoff=0.0,
:param verbose: Show verbose information during ray tracing
:type verbose: bool
'''
- quiet=''
+ quiet = ''
if not verbose:
- quiet='-q'
+ quiet = '-q'
# Render images using sphere raytracer
if method == 'normal':
subprocess.call('cd .. ; ./sphere ' + quiet + \
' --render ' + binary, shell=True)
- else :
+ else:
subprocess.call('cd .. ; ./sphere ' + quiet + \
' --method ' + method + ' {}'.format(max_val) + \
' -l {}'.format(lower_cutoff) + \
t@@ -7556,15 +7401,15 @@ def video(project, out_folder='./', video_format='mp4',
# quiet, panic, fatal, error, warning, info, verbose, debug
loglevel = 'info'
if not verbose:
- loglevel='error'
+ loglevel = 'error'
outfile = out_folder + '/' + project + '.' + video_format
subprocess.call('ffmpeg -loglevel ' + loglevel + ' '
- + '-i ' + graphics_folder + project + '.output%05d.'
- + graphics_format
- + ' -c:v libx264 -profile:v high -pix_fmt yuv420p -g 30'
- + ' -r {} -y '.format(fps)
- + outfile, shell=True)
+ + '-i ' + graphics_folder + project + '.output%05d.'
+ + graphics_format
+ + ' -c:v libx264 -profile:v high -pix_fmt yuv420p -g 30'
+ + ' -r {} -y '.format(fps)
+ + outfile, shell=True)
if verbose:
print('saved to ' + outfile)
t@@ -7594,11 +7439,11 @@ def thinsectionVideo(project, out_folder="./", video_format="mp4", fps=25,
'''
# Render thin section images (png)
- lastfile=status(project)
- sb=sim(fluid=self.fluid)
+ lastfile = status(project)
+ sb = sim(fluid=False)
for i in range(lastfile+1):
- fn="../output/{0}.output{1:0=5}.bin".format(project, i)
- sb.sid=project + ".output{:0=5}".format(i)
+ fn = "../output/{0}.output{1:0=5}.bin".format(project, i)
+ sb.sid = project + ".output{:0=5}".format(i)
sb.readbin(fn, verbose=False)
sb.thinsection_x1x3(cbmax=sb.w_sigma0[0]*4.0)
t@@ -7607,15 +7452,15 @@ def thinsectionVideo(project, out_folder="./", video_format="mp4", fps=25,
# quiet, panic, fatal, error, warning, info, verbose, debug
loglevel = "info"
if not verbose:
- loglevel="error"
+ loglevel = "error"
subprocess.call("ffmpeg -qscale {0} -r {1} -b {2} -y ".format(\
- qscale, fps, bitrate) \
- + "-loglevel " + loglevel + " " \
- + "-i ../img_out/" + project + ".output%05d-ts-x1x3.png " \
- + "-vf 'crop=((in_w/2)*2):((in_h/2)*2)' " \
- + out_folder + "/" + project + "-ts-x1x3." + video_format, \
- shell=True)
+ qscale, fps, bitrate)
+ + "-loglevel " + loglevel + " "
+ + "-i ../img_out/" + project + ".output%05d-ts-x1x3.png "
+ + "-vf 'crop=((in_w/2)*2):((in_h/2)*2)' " \
+ + out_folder + "/" + project + "-ts-x1x3." + video_format,
+ shell=True)
def run(binary, verbose=True, hideinputfile=False):
'''
t@@ -7629,12 +7474,12 @@ def run(binary, verbose=True, hideinputfile=False):
:type hideinputfile: bool
'''
- quiet=''
- stdout=''
+ quiet = ''
+ stdout = ''
if not verbose:
- quiet='-q'
+ quiet = '-q'
if hideinputfile:
- stdout=' > /dev/null'
+ stdout = ' > /dev/null'
subprocess.call('cd ..; ./sphere ' + quiet + ' ' + binary + ' ' + stdout, \
shell=True)
t@@ -7678,11 +7523,11 @@ def torqueScriptParallel3(obj1, obj2, obj3, email='adc@geo.au.dk',
See also :func:`torqueScript()`
'''
- filename=obj1.sid + '_' + obj2.sid + '_' + obj3.sid + '.sh'
+ filename = obj1.sid + '_' + obj2.sid + '_' + obj3.sid + '.sh'
- fh=None
- try :
- fh=open(filename, "w")
+ fh = None
+ try:
+ fh = open(filename, "w")
fh.write('#!/bin/sh\n')
fh.write('#PBS -N ' + obj1.sid + '_' + obj2.sid + '_' + obj3.sid + '\n')
t@@ -7714,7 +7559,7 @@ def torqueScriptParallel3(obj1, obj2, obj3, email='adc@geo.au.dk',
fh.write('echo "End at `date`"\n')
return filename
- finally :
+ finally:
if fh is not None:
fh.close()
t@@ -7730,32 +7575,32 @@ def status(project):
:return type: int
'''
- fh=None
- try :
- filepath="../output/{0}.status.dat".format(project)
- fh=open(filepath)
- data=fh.read()
+ fh = None
+ try:
+ filepath = "../output/{0}.status.dat".format(project)
+ fh = open(filepath)
+ data = fh.read()
return int(data.split()[2]) # Return last file number
- finally :
+ finally:
if fh is not None:
fh.close()
-def cleanup(sim):
+def cleanup(sb):
'''
Removes the input/output files and images belonging to the object simulation
ID from the ``input/``, ``output/`` and ``img_out/`` folders.
- :param spherebin: A sim object
- :type spherebin: sim
+ :param sb: A sphere.sim object
+ :type sb: sim
'''
- subprocess.call("rm -f ../input/" + sim.sid + ".bin", shell=True)
- subprocess.call("rm -f ../output/" + sim.sid + ".*.bin", shell=True)
- subprocess.call("rm -f ../img_out/" + sim.sid + ".*", shell=True)
- subprocess.call("rm -f ../output/" + sim.sid + ".status.dat", shell=True)
- subprocess.call("rm -f ../output/" + sim.sid + ".*.vtu", shell=True)
- subprocess.call("rm -f ../output/fluid-" + sim.sid + ".*.vti", shell=True)
- subprocess.call("rm -f ../output/" + sim.sid + "-conv.png", shell=True)
- subprocess.call("rm -f ../output/" + sim.sid + "-conv.log", shell=True)
+ subprocess.call("rm -f ../input/" + sb.sid + ".bin", shell=True)
+ subprocess.call("rm -f ../output/" + sb.sid + ".*.bin", shell=True)
+ subprocess.call("rm -f ../img_out/" + sb.sid + ".*", shell=True)
+ subprocess.call("rm -f ../output/" + sb.sid + ".status.dat", shell=True)
+ subprocess.call("rm -f ../output/" + sb.sid + ".*.vtu", shell=True)
+ subprocess.call("rm -f ../output/fluid-" + sb.sid + ".*.vti", shell=True)
+ subprocess.call("rm -f ../output/" + sb.sid + "-conv.png", shell=True)
+ subprocess.call("rm -f ../output/" + sb.sid + "-conv.log", shell=True)
def V_sphere(r):
'''
t@@ -7765,5 +7610,3 @@ def V_sphere(r):
:return type: float
'''
return 4.0/3.0 * math.pi * r**3.0
-
-# vim: tabstop=8 expandtab shiftwidth=4 softtabstop=4