tstill improving plot, added starter for half the particle assemblage - 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 c53a170cfe0eea60355d61b90f6b0329a9f18a12
(DIR) parent 23718477440c1df19d0f2da9b54e78177cac7c19
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Tue, 30 Sep 2014 09:02:07 +0200
still improving plot, added starter for half the particle assemblage
Diffstat:
A python/halfshear-starter.py | 63 +++++++++++++++++++++++++++++++
M python/shear-results-forces.py | 27 ++++++++++++++++++---------
2 files changed, 81 insertions(+), 9 deletions(-)
---
(DIR) diff --git a/python/halfshear-starter.py b/python/halfshear-starter.py
t@@ -0,0 +1,63 @@
+#!/usr/bin/env python
+import sphere
+import numpy
+import sys
+
+# launch with:
+# $ python shear-starter.py <DEVICE> <FLUID> <C_PHI> <C_GRAD_P> <SIGMA_0>
+
+device = int(sys.argv[1])
+wet = int(sys.argv[2])
+c_phi = float(sys.argv[3])
+c_grad_p = float(sys.argv[4])
+sigma0 = float(sys.argv[5])
+
+#sim = sphere.sim('diffusivity-sigma0=' + str(sigma0) + '-c_phi=' + \
+# str(c_phi) + '-c_grad_p=' + str(c_grad_p), fluid=True)
+if wet == 1:
+ fluid = True
+else:
+ fluid = False
+
+#sim = sphere.sim('diffusivity-sigma0=' + str(sigma0) +'-c_phi=1.0-c_grad_p=1.0',
+# fluid=True)
+sim = sphere.sim('halfshear-sigma0=' + str(sigma0) + '-hw-relax', fluid=False)
+sim.readlast()
+
+print(sim.sid)
+sim.fluid = fluid
+
+sim.checkerboardColors(nx=6,ny=6,nz=6)
+sim.cleanup()
+sim.adjustUpperWall()
+sim.zeroKinematics()
+
+sim.shear(1.0/20.0)
+
+if fluid:
+ #sim.num[2] *= 2
+ #sim.L[2] *= 2.0
+ sim.initFluid(mu = 1.787e-6, p = 600.0e3, hydrostatic = True)
+ #sim.initFluid(mu = 17.87e-4, p = 1.0e5, hydrostatic = True)
+sim.setFluidBottomNoFlow()
+sim.setFluidTopFixedPressure()
+sim.setDEMstepsPerCFDstep(10)
+sim.setMaxIterations(2e5)
+sim.initTemporal(total = 20.0, file_dt = 0.01, epsilon=0.07)
+sim.c_phi[0] = c_phi
+sim.c_grad_p[0] = c_grad_p
+sim.w_devs[0] = sigma0
+#sim.w_m[0] = numpy.abs(sigma0*sim.L[0]*sim.L[1]/sim.g[2])
+sim.mu_s[0] = 0.5
+sim.mu_d[0] = 0.5
+
+# Fix lowermost particles
+dz = sim.L[2]/sim.num[2]
+I = numpy.nonzero(sim.x[:,2] < 1.5*dz)
+sim.fixvel[I] = 1
+
+sim.run(dry=True)
+sim.run(device=device)
+#sim.writeVTKall()
+#sim.visualize('walls')
+#sim.visualize('fluid-pressure')
(DIR) diff --git a/python/shear-results-forces.py b/python/shear-results-forces.py
t@@ -44,11 +44,14 @@ xdisp = numpy.zeros((len(steps), sim.np))
f_pf = numpy.zeros_like(xdisp)
# pressure - hydrostatic pressure
-dev_p = numpy.zeros((len(steps), sim.num[2]))
+#dev_p = numpy.zeros((len(steps), sim.num[2]))
# mean porosity
phi_bar = numpy.zeros((len(steps), sim.num[2]))
+# mean porosity change
+dphi_bar = numpy.zeros((len(steps), sim.num[2]))
+
# mean per-particle values
xdisp_mean = numpy.zeros((len(steps), sim.num[2]))
f_pf_mean = numpy.zeros((len(steps), sim.num[2]))
t@@ -82,14 +85,18 @@ for step_str in steps:
'''
f_pf[s,:] += sim.f_sum[:,2]
- dev_p[s,:] += \
- numpy.average(numpy.average(sim.p_f, axis=0), axis=0)\
- /nsteps_avg
+ #dev_p[s,:] += \
+ #numpy.average(numpy.average(sim.p_f, axis=0), axis=0)\
+ #/nsteps_avg
phi_bar[s,:] += \
numpy.average(numpy.average(sim.phi, axis=0), axis=0)\
/nsteps_avg
+ dphi_bar[s,:] += \
+ numpy.average(numpy.average(sim.dphi, axis=0), axis=0)\
+ /nsteps_avg/sim.time_dt
+
shear_strain[s] += sim.shearStrain()/nsteps_avg
# calculate mean values of xdisp and f_pf
t@@ -122,14 +129,15 @@ for s in numpy.arange(len(steps)):
ax[s*4+1].plot(f_pf_mean[s], zpos_c[s], color = 'k')
ax[s*4+1].plot([0.0, 0.0], [0.0, sim.L[2]], '--', color='k')
- ax[s*4+2].plot(dev_p[s]/1000.0, zpos_c[s], 'k')
+ #ax[s*4+2].plot(dev_p[s]/1000.0, zpos_c[s], 'k')
+ ax[s*4+2].plot(phi_bar[s,1:], zpos_c[s,1:], '-k', linewidth=3)
#phicolor = '#888888'
#ax[s*4+3].plot(phi_bar[s], zpos_c[s], '-', color = phicolor)
#for tl in ax[s*4+3].get_xticklabels():
#tl.set_color(phicolor)
- ax[s*4+3].plot(phi_bar[s,1:], zpos_c[s,1:], '-k', linewidth=3)
- ax[s*4+3].plot(phi_bar[s,1:], zpos_c[s,1:], '-w', linewidth=2)
+ ax[s*4+3].plot(dphi_bar[s,1:], zpos_c[s,1:], '-k', linewidth=3)
+ ax[s*4+3].plot(dphi_bar[s,1:], zpos_c[s,1:], '-w', linewidth=2)
max_z = numpy.max(zpos_p)
ax[s*4+0].set_ylim([0, max_z])
t@@ -143,9 +151,10 @@ for s in numpy.arange(len(steps)):
ax[s*4+0].set_ylabel('Vertical position $z$ [m]')
ax[s*4+0].set_xlabel('$\\boldsymbol{x}^x_\\text{p}$ [m]')
ax[s*4+1].set_xlabel('$\\boldsymbol{f}^z_\\text{pf}$ [N]')
- ax[s*4+2].set_xlabel('$\\bar{p_\\text{f}}$ [kPa]')
+ #ax[s*4+2].set_xlabel('$\\bar{p_\\text{f}}$ [kPa]')
#ax[s*4+3].set_xlabel('$\\bar{\\phi}$ [-]', color=phicolor)
- ax[s*4+3].set_xlabel('$\\bar{\\phi}$ [-]')
+ ax[s*4+2].set_xlabel('$\\bar{\\phi}$ [-]')
+ ax[s*4+3].set_xlabel('$\\delta \\bar{\\phi}/\\delta t$ [-]')
plt.setp(ax[s*4+1].get_yticklabels(), visible=False)
plt.setp(ax[s*4+2].get_yticklabels(), visible=False)