timprove plots, relabel from ultimate to peak shear strength - 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 f6e31138270bf4cb5bd9c7a1682cfcf9d027bcd1
(DIR) parent 302d661d1a7a2332591837095b77e5533e3130ed
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Wed, 8 Oct 2014 11:59:23 +0200
improve plots, relabel from ultimate to peak shear strength
Diffstat:
M python/permeability-results.py | 13 ++++++++-----
M python/shear-results-internals.py | 33 +++++++++++++++++++++----------
M python/shear-results-pressures.py | 114 +++++++++++++++++--------------
M python/shear-results-strain.py | 11 +++++++----
M python/shear-results.py | 14 ++++++++------
M python/sphere.py | 38 ++++++++++++++++++++++++-------
6 files changed, 140 insertions(+), 83 deletions(-)
---
(DIR) diff --git a/python/permeability-results.py b/python/permeability-results.py
t@@ -93,19 +93,22 @@ ax1 = plt.subplot(3,1,1)
ax2 = plt.subplot(3,1,2, sharex=ax1)
ax3 = plt.subplot(3,1,3, sharex=ax1)
#ax4 = plt.subplot(4,1,4, sharex=ax1)
-lines = ['-', '--', '-.', ':']
-markers = ['o', 'x', '^', '+']
+colors = ['g', 'r', 'c', 'y']
+#lines = ['-', '--', '-.', ':']
+lines = ['-', '-', '-', '-']
+#markers = ['o', 'x', '^', '+']
+markers = ['x', 'x', 'x', 'x']
for c in range(len(cvals)):
dpdz[c] /= 1000.0
#plt.plot(dpdz[c], K[c], 'o-', label='$c$ = %.2f' % (cvals[c]))
#plt.semilogx(dpdz[c], K[c], 'o-', label='$c$ = %.2f' % (cvals[c]))
#plt.semilogy(dpdz[c], K[c], 'o-', label='$c$ = %.2f' % (cvals[c]))
ax1.loglog(dpdz[c], K[c], label='$c$ = %.2f' % (cvals[c]),
- linestyle=lines[c], marker=markers[c], color='black')
+ linestyle=lines[c], marker=markers[c], color=colors[c])
ax2.semilogx(dpdz[c], phi_bar[c], label='$c$ = %.2f' % (cvals[c]),
- linestyle=lines[c], marker=markers[c], color='black')
+ linestyle=lines[c], marker=markers[c], color=colors[c])
ax3.loglog(dpdz[c], Re[c], label='$c$ = %.2f' % (cvals[c]),
- linestyle=lines[c], marker=markers[c], color='black')
+ linestyle=lines[c], marker=markers[c], color=colors[c])
#ax4.loglog(dpdz[c], fp_fsum[c], label='$c$ = %.2f' % (cvals[c]),
#linestyle=lines[c], marker=markers[c], color='black')
(DIR) diff --git a/python/shear-results-internals.py b/python/shear-results-internals.py
t@@ -16,7 +16,8 @@ from matplotlib.ticker import MaxNLocator
#steps = [5, 10, 100]
#steps = [5, 10]
steps = sys.argv[3:]
-nsteps_avg = 5 # no. of steps to average over
+#nsteps_avg = 5 # no. of steps to average over
+nsteps_avg = 100 # no. of steps to average over
sigma0 = float(sys.argv[1])
#c_grad_p = 1.0
t@@ -145,6 +146,16 @@ for step_str in steps:
#fig = plt.figure(figsize=(8,5*(len(steps))+1))
fig = plt.figure(figsize=(16,5*(len(steps))+1))
+def color(c):
+ if c == 1.0:
+ return 'green'
+ elif c == 0.1:
+ return 'red'
+ elif c == 0.01:
+ return 'cyan'
+ else:
+ return 'blue'
+
ax = []
for s in numpy.arange(len(steps)):
t@@ -181,26 +192,27 @@ for s in numpy.arange(len(steps)):
sharex=ax[6])) # 6: v_z^f
#ax[s*n+0].plot(xdisp[s], zpos_p[s], ',', color = '#888888')
- ax[s*n+0].plot(xdisp_mean[s], zpos_c[s], color = 'k')
+ ax[s*n+0].plot(xdisp_mean[s], zpos_c[s], color=color(c_grad_p))
#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)
- ax[s*n+1].plot(phi_bar[s,1:], zpos_c[s,1:], '-k')
+ ax[s*n+1].plot(phi_bar[s,1:], zpos_c[s,1:], '-', color=color(c_grad_p))
#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*n+2].plot(dphi_bar[s,1:], zpos_c[s,1:], '--k')
+ ax[s*n+2].plot(dphi_bar[s,1:], zpos_c[s,1:], '--', color=color(c_grad_p))
#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)
- ax[s*n+3].plot(v_z_p[s]*100.0, zpos_p[s], ',', color = '#888888')
- ax[s*n+3].plot(v_z_p_bar[s]*100.0, zpos_c[s], color = 'k')
+ ax[s*n+3].plot(v_z_p[s]*100.0, zpos_p[s], ',', alpha=0.5,
+ color=color(c_grad_p))
+ ax[s*n+3].plot(v_z_p_bar[s]*100.0, zpos_c[s], color=color(c_grad_p))
#ax[s*n+0].plot([0.0,0.0], [0.0, sim.L[2]], '--', color='k')
# hydrostatic pressure distribution
- ax[s*n+4].plot(dev_p[s]/1000.0, zpos_c[s], 'k')
+ ax[s*n+4].plot(dev_p[s]/1000.0, zpos_c[s], color=color(c_grad_p))
#dz = sim.L[2]/sim.num[2]
#wall0_iz = int(sim.w_x[0]/dz)
#y_top = wall0_iz*dz + 0.5*dz
t@@ -219,12 +231,13 @@ for s in numpy.arange(len(steps)):
f_pf_mean_nonzero = f_pf_mean[s][I]
zpos_c_nonzero = zpos_c[s][I]
- ax[s*n+5].plot(f_pf_nonzero, zpos_p_nonzero, ',', color = '#888888')
+ ax[s*n+5].plot(f_pf_nonzero, zpos_p_nonzero, ',', alpha=0.5,
+ color=color(c_grad_p))
#ax[s*4+1].plot(f_pf_mean[s][1:-2], zpos_c[s][1:-2], color = 'k')
- ax[s*n+5].plot(f_pf_mean_nonzero, zpos_c_nonzero, color = 'k')
+ ax[s*n+5].plot(f_pf_mean_nonzero, zpos_c_nonzero, color=color(c_grad_p))
#ax[s*4+1].plot([0.0, 0.0], [0.0, sim.L[2]], '--', color='k')
- ax[s*n+6].plot(v_z_f_bar[s]*100.0, zpos_c[s], color = 'k')
+ ax[s*n+6].plot(v_z_f_bar[s]*100.0, zpos_c[s], color=color(c_grad_p))
#ax[s*n+2].plot([0.0,0.0], [0.0, sim.L[2]], '--', color='k')
(DIR) diff --git a/python/shear-results-pressures.py b/python/shear-results-pressures.py
t@@ -17,51 +17,56 @@ matplotlib.rcParams['image.cmap'] = 'bwr'
sigma0 = float(sys.argv[1])
#c_grad_p = 1.0
-c_grad_p = float(sys.argv[2])
+c_grad_p = [1.0, 0.1]
c_phi = 1.0
+
#sid = 'shear-sigma0=' + str(sigma0) + '-c_phi=' + \
# str(c_phi) + '-c_grad_p=' + str(c_grad_p) + '-hi_mu-lo_visc'
-sid = 'halfshear-sigma0=' + str(sigma0) + '-c=' + str(c_grad_p) + '-shear'
+sid = 'halfshear-sigma0=' + str(sigma0) + '-c=' + str(c_grad_p[0]) + '-shear'
sim = sphere.sim(sid, fluid=True)
sim.readfirst(verbose=False)
# cell midpoint cell positions
-zpos_c = numpy.zeros(sim.num[2])
+zpos_c = numpy.zeros((len(c_grad_p), sim.num[2]))
dz = sim.L[2]/sim.num[2]
-for i in numpy.arange(sim.num[2]):
- zpos_c[i] = i*dz + 0.5*dz
-
-shear_strain = numpy.zeros(sim.status())
+for c in numpy.arange(len(c_grad_p)):
+ for i in numpy.arange(sim.num[2]):
+ zpos_c[c,i] = i*dz + 0.5*dz
-dev_pres = numpy.zeros((sim.num[2], sim.status()))
-pres_static = numpy.ones_like(dev_pres)*600.0e3
+shear_strain = numpy.zeros((len(c_grad_p), sim.status()))
+dev_pres = numpy.zeros((len(c_grad_p), sim.num[2], sim.status()))
+pres_static = numpy.ones_like(dev_pres)*sim.p_f[0,0,-1]
pres = numpy.zeros_like(dev_pres)
-for i in numpy.arange(sim.status()):
- sim.readstep(i, verbose=False)
+for c in numpy.arange(len(c_grad_p)):
+ sim.sid = 'halfshear-sigma0=' + str(sigma0) + '-c=' + str(c_grad_p[c]) \
+ + '-shear'
+
+ for i in numpy.arange(sim.status()):
+
+ sim.readstep(i, verbose=False)
- pres[:,i] = numpy.average(numpy.average(sim.p_f, axis=0), axis=0)
+ pres[c,:,i] = numpy.average(numpy.average(sim.p_f, axis=0), axis=0)
- dz = sim.L[2]/sim.num[2]
- wall0_iz = int(sim.w_x[0]/dz)
- for z in numpy.arange(0, wall0_iz+1):
- pres_static[z,i] = \
- (wall0_iz*dz - zpos_c[z] + 0.5*dz)*sim.rho_f*numpy.abs(sim.g[2])\
- + sim.p_f[0,0,-1]
- #pres_static[z,i] = zpos_c[z]
- #pres_static[z,i] = z
+ dz = sim.L[2]/sim.num[2]
+ wall0_iz = int(sim.w_x[0]/dz)
+ for z in numpy.arange(0, wall0_iz+1):
+ pres_static[c,z,i] = \
+ (wall0_iz*dz - zpos_c[c,z] + 0.5*dz)\
+ *sim.rho_f*numpy.abs(sim.g[2])\
+ + sim.p_f[0,0,-1]
- shear_strain[i] = sim.shearStrain()
+ shear_strain[c,i] = sim.shearStrain()
-dev_pres = pres - pres_static
+ dev_pres[c] = pres[c] - pres_static[c]
#fig = plt.figure(figsize=(8,6))
#fig = plt.figure(figsize=(8,12))
fig = plt.figure(figsize=(8,15))
-min_p = numpy.min(dev_pres)/1000.0
+min_p = numpy.min(dev_pres[0])/1000.0
#max_p = numpy.min(dev_pres)
max_p = numpy.abs(min_p)
t@@ -69,31 +74,41 @@ max_p = numpy.abs(min_p)
#bounds = [min_p, 0, max_p]
#norm = matplotlib.colors.BoundaryNorm(bounds, cmap.N)
-ax1 = plt.subplot(311)
-#ax1 = plt.subplot(111)
-#ax1 = plt.subplot(211)
-#im1 = ax1.pcolormesh(shear_strain, zpos_c, dev_pres/1000.0, rasterized=True,
-# cmap=cmap, norm=norm)
-im1 = ax1.pcolormesh(shear_strain, zpos_c, dev_pres/1000.0, vmin=min_p,
- vmax=max_p, rasterized=True)
-#im1 = ax1.pcolormesh(shear_strain, zpos_c, dev_pres/1000.0, rasterized=True)
-#ax1.set_xlim([0, shear_strain[-1]])
-#ax1.set_ylim([zpos_c[0], sim.w_x[0]])
-ax1.set_xlabel('Shear strain $\\gamma$ [-]')
-ax1.set_ylabel('Vertical position $z$ [m]')
-cb1 = plt.colorbar(im1)
-#cb1 = plt.colorbar(im1, cmap=cmap, norm=norm)
-cb1.set_label('$p_\\text{f} - p^\\text{hyd}_\\text{f}$ [kPa]')
-cb1.solids.set_rasterized(True)
-
-# annotate plot
-#ax1.text(0.02, 0.15, 'compressive',
- #bbox={'facecolor':'white', 'alpha':0.5, 'pad':10})
-
-#ax1.text(0.12, 0.25, 'dilative',
- #bbox={'facecolor':'white', 'alpha':0.5, 'pad':10})
-
-#'''
+ax = []
+for c in numpy.arange(len(c_grad_p)):
+
+ ax.append(plt.subplot(len(c_grad_p), 1, c+1))
+
+ #im1 = ax[c].pcolormesh(shear_strain[c], zpos_c[c], dev_pres[c]/1000.0,
+ #vmin=min_p, vmax=max_p, rasterized=True)
+ im1 = ax[c].pcolormesh(shear_strain[c], zpos_c[c], dev_pres[c]/1000.0,
+ rasterized=True)
+ ax[c].set_xlim([0, shear_strain[c,-1]])
+ ax[c].set_ylim([zpos_c[0,0], sim.w_x[0]])
+ ax[c].set_xlabel('Shear strain $\\gamma$ [-]')
+ ax[c].set_ylabel('Vertical position $z$ [m]')
+
+ #plt.text(0.0, 0.15, '$c = %.2f$' % (c_grad_p[c]),\
+ # horizontalalignment='left', fontsize=22,
+ # transform=ax[c].transAxes)
+ ax[c].set_title('$c = %.2f$' % (c_grad_p[c]))
+
+ #cb = plt.colorbar(im1, orientation='horizontal')
+ cb = plt.colorbar(im1)
+ cb.set_label('$p_\\text{f} - p^\\text{hyd}_\\text{f}$ [kPa]')
+ cb.solids.set_rasterized(True)
+
+ # annotate plot
+ #ax1.text(0.02, 0.15, 'compressive',
+ #bbox={'facecolor':'white', 'alpha':0.5, 'pad':10})
+
+ #ax1.text(0.12, 0.25, 'dilative',
+ #bbox={'facecolor':'white', 'alpha':0.5, 'pad':10})
+
+#cb = plt.colorbar(im1, orientation='horizontal')
+#cb.set_label('$p_\\text{f} - p^\\text{hyd}_\\text{f}$ [kPa]')
+#cb.solids.set_rasterized(True)
+'''
ax2 = plt.subplot(312)
im2 = ax2.pcolormesh(shear_strain, zpos_c, pres/1000.0, rasterized=True)
#ax2.set_xlim([0, shear_strain[-1]])
t@@ -104,7 +119,6 @@ cb2 = plt.colorbar(im2)
cb2.set_label('Pressure $p_\\text{f}$ [kPa]')
cb2.solids.set_rasterized(True)
-#'''
ax3 = plt.subplot(313)
im3 = ax3.pcolormesh(shear_strain, zpos_c, pres_static/1000.0, rasterized=True)
#ax3.set_xlim([0, shear_strain[-1]])
t@@ -114,7 +128,7 @@ ax3.set_ylabel('Vertical position $z$ [m]')
cb3 = plt.colorbar(im3)
cb3.set_label('Static Pressure $p_\\text{f}$ [kPa]')
cb3.solids.set_rasterized(True)
-#'''
+'''
#plt.MaxNLocator(nbins=4)
(DIR) diff --git a/python/shear-results-strain.py b/python/shear-results-strain.py
t@@ -84,7 +84,9 @@ for c in cvals:
fig = plt.figure(figsize=(8,6))
ax = []
-linetype = ['-', '--', '-.']
+#linetype = ['-', '--', '-.']
+linetype = ['-', '-', '-']
+color = ['b','g','r','c']
for s in numpy.arange(len(cvals)):
ax.append(plt.subplot(111))
t@@ -98,9 +100,10 @@ for s in numpy.arange(len(cvals)):
else:
legend = 'wet, c = ' + str(cvals[s])
- ax[0].plot(xdisp[s], zpos_p[s], ',', color = '#888888')
- ax[0].plot(xdisp_mean[s], zpos_c[s], linetype[s], color='k', label = legend,
- linewidth=2)
+ #ax[0].plot(xdisp[s], zpos_p[s], ',', color = '#888888')
+ ax[0].plot(xdisp[s], zpos_p[s], ',', color=color[s], alpha=0.5)
+ ax[0].plot(xdisp_mean[s], zpos_c[s], linetype[s], color=color[s],
+ label=legend, linewidth=2)
ax[0].set_ylabel('Vertical position $z$ [m]')
ax[0].set_xlabel('$\\boldsymbol{x}^x_\\text{p}$ [m]')
(DIR) diff --git a/python/shear-results.py b/python/shear-results.py
t@@ -87,7 +87,7 @@ def smooth(x, window_len=10, window='hanning'):
return y[window_len-1:-window_len+1]
-smooth_window = 40
+smooth_window = 10
shear_strain = [[], [], []]
friction = [[], [], []]
t@@ -188,18 +188,19 @@ ax2 = plt.subplot(212, sharex=ax1)
#ax4 = plt.subplot(414, sharex=ax1)
alpha = 0.5
#ax1.plot(shear_strain[0], friction[0], label='dry', alpha = 0.5)
-ax1.plot(shear_strain[0], friction_smooth[0], label='dry')
-ax2.plot(shear_strain[0], dilation[0], label='dry')
+ax1.plot(shear_strain[0], friction_smooth[0], label='dry', linewidth=1,
+ alpha=0.5)
+ax2.plot(shear_strain[0], dilation[0], label='dry', linewidth=2)
#ax4.plot(shear_strain[0], f_n_mean[0], '-', label='dry', color='blue')
#ax4.plot(shear_strain[0], f_n_max[0], '--', color='blue')
-color = ['b','g','r']
+color = ['b','g','r','c']
for c in numpy.arange(1,len(cvals)+1):
#ax1.plot(shear_strain[c][1:], friction[c][1:], \
#label='$c$ = %.2f' % (cvals[c-1]))
ax1.plot(shear_strain[c][1:], friction_smooth[c][1:], \
- label='$c$ = %.2f' % (cvals[c-1]))
+ label='$c$ = %.2f' % (cvals[c-1]), linewidth=1, alpha=0.3)
ax2.plot(shear_strain[c][1:], dilation[c][1:], \
label='$c$ = %.2f' % (cvals[c-1]), linewidth=2)
t@@ -241,7 +242,8 @@ ax2.grid()
#ax4.grid()
legend_alpha=0.5
-ax1.legend(loc='best', prop={'size':18}, fancybox=True, framealpha=legend_alpha)
+ax1.legend(loc='lower right', prop={'size':18}, fancybox=True,
+ framealpha=legend_alpha)
ax2.legend(loc='lower right', prop={'size':18}, fancybox=True,
framealpha=legend_alpha)
#ax3.legend(loc='lower right', prop={'size':18}, fancybox=True,
(DIR) diff --git a/python/sphere.py b/python/sphere.py
t@@ -4789,6 +4789,28 @@ class sim:
'''
self.ndem = numpy.asarray(ndem)
+ def shearStress(self):
+ '''
+ Calculates the sum of shear stress values measured on any moving
+ particles.
+
+ :returns: The shear stress in Pa
+ :return type: numpy.array
+ '''
+
+ fixvel = numpy.nonzero(self.fixvel > 0.0)
+ shearvel = self.shearVel()
+
+ force = numpy.zeros(3)
+
+ # Summation of shear stress contributions
+ for i in fixvel[0]:
+ if (sb.vel[i,0] > 0.0):
+ force += -sb.force[i,:]
+
+ return force/(sim.L[0]*sim.L[1])
+
+
def visualize(self, method = 'energy', savefig = True, outformat = 'png'):
'''
Visualize output from the simulation, where the temporal progress is
t@@ -5069,9 +5091,9 @@ class sim:
self.dilation = numpy.zeros(lastfile+1, dtype=numpy.float64)
# Upper wall position
- self.tau_u = 0.0 # Peak shear stress
+ self.tau_p = 0.0 # Peak shear stress
# Shear strain value of peak sh. stress
- self.tau_u_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)
t@@ -5082,7 +5104,7 @@ class sim:
# Summation of shear stress contributions
for j in fixvel[0]:
if (sb.vel[j,0] > 0.0):
- self.tau[i] += -sb.force[j,0]
+ self.tau[i] += -sb.force[j,0]/A
if (i > 0):
self.xdisp[i] = self.xdisp[i-1] +sb.time_file_dt[0]*shearvel
t@@ -5103,17 +5125,17 @@ class sim:
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_u):
- self.tau_u = self.tau[i]
- self.tau_u_shearstrain = self.xdisp[i]/w_x0
+ if (self.tau[i] > self.tau_p):
+ self.tau_p = self.tau[i]
+ self.tau_p_shearstrain = self.xdisp[i]/w_x0
self.shear_strain = self.xdisp/w_x0
# Plot stresses
if (outformat != 'txt'):
- shearinfo = "$\\tau_u$ = {:.3} Pa at $\gamma$ = {:.3}".format(\
- self.tau_u, self.tau_u_shearstrain)
+ 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))