timproved plots - sphere - GPU-based 3D discrete element method algorithm with … | |
git clone git://src.adamsgaard.dk/sphere | |
Log | |
Files | |
Refs | |
LICENSE | |
--- | |
commit 09657d844261b3e1730f2fb2132d4dd66dba5e5f | |
parent f747e43ef3abf99c7365e9eb6451f7bbd8828304 | |
Author: Anders Damsgaard <[email protected]> | |
Date: Tue, 7 Oct 2014 14:33:19 +0200 | |
improved plots | |
Diffstat: | |
A python/shear-results-internals.py | 339 +++++++++++++++++++++++++++++… | |
M python/shear-results-strain.py | 4 ++-- | |
M python/shear-results.py | 123 +++++++++++++++++++++++++----… | |
3 files changed, 443 insertions(+), 23 deletions(-) | |
--- | |
diff --git a/python/shear-results-internals.py b/python/shear-results-internals… | |
t@@ -0,0 +1,339 @@ | |
+#!/usr/bin/env python | |
+import matplotlib | |
+matplotlib.use('Agg') | |
+matplotlib.rcParams.update({'font.size': 18, 'font.family': 'serif'}) | |
+matplotlib.rc('text', usetex=True) | |
+matplotlib.rcParams['text.latex.preamble']=[r"\usepackage{amsmath}"] | |
+import shutil | |
+ | |
+import os | |
+import numpy | |
+import sphere | |
+from permeabilitycalculator import * | |
+import matplotlib.pyplot as plt | |
+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 | |
+ | |
+sigma0 = float(sys.argv[1]) | |
+#c_grad_p = 1.0 | |
+c_grad_p = float(sys.argv[2]) | |
+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' | |
+sim = sphere.sim(sid, fluid=True) | |
+sim.readfirst(verbose=False) | |
+ | |
+# particle z positions | |
+zpos_p = numpy.zeros((len(steps), sim.np)) | |
+ | |
+# cell midpoint cell positions | |
+zpos_c = numpy.zeros((len(steps), 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 | |
+ | |
+# particle x displacements | |
+xdisp = numpy.zeros((len(steps), sim.np)) | |
+ | |
+# particle z velocity | |
+v_z_p = numpy.zeros((len(steps), sim.np)) | |
+ | |
+# fluid z velocity | |
+v_z_f = numpy.zeros((len(steps), sim.num[0], sim.num[1], sim.num[2])) | |
+ | |
+# pressure - hydrostatic pressure | |
+dev_p = numpy.zeros((len(steps), sim.num[2])) | |
+ | |
+# mean per-particle values | |
+v_z_p_bar = numpy.zeros((len(steps), sim.num[2])) | |
+v_z_f_bar = numpy.zeros((len(steps), sim.num[2])) | |
+ | |
+# particle-fluid force per particle | |
+f_pf = numpy.zeros_like(xdisp) | |
+ | |
+# pressure - hydrostatic pressure | |
+#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])) | |
+ | |
+shear_strain = numpy.zeros(len(steps)) | |
+ | |
+s = 0 | |
+for step_str in steps: | |
+ | |
+ step = int(step_str) | |
+ | |
+ if os.path.isfile('../output/' + sid + '.status.dat'): | |
+ | |
+ for substep in numpy.arange(nsteps_avg): | |
+ | |
+ if step + substep > sim.status(): | |
+ raise Exception( | |
+ 'Simulation step %d not available (sim.status = %d).' | |
+ % (step, sim.status())) | |
+ | |
+ sim.readstep(step + substep, verbose=False) | |
+ | |
+ zpos_p[s,:] += sim.x[:,2]/nsteps_avg | |
+ | |
+ xdisp[s,:] += sim.xyzsum[:,0]/nsteps_avg | |
+ v_z_p[s,:] += sim.vel[:,2]/nsteps_avg | |
+ | |
+ ''' | |
+ for i in numpy.arange(sim.np): | |
+ f_pf[s,i] += \ | |
+ sim.f_sum[i].dot(sim.f_sum[i])/nsteps_avg | |
+ ''' | |
+ f_pf[s,:] += sim.f_sum[:,2] | |
+ | |
+ dev_p[s,:] += \ | |
+ numpy.average(numpy.average(sim.p_f, axis=0), axis=0)\ | |
+ /nsteps_avg | |
+ | |
+ v_z_f[s,:] += sim.v_f[:,:,:,2]/nsteps_avg | |
+ | |
+ v_z_f_bar[s,:] += \ | |
+ numpy.average(numpy.average(sim.v_f[:,:,:,2], axis=0), axi… | |
+ /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 | |
+ for iz in numpy.arange(sim.num[2]): | |
+ z_bot = iz*dz | |
+ z_top = (iz+1)*dz | |
+ I = numpy.nonzero((zpos_p[s,:] >= z_bot) & (zpos_p[s,:] < z_top)) | |
+ if len(I) > 0: | |
+ xdisp_mean[s,iz] = numpy.mean(xdisp[s,I]) | |
+ v_z_p_bar[s,iz] = numpy.mean(v_z_p[s,I]) | |
+ f_pf_mean[s,iz] = numpy.mean(f_pf[s,I]) | |
+ | |
+ else: | |
+ print(sid + ' not found') | |
+ s += 1 | |
+ | |
+#fig = plt.figure(figsize=(8,4*(len(steps))+1)) | |
+#fig = plt.figure(figsize=(8,5*(len(steps))+1)) | |
+fig = plt.figure(figsize=(16,5*(len(steps))+1)) | |
+ | |
+ax = [] | |
+for s in numpy.arange(len(steps)): | |
+ | |
+ strain_str = 'Shear strain $\\gamma = %.3f$' % (shear_strain[s]) | |
+ | |
+ n = 7 | |
+ if s == 0: | |
+ ax.append(plt.subplot(len(steps), n-1, s*(n-1)+1)) # 0: xdisp | |
+ ax.append(plt.subplot(len(steps), n-1, s*(n-1)+2, sharey=ax[s*n+0])) #… | |
+ ax.append(ax[s*n+1].twiny()) # 2: dphi/dt | |
+ ax.append(plt.subplot(len(steps), n-1, s*(n-1)+3, sharey=ax[s*n+0])) #… | |
+ ax.append(plt.subplot(len(steps), n-1, s*(n-1)+4, sharey=ax[s*n+0])) #… | |
+ ax.append(plt.subplot(len(steps), n-1, s*(n-1)+5, sharey=ax[s*n+0])) #… | |
+ ax.append(plt.subplot(len(steps), n-1, s*(n-1)+6, sharey=ax[s*n+0])) #… | |
+ else: | |
+ ax.append(plt.subplot(len(steps), n-1, s*(n-1)+1, sharex=ax[0])) # 0: … | |
+ ax.append(plt.subplot(len(steps), n-1, s*(n-1)+2, sharey=ax[s*n+0], | |
+ sharex=ax[1])) # 1: phi | |
+ ax.append(ax[s*n+1].twiny()) # 2: dphi/dt | |
+ ax.append(plt.subplot(len(steps), n-1, s*(n-1)+3, sharey=ax[s*n+0], | |
+ sharex=ax[3])) # 3: v_z^p | |
+ ax.append(plt.subplot(len(steps), n-1, s*(n-1)+4, sharey=ax[s*n+0], | |
+ sharex=ax[4])) # 4: p_f | |
+ ax.append(plt.subplot(len(steps), n-1, s*(n-1)+5, sharey=ax[s*n+0], | |
+ sharex=ax[5])) # 5: f_pf_z | |
+ ax.append(plt.subplot(len(steps), n-1, s*(n-1)+6, sharey=ax[s*n+0], | |
+ 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*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') | |
+ | |
+ #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*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+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') | |
+ y_top = sim.w_x[0] | |
+ x_top = sim.p_f[0,0,-1] | |
+ y_bot = 0.0 | |
+ x_bot = x_top + (y_top - y_bot)*sim.rho*numpy.abs(sim.g[2]) | |
+ ax[s*n+4].plot([x_top/1000.0, x_bot/1000.0], [y_top, y_bot], '--', color='… | |
+ #ax[s*n+1].set_title(strain_str) | |
+ #ax[s*n+1].set_title(' ') | |
+ | |
+ # remove particles with 0.0 pressure force | |
+ I = numpy.nonzero(numpy.abs(f_pf[s]) > .01) | |
+ f_pf_nonzero = f_pf[s][I] | |
+ zpos_p_nonzero = zpos_p[s][I] | |
+ I = numpy.nonzero(numpy.abs(f_pf_mean[s]) > .01) | |
+ 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*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*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+2].plot([0.0,0.0], [0.0, sim.L[2]], '--', color='k') | |
+ | |
+ | |
+ #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])) | |
+ #plt.loglog(dpdz[c], K[c], 'o-', label='$c$ = %.2f' % (cvals[c])) | |
+ | |
+ #for x in numpy.arange(sim.num[0]): | |
+ #for y in numpy.arange(sim.num[1]): | |
+ #ax[s*n+2].plot(v_z_f[s,x,y,:], zpos_c[s], ',', color = '#888888') | |
+ | |
+ | |
+ #phicolor = '#888888' | |
+ #ax[s*n+3].plot(phi_bar[s], zpos_c[s], '-', color = phicolor) | |
+ #for tl in ax[s*n+3].get_xticklabels(): | |
+ #tl.set_color(phicolor) | |
+ | |
+ max_z = numpy.max(zpos_p) | |
+ ax[s*n+0].set_ylim([0, max_z]) | |
+ #ax[s*n+1].set_ylim([0, max_z]) | |
+ #ax[s*n+2].set_ylim([0, max_z]) | |
+ | |
+ #ax[s*n+0].set_xlim([-0.01,0.01]) | |
+ #ax[s*n+0].set_xlim([-0.005,0.005]) | |
+ #ax[s*n+0].set_xlim([-0.25,0.75]) | |
+ ax[s*n+4].set_xlim([595,625]) # p_f | |
+ #ax[s*n+2].set_xlim([-0.0005,0.0005]) | |
+ #ax[s*n+2].set_xlim([-0.08,0.08]) | |
+ | |
+ #ax[s*4+1].set_xlim([0.15, 0.46]) # f_pf | |
+ #ax[s*n+1].set_xlim([0.235, 0.409]) # f_pf | |
+ | |
+ ax[s*n+1].set_xlim([0.33, 0.6]) # phi | |
+ ax[s*n+2].set_xlim([-0.09, 0.035]) # dphi/dt | |
+ | |
+ | |
+ #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])) | |
+ #plt.loglog(dpdz[c], K[c], 'o-', label='$c$ = %.2f' % (cvals[c])) | |
+ | |
+ ax[s*n+0].set_ylabel('Vertical position $z$ [m]') | |
+ ax[s*n+0].set_xlabel('$\\boldsymbol{x}^x_\\text{p}$ [m]') | |
+ ax[s*n+1].set_xlabel('$\\bar{\\phi}$ [-] (solid)') | |
+ ax[s*n+2].set_xlabel('$\\delta \\bar{\\phi}/\\delta t$ [-] (dashed)') | |
+ ax[s*n+3].set_xlabel('$\\boldsymbol{v}^z_\\text{p}$ [cms$^-1$]') | |
+ ax[s*n+4].set_xlabel('$\\bar{p_\\text{f}}$ [kPa]') | |
+ ax[s*n+5].set_xlabel('$\\boldsymbol{f}^z_\\text{pf}$ [N]') | |
+ ax[s*n+6].set_xlabel('$\\bar{\\boldsymbol{v}}^z_\\text{f}$ [cms$^-1$]') | |
+ | |
+ plt.setp(ax[s*n+1].get_yticklabels(), visible=False) | |
+ plt.setp(ax[s*n+2].get_yticklabels(), visible=False) | |
+ plt.setp(ax[s*n+3].get_yticklabels(), visible=False) | |
+ plt.setp(ax[s*n+4].get_yticklabels(), visible=False) | |
+ plt.setp(ax[s*n+5].get_yticklabels(), visible=False) | |
+ plt.setp(ax[s*n+6].get_yticklabels(), visible=False) | |
+ | |
+ #nbins = 4 | |
+ #ax[s*n+0].get_xaxis().set_major_locator(MaxNLocator(nbins=nbins)) | |
+ #ax[s*n+1].get_xaxis().set_major_locator(MaxNLocator(nbins=nbins)) | |
+ #ax[s*n+2].get_xaxis().set_major_locator(MaxNLocator(nbins=nbins)) | |
+ | |
+ plt.setp(ax[s*n+0].xaxis.get_majorticklabels(), rotation=90) | |
+ plt.setp(ax[s*n+1].xaxis.get_majorticklabels(), rotation=90) | |
+ plt.setp(ax[s*n+2].xaxis.get_majorticklabels(), rotation=90) | |
+ plt.setp(ax[s*n+3].xaxis.get_majorticklabels(), rotation=90) | |
+ plt.setp(ax[s*n+4].xaxis.get_majorticklabels(), rotation=90) | |
+ plt.setp(ax[s*n+5].xaxis.get_majorticklabels(), rotation=90) | |
+ plt.setp(ax[s*n+6].xaxis.get_majorticklabels(), rotation=90) | |
+ | |
+ #if s == 0: | |
+ #y = 0.95 | |
+ #if s == 1: | |
+ #y = 0.55 | |
+ | |
+ #plt.ticklabel_format(style='sci', axis='x', scilimits=(0,0)) | |
+ #ax[s*n+0].ticklabel_format(style='sci', axis='x', scilimits=(-3,3)) | |
+ #ax[s*n+1].ticklabel_format(style='sci', axis='x', scilimits=(-3,3)) | |
+ #ax[s*n+2].ticklabel_format(style='sci', axis='x', scilimits=(-3,3)) | |
+ | |
+ #fig.text(0.1, y, strain_str, horizontalalignment='left', fontsize=22) | |
+ #ax[s*4+0].annotate(strain_str, xytext=(0,1.1), textcoords='figure fractio… | |
+ #horizontalalignment='left', fontsize=22) | |
+ #plt.text(0.05, 1.06, strain_str, horizontalalignment='left', fontsize=22, | |
+ #transform=ax[s*n+0].transAxes) | |
+ #ax[s*4+0].set_title(strain_str) | |
+ ax[s*n+0].set_title('a') | |
+ ax[s*n+1].set_title('b') | |
+ ax[s*n+3].set_title('c') | |
+ ax[s*n+4].set_title('d') | |
+ ax[s*n+5].set_title('e') | |
+ ax[s*n+6].set_title('f') | |
+ | |
+ ax[s*n+0].grid() | |
+ ax[s*n+1].grid() | |
+ #ax[s*n+2].grid() | |
+ ax[s*n+3].grid() | |
+ ax[s*n+4].grid() | |
+ ax[s*n+5].grid() | |
+ ax[s*n+6].grid() | |
+ #ax1.legend(loc='lower right', prop={'size':18}) | |
+ #ax2.legend(loc='lower right', prop={'size':18}) | |
+ | |
+ strain_str = 'Shear strain $\\gamma = %.3f$' % (shear_strain[s]) | |
+ #fig.text(0.1, y, strain_str, horizontalalignment='left', fontsize=22) | |
+ #ax[s*4+0].annotate(strain_str, xytext=(0,1.1), textcoords='figure fractio… | |
+ #horizontalalignment='left', fontsize=22) | |
+ plt.text(0.05, 1.06, strain_str, horizontalalignment='left', fontsize=22, | |
+ transform=ax[s*n+0].transAxes) | |
+ | |
+#plt.title(' ') | |
+plt.MaxNLocator(nbins=4) | |
+ #ax1.legend(loc='lower right', prop={'size':18}) | |
+ #ax2.legend(loc='lower right', prop={'size':18}) | |
+ | |
+#plt.title(' ') | |
+#plt.MaxNLocator(nbins=4) | |
+#plt.subplots_adjust(wspace = .05) | |
+#plt.subplots_adjust(hspace = 1.05) | |
+plt.tight_layout() | |
+#plt.MaxNLocator(nbins=4) | |
+ | |
+filename = 'shear-' + str(int(sigma0/1000.0)) + 'kPa-internals.pdf' | |
+plt.savefig(filename) | |
+shutil.copyfile(filename, '/home/adc/articles/own/2-org/' + filename) | |
+print(filename) | |
diff --git a/python/shear-results-strain.py b/python/shear-results-strain.py | |
t@@ -96,7 +96,7 @@ for s in numpy.arange(len(cvals)): | |
if cvals[s] == 'dry': | |
legend = 'dry' | |
else: | |
- legend = 'c = ' + str(cvals[s]) | |
+ 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 = legen… | |
t@@ -134,7 +134,7 @@ for s in numpy.arange(len(cvals)): | |
#ax2.legend(loc='lower right', prop={'size':18}) | |
legend_alpha=0.5 | |
-ax[0].legend(loc='best', prop={'size':18}, fancybox=True, framealpha=legend_al… | |
+ax[0].legend(loc='lower right', prop={'size':18}, fancybox=True, framealpha=le… | |
ax[0].grid() | |
plt.tight_layout() | |
plt.subplots_adjust(wspace = .05) | |
diff --git a/python/shear-results.py b/python/shear-results.py | |
t@@ -19,8 +19,79 @@ sigma0 = float(sys.argv[1]) | |
cvals = [1.0, 0.1] | |
#cvals = [1.0] | |
+# return a smoothed version of in. The returned array is smaller than the | |
+# original input array | |
+''' | |
+def smooth(in_arr, plus_minus_steps): | |
+ out_arr = numpy.zeros(in_arr.size - 2*plus_minus_steps + 1) | |
+ s = 0 | |
+ for i in numpy.arange(in_arr.size): | |
+ if i >= plus_minus_steps and i < plus_minus_steps: | |
+ for i in numpy.arange(-plus_minus_steps, plus_minus_steps+1): | |
+ out_arr[s] += in_arr[s+i]/(2.0*plus_minus_steps) | |
+ s += 1 | |
+''' | |
+ | |
+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 = 11 | |
+ | |
shear_strain = [[], [], []] | |
friction = [[], [], []] | |
+friction_smooth = [[], [], []] | |
dilation = [[], [], []] | |
p_min = [[], [], []] | |
p_mean = [[], [], []] | |
t@@ -39,6 +110,7 @@ sim.visualize('shear') | |
shear_strain[0] = sim.shear_strain | |
#shear_strain[0] = numpy.arange(sim.status()+1) | |
friction[0] = sim.tau/sim.sigma_eff | |
+friction_smooth[0] = smooth(friction[0], smooth_window) | |
dilation[0] = sim.dilation | |
f_n_mean[0] = numpy.zeros_like(shear_strain[0]) | |
t@@ -64,6 +136,7 @@ for c in numpy.arange(1,len(cvals)+1): | |
shear_strain[c] = numpy.zeros(sim.status()) | |
friction[c] = numpy.zeros_like(shear_strain[c]) | |
dilation[c] = numpy.zeros_like(shear_strain[c]) | |
+ friction_smooth[c] = numpy.zeros_like(shear_strain[c]) | |
sim.readlast(verbose=False) | |
sim.visualize('shear') | |
t@@ -71,6 +144,7 @@ for c in numpy.arange(1,len(cvals)+1): | |
#shear_strain[c] = numpy.arange(sim.status()+1) | |
friction[c] = sim.tau/sim.sigma_eff | |
dilation[c] = sim.dilation | |
+ friction_smooth[c] = smooth(friction[c], smooth_window) | |
# fluid pressures and particle forces | |
p_mean[c] = numpy.zeros_like(shear_strain[c]) | |
t@@ -100,31 +174,37 @@ for c in numpy.arange(1,len(cvals)+1): | |
#fig = plt.figure(figsize=(8,8)) # (w,h) | |
+fig = plt.figure(figsize=(8,10)) | |
#fig = plt.figure(figsize=(8,12)) | |
-fig = plt.figure(figsize=(8,16)) | |
+#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)) | |
-ax1 = plt.subplot(411) | |
-ax2 = plt.subplot(412, sharex=ax1) | |
-ax3 = plt.subplot(413, sharex=ax1) | |
-ax4 = plt.subplot(414, sharex=ax1) | |
-ax1.plot(shear_strain[0], friction[0], label='dry') | |
+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 | |
+#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') | |
-ax4.plot(shear_strain[0], f_n_mean[0], '-', label='dry', color='blue') | |
-ax4.plot(shear_strain[0], f_n_max[0], '--', color='blue') | |
+#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'] | |
for c in numpy.arange(1,len(cvals)+1): | |
- ax1.plot(shear_strain[c][1:], friction[c][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])) | |
ax2.plot(shear_strain[c][1:], dilation[c][1:], \ | |
label='$c$ = %.2f' % (cvals[c-1]), linewidth=2) | |
+ ''' | |
alpha = 0.5 | |
ax3.plot(shear_strain[c][1:], p_max[c][1:], '-' + color[c], alpha=alpha) | |
ax3.plot(shear_strain[c][1:], p_mean[c][1:], '-' + color[c], \ | |
t@@ -139,34 +219,35 @@ for c in numpy.arange(1,len(cvals)+1): | |
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$ [-]') | |
+#ax4.set_xlabel('Shear strain $\\gamma$ [-]') | |
ax1.set_ylabel('Shear friction $\\tau/\\sigma\'$ [-]') | |
ax2.set_ylabel('Dilation $\\Delta h/(2r)$ [-]') | |
-ax3.set_ylabel('Fluid pressure $p_\\text{f}$ [kPa]') | |
-ax4.set_ylabel('Particle contact force $||\\boldsymbol{f}_\\text{p}||$ [N]') | |
+#ax3.set_ylabel('Fluid pressure $p_\\text{f}$ [kPa]') | |
+#ax4.set_ylabel('Particle contact force $||\\boldsymbol{f}_\\text{p}||$ [N]') | |
#ax1.set_xlim([200,300]) | |
-ax3.set_ylim([595,608]) | |
+#ax3.set_ylim([595,608]) | |
plt.setp(ax1.get_xticklabels(), visible=False) | |
-plt.setp(ax2.get_xticklabels(), visible=False) | |
-plt.setp(ax3.get_xticklabels(), visible=False) | |
+#plt.setp(ax2.get_xticklabels(), visible=False) | |
+#plt.setp(ax3.get_xticklabels(), visible=False) | |
ax1.grid() | |
ax2.grid() | |
-ax3.grid() | |
-ax4.grid() | |
+#ax3.grid() | |
+#ax4.grid() | |
legend_alpha=0.5 | |
ax1.legend(loc='best', prop={'size':18}, fancybox=True, framealpha=legend_alph… | |
ax2.legend(loc='lower right', prop={'size':18}, fancybox=True, | |
framealpha=legend_alpha) | |
-ax3.legend(loc='lower right', prop={'size':18}, fancybox=True, | |
- framealpha=legend_alpha) | |
-ax4.legend(loc='best', prop={'size':18}, fancybox=True, | |
- framealpha=legend_alpha) | |
+#ax3.legend(loc='lower right', prop={'size':18}, fancybox=True, | |
+ #framealpha=legend_alpha) | |
+#ax4.legend(loc='best', prop={'size':18}, fancybox=True, | |
+ #framealpha=legend_alpha) | |
plt.tight_layout() | |
plt.subplots_adjust(hspace=0.05) |