tadd plot with shear velocities - 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 cb034e35510563d64958f7a01f84ea6c5f932e90
(DIR) parent a06ced8a9b36a17e0c4addf19f7200031ec27713
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Tue, 10 Feb 2015 20:36:50 +0100
add plot with shear velocities
Diffstat:
M python/halfshear-darcy-strength-di… | 19 ++++++++++++-------
M python/sphere.py | 15 +++++++++++++--
2 files changed, 25 insertions(+), 9 deletions(-)
---
(DIR) diff --git a/python/halfshear-darcy-strength-dilation-rate.py b/python/halfshear-darcy-strength-dilation-rate.py
t@@ -22,8 +22,11 @@ sigma0 = 20000.0
#k_c_vals = [3.5e-13, 3.5e-15]
k_c = 3.5e-15
#k_c = 3.5e-13
-mu_f_vals = [1.797e-06, 1.204e-06, 5.0e-8, 1.797e-08]
-#velfac_vals = [0.5, 1.0, 2.0]
+
+# 5.0e-8 results present
+#mu_f_vals = [1.797e-06, 1.204e-06, 5.0e-8, 1.797e-08]
+#mu_f_vals = [1.797e-06, 1.204e-06, 3.594e-07, 1.797e-08]
+mu_f_vals = [1.797e-06, 1.204e-06, 1.797e-08]
velfac = 1.0
t@@ -130,12 +133,13 @@ for c, mu_f in enumerate(mu_f_vals):
if numpy.isclose(mu_f, 1.797e-6):
label = 'ref. shear velocity'
- elif numpy.isclose(mu_f, 1.204e-6):
- label = 'ref. shear velocity$\\times$0.67'
- elif numpy.isclose(mu_f, 1.797e-8):
- label = 'ref. shear velocity$\\times$0.01'
+ #elif numpy.isclose(mu_f, 1.204e-6):
+ #label = 'ref. shear velocity$\\times$0.67'
+ #elif numpy.isclose(mu_f, 1.797e-8):
+ #label = 'ref. shear velocity$\\times$0.01'
else:
- label = '$\\mu_\\text{{f}}$ = {:.3e} Pa s'.format(mu_f)
+ #label = '$\\mu_\\text{{f}}$ = {:.3e} Pa s'.format(mu_f)
+ label = 'ref. shear velocity$\\times${:.2}'.format(mu_f/mu_f_vals[0])
ax1.plot(shear_strain[c][1:], friction[c][1:], \
label=label, linewidth=1,
t@@ -206,6 +210,7 @@ ax1.legend(loc='upper right', prop={'size':18}, fancybox=True,
#framealpha=legend_alpha)
ax1.set_xlim([0.0, 0.2])
+ax1.set_ylim([-7, 45])
ax2.set_xlim([0.0, 0.2])
#ax1.set_ylim([0.0, 1.0])
if pressures:
(DIR) diff --git a/python/sphere.py b/python/sphere.py
t@@ -5930,6 +5930,9 @@ class sim:
d_bar = numpy.mean(sb.radius)*2.0
+ # Shear velocity
+ self.v = numpy.zeros(lastfile+1, dtype=numpy.float64)
+
time[i] = sb.time_current[0]
if (i == 1):
t@@ -5942,6 +5945,7 @@ class sim:
if (i > 0):
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()
t@@ -5975,7 +5979,7 @@ class sim:
#fontproperties=FontProperties(size=14))
# Upper plot
- ax1 = plt.subplot(2, 1, 1)
+ ax1 = plt.subplot(3,1,1)
ax1.plot(time, self.xdisp, 'k', label='Displacement')
ax1.set_ylabel('Horizontal displacement [m]')
t@@ -5994,8 +5998,13 @@ class sim:
for tl in ax2.get_yticklabels():
tl.set_color(ax2color)
+ # Middle plot
+ ax5 = plt.subplot(3, 1, 2, sharex=ax1)
+ ax5.semilogy(time[1:], self.vel[1:], label='Shear velocity')
+ ax5.set_ylabel('Shear velocity [ms$^{-1}$]')
+
# Lower plot
- ax3 = plt.subplot(2, 1, 2, sharex=ax1)
+ ax3 = plt.subplot(3, 1, 3, sharex=ax1)
if sb.w_sigma0_A > 1.0e-3:
lns0 = ax3.plot(time, self.sigma_def/1000.0,
'-k', label="$\\sigma_0$")
t@@ -6044,6 +6053,7 @@ class sim:
if xlim:
ax4.set_xlim(xlim)
+
# aesthetics
ax3.set_xlabel('Time [s]')
#ax1.grid()
t@@ -6053,6 +6063,7 @@ class sim:
ax1.set_xlim(xlim)
ax2.set_xlim(xlim)
ax3.set_xlim(xlim)
+ ax5.set_xlim(xlim)
plt.setp(ax1.get_xticklabels(), visible=False)
fig.tight_layout()