timprove friction,dilation,pressure vs rate plot - sphere - GPU-based 3D discre… | |
git clone git://src.adamsgaard.dk/sphere | |
Log | |
Files | |
Refs | |
LICENSE | |
--- | |
commit e1ca057c3a2749911d17d6f9356f67e8470be834 | |
parent e05f1f05655a65d9f3920854b5e072c227285260 | |
Author: Anders Damsgaard <[email protected]> | |
Date: Tue, 14 Apr 2015 16:02:19 +0200 | |
improve friction,dilation,pressure vs rate plot | |
Diffstat: | |
M python/halfshear-darcy-strength-di… | 393 +++++++++++++++++++++------… | |
1 file changed, 275 insertions(+), 118 deletions(-) | |
--- | |
diff --git a/python/halfshear-darcy-strength-dilation-rate.py b/python/halfshea… | |
t@@ -23,6 +23,8 @@ sns.despine() # remove right and top spines | |
pressures = True | |
zflow = False | |
contact_forces = False | |
+smooth_friction = True | |
+smooth_window = 30 | |
#sigma0_list = numpy.array([1.0e3, 2.0e3, 4.0e3, 10.0e3, 20.0e3, 40.0e3]) | |
sigma0 = 20000.0 | |
t@@ -32,10 +34,72 @@ k_c = 3.5e-15 | |
# 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, 3.594e-07, 1.797e-08] | |
+#mu_f_vals = [1.797e-06, 1.797-07, 1.797e-08] | |
+#mu_f_vals = [1.797e-06, 1.797-07, 1.797e-08] | |
#mu_f_vals = [1.797e-06, 1.204e-06, 1.797e-08] | |
+mu_f_vals = [1.797e-06, 3.594e-07, 1.797e-08] | |
velfac = 1.0 | |
+# return a smoothed version of in. The returned array is smaller than the | |
+# original input array | |
+def smooth(x, window_len=10, window='hanning'): | |
+ """smooth the data using a window with requested size. | |
+ | |
+ This method is based on the convolution of a scaled window with the signal. | |
+ The signal is prepared by introducing reflected copies of the signal | |
+ (with the window size) in both ends so that transient parts are minimized | |
+ in the begining and end part of the output signal. | |
+ | |
+ input: | |
+ x: the input signal | |
+ window_len: the dimension of the smoothing window | |
+ window: the type of window from 'flat', 'hanning', 'hamming', 'bartlet… | |
+ flat window will produce a moving average smoothing. | |
+ | |
+ output: | |
+ the smoothed signal | |
+ | |
+ example: | |
+ | |
+ import numpy as np | |
+ t = np.linspace(-2,2,0.1) | |
+ x = np.sin(t)+np.random.randn(len(t))*0.1 | |
+ y = smooth(x) | |
+ | |
+ see also: | |
+ | |
+ numpy.hanning, numpy.hamming, numpy.bartlett, numpy.blackman, numpy.convol… | |
+ scipy.signal.lfilter | |
+ | |
+ TODO: the window parameter could be the window itself if an array instead … | |
+ """ | |
+ | |
+ if x.ndim != 1: | |
+ raise ValueError, "smooth only accepts 1 dimension arrays." | |
+ | |
+ if x.size < window_len: | |
+ raise ValueError, "Input vector needs to be bigger than window size." | |
+ | |
+ if window_len < 3: | |
+ return x | |
+ | |
+ if not window in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']: | |
+ raise ValueError, "Window is on of 'flat', 'hanning', 'hamming', 'bart… | |
+ | |
+ s=numpy.r_[2*x[0]-x[window_len:1:-1], x, 2*x[-1]-x[-1:-window_len:-1]] | |
+ #print(len(s)) | |
+ | |
+ if window == 'flat': #moving average | |
+ w = numpy.ones(window_len,'d') | |
+ else: | |
+ w = getattr(numpy, window)(window_len) | |
+ y = numpy.convolve(w/w.sum(), s, mode='same') | |
+ return y[window_len-1:-window_len+1] | |
+ | |
+ | |
+smooth_window = 10 | |
+ | |
shear_strain = [[], [], [], []] | |
friction = [[], [], [], []] | |
t@@ -123,37 +187,39 @@ for c in numpy.arange(0,len(mu_f_vals)): | |
if zflow or pressures: | |
#fig = plt.figure(figsize=(8,10)) | |
- fig = plt.figure(figsize=(3.74, 2*3.74)) | |
+ #fig = plt.figure(figsize=(3.74, 2*3.74)) | |
+ fig = plt.figure(figsize=(2*3.74, 2*3.74)) | |
else: | |
fig = plt.figure(figsize=(8,8)) # (w,h) | |
#fig = plt.figure(figsize=(8,12)) | |
#fig = plt.figure(figsize=(8,16)) | |
-fig.subplots_adjust(hspace=0.0) | |
#plt.subplot(3,1,1) | |
#plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0)) | |
-if zflow or pressures: | |
- ax1 = plt.subplot(311) | |
- ax2 = plt.subplot(312, sharex=ax1) | |
- ax3 = plt.subplot(313, sharex=ax1) | |
-else: | |
- ax1 = plt.subplot(211) | |
- ax2 = plt.subplot(212, sharex=ax1) | |
-#ax3 = plt.subplot(413, sharex=ax1) | |
-#ax4 = plt.subplot(414, sharex=ax1) | |
-#alpha = 0.5 | |
-alpha = 1.0 | |
-#ax1.plot(shear_strain[0], friction[0], label='dry', linewidth=1, alpha=alpha) | |
-#ax2.plot(shear_strain[0], dilation[0], label='dry', linewidth=1) | |
-#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','c'] | |
-#color = ['g','r','c'] | |
-color = sns.color_palette() | |
-#for c, mu_f in enumerate(mu_f_vals): | |
-for c in numpy.arange(len(mu_f_vals)-1, -1, -1): | |
+for c in numpy.arange(0,len(mu_f_vals)): | |
+ | |
+ if zflow or pressures: | |
+ ax1 = plt.subplot(3, len(mu_f_vals), 1+c) | |
+ ax2 = plt.subplot(3, len(mu_f_vals), 4+c, sharex=ax1) | |
+ ax3 = plt.subplot(3, len(mu_f_vals), 7+c, sharex=ax1) | |
+ else: | |
+ ax1 = plt.subplot(211) | |
+ ax2 = plt.subplot(212, sharex=ax1) | |
+ #ax3 = plt.subplot(413, sharex=ax1) | |
+ #ax4 = plt.subplot(414, sharex=ax1) | |
+ #alpha = 0.5 | |
+ alpha = 1.0 | |
+ #ax1.plot(shear_strain[0], friction[0], label='dry', linewidth=1, alpha=al… | |
+ #ax2.plot(shear_strain[0], dilation[0], label='dry', linewidth=1) | |
+ #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','c'] | |
+ #color = ['g','r','c'] | |
+ color = sns.color_palette() | |
+ #for c, mu_f in enumerate(mu_f_vals): | |
+ #for c in numpy.arange(len(mu_f_vals)-1, -1, -1): | |
mu_f = mu_f_vals[c] | |
if numpy.isclose(mu_f, 1.797e-6): | |
t@@ -164,12 +230,20 @@ for c in numpy.arange(len(mu_f_vals)-1, -1, -1): | |
#label = 'ref. shear velocity$\\times$0.01' | |
else: | |
#label = '$\\mu_\\text{{f}}$ = {:.3e} Pa s'.format(mu_f) | |
- label = 'ref. shear velocity$\\times${:.2}'.format(mu_f/mu_f_vals[0]) | |
+ label = 'shear velocity$\\times${:.2}'.format(mu_f/mu_f_vals[0]) | |
+ | |
+ # unsmoothed | |
+ ax1.plot(shear_strain[c][1:], friction[c][1:], \ | |
+ label=label, linewidth=1, | |
+ alpha=0.2, color='gray') | |
+ #alpha=alpha, color=color[c]) | |
- ax1.plot(shear_strain[c], friction[c], \ | |
+ # smoothed | |
+ ax1.plot(shear_strain[c][1:], smooth(friction[c], smooth_window)[1:], \ | |
label=label, linewidth=1, | |
alpha=alpha, color=color[c]) | |
+ | |
ax2.plot(shear_strain[c], dilation[c], \ | |
label=label, linewidth=1, | |
color=color[c]) | |
t@@ -179,111 +253,194 @@ for c in numpy.arange(len(mu_f_vals)-1, -1, -1): | |
label=label, linewidth=1) | |
if pressures: | |
- #ax3.plot(shear_strain[c][1:], p_max[c][1:], '-', color=color[c], | |
- #alpha=0.5) | |
+ #ax3.plot(shear_strain[c], p_max[c], '-', color=color[c], alpha=0.5) | |
+ | |
ax3.plot(shear_strain[c], p_mean[c], '-', color=color[c], \ | |
label=label, linewidth=1) | |
- #ax3.plot(shear_strain[c][1:], p_min[c][1:], '-', color=color[c], | |
- #alpha=0.5) | |
- #ax3.fill_between(shear_strain[c][1:], p_min[c][1:], p_max[c][1:], | |
- #where=p_min[c][1:]<=p_max[c][1:], facecolor=color[c], | |
- #interpolate=True, alpha=0.5) | |
+ #ax3.plot(shear_strain[c], p_min[c], '-', color=color[c], alpha=0.5) | |
+ | |
+ | |
+ ax3.fill_between(shear_strain[c], p_min[c], p_max[c], | |
+ where=p_min[c]<=p_max[c], facecolor=color[c], | |
+ interpolate=True, alpha=0.5) | |
#ax4.plot(shear_strain[c][1:], f_n_mean[c][1:], '-' + color[c], | |
#label='$c$ = %.2f' % (cvals[c-1]), linewidth=2) | |
#ax4.plot(shear_strain[c][1:], f_n_max[c][1:], '--' + color[c]) | |
#label='$c$ = %.2f' % (cvals[c-1]), linewidth=2) | |
-#ax4.set_xlabel('Shear strain $\\gamma$ [-]') | |
-if zflow or pressures: | |
- ax3.set_xlabel('Shear strain $\\gamma$ [-]') | |
-else: | |
- ax2.set_xlabel('Shear strain $\\gamma$ [-]') | |
+ #ax4.set_xlabel('Shear strain $\\gamma$ [-]') | |
+ if zflow or pressures: | |
+ ax3.set_xlabel('Shear strain $\\gamma$ [-]') | |
+ else: | |
+ ax2.set_xlabel('Shear strain $\\gamma$ [-]') | |
-ax1.set_ylabel('Shear friction $\\tau/\\sigma_0$ [-]') | |
-#ax1.set_ylabel('Shear stress $\\tau$ [kPa]') | |
-ax2.set_ylabel('Dilation $\\Delta h/(2r)$ [-]') | |
-if zflow: | |
- ax3.set_ylabel('$\\boldsymbol{v}_\\text{f}^z h$ [ms$^{-1}$]') | |
-if pressures: | |
- ax3.set_ylabel('Mean fluid pressure $\\bar{p}_\\text{f}$ [kPa]') | |
-#ax4.set_ylabel('Particle contact force $||\\boldsymbol{f}_\\text{p}||$ [N]') | |
+ if c == 0: | |
+ ax1.set_ylabel('Shear friction $\\tau/\\sigma_0$ [-]') | |
+ #ax1.set_ylabel('Shear stress $\\tau$ [kPa]') | |
+ ax2.set_ylabel('Dilation $\\Delta h/(2r)$ [-]') | |
+ if zflow: | |
+ ax3.set_ylabel('$\\boldsymbol{v}_\\text{f}^z h$ [ms$^{-1}$]') | |
+ if pressures: | |
+ ax3.set_ylabel('Mean fluid pressure $\\bar{p}_\\text{f}$ [kPa]') | |
+ #ax4.set_ylabel('Particle contact force $||\\boldsymbol{f}_\\text{p}||… | |
+ | |
+ #ax1.set_xlim([200,300]) | |
+ #ax3.set_ylim([595,608]) | |
+ | |
+ plt.setp(ax1.get_xticklabels(), visible=False) | |
+ if zflow or pressures: | |
+ plt.setp(ax2.get_xticklabels(), visible=False) | |
+ #plt.setp(ax2.get_xticklabels(), visible=False) | |
+ #plt.setp(ax3.get_xticklabels(), visible=False) | |
+ | |
+ ''' | |
+ ax1.grid() | |
+ ax2.grid() | |
+ if zflow or pressures: | |
+ ax3.grid() | |
+ #ax4.grid() | |
+ ''' | |
+ | |
+ if c == 0: # left | |
+ # remove box at top and right | |
+ ax1.spines['top'].set_visible(False) | |
+ ax1.spines['bottom'].set_visible(False) | |
+ ax1.spines['right'].set_visible(False) | |
+ #ax1.spines['left'].set_visible(True) | |
+ # remove ticks at top and right | |
+ ax1.get_xaxis().set_ticks_position('none') | |
+ ax1.get_yaxis().set_ticks_position('none') | |
+ ax1.get_yaxis().tick_left() | |
+ | |
+ # remove box at top and right | |
+ ax2.spines['top'].set_visible(False) | |
+ ax2.spines['right'].set_visible(False) | |
+ ax2.spines['bottom'].set_visible(False) | |
+ # remove ticks at top and right | |
+ ax2.get_xaxis().set_ticks_position('none') | |
+ ax2.get_yaxis().set_ticks_position('none') | |
+ ax2.get_yaxis().tick_left() | |
+ | |
+ # remove box at top and right | |
+ ax3.spines['top'].set_visible(False) | |
+ ax3.spines['right'].set_visible(False) | |
+ # remove ticks at top and right | |
+ ax3.get_xaxis().set_ticks_position('none') | |
+ ax3.get_yaxis().set_ticks_position('none') | |
+ ax3.get_xaxis().tick_bottom() | |
+ ax3.get_yaxis().tick_left() | |
+ | |
+ elif c == len(mu_f_vals)-1: # right | |
+ # remove box at top and right | |
+ ax1.spines['top'].set_visible(False) | |
+ ax1.spines['bottom'].set_visible(False) | |
+ ax1.spines['right'].set_visible(True) | |
+ ax1.spines['left'].set_visible(False) | |
+ # remove ticks at top and right | |
+ ax1.get_xaxis().set_ticks_position('none') | |
+ ax1.get_yaxis().set_ticks_position('none') | |
+ ax1.get_yaxis().tick_right() | |
+ | |
+ # remove box at top and right | |
+ ax2.spines['top'].set_visible(False) | |
+ ax2.spines['right'].set_visible(True) | |
+ ax2.spines['bottom'].set_visible(False) | |
+ ax2.spines['left'].set_visible(False) | |
+ # remove ticks at top and right | |
+ ax2.get_xaxis().set_ticks_position('none') | |
+ ax2.get_yaxis().set_ticks_position('none') | |
+ #ax2.get_yaxis().tick_left() | |
+ ax2.get_yaxis().tick_right() | |
+ | |
+ # remove box at top and right | |
+ ax3.spines['top'].set_visible(False) | |
+ ax3.spines['right'].set_visible(True) | |
+ ax3.spines['left'].set_visible(False) | |
+ # remove ticks at top and right | |
+ ax3.get_xaxis().set_ticks_position('none') | |
+ ax3.get_yaxis().set_ticks_position('none') | |
+ ax3.get_xaxis().tick_bottom() | |
+ #ax3.get_yaxis().tick_left() | |
+ ax3.get_yaxis().tick_right() | |
+ | |
+ else: # middle | |
+ # remove box at top and right | |
+ ax1.spines['top'].set_visible(False) | |
+ ax1.spines['bottom'].set_visible(False) | |
+ ax1.spines['right'].set_visible(False) | |
+ ax1.spines['left'].set_visible(False) | |
+ # remove ticks at top and right | |
+ ax1.get_xaxis().set_ticks_position('none') | |
+ ax1.get_yaxis().set_ticks_position('none') | |
+ #ax1.get_yaxis().tick_left() | |
+ plt.setp(ax1.get_yticklabels(), visible=False) | |
+ | |
+ # remove box at top and right | |
+ ax2.spines['top'].set_visible(False) | |
+ ax2.spines['right'].set_visible(False) | |
+ ax2.spines['bottom'].set_visible(False) | |
+ ax2.spines['left'].set_visible(False) | |
+ # remove ticks at top and right | |
+ ax2.get_xaxis().set_ticks_position('none') | |
+ ax2.get_yaxis().set_ticks_position('none') | |
+ #ax2.get_yaxis().tick_left() | |
+ plt.setp(ax2.get_yticklabels(), visible=False) | |
+ | |
+ # remove box at top and right | |
+ ax3.spines['top'].set_visible(False) | |
+ ax3.spines['right'].set_visible(False) | |
+ ax3.spines['left'].set_visible(False) | |
+ # remove ticks at top and right | |
+ ax3.get_xaxis().set_ticks_position('none') | |
+ ax3.get_yaxis().set_ticks_position('none') | |
+ ax3.get_xaxis().tick_bottom() | |
+ #ax3.get_yaxis().tick_left() | |
+ plt.setp(ax3.get_yticklabels(), visible=False) | |
+ | |
+ | |
+ # vertical grid lines | |
+ ax1.get_xaxis().grid(True, linestyle=':', linewidth=0.5) | |
+ ax2.get_xaxis().grid(True, linestyle=':', linewidth=0.5) | |
+ ax3.get_xaxis().grid(True, linestyle=':', linewidth=0.5) | |
+ | |
+ | |
+ # horizontal grid lines | |
+ ax1.get_yaxis().grid(True, linestyle=':', linewidth=0.5) | |
+ ax2.get_yaxis().grid(True, linestyle=':', linewidth=0.5) | |
+ ax3.get_yaxis().grid(True, linestyle=':', linewidth=0.5) | |
+ | |
+ ax1.set_title(label) | |
+ #ax1.legend(loc='best') | |
+ #legend_alpha=0.5 | |
+ #ax1.legend(loc='upper right', prop={'size':18}, fancybox=True, | |
+ #framealpha=legend_alpha) | |
+ #ax2.legend(loc='lower right', prop={'size':18}, fancybox=True, | |
+ #framealpha=legend_alpha) | |
+ #if zflow or pressures: | |
+ #ax3.legend(loc='upper right', prop={'size':18}, fancybox=True, | |
+ #framealpha=legend_alpha) | |
+ #ax4.legend(loc='best', prop={'size':18}, fancybox=True, | |
+ #framealpha=legend_alpha) | |
-#ax1.set_xlim([200,300]) | |
-#ax3.set_ylim([595,608]) | |
+ #ax1.set_xlim([0.0, 0.09]) | |
+ #ax2.set_xlim([0.0, 0.09]) | |
+ #ax2.set_xlim([0.0, 0.2]) | |
-plt.setp(ax1.get_xticklabels(), visible=False) | |
-if zflow or pressures: | |
- plt.setp(ax2.get_xticklabels(), visible=False) | |
-#plt.setp(ax2.get_xticklabels(), visible=False) | |
-#plt.setp(ax3.get_xticklabels(), visible=False) | |
+ #ax1.set_ylim([-7, 45]) | |
+ ax1.set_xlim([0.0, 1.0]) | |
+ ax1.set_ylim([0.0, 1.0]) | |
+ ax2.set_ylim([0.0, 1.0]) | |
+ ax3.set_ylim([-200., 200.]) | |
-''' | |
-ax1.grid() | |
-ax2.grid() | |
-if zflow or pressures: | |
- ax3.grid() | |
-#ax4.grid() | |
-''' | |
- | |
- | |
-# remove box at top and right | |
-ax1.spines['top'].set_visible(False) | |
-ax1.spines['bottom'].set_visible(False) | |
-ax1.spines['right'].set_visible(False) | |
-#ax1.spines['left'].set_visible(True) | |
-# remove ticks at top and right | |
-ax1.get_xaxis().set_ticks_position('none') | |
-ax1.get_yaxis().set_ticks_position('none') | |
-ax1.get_yaxis().tick_left() | |
-ax1.get_xaxis().grid(True, linestyle='--', linewidth=0.5) | |
- | |
-# remove box at top and right | |
-ax2.spines['top'].set_visible(False) | |
-ax2.spines['right'].set_visible(False) | |
-ax2.spines['bottom'].set_visible(False) | |
-# remove ticks at top and right | |
-ax2.get_xaxis().set_ticks_position('none') | |
-ax2.get_yaxis().set_ticks_position('none') | |
-ax2.get_yaxis().tick_left() | |
-ax2.get_xaxis().grid(True, linestyle='--', linewidth=0.5) | |
- | |
-# remove box at top and right | |
-ax3.spines['top'].set_visible(False) | |
-ax3.spines['right'].set_visible(False) | |
-# remove ticks at top and right | |
-ax3.get_xaxis().set_ticks_position('none') | |
-ax3.get_yaxis().set_ticks_position('none') | |
-ax3.get_xaxis().tick_bottom() | |
-ax3.get_yaxis().tick_left() | |
-ax3.get_xaxis().grid(True, linestyle='--', linewidth=0.5) | |
- | |
-ax1.legend(loc='best') | |
-#legend_alpha=0.5 | |
-#ax1.legend(loc='upper right', prop={'size':18}, fancybox=True, | |
- #framealpha=legend_alpha) | |
-#ax2.legend(loc='lower right', prop={'size':18}, fancybox=True, | |
- #framealpha=legend_alpha) | |
-#if zflow or pressures: | |
- #ax3.legend(loc='upper right', prop={'size':18}, fancybox=True, | |
- #framealpha=legend_alpha) | |
-#ax4.legend(loc='best', prop={'size':18}, fancybox=True, | |
- #framealpha=legend_alpha) | |
- | |
-#ax1.set_xlim([0.0, 0.09]) | |
-#ax2.set_xlim([0.0, 0.09]) | |
-#ax2.set_xlim([0.0, 0.2]) | |
- | |
-#ax1.set_ylim([-7, 45]) | |
-ax2.set_ylim([0.0, 0.8]) | |
-#ax1.set_ylim([0.0, 1.0]) | |
-#if pressures: | |
- #ax3.set_ylim([-1400, 900]) | |
- #ax3.set_ylim([-200, 200]) | |
- #ax3.set_xlim([0.0, 0.09]) | |
- | |
-plt.tight_layout() | |
+ #ax1.set_ylim([0.0, 1.0]) | |
+ #if pressures: | |
+ #ax3.set_ylim([-1400, 900]) | |
+ #ax3.set_ylim([-200, 200]) | |
+ #ax3.set_xlim([0.0, 0.09]) | |
+ | |
+#plt.tight_layout() | |
#plt.subplots_adjust(hspace=0.05) | |
plt.subplots_adjust(hspace=0.15) | |
#filename = 'shear-' + str(int(sigma0/1000.0)) + 'kPa-stress-dilation.pdf' |