tUse consistent format for function definitions and default argument values - 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 346ae8299917811199c587382cb5c683d2b37995
(DIR) parent 3b2239705ebb3dff84425256a5d90456ba28fe5a
(HTM) Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Mon, 2 Sep 2019 06:19:34 +0200
Use consistent format for function definitions and default argument values
Diffstat:
M python/sphere.py | 3186 +++++++++++++++----------------
1 file changed, 1527 insertions(+), 1659 deletions(-)
---
(DIR) diff --git a/python/sphere.py b/python/sphere.py
t@@ -8,32 +8,32 @@ 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
+ py_mpl=False
import subprocess
import pickle as pl
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@@ -43,17 +43,17 @@ class sim:
Contains functions for reading and writing binaries, as well as simulation
setup and data analysis. Most arrays are initialized to default values.
- :param np: The number of particles to allocate memory for (default = 1)
+ :param np: The number of particles to allocate memory for (default=1)
:type np: int
- :param nd: The number of spatial dimensions (default = 3). Note that 2D and
+ :param nd: The number of spatial dimensions (default=3). Note that 2D and
1D simulations currently are not possible.
:type nd: int
- :param nw: The number of dynamic walls (default = 1)
+ :param nw: The number of dynamic walls (default=1)
:type nw: int
- :param sid: The simulation id (default = 'unnamed'). The simulation files
+ :param sid: The simulation id (default='unnamed'). The simulation files
will be written with this base name.
:type sid: str
- :param fluid: Setup fluid simulation (default = False)
+ :param fluid: Setup fluid simulation (default=False)
:type fluid: bool
:param cfd_solver: Fluid solver to use if fluid == True. 0: Navier-Stokes
(default), 1: Darcy.
t@@ -69,286 +69,286 @@ class sim:
cfd_solver=0):
# 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)
+ 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),
+ 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),
+ 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),
+ self.force =numpy.zeros((self.np, self.nd),
dtype=numpy.float64)
# The angular position vectors [rad]
- self.angpos = numpy.zeros((self.np, self.nd),
+ 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),
+ self.angvel =numpy.zeros((self.np, self.nd),
dtype=numpy.float64)
# The torque vectors [N*m]
- self.torque = numpy.zeros((self.np, self.nd),
+ 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)
+ # nw=1: Uniaxial (also used for shear experiments)
+ # nw=2: Biaxial
+ # nw=5: Triaxial
+ 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),
+ 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),
+ 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),
+ 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.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]),
+ 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]),
+ 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]),
+ 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@@ -357,10 +357,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@@ -368,18 +368,18 @@ 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)
+ # 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)
# Hold pressures constant in fluid cell (0: True, 1: False)
- self.p_f_constant = numpy.zeros(
+ self.p_f_constant=numpy.zeros(
(self.num[0], self.num[1], self.num[2]),
dtype=numpy.int32)
t@@ -387,77 +387,77 @@ class sim:
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)
+ # 0.0=no smoothing.
+ 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)
+ # 1.0=no under-relaxation
+ 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),
+ self.f_d=numpy.zeros((self.np, self.nd),
dtype=numpy.float64)
- self.f_p = numpy.zeros((self.np, self.nd),
+ self.f_p=numpy.zeros((self.np, self.nd),
dtype=numpy.float64)
- self.f_v = numpy.zeros((self.np, self.nd),
+ self.f_v=numpy.zeros((self.np, self.nd),
dtype=numpy.float64)
- self.f_sum = numpy.zeros((self.np, self.nd),
+ 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),
+ 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
+ # Fluid bulk modulus=1/self.beta_f
+ 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]) + ')')
# 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@@ -812,7 +812,7 @@ class sim:
if sid == '':
return self.sid
else:
- self.sid = sid
+ self.sid=sid
def idAppend(self, string):
'''
t@@ -827,19 +827,19 @@ class sim:
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):
+ 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@@ -848,52 +848,52 @@ 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)
+ :param es_dot: The particle shear energy loss rate (default=0)
:type es_dot: float
- :param es: The particle shear energy loss (default = 0)
+ :param es: The particle shear energy loss (default=0)
:type es: float
- :param ev_dot: The particle viscous energy rate loss (default = 0)
+ :param ev_dot: The particle viscous energy rate loss (default=0)
:type ev_dot: float
- :param ev: The particle viscous energy loss (default = 0)
+ :param ev: The particle viscous energy loss (default=0)
:type ev: float
- :param p: The particle pressure (default = 0)
+ :param p: The particle pressure (default=0)
: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=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@@ -920,56 +920,56 @@ class sim:
'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)
+ 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)
-
- def readbin(self, targetbin, verbose = True, bonds = True, sigma0mod = True,
- esysparticle = False):
+ 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):
'''
Reads a target ``sphere`` binary file.
t@@ -978,36 +978,36 @@ class sim:
:param targetbin: The path to the binary ``sphere`` file
:type targetbin: str
- :param verbose: Show diagnostic information (default = True)
+ :param verbose: Show diagnostic information (default=True)
:type verbose: bool
- :param bonds: The input file contains bond information (default = True).
+ :param bonds: The input file contains bond information (default=True).
This parameter should be true for all recent ``sphere`` versions.
:type bonds: bool
:param sigma0mod: The input file contains information about modulating
- stresses at the top wall (default = True). This parameter should be
+ stresses at the top wall (default=True). This parameter should be
true for all recent ``sphere`` versions.
:type sigma0mod: bool
:param esysparticle: Stop reading the file after reading the kinematics,
which is useful for reading output files from other DEM programs.
- (default = False)
+ (default=False)
: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 = \
+ self.time_dt=\
numpy.fromfile(fh, dtype=numpy.float64, count=1)
self.time_current =\
numpy.fromfile(fh, dtype=numpy.float64, count=1)
t@@ -1019,27 +1019,27 @@ class sim:
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):
t@@ -1049,10 +1049,10 @@ class sim:
numpy.fromfile(fh, dtype=numpy.float64, count=1)
if self.version >= 1.03:
- self.xyzsum = numpy.fromfile(fh, dtype=numpy.float64,\
+ 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,\
+ self.xyzsum=numpy.fromfile(fh, dtype=numpy.float64,\
count=self.np*2).reshape(self.np,2)
for i in numpy.arange(self.np):
t@@ -1061,147 +1061,147 @@ class sim:
self.fixvel[i] =\
numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.force = numpy.fromfile(fh, dtype=numpy.float64,\
+ 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,\
+ 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,\
+ 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,\
+ 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)\
+ 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.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.wmode =numpy.fromfile(fh, dtype=numpy.int32, count=self.nw)
for i in numpy.arange(self.nw):
self.w_n[i,:] =\
numpy.fromfile(fh, dtype=numpy.float64, count=self.nd)
- 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_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.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,
+ self.bonds[i,0]=numpy.fromfile(fh, dtype=numpy.uint32,
count=1)
- self.bonds[i,1] = numpy.fromfile(fh, dtype=numpy.uint32,
+ self.bonds[i,1]=numpy.fromfile(fh, dtype=numpy.uint32,
count=1)
- self.bonds_delta_n = numpy.fromfile(fh, dtype=numpy.float64,
+ self.bonds_delta_n=numpy.fromfile(fh, dtype=numpy.float64,
count=self.nb0)
- self.bonds_delta_t = numpy.fromfile(fh, dtype=numpy.float64,
+ 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,
+ self.bonds_omega_n=numpy.fromfile(fh, dtype=numpy.float64,
count=self.nb0)
- self.bonds_omega_t = numpy.fromfile(fh, dtype=numpy.float64,
+ 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,
+ self.cfd_solver=numpy.fromfile(fh, dtype=numpy.int32,
count=1)
else:
- self.cfd_solver = numpy.zeros(1, dtype=numpy.int32)
+ self.cfd_solver=numpy.zeros(1, dtype=numpy.int32)
- self.mu = numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ self.mu=numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.v_f = numpy.empty(
+ self.v_f=numpy.empty(
(self.num[0], self.num[1], self.num[2], self.nd),
dtype=numpy.float64)
- self.p_f = \
+ self.p_f=\
numpy.empty((self.num[0],self.num[1],self.num[2]),
dtype=numpy.float64)
- self.phi = \
+ self.phi=\
numpy.empty((self.num[0],self.num[1],self.num[2]),
dtype=numpy.float64)
- self.dphi = \
+ 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] = \
+ self.v_f[x,y,z,0]=\
numpy.fromfile(fh, dtype=numpy.float64,\
count=1)
- self.v_f[x,y,z,1] = \
+ self.v_f[x,y,z,1]=\
numpy.fromfile(fh, dtype=numpy.float64,\
count=1)
- self.v_f[x,y,z,2] = \
+ self.v_f[x,y,z,2]=\
numpy.fromfile(fh, dtype=numpy.float64,\
count=1)
- self.p_f[x,y,z] = \
+ self.p_f[x,y,z]=\
numpy.fromfile(fh, dtype=numpy.float64,\
count=1)
- self.phi[x,y,z] = \
+ self.phi[x,y,z]=\
numpy.fromfile(fh, dtype=numpy.float64,\
count=1)
- self.dphi[x,y,z] = \
+ self.dphi[x,y,z]=\
numpy.fromfile(fh, dtype=numpy.float64,\
count=1)/(self.time_dt*self.ndem)
t@@ -1239,45 +1239,45 @@ class sim:
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 = \
+ 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] = \
+ self.p_f_constant[x,y,z]=\
numpy.fromfile(fh,
dtype=numpy.int32,
count=1)
else:
- self.p_f_constant = numpy.zeros(
+ 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 = \
+ self.gamma=\
numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.theta = \
+ self.theta=\
numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.beta = \
+ self.beta =\
numpy.fromfile(fh, dtype=numpy.float64, count=1)
self.tolerance =\
numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.maxiter = \
+ self.maxiter=\
numpy.fromfile(fh, dtype=numpy.uint32, count=1)
if self.version >= 1.01:
- self.ndem = \
+ 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 = \
+ self.c_phi=\
numpy.fromfile(fh, dtype=numpy.float64, count=1)
self.c_v =\
numpy.fromfile(fh, dtype=numpy.float64, count=1)
t@@ -1290,75 +1290,75 @@ class sim:
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,:] = \
+ self.f_d[i,:]=\
numpy.fromfile(fh, dtype=numpy.float64,
count=self.nd)
for i in numpy.arange(self.np):
- self.f_p[i,:] = \
+ self.f_p[i,:]=\
numpy.fromfile(fh, dtype=numpy.float64,
count=self.nd)
for i in numpy.arange(self.np):
- self.f_v[i,:] = \
+ self.f_v[i,:]=\
numpy.fromfile(fh, dtype=numpy.float64,
count=self.nd)
for i in numpy.arange(self.np):
- self.f_sum[i,:] = \
+ self.f_sum[i,:]=\
numpy.fromfile(fh, dtype=numpy.float64,
count=self.nd)
else:
- self.f_d = numpy.zeros((self.np, self.nd),
+ self.f_d=numpy.zeros((self.np, self.nd),
dtype=numpy.float64)
- self.f_p = numpy.zeros((self.np, self.nd),
+ self.f_p=numpy.zeros((self.np, self.nd),
dtype=numpy.float64)
- self.f_v = numpy.zeros((self.np, self.nd),
+ self.f_v=numpy.zeros((self.np, self.nd),
dtype=numpy.float64)
- self.f_sum = numpy.zeros((self.np, self.nd),
+ self.f_sum=numpy.zeros((self.np, self.nd),
dtype=numpy.float64)
elif self.version >= 2.0 and self.cfd_solver == 1:
- self.tolerance = \
+ self.tolerance=\
numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.maxiter = \
+ self.maxiter=\
numpy.fromfile(fh, dtype=numpy.uint32, count=1)
- self.ndem = numpy.fromfile(fh, dtype=numpy.uint32, count=1)
- self.c_phi = \
+ 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.f_p=numpy.empty_like(self.x)
for i in numpy.arange(self.np):
- self.f_p[i,:] = \
+ self.f_p[i,:]=\
numpy.fromfile(fh, dtype=numpy.float64,
count=self.nd)
- self.beta_f = \
+ self.beta_f=\
numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.k_c = \
+ 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)
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()
- def writebin(self, folder = "../input/", verbose = True):
+ def writebin(self, folder="../input/", verbose=True):
'''
Writes a ``sphere`` binary file to the ``../input/`` folder by default.
The file name will be in the format ``<self.sid>.bin``.
t@@ -1367,16 +1367,16 @@ class sim:
:param folder: The folder where to place the output binary file
:type folder: str
- :param verbose: Show diagnostic information (default = True)
+ :param verbose: Show diagnostic information (default=True)
:type verbose: bool
'''
- fh = None
+ fh=None
try :
- targetbin = folder + "/" + self.sid + ".bin"
+ 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@@ -1558,7 +1558,7 @@ class sim:
if fh is not None:
fh.close()
- def writeVTKall(self, cell_centered = True, verbose = True, forces = False):
+ def writeVTKall(self, cell_centered=True, verbose=True, forces=False):
'''
Writes a VTK file for each simulation output file with particle
information and the fluid grid to the ``../output/`` folder by default.
t@@ -1605,20 +1605,20 @@ class sim:
the :func:`writeVTKall()` function), it is able to step the
visualization through time by using the ParaView controls.
- :param verbose: Show diagnostic information (default = True)
+ :param verbose: Show diagnostic information (default=True)
:type verbose: bool
:param cell_centered: Write fluid values to cell centered positions
- (default = true)
+ (default=true)
:type cell_centered: bool
- :param forces: Write contact force files (slow) (default = False)
+ :param forces: Write contact force files (slow) (default=False)
:type forces: bool
'''
- lastfile = status(self.sid)
- sb = sim(fluid = self.fluid)
+ lastfile=status(self.sid)
+ sb=sim(fluid=self.fluid)
for i in range(lastfile+1):
- fn = "../output/{0}.output{1:0=5}.bin".format(self.sid, i)
- sb.sid = self.sid + ".{:0=5}".format(i)
- sb.readbin(fn, verbose = False)
+ fn="../output/{0}.output{1:0=5}.bin".format(self.sid, i)
+ sb.sid=self.sid + ".{:0=5}".format(i)
+ sb.readbin(fn, verbose=False)
if sb.np > 0:
if i == 0:
sb.writeVTK(verbose=verbose)
t@@ -1649,7 +1649,7 @@ class sim:
else:
sb.writeFluidVTK(verbose=False, cell_centered=cell_centered)
- def writeVTK(self, folder = '../output/', verbose = True):
+ def writeVTK(self, folder='../output/', verbose=True):
'''
Writes a VTK file with particle information to the ``../output/`` folder
by default. The file name will be in the format ``<self.sid>.vtu``.
t@@ -1676,19 +1676,19 @@ class sim:
visualization through time by using the ParaView controls.
:param folder: The folder where to place the output binary file (default
- (default = '../output/')
+ (default='../output/')
:type folder: str
- :param verbose: Show diagnostic information (default = True)
+ :param verbose: Show diagnostic information (default=True)
:type verbose: bool
'''
- fh = None
+ fh=None
try :
- targetbin = folder + '/' + self.sid + '.vtu' # unstructured grid
+ 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
t@@ -1920,7 +1920,7 @@ class sim:
if fh is not None:
fh.close()
- def writeVTKforces(self, folder = '../output/', verbose = True):
+ def writeVTKforces(self, folder='../output/', verbose=True):
'''
Writes a VTK file with particle-interaction information to the
``../output/`` folder by default. The file name will be in the format
t@@ -1930,9 +1930,9 @@ class sim:
filter.
:param folder: The folder where to place the output file (default
- (default = '../output/')
+ (default='../output/')
:type folder: str
- :param verbose: Show diagnostic information (default = True)
+ :param verbose: Show diagnostic information (default=True)
:type verbose: bool
'''
t@@ -1940,28 +1940,28 @@ 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)
t@@ -1969,7 +1969,7 @@ class sim:
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()
+ 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)
t@@ -1978,7 +1978,7 @@ class sim:
stresses.SetValue(i, self.sigma_contacts[i])
# initalize VTK data structure
- polydata = vtk.vtkPolyData()
+ polydata=vtk.vtkPolyData()
polydata.SetPoints(points)
polydata.SetLines(lines)
t@@ -1991,7 +1991,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@@ -2003,8 +2003,8 @@ class sim:
print('Output file: ' + filename)
- def writeFluidVTK(self, folder = '../output/', cell_centered = True,
- verbose = True):
+ def writeFluidVTK(self, folder='../output/', cell_centered=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@@ -2012,11 +2012,11 @@ class sim:
The scalars (pressure, porosity, porosity change) and the velocity
vectors are either placed in a grid where the grid corners correspond to
- the computational grid center (cell_centered = False). This results in a
+ the computational grid center (cell_centered=False). This results in a
grid that doesn't appears to span the simulation domain, and values are
smoothly interpolated on the cell faces. Alternatively, the
visualization grid is equal to the computational grid, and cells face
- colors are not interpolated (cell_centered = True, default behavior).
+ colors are not interpolated (cell_centered=True, default behavior).
The fluid grid is visualized by opening the vti files, and pressing
"Apply" to import all fluid field properties. To visualize the scalar
t@@ -2048,23 +2048,23 @@ class sim:
visualization through time by using the ParaView controls.
:param folder: The folder where to place the output binary file (default
- (default = '../output/')
+ (default='../output/')
:type folder: str
:param cell_centered: put scalars and vectors at cell centers (True) or
- cell corners (False), (default = True)
+ cell corners (False), (default=True)
:type cell_centered: bool
- :param verbose: Show diagnostic information (default = True)
+ :param verbose: Show diagnostic information (default=True)
:type verbose: bool
'''
if not py_vtk:
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@@ -2076,7 +2076,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@@ -2085,7 +2085,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@@ -2094,7 +2094,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@@ -2103,7 +2103,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@@ -2113,7 +2113,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@@ -2124,7 +2124,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@@ -2133,7 +2133,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@@ -2141,7 +2141,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@@ -2153,7 +2153,7 @@ 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;
+ 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])
t@@ -2187,7 +2187,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@@ -2212,19 +2212,19 @@ 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)
t@@ -2239,27 +2239,27 @@ class sim:
for i in numpy.arange(self.np):
# create source
- source = vtk.vtkSphereSource()
+ 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)
+ 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
t@@ -2277,14 +2277,14 @@ class sim:
Read the first output file from the ``../output/`` folder, corresponding
to the object simulation id (``self.sid``).
- :param verbose: Display diagnostic information (default = True)
+ :param verbose: Display diagnostic information (default=True)
:type verbose: bool
See also :func:`readbin()`, :func:`readlast()`, :func:`readsecond`, and
:func:`readstep`.
'''
- fn = '../output/' + self.sid + '.output00000.bin'
+ fn='../output/' + self.sid + '.output00000.bin'
self.readbin(fn, verbose)
def readsecond(self, verbose=True):
t@@ -2292,13 +2292,13 @@ class sim:
Read the second output file from the ``../output/`` folder,
corresponding to the object simulation id (``self.sid``).
- :param verbose: Display diagnostic information (default = True)
+ :param verbose: Display diagnostic information (default=True)
:type verbose: bool
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@@ -2308,13 +2308,13 @@ class sim:
:param step: The output file number to read, starting from 0.
:type step: int
- :param verbose: Display diagnostic information (default = True)
+ :param verbose: Display diagnostic information (default=True)
:type verbose: bool
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@@ -2322,14 +2322,14 @@ class sim:
Read the last output file from the ``../output/`` folder, corresponding
to the object simulation id (``self.sid``).
- :param verbose: Display diagnostic information (default = True)
+ :param verbose: Display diagnostic information (default=True)
:type verbose: bool
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@@ -2345,26 +2345,24 @@ 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))
- 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,
- histogram = False):
+ def generateRadii(self, psd='logn', mean=440e-6, variance=8.8e-9,
+ histogram=False):
'''
Draw random particle radii from a selected probability distribution.
The larger the variance of radii is, the slower the computations will
t@@ -2380,7 +2378,7 @@ class sim:
value is ``uni``, which is a uniform distribution from
``mean - variance`` to ``mean + variance``.
:type psd: str
- :param mean: The mean radius [m] (default = 440e-6 m)
+ :param mean: The mean radius [m] (default=440e-6 m)
:type mean: float
:param variance: The variance in the probability distribution
[m].
t@@ -2390,26 +2388,26 @@ class sim:
'''
if psd == 'logn': # Log-normal probability distribution
- mu = math.log(\
+ 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)
+ 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\'')
# Show radii as histogram
if histogram and py_mpl:
- fig = plt.figure(figsize=(8,8))
- figtitle = 'Particle size distribution, {0} particles'.format(\
+ 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
+ bins=20
# Create histogram
plt.hist(self.radius, bins)
t@@ -2421,11 +2419,8 @@ class sim:
fig.savefig(self.sid + '-psd.png')
fig.clf()
- def generateBimodalRadii(self,
- r_small = 0.005,
- r_large = 0.05,
- ratio = 0.2,
- verbose = True):
+ def generateBimodalRadii(self, r_small=0.005, r_large=0.05, ratio=0.2,
+ verbose=True):
'''
Draw random radii from two distinct sizes.
t@@ -2442,17 +2437,17 @@ 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")
t@@ -2461,7 +2456,7 @@ class sim:
+ " large particles, and " + str(self.np - nlarge)
+ " small")
- def checkerboardColors(self, nx = 6, ny = 6, nz = 6):
+ def checkerboardColors(self, nx=6, ny=6, nz=6):
'''
Assign checkerboard color values to the particles in an orthogonal grid.
t@@ -2472,17 +2467,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@@ -2492,10 +2487,10 @@ class sim:
the viscous-frictional contact model is significantly faster.
:param contactmodel: The type of tangential contact model to use
- (visco-frictional = 1, elasto-visco-frictional = 2)
+ (visco-frictional=1, elasto-visco-frictional=2)
:type contactmodel: int
'''
- self.contactmodel[0] = contactmodel
+ self.contactmodel[0]=contactmodel
def wall0iz(self):
'''
t@@ -2516,7 +2511,7 @@ class sim:
See also :func:`periodicBoundariesXY()` and
:func:`periodicBoundariesX()`
'''
- self.periodic[0] = 0
+ self.periodic[0]=0
def periodicBoundariesXY(self):
'''
t@@ -2525,7 +2520,7 @@ class sim:
See also :func:`normalBoundariesXY()` and
:func:`periodicBoundariesX()`
'''
- self.periodic[0] = 1
+ self.periodic[0]=1
def periodicBoundariesX(self):
'''
t@@ -2534,7 +2529,7 @@ class sim:
See also :func:`normalBoundariesXY()` and
:func:`periodicBoundariesXY()`
'''
- self.periodic[0] = 2
+ self.periodic[0]=2
def adaptiveGrid(self):
'''
t@@ -2545,7 +2540,7 @@ class sim:
See also :func:`staticGrid()`
'''
- self.adaptive[0] = 1
+ self.adaptive[0]=1
def staticGrid(self):
'''
t@@ -2553,9 +2548,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])):
'''
Initialize particle positions in completely random configuration. Radii
*must* be set beforehand. If the x and y boundaries are set as periodic,
t@@ -2564,7 +2559,7 @@ class sim:
make space for their radii within the bounding box.
:param gridnum: The number of sorting cells in each spatial direction
- (default = [12, 12, 36])
+ (default=[12, 12, 36])
:type gridnum: numpy.array
:param dx: The cell width in any direction. If the default value is used
(-1), the cell width is calculated to fit the largest particle.
t@@ -2572,37 +2567,37 @@ class sim:
'''
# Calculate cells in grid
- self.num = gridnum
+ self.num=gridnum
# 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)) \
+ 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)))
t@@ -2620,7 +2615,7 @@ class sim:
:param L: The upper boundary of the domain [m]
:type L: numpy.array
:param origo: The lower boundary of the domain [m]. Negative values
- won't work. Default = [0.0, 0.0, 0.0].
+ won't work. Default=[0.0, 0.0, 0.0].
:type origo: numpy.array
:param dx: The cell width in any direction. If the default value is used
(-1), the cell width is calculated to fit the largest particle.
t@@ -2629,25 +2624,25 @@ 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)
+ 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):
t@@ -2673,12 +2668,12 @@ 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 "
t@@ -2687,10 +2682,10 @@ class sim:
# 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):
+ def initGridAndWorldsize(self, margin=2.0):
'''
Initialize grid suitable for the particle positions set previously.
The margin parameter adjusts the distance (in no. of max. radii)
t@@ -2702,33 +2697,33 @@ 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[:]),
+ 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[:]),
+ 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)
+ 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])):
+ def initGridPos(self, gridnum=numpy.array([12, 12, 36])):
'''
Initialize particle positions in loose, cubic configuration.
``gridnum`` is the number of cells in the x, y and z directions.
t@@ -2740,19 +2735,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@@ -2769,13 +2764,13 @@ class sim:
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:
t@@ -2790,8 +2785,7 @@ class sim:
self.num[1] += 1
- def initRandomGridPos(self,
- gridnum = numpy.array([12, 12, 32]),
+ def initRandomGridPos(self, gridnum=numpy.array([12, 12, 32]),
padding=2.1):
'''
Initialize particle positions in loose, cubic configuration with some
t@@ -2810,46 +2804,46 @@ 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 \
+ 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@@ -2866,21 +2860,21 @@ 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]
t@@ -2899,8 +2893,8 @@ class sim:
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@@ -2922,11 +2916,11 @@ class sim:
:type spacing: float
'''
- bondparticles = numpy.unique(\
+ 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[:-1].copy()
bondparticles =\
bondparticles.reshape(int(bondparticles.size/2), 2).copy()
t@@ -2939,20 +2933,20 @@ 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)\
+ 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)\
+ 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)\
+ 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)\
+ 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):
+ def adjustUpperWall(self, z_adjust=1.1):
'''
Included for legacy purposes, calls :func:`adjustWall()` with ``idx=0``.
t@@ -2962,20 +2956,20 @@ 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=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_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)
+ 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, wall_mass_factor = 1.0):
+ def adjustWall(self, idx, adjust=1.1, wall_mass_factor=1.0):
'''
Adjust grid and dynamic wall to max. particle position. The wall
thickness will by standard equal the maximum particle diameter. The
t@@ -2984,8 +2978,8 @@ class sim:
The total wall mass can be linearly scaled by the `wall_mass_factor`
parameter.
- :param: idx: The wall to adjust. 0 = +z, upper wall (default), 1 = -x,
- left wall, 2 = +x, right wall, 3 = -y, front wall, 4 = +y, back
+ :param: idx: The wall to adjust. 0=+z, upper wall (default), 1=-x,
+ left wall, 2=+x, right wall, 3=-y, front wall, 4=+y, back
wall.
:type idx: int
:param z_adjust: Increase the world and grid size by this amount to
t@@ -2995,34 +2989,34 @@ 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)
+ 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
+ cellsize=self.L[0] / self.num[0]
+ d_max=numpy.max(self.radius)*2.0
- self.num[dim] = numpy.ceil(((xmax-xmin)*adjust + xmin)/cellsize)
- self.L[dim] = (xmax-xmin)*adjust + xmin
+ 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 \
+ 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_m[idx]=numpy.array([self.rho*self.L[0]*self.L[1]*d_max])
+ self.w_m[idx]=self.totalMass()
- def consolidate(self, normal_stress = 10e3):
+ def consolidate(self, normal_stress=10e3):
'''
Setup consolidation experiment. Specify the upper wall normal stress in
Pascal, default value is 10 kPa.
t@@ -3031,7 +3025,7 @@ 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 '
t@@ -3044,16 +3038,16 @@ 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):
+ def uniaxialStrainRate(self, wvel=-0.001):
'''
Setup consolidation experiment. Specify the upper wall velocity in m/s,
default value is -0.001 m/s (i.e. downwards).
t@@ -3068,10 +3062,10 @@ 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):
+ def triaxial(self, wvel=-0.001, normal_stress=10.0e3):
'''
Setup triaxial experiment. The upper wall is moved at a fixed velocity
in m/s, default values is -0.001 m/s (i.e. downwards). The side walls
t@@ -3088,19 +3082,19 @@ 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]),
+ 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.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)
- def shear(self, shear_strain_rate = 1.0, shear_stress = False):
+ def shear(self, shear_strain_rate=1.0, shear_stress=False):
'''
Setup shear experiment either by a constant shear rate or a constant
shear stress. The shear strain rate is the shear velocity divided by
t@@ -3115,60 +3109,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] <
+ 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
+ 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] >
+ 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
+ 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 shear_stress == False:
- 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@@ -3195,27 +3189,27 @@ 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(\
+ 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
t@@ -3230,10 +3224,10 @@ class sim:
# 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@@ -3243,10 +3237,10 @@ 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()
t@@ -3262,7 +3256,7 @@ class sim:
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
+ self.dx=self.L/self.num
def hydraulicConductivity(self, phi=0.35):
'''
t@@ -3276,13 +3270,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, phi_min=0.3, phi_max=0.99):
'''
Determine the hydraulic permeability (k) [m*m] from the Kozeny-Carman
relationship, using the permeability prefactor (`self.k_c`), and the
t@@ -3305,17 +3299,13 @@ class sim:
'''
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))
+ phi_bar=numpy.mean(self.phi)
+ self.D=self.K_c/(self.rho_f*g*(self.k_n[0] + phi_bar*K))
else:
raise Exception('This function only works for the Darcy solver')
- def initTemporal(self, total,
- current = 0.0,
- file_dt = 0.05,
- step_count = 0,
- dt = -1,
- epsilon = 0.01):
+ def initTemporal(self, total, current=0.0, file_dt=0.05, step_count=0,
+ dt=-1, epsilon=0.01):
'''
Set temporal parameters for the simulation. *Important*: Particle radii,
physical parameters, and the optional fluid grid need to be set prior to
t@@ -3327,42 +3317,42 @@ class sim:
:param total: The time at which to end the simulation [s]
:type total: float
- :param current: The current time [s] (default = 0.0 s)
+ :param current: The current time [s] (default=0.0 s)
:type total: float
- :param file_dt: The interval between output files [s] (default = 0.05 s)
+ :param file_dt: The interval between output files [s] (default=0.05 s)
:type total: float
- :step_count: The number of the first output file (default = 0)
+ :step_count: The number of the first output file (default=0)
:type step_count: int
:param dt: The computational time step length [s]
:type total: float
- :param epsilon: Time step multiplier (default = 0.01)
+ :param epsilon: Time step multiplier (default=0.01)
:type epsilon: float
'''
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.")
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 self.fluid == False:
raise Exception('Error: Could not automatically set a time step.')
t@@ -3371,15 +3361,15 @@ 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]=0
def dry(self):
'''
t@@ -3387,7 +3377,7 @@ class sim:
See also :func:`wet()`
'''
- self.fluid = False
+ self.fluid=False
def wet(self):
'''
t@@ -3395,11 +3385,11 @@ 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, cfd_solver = 0):
+ def initFluid(self, mu=8.9e-4, rho=1.0e3, p=0.0, hydrostatic=False,
+ cfd_solver=0):
'''
Initialize the fluid arrays and the fluid viscosity. The default value
of ``mu`` equals the dynamic viscosity of water at 25 degrees Celcius.
t@@ -3421,93 +3411,93 @@ 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]),
+ 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),
+ 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]),
+ 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]),
+ 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.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@@ -3526,7 +3516,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@@ -3538,7 +3528,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@@ -3548,7 +3538,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@@ -3558,7 +3548,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@@ -3570,8 +3560,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@@ -3581,7 +3571,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@@ -3591,7 +3581,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@@ -3601,7 +3591,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@@ -3613,8 +3603,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@@ -3626,8 +3616,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@@ -3639,8 +3629,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@@ -3651,8 +3641,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@@ -3663,8 +3653,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@@ -3675,8 +3665,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@@ -3687,8 +3677,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@@ -3718,67 +3708,56 @@ 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
- print('Hydraulic permeability limits for porosity phi = ' + \
+ 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')
- print('Hydraulic conductivity limits for porosity phi = ' + \
+ print('\tk=' + str(k) + ' m*m')
+ print('Hydraulic conductivity limits for porosity phi=' + \
str(phi) + ':')
- print('\tK = ' + str(K) + ' m/s')
+ 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):
'''
Calculates the hydrological permeabilities from the Kozeny-Carman
relationship. These values are only relevant when the Darcy solver is
- used (`self.cfd_solver = 1`). The permeability pre-factor `self.k_c`
+ used (`self.cfd_solver=1`). The permeability pre-factor `self.k_c`
and the assemblage porosities must be set beforehand. The former values
are set if a file from the `output/` folder is read using
`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):
'''
Calculates the hydrological conductivities from the Kozeny-Carman
relationship. These values are only relevant when the Darcy solver is
- used (`self.cfd_solver = 1`). The permeability pre-factor `self.k_c`
+ used (`self.cfd_solver=1`). The permeability pre-factor `self.k_c`
and the assemblage porosities must be set beforehand. The former values
are set if a file from the `output/` folder is read using
`self.readbin`.
'''
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)')
-
- 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,
- gamma_wn = 0.0,
- gamma_wt = 0.0,
- capillaryCohesion = 0):
+ + '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,
+ gamma_wn=0.0, gamma_wt=0.0, capillaryCohesion=0):
'''
Initialize particle parameters to default values.
t@@ -3807,76 +3786,76 @@ class sim:
:param gamma_wt: Wall-particle contact tangential viscosity [Ns/m]
:type gamma_wt: float
:param capillaryCohesion: Enable particle-particle capillary cohesion
- interaction model (0 = no (default), 1 = yes)
+ interaction model (0=no (default), 1=yes)
:type capillaryCohesion: int
'''
# 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
+ # Contact shear elastic stiffness (for contactmodel=2), N/m
+ 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
+ # Contact rolling elastic stiffness (for contactmodel=2), N/m
+ 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.
- #self.gamma_n = numpy.ones(self.np, dtype=numpy.float64) \
+ # Normal force component elastic if nu=0.0.
+ #self.gamma_n=numpy.ones(self.np, dtype=numpy.float64) \
# * nu_frac * 2.0 * math.sqrt(4.0/3.0 * math.pi \
# * numpy.amin(self.radius)**3 \
# * self.rho[0] * self.k_n[0])
- 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
+ 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]=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@@ -3886,7 +3865,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@@ -3896,7 +3875,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@@ -3911,7 +3890,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@@ -3926,29 +3905,29 @@ 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 = '
+ print('Info: The system is under-dampened (ratio='
+ str(damping_ratio)
- + ') in the normal component. \nCritical damping = '
+ + ') in the normal component. \nCritical damping='
+ str(critical_gamma) + '. This is ok.')
elif damping_ratio > 1.0:
if over_damping:
- print('Warning: The system is over-dampened (ratio = '
+ print('Warning: The system is over-dampened (ratio='
+ str(damping_ratio) + ') in the normal component. '
- + '\nCritical damping = ' + str(critical_gamma) + '.')
+ + '\nCritical damping=' + str(critical_gamma) + '.')
else:
- raise Exception('Warning: The system is over-dampened (ratio = '
+ 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 = '
+ + 'if this is what you want. \nCritical damping='
+ str(critical_gamma) + '.')
else:
- print('Warning: The system is critically dampened (ratio = '
+ print('Warning: The system is critically dampened (ratio='
+ str(damping_ratio) + ') in the normal component. '
- + '\nCritical damping = ' + str(critical_gamma) + '.')
+ + '\nCritical damping=' + str(critical_gamma) + '.')
def setDampingTangential(self, gamma, over_damping=False):
'''
t@@ -3963,23 +3942,23 @@ 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 = '
+ print('Info: The system is under-dampened (ratio='
+ str(damping_ratio)
+ ') in the tangential component. This is ok.')
elif damping_ratio > 1.0:
if over_damping:
- print('Warning: The system is over-dampened (ratio = '
+ print('Warning: The system is over-dampened (ratio='
+ str(damping_ratio) + ') in the tangential component.')
else:
- raise Exception('Warning: The system is over-dampened (ratio = '
+ 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.')
else:
- print('Warning: The system is critically dampened (ratio = '
+ print('Warning: The system is critically dampened (ratio='
+ str(damping_ratio) + ') in the tangential component.')
def setStaticFriction(self, mu_s):
t@@ -3994,7 +3973,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@@ -4011,7 +3990,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@@ -4026,10 +4005,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@@ -4043,7 +4022,7 @@ class sim:
See also: :func:`setFluidDensity()` and
:func:`setFluidCompressibility()`
'''
- self.mu[0] = mu
+ self.mu[0]=mu
def setFluidDensity(self, rho_f):
'''
t@@ -4056,7 +4035,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@@ -4088,37 +4067,37 @@ 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 hasattr(self, 'bonds') == False:
- self.bonds = numpy.array([[i,j]], dtype=numpy.uint32)
+ self.bonds=numpy.array([[i,j]], dtype=numpy.uint32)
else :
- self.bonds = numpy.vstack((self.bonds, [i,j]))
+ self.bonds=numpy.vstack((self.bonds, [i,j]))
if hasattr(self, 'bonds_delta_n') == False:
- self.bonds_delta_n = numpy.array([0.0], dtype=numpy.uint32)
+ 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.vstack((self.bonds_delta_n, [0.0]))
+ self.bonds_delta_n=numpy.append(self.bonds_delta_n, [0.0])
if hasattr(self, 'bonds_delta_t') == False:
- self.bonds_delta_t = numpy.array([[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,\
+ self.bonds_delta_t=numpy.vstack((self.bonds_delta_t,\
[0.0, 0.0, 0.0]))
if hasattr(self, 'bonds_omega_n') == False:
- self.bonds_omega_n = numpy.array([0.0], dtype=numpy.uint32)
+ 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.vstack((self.bonds_omega_n, [0.0]))
+ self.bonds_omega_n=numpy.append(self.bonds_omega_n, [0.0])
if hasattr(self, 'bonds_omega_t') == False:
- self.bonds_omega_t = numpy.array([[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,\
+ self.bonds_omega_t=numpy.vstack((self.bonds_omega_t,\
[0.0, 0.0, 0.0]))
# Increment the number of bonds with one
t@@ -4184,7 +4163,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@@ -4240,7 +4219,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@@ -4263,7 +4242,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@@ -4286,7 +4265,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@@ -4309,7 +4288,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@@ -4331,20 +4310,20 @@ class sim:
'''
if method == 'pot':
- m = numpy.ones(self.np)*4.0/3.0*math.pi*self.radius**3*self.rho
+ 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
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(\
t@@ -4365,19 +4344,19 @@ class sim:
elif method == 'bondpot':
if self.nb0 > 0:
- R_bar = self.lambda_bar*numpy.minimum(\
+ 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(\
+ 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 :
t@@ -4394,15 +4373,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@@ -4417,32 +4396,30 @@ 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):
+ def porosity(self, slices=10, verbose=False):
'''
Calculates the porosity as a function of depth, by averaging values in
horizontal slabs. Returns porosity values and their corresponding depth.
t@@ -4461,24 +4438,24 @@ class sim:
self.writebin(verbose=False)
# Run porosity program on binary
- pipe = subprocess.Popen(\
+ pipe=subprocess.Popen(\
["../porosity",\
"-s","{}".format(slices),\
"../input/" + self.sid + ".bin"],\
stdout=subprocess.PIPE)
- output, err = pipe.communicate()
+ 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@@ -4513,34 +4490,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 verbose == False:
- 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@@ -4552,15 +4529,11 @@ class sim:
'''
cleanup(self)
- def torqueScript(self,
- email='adc@geo.au.dk',
- email_alerts='ae',
- walltime='24:00:00',
- queue='qfermi',
- cudapath='/com/cuda/4.0.17/cuda',
- spheredir='/home/adc/code/sphere',
- use_workdir=False,
- workdir='/scratch'):
+ def torqueScript(self, email='adc@geo.au.dk', email_alerts='ae',
+ walltime='24:00:00', queue='qfermi',
+ cudapath='/com/cuda/4.0.17/cuda',
+ spheredir='/home/adc/code/sphere',
+ use_workdir=False, workdir='/scratch'):
'''
Creates a job script for the Torque queue manager for the simulation
object.
t@@ -4594,10 +4567,10 @@ class sim:
'''
- filename = self.sid + ".sh"
- fh = None
+ filename=self.sid + ".sh"
+ fh=None
try :
- fh = open(filename, "w")
+ fh=open(filename, "w")
fh.write('#!/bin/sh\n')
fh.write('#PBS -N ' + self.sid + '\n')
t@@ -4630,78 +4603,8 @@ class sim:
if fh is not None:
fh.close()
- def torqueScriptPenguin(self,
- email='adc@geo.au.dk',
- email_alerts='ae',
- walltime='1920:00:00',
- queue='H30G',
- spheredir='/home/adc/code/sphere'):
- '''
- Creates a job script for the Torque queue manager for the simulation
- object.
-
- :param email: The e-mail address that Torque messages should be sent to
- :type email: str
- :param email_alerts: The type of Torque messages to send to the e-mail
- address. The character 'b' causes a mail to be sent when the
- execution begins. The character 'e' causes a mail to be sent when
- the execution ends normally. The character 'a' causes a mail to be
- sent if the execution ends abnormally. The characters can be written
- in any order.
- :type email_alerts: str
- :param walltime: The maximal allowed time for the job, in the format
- 'HH:MM:SS'.
- :type walltime: str
- :param queue: The Torque queue to schedule the job for
- :type queue: str
- :param spheredir: The path to the root directory of sphere on the
- cluster
- :type spheredir: str
- '''
-
- self.writebin()
- filename = self.sid + ".sh"
- fh = None
- try :
- fh = open(filename, "w")
-
- fh.write('#!/bin/sh\n')
- fh.write('#PBS -S /bin/bash\n')
- fh.write('#PBS -N ' + self.sid + '\n')
- fh.write('#PBS -l nodes=1:ppn=8\n')
- fh.write('#PBS -l walltime=' + walltime + '\n')
- fh.write('#PBS -q ' + queue + '\n')
- fh.write('#PBS -M ' + email + '\n')
- fh.write('#PBS -m ' + email_alerts + '\n')
- fh.write('module load cmake/2.8.11.2\n')
- fh.write('module load cuda/6.0\n')
- fh.write('module load python/2.7.4\n')
- fh.write('module load numpy/1.7.1/python.2.7.4\n')
- #fh.write('module load matplotlib/1.7.1/python.2.7.4\n')
- fh.write('echo "`whoami`@`hostname`\n')
- fh.write('echo "Start at `date`\n')
- fh.write('nvidia-smi\n')
- fh.write('cd ' + spheredir + '\n')
- fh.write('rm CMakeCache.txt\n')
- fh.write('cmake . && make\n')
- fh.write('./sphere input/' + self.sid + '.bin > /dev/null &\n')
- fh.write('wait\n')
- if use_workdir:
- fh.write('cp $WORKDIR/output/* $ORIGDIR/output/\n')
- fh.write('echo "End at `date`"\n')
-
- print('submit with `qsub ' + filename + '`')
-
- finally :
- if fh is not None:
- fh.close()
-
- def render(self,
- method = "pres",
- max_val = 1e3,
- lower_cutoff = 0.0,
- graphics_format = "png",
- verbose=True):
+ def render(self, method="pres", max_val=1e3, lower_cutoff=0.0,
+ graphics_format="png", verbose=True):
'''
Using the built-in ray tracer, render all output files that belong to
the simulation, determined by the simulation id (``sid``).
t@@ -4727,9 +4630,9 @@ class sim:
print("Rendering {} images with the raytracer".format(self.sid))
- quiet = ""
+ quiet=""
if verbose == False:
- quiet = "-q"
+ quiet="-q"
# Render images using sphere raytracer
if method == "normal":
t@@ -4746,15 +4649,9 @@ class sim:
# Convert images to compressed format
convert(graphics_format=graphics_format)
- def video(self,
- out_folder = "./",
- video_format = "mp4",
- graphics_folder = "../img_out/",
- graphics_format = "png",
- fps = 25,
- qscale = 1,
- bitrate = 1800,
- verbose = False):
+ def video(self, out_folder="./", video_format="mp4",
+ graphics_folder="../img_out/", graphics_format="png", fps=25,
+ qscale=1, bitrate=1800, verbose=False):
'''
Uses ffmpeg to combine images to animation. All images should be
rendered beforehand using :func:`render()`.
t@@ -4793,8 +4690,8 @@ 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)
+ #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@@ -4809,8 +4706,8 @@ 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)
+ #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@@ -4832,10 +4729,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@@ -4852,8 +4749,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@@ -4883,10 +4780,10 @@ 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]),
+ 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])
+ self.overlaps=numpy.array(contactdata[:,2])
def findCoordinationNumber(self):
'''
t@@ -4894,7 +4791,7 @@ 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
t@@ -4920,7 +4817,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@@ -4935,10 +4832,10 @@ 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)*(
+ 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
t@@ -4957,7 +4854,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@@ -4999,13 +4896,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@@ -5042,9 +4939,9 @@ 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 \
t@@ -5070,42 +4967,42 @@ 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
+ 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=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
+ dz=zlower - zupper
+ dhoriz=numpy.sqrt(dx**2 + dy**2)
# Find dip angle
diplist.append(math.degrees(math.atan((zupper - zlower)/dhoriz)))
t@@ -5118,7 +5015,7 @@ class sim:
plt.figure(figsize=[4,4])
- ax = plt.subplot(111, polar=True)
+ ax=plt.subplot(111, polar=True)
ax.scatter(strikelist, diplist, c='k', marker='+')
ax.set_rmax(90)
ax.set_rticks([])
t@@ -5141,32 +5038,32 @@ 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
+ 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=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
+ dz=zlower - zupper
+ dhoriz=numpy.sqrt(dx**2 + dy**2)
# Find dip angle
diplist.append(math.degrees(math.atan((zupper - zlower)/dhoriz)))
t@@ -5178,7 +5075,7 @@ class sim:
strikelist.append(2.0*numpy.pi - math.acos(dx/dhoriz))
plt.figure(figsize=[4,4])
- ax = plt.subplot(111, polar=True)
+ ax=plt.subplot(111, polar=True)
ax.scatter(strikelist, diplist, c='k', marker='+')
ax.set_rmax(90)
ax.set_rticks([])
t@@ -5213,7 +5110,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@@ -5232,33 +5129,33 @@ 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=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)
t@@ -5282,10 +5179,10 @@ 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=plt.subplot(111)
ax.plot(porosity, depth,
c='black', linestyle='-', linewidth=1.4)
ax.set_xlabel('Horizontally averaged porosity, [-]')
t@@ -5293,14 +5190,9 @@ class sim:
plt.savefig(self.sid + '-porositiy.' + graphics_format,
transparent=True)
- def thinsection_x1x3(self,
- x2 = 'center',
- graphics_format = 'png',
- cbmax = None,
- arrowscale = 0.01,
- velarrowscale = 1.0,
- slipscale = 1.0,
- verbose = False):
+ def thinsection_x1x3(self, x2='center', graphics_format='png', cbmax=None,
+ arrowscale=0.01, velarrowscale=1.0, slipscale=1.0,
+ verbose=False):
'''
Produce a 2D image of particles on a x1,x3 plane, intersecting the
second axis at x2. Output is saved as '<sid>-ts-x1x3.txt' in the
t@@ -5338,31 +5230,31 @@ class sim:
return
if x2 == 'center' :
- x2 = (self.L[1] - self.origo[1]) / 2.0
+ 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 = []
+ 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 = []
+ 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@@ -5374,9 +5266,9 @@ class sim:
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
t@@ -5386,10 +5278,10 @@ class sim:
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
t@@ -5423,7 +5315,7 @@ class sim:
raise Exception("Error, circle radius is larger than the "
+ "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@@ -5431,10 +5323,10 @@ class sim:
print("Value limited to: " + str(cbmax) + " Pa")
# Save circle data
- filename = '../gnuplot/data/' + self.sid + '-ts-x1x3.txt'
- fh = None
+ filename='../gnuplot/data/' + self.sid + '-ts-x1x3.txt'
+ fh=None
try :
- fh = open(filename, 'w')
+ 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))
t@@ -5444,10 +5336,10 @@ class sim:
fh.close()
# Save circle data for articles spinning with pos. y
- filename = '../gnuplot/data/' + self.sid + '-ts-x1x3-circ.txt'
- fh = None
+ filename='../gnuplot/data/' + self.sid + '-ts-x1x3-circ.txt'
+ fh=None
try :
- fh = open(filename, 'w')
+ fh=open(filename, 'w')
for (x, y, r) in zip(cxlist, cylist, crlist):
fh.write("{}\t{}\t{}\n".format(x, y, r))
t@@ -5460,10 +5352,10 @@ 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
+ filename='../gnuplot/data/' + self.sid + '-ts-x1x3-arrows.txt'
+ fh=None
try :
- fh = open(filename, 'w')
+ 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))
t@@ -5475,10 +5367,10 @@ class sim:
# 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
+ filename='../gnuplot/data/' + self.sid + '-ts-x1x3-velarrows.txt'
+ fh=None
try :
- fh = open(filename, 'w')
+ 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))
t@@ -5489,12 +5381,12 @@ class sim:
# 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@@ -5504,43 +5396,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@@ -5555,10 +5447,10 @@ class sim:
# Write slip lines to text file
- filename = '../gnuplot/data/' + self.sid + '-ts-x1x3-slips.txt'
- fh = None
+ filename='../gnuplot/data/' + self.sid + '-ts-x1x3-slips.txt'
+ fh=None
try :
- fh = open(filename, 'w')
+ 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))
t@@ -5568,7 +5460,7 @@ class sim:
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(\
t@@ -5576,33 +5468,33 @@ class sim:
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)
fig.clf()
- def plotContacts(self, graphics_format = 'png', figsize=[4,4], title=None,
- lower_limit = 0.0, upper_limit = 1.0, alpha=1.0,
+ 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,
+ f_min=None, f_max=None, histogram=True,
forcechains=True):
'''
Plot current contact orientations on polar plot
t@@ -5627,66 +5519,66 @@ 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
+ 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=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
+ dz=zlower - zupper
+ 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)
+ strikelist[j]=numpy.arccos(dx/dhoriz)
else :
- strikelist[j] = 2.0*numpy.pi - numpy.arccos(dx/dhoriz)
+ 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',
+ 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,
t@@ -5714,7 +5606,7 @@ class sim:
if title:
plt.title(title)
else:
- plt.title('t = {:.2f} s'.format(self.currentTime()))
+ plt.title('t={:.2f} s'.format(self.currentTime()))
#plt.tight_layout()
plt.savefig(outfolder + '/contacts-' + self.sid + '-' + \
t@@ -5726,8 +5618,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)
+ n, bins, patches=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@@ -5740,10 +5632,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)
+ n, bins, patches=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@@ -5762,27 +5654,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)
- 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@@ -5807,8 +5699,7 @@ class sim:
if return_data:
return data, strikelist, diplist, forcemagnitude, alpha, f_n_max
- def plotFluidPressuresY(self, y = -1, graphics_format = 'png',
- verbose = True):
+ def plotFluidPressuresY(self, y=-1, graphics_format='png', verbose=True):
'''
Plot fluid pressures in a plane normal to the second axis.
The plot is saved in the current folder with the format
t@@ -5830,26 +5721,25 @@ class sim:
return
if y == -1:
- y = self.num[1]/2
+ y=self.num[1]/2
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)
plt.clf()
plt.close()
- def plotFluidPressuresZ(self, z = -1, graphics_format = 'png',
- verbose = True):
+ def plotFluidPressuresZ(self, z=-1, graphics_format='png', verbose=True):
'''
Plot fluid pressures in a plane normal to the third axis.
The plot is saved in the current folder with the format
t@@ -5871,26 +5761,25 @@ class sim:
return
if z == -1:
- z = self.num[2]/2
+ z=self.num[2]/2
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)
plt.clf()
plt.close()
- def plotFluidVelocitiesY(self, y = -1, graphics_format = 'png',
- verbose = True):
+ def plotFluidVelocitiesY(self, y=-1, graphics_format='png', verbose=True):
'''
Plot fluid velocities in a plane normal to the second axis.
The plot is saved in the current folder with the format
t@@ -5912,50 +5801,49 @@ 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.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')
plt.title("$v_1$")
plt.xlabel('$x_1$')
plt.ylabel('$x_3$')
- plt.colorbar(orientation = 'horizontal')
+ 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')
plt.title("$v_2$")
plt.xlabel('$x_1$')
plt.ylabel('$x_3$')
- plt.colorbar(orientation = 'horizontal')
+ 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')
plt.title("$v_3$")
plt.xlabel('$x_1$')
plt.ylabel('$x_3$')
- plt.colorbar(orientation = 'horizontal')
+ 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)
plt.clf()
plt.close()
- def plotFluidVelocitiesZ(self, z = -1, graphics_format = 'png',
- verbose = True):
+ def plotFluidVelocitiesZ(self, z=-1, graphics_format='png', verbose=True):
'''
Plot fluid velocities in a plane normal to the third axis.
The plot is saved in the current folder with the format
t@@ -5976,52 +5864,52 @@ 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.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')
plt.title("$v_1$")
plt.xlabel('$x_1$')
plt.ylabel('$x_2$')
- plt.colorbar(orientation = 'horizontal')
+ 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')
plt.title("$v_2$")
plt.xlabel('$x_1$')
plt.ylabel('$x_2$')
- plt.colorbar(orientation = 'horizontal')
+ 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')
plt.title("$v_3$")
plt.xlabel('$x_1$')
plt.ylabel('$x_2$')
- plt.colorbar(orientation = 'horizontal')
+ 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)
plt.clf()
plt.close()
- def plotFluidDiffAdvPresZ(self, graphics_format = 'png', verbose = True):
+ def plotFluidDiffAdvPresZ(self, graphics_format='png', verbose=True):
'''
Compare contributions to the velocity from diffusion and advection,
- assuming the flow is 1D along the z-axis, phi = 1, and dphi = 0. This
+ assuming the flow is 1D along the z-axis, phi=1, and dphi=0. This
solution is analog to the predicted velocity and not constrained by the
conservation of mass. The plot is saved in the output folder with the
name format '<simulation id>-diff_adv-t=<current time>s-mu=<dynamic
t@@ -6037,25 +5925,25 @@ class sim:
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
+ 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()
+ fig=plt.figure()
titlesize=12
plt.subplot(1,3,1)
t@@ -6079,12 +5967,12 @@ class sim:
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(
+ 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:
t@@ -6094,7 +5982,7 @@ class sim:
def ReynoldsNumber(self):
'''
- Estimate the per-cell Reynolds number by: Re = rho * ||v_f|| * dx/mu.
+ Estimate the per-cell Reynolds number by: Re=rho * ||v_f|| * dx/mu.
This value is returned and also stored in `self.Re`.
:returns: Reynolds number
t@@ -6102,14 +5990,14 @@ 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_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 + \
+ self.Re=self.rho_f*self.v_f_magn*self.L[0]/self.num[0]/(self.mu + \
1.0e-16)
return self.Re
t@@ -6132,57 +6020,57 @@ 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
+ 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
+ self.H0=H[0]
+ self.H100=H[-1]
+ self.H50=(self.H0 + self.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
+ i_lower=0
+ i_upper=self.status()-1
while (i_upper - i_lower > 1):
- i_midpoint = int((i_upper + i_lower)/2)
+ i_midpoint=int((i_upper + i_lower)/2)
if self.H50 < H[i_midpoint]:
- i_lower = i_midpoint
+ i_lower=i_midpoint
else:
- i_upper = i_midpoint
- self.t50 = t[i_lower] + (t[i_upper] - t[i_lower]) * \
+ 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])
- self.c_coeff = T50*self.H50**2.0/(self.t50)
+ self.c_coeff=T50*self.H50**2.0/(self.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()
+ self.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' \
+ 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))
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=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.grid()
- filename = self.sid + '-loadcurve.' + graphics_format
+ filename=self.sid + '-loadcurve.' + graphics_format
plt.savefig(filename)
if verbose:
print('saved to ' + filename)
t@@ -6197,7 +6085,7 @@ class sim:
See also: :func:`plotConvergence()`
'''
- self.conv = numpy.loadtxt('../output/' + self.sid + '-conv.log',
+ self.conv=numpy.loadtxt('../output/' + self.sid + '-conv.log',
dtype=numpy.int32)
def plotConvergence(self, graphics_format='png', verbose=True):
t@@ -6217,7 +6105,7 @@ class sim:
print('Error: matplotlib module not found (plotConvergence).')
return
- fig = plt.figure()
+ fig=plt.figure()
self.convergence()
plt.title('Convergence evolution in CFD solver in "' + self.sid + '"')
t@@ -6225,18 +6113,15 @@ class sim:
plt.ylabel('Jacobi iterations')
plt.plot(self.conv[:,0], self.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)
plt.clf()
plt.close(fig)
- def plotSinFunction(self,
- baseval, A, f, phi = 0.0,
- xlabel = '$t$ [s]', ylabel = '$y$',
- plotstyle = '.',
- outformat = 'png'):
+ def plotSinFunction(self, baseval, A, f, phi=0.0, xlabel='$t$ [s]',
+ ylabel='$y$', plotstyle='.', outformat='png'):
'''
Plot the values of a sinusoidal modulated base value. Saves the output
as a plot in the current folder.
t@@ -6265,17 +6150,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@@ -6297,8 +6182,8 @@ 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,
t@@ -6311,7 +6196,7 @@ class sim:
See also: :func:`setTopWallNormalStressModulation()`
'''
- self.setTopWallNormalStressModulation(A = 0.0, f = 0.0)
+ self.setTopWallNormalStressModulation(A=0.0, f=0.0)
def setFluidPressureModulation(self, A, f, phi=0.0, plot=False):
'''
t@@ -6330,9 +6215,9 @@ 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,
t@@ -6345,7 +6230,7 @@ class sim:
See also: :func:`setFluidPressureModulation()`
'''
- self.setFluidPressureModulation(A = 0.0, f = 0.0)
+ self.setFluidPressureModulation(A=0.0, f=0.0)
def plotPrescribedFluidPressures(self, graphics_format='png',
verbose=True):
t@@ -6360,20 +6245,20 @@ class sim:
'(plotPrescribedFluidPressures).')
return
- fig = plt.figure()
- conv = numpy.loadtxt('../output/' + self.sid + '-conv.log')
+ fig=plt.figure()
+ conv=numpy.loadtxt('../output/' + self.sid + '-conv.log')
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,\
+ t=numpy.linspace(0,self.time_total,\
self.time_total/self.time_file_dt)
- p = self.p_f[0,0,-1] + \
+ 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@@ -6392,7 +6277,7 @@ class sim:
:return type: numpy.array
'''
if idx == -1:
- idx = range(self.np)
+ idx=range(self.np)
return self.force[idx,:]/(V_sphere(self.radius[idx])*self.rho[0]) + \
self.g
t@@ -6413,7 +6298,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@@ -6433,7 +6318,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@@ -6453,7 +6338,7 @@ class sim:
:func:`setDEMstepsPerCFDstep()` and
:func:`setMaxIterations()`
'''
- self.beta = numpy.asarray(beta)
+ self.beta=numpy.asarray(beta)
def setTolerance(self, tolerance):
'''
t@@ -6471,7 +6356,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@@ -6491,7 +6376,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@@ -6507,7 +6392,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@@ -6526,8 +6411,8 @@ 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]:
t@@ -6541,8 +6426,9 @@ class sim:
def visualize(self, method='energy', savefig=True, outformat='png',
- figsize=False, pickle=False, xlim=False, firststep=0, f_min=None,
- f_max=None, cmap=None, smoothing=0, smoothing_window='hanning'):
+ figsize=False, pickle=False, xlim=False, firststep=0,
+ f_min=None, f_max=None, cmap=None, smoothing=0,
+ smoothing_window='hanning'):
'''
Visualize output from the simulation, where the temporal progress is
of interest. The output will be saved in the current folder with a name
t@@ -6561,7 +6447,7 @@ class sim:
:param figsize: Specify output figure size in inches
:type figsize: array
:param pickle: Save all figure content as a Python pickle file. It can
- be opened later using `fig = pickle.load(open('file.pickle','rb'))`.
+ be opened later using `fig=pickle.load(open('file.pickle','rb'))`.
:type pickle: bool
:param xlim: Set custom limits to the x axis. If not specified, the x
range will correspond to the entire data interval.
t@@ -6581,8 +6467,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@@ -6591,111 +6477,111 @@ 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] +\
+ 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] +\
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')
t@@ -6727,33 +6613,33 @@ 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)
-
- 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)
+ 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()
+ 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@@ -6761,27 +6647,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@@ -6799,28 +6685,28 @@ class sim:
# Read energy values from simulation binaries
for i in numpy.arange(firststep, lastfile+1):
- sb.readstep(i, verbose = False)
+ 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@@ -6831,20 +6717,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@@ -6858,12 +6744,12 @@ class sim:
# Read stress values from simulation binaries
for i in numpy.arange(firststep, lastfile+1):
- sb.readstep(i, verbose = False)
+ sb.readstep(i, verbose=False)
# First iteration: Allocate arrays and find constant values
if i == firststep:
# Shear displacement
- self.xdisp = numpy.zeros(lastfile+1, dtype=numpy.float64)
+ self.xdisp =numpy.zeros(lastfile+1, dtype=numpy.float64)
# Normal stress
self.sigma_eff= numpy.zeros(lastfile+1, dtype=numpy.float64)
t@@ -6872,24 +6758,24 @@ class sim:
self.sigma_def= numpy.zeros(lastfile+1, dtype=numpy.float64)
# Shear stress
- self.tau = numpy.zeros(lastfile+1, dtype=numpy.float64)
+ self.tau =numpy.zeros(lastfile+1, dtype=numpy.float64)
# Upper wall position
- self.dilation = numpy.zeros(lastfile+1, dtype=numpy.float64)
+ self.dilation=numpy.zeros(lastfile+1, dtype=numpy.float64)
# Upper wall position
- self.tau_p = 0.0 # Peak shear stress
+ self.tau_p=0.0 # Peak shear stress
# Shear strain value of peak sh. stress
- self.tau_p_shearstrain = 0.0
+ self.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]:
t@@ -6897,38 +6783,38 @@ class sim:
self.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]
+ 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]
# 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')
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
+ self.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
+ self.tau_p=self.tau[i]
+ self.tau_p_shearstrain=self.xdisp[i]/w_x0
- self.shear_strain = self.xdisp/w_x0
+ self.shear_strain=self.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=self.shear_strain
+ tau=self.tau
+ sigma_def=self.sigma_def
# Optionally smooth the shear stress
if smoothing > 2:
t@@ -6937,23 +6823,23 @@ class sim:
'bartlett', 'blackman']:
raise ValueError
- s = numpy.r_[2*tau[0]-tau[smoothing:1:-1], tau,
+ 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(numpy, smoothing_window)(smoothing)
+ y=numpy.convolve(w/w.sum(), s, mode='same')
+ tau=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\'$ [-]')
t@@ -6976,7 +6862,7 @@ class sim:
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 [%]')
t@@ -6997,12 +6883,12 @@ class sim:
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
+ fh=None
try :
- fh = open(filename, "w")
- L = sb.L[2] - sb.origo[2] # Initial height
+ fh=open(filename, "w")
+ L=sb.L[2] - sb.origo[2] # Initial height
for i in numpy.arange(firststep, lastfile+1):
# format: shear distance [mm], sigma [kPa], tau [kPa],
# Dilation [%]
t@@ -7016,16 +6902,16 @@ class sim:
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)
+ sb.readstep(i, verbose=False)
# First iteration: Allocate arrays and find constant values
if i == firststep:
# Shear displacement
- self.xdisp = numpy.zeros(lastfile+1, dtype=numpy.float64)
+ self.xdisp =numpy.zeros(lastfile+1, dtype=numpy.float64)
# Normal stress
self.sigma_eff= numpy.zeros(lastfile+1, dtype=numpy.float64)
t@@ -7034,37 +6920,37 @@ class sim:
self.sigma_def= numpy.zeros(lastfile+1, dtype=numpy.float64)
# Shear stress
- self.tau_eff = numpy.zeros(lastfile+1, dtype=numpy.float64)
+ self.tau_eff =numpy.zeros(lastfile+1, dtype=numpy.float64)
# Upper wall position
- self.dilation = numpy.zeros(lastfile+1, dtype=numpy.float64)
+ self.dilation=numpy.zeros(lastfile+1, dtype=numpy.float64)
# Mean porosity
- self.phi_bar = numpy.zeros(lastfile+1, dtype=numpy.float64)
+ self.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)
+ self.p_f_bar=numpy.zeros(lastfile+1, dtype=numpy.float64)
+ self.p_f_top=numpy.zeros(lastfile+1, dtype=numpy.float64)
# Upper wall position
- self.tau_p = 0.0 # Peak shear stress
+ self.tau_p=0.0 # Peak shear stress
# Shear strain value of peak sh. stress
- self.tau_p_shearstrain = 0.0
+ self.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)
+ self.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]:
t@@ -7072,53 +6958,53 @@ class sim:
self.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()
+ self.xdisp[i]=sb.xyzsum[fixvel,0].max()
+ self.v[i] =sb.vel[fixvel,0].max()
- self.sigma_eff[i] = sb.w_force[0]/A
- self.sigma_def[i] = sb.currentNormalStress()
+ self.sigma_eff[i]=sb.w_force[0]/A
+ self.sigma_def[i]=sb.currentNormalStress()
# dilation in number of mean particle diameters
- self.dilation[i] = (sb.w_x[0] - w_x0)/d_bar
+ self.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])
+ self.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]
+ 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]
# 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
+ self.tau_p=self.tau_eff[i]
+ self.tau_p_shearstrain=self.xdisp[i]/w_x0
- self.shear_strain = self.xdisp/w_x0
+ self.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))
+ fig=plt.figure(figsize=(8,12))
- #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))
# Upper plot
- ax1 = plt.subplot(3,1,1)
+ ax1=plt.subplot(3,1,1)
ax1.plot(time, self.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')
t@@ -7131,12 +7017,12 @@ class sim:
tl.set_color(ax2color)
# Middle plot
- ax5 = plt.subplot(3, 1, 2, sharex=ax1)
+ ax5=plt.subplot(3, 1, 2, sharex=ax1)
ax5.semilogy(time[1:], self.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),
t@@ -7145,16 +7031,16 @@ class sim:
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,
+ lns0=ax3.plot(time, self.sigma_def/1000.0,
'-k', label="$\\sigma_0$")
- lns1 = ax3.plot(time, self.sigma_eff/1000.0,
+ lns1=ax3.plot(time, self.sigma_eff/1000.0,
'--k', label="$\\sigma'$")
- lns2 = ax3.plot(time,
+ 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,
+ lns3=ax3.plot(time, self.tau_eff/1000.0,
'--r', label="$\\tau'$")
ax3.set_ylabel('Stress [kPa]')
#ax3.legend(loc='upper left')
t@@ -7172,13 +7058,13 @@ class sim:
if self.fluid:
- ax4 = ax3.twinx()
- #ax4color = '#666666'
- ax4color = ax2color
- lns4 = ax4.plot(time, self.p_f_top/1000.0, '-',
+ 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, '--',
+ lns5=ax4.plot(time, self.p_f_bar/1000.0, '--',
color=ax4color,
label='$\\bar{p}_\\text{f}$')
ax4.set_ylabel('Mean fluid pressure '
t@@ -7187,8 +7073,8 @@ class sim:
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)
if xlim:
t@@ -7216,60 +7102,60 @@ 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
+ #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()
# 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(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)
#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,
+ CS=ax1.scatter(friction, shearstrainrate_nonzero,
c=shearstrain_nonzero, linewidth=0.1,
cmap=cmap)
else:
- CS = ax1.scatter(friction, shearstrainrate_nonzero,
+ 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(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
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]')
t@@ -7278,7 +7164,7 @@ class sim:
ax1.set_ylabel('Shear strain rate $\\dot{\\gamma}$ [s$^{-1}$]')
'''
- ax2 = plt.subplot(212)
+ ax2=plt.subplot(212)
ax2.plot(tau/N, v, '.')
ax1.set_xlabel('Friction $\\tau/N$ [-]')
ax1.set_ylabel('Shear velocity [m/s]')
t@@ -7286,13 +7172,13 @@ class sim:
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()
+ sb.readstep(i, verbose=False)
+ t[i]=sb.currentTime()
+ I[i]=sb.inertiaParameterPlanarShear()
# Plotting
if outformat != 'txt':
t@@ -7301,7 +7187,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@@ -7312,15 +7198,15 @@ class sim:
# Read pressure values from simulation binaries
for i in numpy.arange(firststep, lastfile+1):
- sb.readstep(i, verbose = False)
+ sb.readstep(i, verbose=False)
# 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@@ -7329,7 +7215,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@@ -7339,64 +7225,64 @@ 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)
+ 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)
# 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]
+ 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))
+ p_ext=numpy.max(numpy.abs(pres[0:9,:])) # for article2
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(
+ im1=ax.pcolormesh(
x, zpos_c, pres,
#cmap=matplotlib.cm.get_cmap('bwr'),
cmap=cmap,
vmin=-p_ext, vmax=p_ext,
rasterized=True)
else:
- im1 = ax.pcolormesh(
+ im1=ax.pcolormesh(
x, zpos_c, pres,
#cmap=matplotlib.cm.get_cmap('bwr'),
cmap=matplotlib.cm.get_cmap('RdBu_r'),
t@@ -7423,7 +7309,7 @@ class sim:
# 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@@ -7435,44 +7321,44 @@ class sim:
raise Exception('Porosities can only be visualized in wet ' +
'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)
+ 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)
# 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(
+ im1=ax.pcolormesh(
x, zpos_c, poros,
cmap=cmap,
vmin=poros_min, vmax=poros_max,
rasterized=True)
else:
- im1 = ax.pcolormesh(
+ im1=ax.pcolormesh(
x, zpos_c, poros,
cmap=matplotlib.cm.get_cmap('Blues_r'),
#cmap=matplotlib.cm.get_cmap('bwr'),
t@@ -7497,32 +7383,32 @@ class sim:
#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()
- plt.subplots_adjust(wspace = .05)
+ plt.subplots_adjust(wspace=.05)
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(
+ 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(
+ title="t={:.2f} s, $N$={:.0f} kPa".format(
sb.currentTime(),
sb.currentNormalStress('defined')/1000.),
- outfolder = '../img_out/')
+ outfolder='../img_out/')
# render images to movie
subprocess.call('cd ../img_out/ && ' +
t@@ -7538,10 +7424,10 @@ class sim:
# Optional save of figure content
if xlim:
- filename = '{0}-{1}-{3}.{2}'.format(self.sid, method, outformat,
+ 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@@ -7556,7 +7442,7 @@ class sim:
plt.show()
-def convert(graphics_format = 'png', folder = '../img_out'):
+def convert(graphics_format='png', folder='../img_out'):
'''
Converts all PPM images in img_out to graphics_format using Imagemagick. All
PPM images are subsequently removed.
t@@ -7567,8 +7453,8 @@ def convert(graphics_format = 'png', folder = '../img_out'):
:type folder: str
'''
- #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@@ -7578,12 +7464,8 @@ def convert(graphics_format = 'png', folder = '../img_out'):
# Remove PPM files
subprocess.call('rm ' + folder + '/*.ppm', shell=True)
-def render(binary,
- method = 'pres',
- max_val = 1e3,
- lower_cutoff = 0.0,
- graphics_format = 'png',
- verbose=True):
+def render(binary, method='pres', max_val=1e3, lower_cutoff=0.0,
+ graphics_format='png', verbose=True):
'''
Render target binary using the ``sphere`` raytracer.
t@@ -7605,9 +7487,9 @@ def render(binary,
:param verbose: Show verbose information during ray tracing
:type verbose: bool
'''
- quiet = ''
+ quiet=''
if verbose == False:
- quiet = '-q'
+ quiet='-q'
# Render images using sphere raytracer
if method == 'normal':
t@@ -7622,15 +7504,9 @@ def render(binary,
# Convert images to compressed format
convert(graphics_format)
-def video(project,
- out_folder = './',
- video_format = 'mp4',
- graphics_folder = '../img_out/',
- graphics_format = 'png',
- fps = 25,
- qscale = 1,
- bitrate = 1800,
- verbose = False):
+def video(project, out_folder='./', video_format='mp4',
+ graphics_folder='../img_out/', graphics_format='png', fps=25,
+ qscale=1, bitrate=1800, verbose=False):
'''
Uses ffmpeg to combine images to animation. All images should be
rendered beforehand using :func:`render()`.
t@@ -7656,9 +7532,9 @@ def video(project,
'''
# Possible loglevels:
# quiet, panic, fatal, error, warning, info, verbose, debug
- loglevel = 'info' # verbose = True
+ loglevel='info' # verbose=True
if verbose == False:
- loglevel = 'error'
+ loglevel='error'
subprocess.call(\
'ffmpeg -qscale {0} -r {1} -b {2} -y '.format(\
t@@ -7668,13 +7544,8 @@ def video(project,
+ graphics_format + ' ' \
+ out_folder + '/' + project + '.' + video_format, shell=True)
-def thinsectionVideo(project,
- out_folder = "./",
- video_format = "mp4",
- fps = 25,
- qscale = 1,
- bitrate = 1800,
- verbose = False):
+def thinsectionVideo(project, out_folder="./", video_format="mp4", fps=25,
+ qscale=1, bitrate=1800, verbose=False):
'''
Uses ffmpeg to combine thin section images to an animation. This function
will implicity render the thin section images beforehand.
t@@ -7699,20 +7570,20 @@ def thinsectionVideo(project,
'''
# Render thin section images (png)
- lastfile = status(project)
- sb = sim(fluid = self.fluid)
+ lastfile=status(project)
+ sb=sim(fluid=self.fluid)
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)
- sb.readbin(fn, verbose = False)
- sb.thinsection_x1x3(cbmax = sb.w_sigma0[0]*4.0)
+ 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)
# Combine images to animation
# Possible loglevels:
# quiet, panic, fatal, error, warning, info, verbose, debug
- loglevel = "info" # verbose = True
+ loglevel="info" # verbose=True
if verbose == False:
- loglevel = "error"
+ loglevel="error"
subprocess.call(\
"ffmpeg -qscale {0} -r {1} -b {2} -y ".format(\
t@@ -7735,24 +7606,21 @@ def run(binary, verbose=True, hideinputfile=False):
:type hideinputfile: bool
'''
- quiet = ''
- stdout = ''
+ quiet=''
+ stdout=''
if verbose == False:
- quiet = '-q'
+ quiet='-q'
if hideinputfile:
- stdout = ' > /dev/null'
+ stdout=' > /dev/null'
subprocess.call('cd ..; ./sphere ' + quiet + ' ' + binary + ' ' + stdout, \
shell=True)
-def torqueScriptParallel3(obj1, obj2, obj3,
- email='adc@geo.au.dk',
- email_alerts='ae',
- walltime='24:00:00',
- queue='qfermi',
- cudapath='/com/cuda/4.0.17/cuda',
- spheredir='/home/adc/code/sphere',
- use_workdir=False,
- workdir='/scratch'):
+def torqueScriptParallel3(obj1, obj2, obj3, email='adc@geo.au.dk',
+ email_alerts='ae', walltime='24:00:00',
+ queue='qfermi', cudapath='/com/cuda/4.0.17/cuda',
+ spheredir='/home/adc/code/sphere',
+ use_workdir=False,
+ workdir='/scratch'):
'''
Create job script for the Torque queue manager for three binaries,
executed in parallel, ideally on three GPUs.
t@@ -7787,11 +7655,11 @@ def torqueScriptParallel3(obj1, obj2, obj3,
See also :func:`torqueScript()`
'''
- filename = obj1.sid + '_' + obj2.sid + '_' + obj3.sid + '.sh'
+ filename=obj1.sid + '_' + obj2.sid + '_' + obj3.sid + '.sh'
- fh = None
+ fh=None
try :
- fh = open(filename, "w")
+ fh=open(filename, "w")
fh.write('#!/bin/sh\n')
fh.write('#PBS -N ' + obj1.sid + '_' + obj2.sid + '_' + obj3.sid + '\n')
t@@ -7839,11 +7707,11 @@ def status(project):
:return type: int
'''
- fh = None
+ fh=None
try :
- filepath = "../output/{0}.status.dat".format(project)
- fh = open(filepath)
- data = fh.read()
+ filepath="../output/{0}.status.dat".format(project)
+ fh=open(filepath)
+ data=fh.read()
return int(data.split()[2]) # Return last file number
finally :
if fh is not None: