tfix velocity limit - 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 d0149c57c01d10ebbe5943c6c9045fb71db3654d
(DIR) parent 0fb9ef5ab8b31d748fa0c42cb04c13dd3fb349fe
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Thu, 12 Mar 2015 15:31:42 +0100
fix velocity limit
Diffstat:
M python/halfshear-darcy-combined.py | 25 ++++++++++++++++++++-----
M python/halfshear-darcy-rate-depend… | 6 +++---
2 files changed, 23 insertions(+), 8 deletions(-)
---
(DIR) diff --git a/python/halfshear-darcy-combined.py b/python/halfshear-darcy-combined.py
t@@ -53,7 +53,6 @@ v = numpy.empty_like(t)
# displacement and mean porosity plot
xdisp = numpy.empty_like(t)
-xdispint = numpy.zeros_like(t)
phi_bar = numpy.empty_like(t)
# mean horizontal porosity plot
t@@ -136,7 +135,22 @@ if calculateforcechains:
else:
n, nkept, coordinationnumber = numpy.loadtxt(sid + '-fc.txt')
+# Transform time from model time to real time [s]
+t = t/t_DEM_to_t_real
+## integrate velocities to displacement along x (xdispint)
+# Taylor two term expansion
+xdispint = numpy.zeros_like(t)
+v_limit = 1.0e-7
+dt = (t[1] - t[0])
+dt2 = dt*2.
+for i in numpy.arange(t.size):
+ if i > 0 and i < t.size-1:
+ acc = (numpy.min([v[i+1], v_limit]) - numpy.min([v[i-1], v_limit]))/dt2
+ xdispint[i] = xdispint[i-1] +\
+ numpy.min([v[i], v_limit])*dt + 0.5*acc*dt**2
+ elif i == t.size-1:
+ xdispint[i] = xdispint[i-1] + numpy.min([v[i], v_limit])*dt
##################
t@@ -149,7 +163,8 @@ horizontalalignment='left'
fontweight='bold'
bbox={'facecolor':'white', 'alpha':1.0, 'pad':3}
-t = t/t_DEM_to_t_real / (60.*60.*24.)
+# Time in days
+t = t/(60.*60.*24.)
fig = plt.figure(figsize=[3.5,8])
t@@ -215,7 +230,7 @@ ax3.text(bbox_x, bbox_y, 'b',
## ax5: xdisp, ax6: mean(phi)
ax5 = plt.subplot(5, 1, 3, sharex=ax1)
-ax5.plot(t, xdisp, 'k', linewidth=linewidth)
+ax5.plot(t, xdispint, 'k', linewidth=linewidth)
ax5.set_ylabel('Shear displacement [m]')
ax6color='blue'
t@@ -235,14 +250,14 @@ ax6.text(bbox_x, bbox_y, 'c',
## ax7: n_heavy, dn_heavy, ax8: z
ax7 = plt.subplot(5, 1, 4, sharex=ax1)
-ax7.semilogy(t, n, 'k', label='$n_\\text{heavy}$', linewidth=linewidth)
+ax7.semilogy(t[:n.size], n, 'k', label='$n_\\text{heavy}$', linewidth=linewidth)
ax7.set_ylabel('Number of heavily loaded contacts [-]')
#ax7.semilogy(t, n - nkept, 'b', label='$\Delta n_\\text{heavy}$')
ax7.set_ylim([1.0e1, 2.0e4])
ax8 = ax7.twinx()
ax8color='green'
-ax8.plot(t, coordinationnumber, color=ax8color, linewidth=linewidth)
+ax8.plot(t[:n.size], coordinationnumber, color=ax8color, linewidth=linewidth)
ax8.set_ylabel('Contacts per particle [-]')
ax8.yaxis.label.set_color(ax8color)
for tl in ax8.get_yticklabels():
(DIR) diff --git a/python/halfshear-darcy-rate-dependence.py b/python/halfshear-darcy-rate-dependence.py
t@@ -104,7 +104,7 @@ for sid in sids:
popt2, pvoc2 = scipy.optimize.curve_fit(
creep_rheology, tau_nonzero[idxfit2]/N_nonzero[idxfit2],
shearstrainrate_nonzero[idxfit2])
- print '# Consolidatet state'
+ print '# Consolidated state'
print popt2
print pvoc2
n2 = popt2[0] # stress exponent
t@@ -133,7 +133,7 @@ for sid in sids:
#cmap=matplotlib.cm.get_cmap('afmhot'))
## plastic limit
- x = [0.3, 0.3]
+ x = [0.28, 0.28]
y = ax1.get_ylim()
limitcolor = '#333333'
ax1.plot(x, y, '--', linewidth=2, color=limitcolor)
t@@ -145,7 +145,7 @@ for sid in sids:
## Fit
ax1.plot(friction_fit, strainrate_fit)
#ax1.plot(friction_fit2, strainrate_fit2)
- ax1.annotate('$\\dot{\\gamma} = (\\tau/N)^{6.4}$',
+ ax1.annotate('$\\dot{{\\gamma}} = (\\tau/N)^{{{:.1f}}}$'.format(n),
xy = (friction_fit[40], strainrate_fit[40]),
xytext = (0.32+0.05, 2.0e-9),
arrowprops=dict(facecolor='blue', edgecolor='blue', shrink=0.1,