tupdate histogram to log, plot creep snapshots in loop - 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 df5ef21559c8d9f76c721a15950a144290b3e4f8
(DIR) parent 14fb2bf075a6836438322eb5cfb2b70536a5c045
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Mon, 18 May 2015 13:17:31 +0200
update histogram to log, plot creep snapshots in loop
Diffstat:
M python/halfshear-darcy-creep-dynam… | 475 +++++++++----------------------
M python/sphere.py | 1 +
2 files changed, 139 insertions(+), 337 deletions(-)
---
(DIR) diff --git a/python/halfshear-darcy-creep-dynamics.py b/python/halfshear-darcy-creep-dynamics.py
t@@ -29,6 +29,7 @@ matplotlib.rc('grid', linestyle=':', linewidth=0.2)
outformat='pdf'
+plotContacts=False
#sids = ['halfshear-darcy-sigma0=10000.0-k_c=2e-16-mu=2.08e-07-ss=2000.0-A=4000.0-f=0.2']
sids = ['halfshear-darcy-sigma0=10000.0-k_c=2e-16-mu=2.08e-07-ss=2000.0-A=4375.0-f=0.2']
t@@ -59,6 +60,20 @@ f_n_max = 50 # for force chain plots
N = numpy.zeros_like(steps, dtype=numpy.float64)
t = numpy.zeros_like(steps, dtype=numpy.float64)
+scatter=False
+
+# insert plot positions
+Lx=[.17, .37, .65]
+Ly=[.3, .7, .4]
+dx=.23; dy=.23
+xytext_x=[Lx[0]+.15*dx, Lx[1]+.6*dx, Lx[2]+1.15*dx]
+xytext_y=[Ly[0]+dy+.07, Ly[1]+.15, .2]
+fc_x = [Lx[0]-.007, Lx[1]+dx+.08, Lx[2]-.007]
+fc_y = [Ly[0]-.07*dy, Ly[1]+.035, Ly[2]-.7*dy]
+fc_dx = [dx, dx, dx]
+fc_dy = [dy*.7, dy*.7, dy*.7]
+
+
s=0
for sid in sids:
t@@ -77,15 +92,14 @@ for sid in sids:
N[i] = sim.currentNormalStress('defined')
t[i] = sim.currentTime()
- '''
- outfolder='../img_out/' + sim.id() + '/'
- if not os.path.exists(outfolder):
- os.makedirs(outfolder)
- sim.plotContacts(outfolder=outfolder,
- alpha=0.9,
- lower_limit=lower_limit,
- upper_limit=upper_limit)
- '''
+ if plotContacts:
+ outfolder='../img_out/' + sim.id() + '/'
+ if not os.path.exists(outfolder):
+ os.makedirs(outfolder)
+ sim.plotContacts(outfolder=outfolder,
+ alpha=0.9,
+ lower_limit=lower_limit,
+ upper_limit=upper_limit)
if (step == plotsteps).any():
#contactdata.append(sim.plotContacts(return_data=True))
t@@ -150,338 +164,125 @@ for sid in sids:
ax1.locator_params(axis='y', nbins=4)
## Contact scatter plots
- dx=.23; dy=.23
# Scatter plot 1
- sc=0
- Lx=.17; Ly=.3;
- #xytext=(Lx+.5*dx, Ly+.5*dy)
- xytext=(Lx+.15*dx, Ly+dy+.07)
- xy=(ts[0]*scalingfactor, Ns[0]/1000.)
- #print xytext
- #print xy
- ax1.annotate('',
- xytext=xytext, textcoords='axes fraction',
- xy=xy, xycoords='data',
- arrowprops=dict(arrowstyle='->', connectionstyle='arc3'))
-
- axsc1 = fig.add_axes([Lx, Ly, dx, dy], polar=True)
- cs = axsc1.scatter(strikelists[sc], 90. - diplists[sc], marker='o',
- c=forcemagnitudes[sc],
- s=forcemagnitudes[sc]/f_n_maxs[sc]*5.,
- alpha=alphas[sc],
- edgecolors='none',
- vmin=f_n_maxs[sc]*lower_limit,
- vmax=f_n_maxs[sc]*upper_limit,
- cmap=matplotlib.cm.get_cmap('afmhot_r'))#,
- #norm=matplotlib.colors.LogNorm())
- # tick locations
- #thetaticks = numpy.arange(0,360,90)
- # set ticklabels location at 1.3 times the axes' radius
- #ax.set_thetagrids(thetaticks, frac=1.3)
- axsc1.set_xticklabels([])
- axsc1.set_yticklabels([])
-
- axsc1.set_title('\\textbf{1. Slow creep}', fontsize=7)
- if upper_limit < 1.0:
- cbar = plt.colorbar(cs, extend='max', fraction=0.035, pad=0.04)
- else:
- cbar = plt.colorbar(cs, fraction=0.045, pad=0.04)
- #cbar.set_label('$||\\boldsymbol{f}_n||$')
- cbar.set_label('$\\boldsymbol{f}_\\text{n}$ [N]')
- cbar.locator = matplotlib.ticker.MaxNLocator(nbins=4)
- cbar.update_ticks()
-
- # plot defined max compressive stress from tau/N ratio
- axsc1.scatter(0., numpy.degrees(numpy.arctan(taus[sc]/Ns[sc])),
- marker='+', c='none', edgecolor='red', s=100)
- axsc1.scatter(0., numpy.degrees(numpy.arctan(taus[sc]/Ns[sc])),
- marker='o', c='none', edgecolor='red', s=100)
- '''
- ax.sc1scatter(0., # actual stress
- numpy.degrees(numpy.arctan(
- self.shearStress('effective')/
- self.currentNormalStress('effective'))),
- marker='+', color='blue', s=300)
- '''
-
- axsc1.set_rmax(90)
- axsc1.set_rticks([])
- #axsc1.grid(False)
-
- # force chain plot
-
- axfc1 = fig.add_axes([Lx-0.007, Ly-0.7*dy, dx, dy*0.7])
-
- data = datalists[sc]
-
- # find the max. value of the normal force
- #f_n_max = numpy.amax(data[:,6])
-
- # specify the lower limit of force chains to do statistics on
- f_n_lim = lower_limit * f_n_max
-
- # find the indexes of these contacts
- I = numpy.nonzero(data[:,6] >= f_n_lim)
-
- #color = matplotlib.cm.spectral(data[:,6]/f_n_max)
- for i in I[0]:
-
- x1 = data[i,0]
- #y1 = data[i,1]
- z1 = data[i,2]
- x2 = data[i,3]
- #y2 = data[i,4]
- z2 = data[i,5]
- f_n = data[i,6]
-
- lw_max = 1.0
- if f_n >= f_n_max:
- lw = lw_max
- else:
- lw = (f_n - f_n_lim)/(f_n_max - f_n_lim)*lw_max
-
- #print lw
- axfc1.plot([x1,x2], [z1,z2], '-k', linewidth=lw)
- #axfc1.plot([x1,x2], [z1,z2], '-', linewidth=lw, color=forcemagnitudes[sc])
- #axfc1.plot([x1,x2], [z1,z2], '-', linewidth=lw, color=color)
+ for sc in range(len(Lx)):
+ xy=(ts[sc]*scalingfactor, Ns[sc]/1000.)
+ #print xytext
+ #print xy
+ ax1.annotate('',
+ xytext=(xytext_x[sc], xytext_y[sc]), textcoords='axes fraction',
+ xy=xy, xycoords='data',
+ arrowprops=dict(arrowstyle='->', connectionstyle='arc3'))
+
+ if scatter:
+ axsc1 = fig.add_axes([Lx[sc], Ly[sc], dx, dy], polar=True)
+ cs = axsc1.scatter(strikelists[sc], 90. - diplists[sc], marker='o',
+ c=forcemagnitudes[sc],
+ s=forcemagnitudes[sc]/f_n_maxs[sc]*5.,
+ alpha=alphas[sc],
+ edgecolors='none',
+ vmin=f_n_maxs[sc]*lower_limit,
+ vmax=f_n_maxs[sc]*upper_limit,
+ cmap=matplotlib.cm.get_cmap('afmhot_r'))#,
+ #norm=matplotlib.colors.LogNorm())
+ # tick locations
+ #thetaticks = numpy.arange(0,360,90)
+ # set ticklabels location at 1.3 times the axes' radius
+ #ax.set_thetagrids(thetaticks, frac=1.3)
+ axsc1.set_xticklabels([])
+ axsc1.set_yticklabels([])
+
+ if sc == 0:
+ title = '1. Slow creep'
+ elif sc == 1:
+ title = '2. Fast creep'
+ elif sc == 2:
+ title = '3. Slip'
+ else:
+ title = 'Unknown'
+ axsc1.set_title('\\textbf{' + title + '}', fontsize=7)
+ if upper_limit < 1.0:
+ cbar = plt.colorbar(cs, extend='max', fraction=0.035, pad=0.04)
+ else:
+ cbar = plt.colorbar(cs, fraction=0.045, pad=0.04)
+ #cbar.set_label('$||\\boldsymbol{f}_n||$')
+ cbar.set_label('$\\boldsymbol{f}_\\text{n}$ [N]')
+ cbar.locator = matplotlib.ticker.MaxNLocator(nbins=4)
+ cbar.update_ticks()
+
+ # plot defined max compressive stress from tau/N ratio
+ axsc1.scatter(0., numpy.degrees(numpy.arctan(taus[sc]/Ns[sc])),
+ marker='+', c='none', edgecolor='red', s=100)
+ axsc1.scatter(0., numpy.degrees(numpy.arctan(taus[sc]/Ns[sc])),
+ marker='o', c='none', edgecolor='red', s=100)
+ '''
+ ax.sc1scatter(0., # actual stress
+ numpy.degrees(numpy.arctan(
+ self.shearStress('effective')/
+ self.currentNormalStress('effective'))),
+ marker='+', color='blue', s=300)
+ '''
+
+ axsc1.set_rmax(90)
+ axsc1.set_rticks([])
+ #axsc1.grid(False)
- axfc1.spines['right'].set_visible(False)
- axfc1.spines['left'].set_visible(False)
- # Only show ticks on the left and bottom spines
- axfc1.xaxis.set_ticks_position('none')
- axfc1.yaxis.set_ticks_position('none')
- axfc1.set_xticklabels([])
- axfc1.set_yticklabels([])
- axfc1.set_xlim([numpy.min(data[I[0],0]), numpy.max(data[I[0],0])])
- axfc1.set_ylim([numpy.min(data[I[0],2]), numpy.max(data[I[0],2])])
- axfc1.set_aspect('equal')
-
-
-
- # Scatter plot 2
- sc=1
- Lx=.37; Ly=.7;
- #xytext=(Lx+.5*dx, Ly+.5*dy)
- #xytext=(Lx+.5*dx, Ly+.05)
- xytext=(Lx+.6*dx, Ly+.15)
- xy=(ts[sc]*scalingfactor, Ns[sc]/1000.)
- #print xytext
- #print xy
- ax1.annotate('',
- xytext=xytext, textcoords='axes fraction',
- xy=xy, xycoords='data',
- arrowprops=dict(arrowstyle='->', connectionstyle='arc3'))
-
- axsc2 = fig.add_axes([Lx, Ly, dx, dy], polar=True)
-
- cs = axsc2.scatter(strikelists[sc], 90. - diplists[sc], marker='o',
- c=forcemagnitudes[sc],
- s=forcemagnitudes[sc]/f_n_max*5.,
- alpha=alphas[sc],
- edgecolors='none',
- vmin=f_n_maxs[sc]*lower_limit,
- vmax=f_n_maxs[sc]*upper_limit,
- cmap=matplotlib.cm.get_cmap('afmhot_r'))#,
- #norm=matplotlib.colors.LogNorm())
- # tick locations
- #thetaticks = numpy.arange(0,360,90)
- # set ticklabels location at 1.3 times the axes' radius
- #ax.set_thetagrids(thetaticks, frac=1.3)
- axsc2.set_xticklabels([])
- axsc2.set_yticklabels([])
-
- axsc2.set_title('\\textbf{2. Fast creep}', fontsize=7)
- if upper_limit < 1.0:
- cbar = plt.colorbar(cs, extend='max', fraction=0.035, pad=0.04)
- else:
- cbar = plt.colorbar(cs, fraction=0.045, pad=0.04)
- cbar.set_label('$\\boldsymbol{f}_\\text{n}$ [N]')
- #cbar.set_label('Contact force [N]')
- cbar.locator = matplotlib.ticker.MaxNLocator(nbins=4)
- cbar.update_ticks()
-
- # plot defined max compressive stress from tau/N ratio
- axsc2.scatter(0., numpy.degrees(numpy.arctan(taus[sc]/Ns[sc])),
- marker='+', c='none', edgecolor='red', s=100)
- axsc2.scatter(0., numpy.degrees(numpy.arctan(taus[sc]/Ns[sc])),
- marker='o', c='none', edgecolor='red', s=100)
- '''
- ax.sc2scatter(0., # actual stress
- numpy.degrees(numpy.arctan(
- self.shearStress('effective')/
- self.currentNormalStress('effective'))),
- marker='+', color='blue', s=300)
- '''
-
- axsc2.set_rmax(90)
- axsc2.set_rticks([])
- #axsc2.grid(False)
-
- # force chain plot
-
- #axfc2 = fig.add_axes([Lx+dx+0.05, Ly+0.03, dx, dy*0.7])
- axfc2 = fig.add_axes([Lx+dx+0.08, Ly+0.035, dx, dy*0.7])
-
- data = datalists[sc]
-
- # find the max. value of the normal force
- #f_n_max = numpy.amax(data[:,6])
-
- # specify the lower limit of force chains to do statistics on
- f_n_lim = lower_limit * f_n_max
-
- # find the indexes of these contacts
- I = numpy.nonzero(data[:,6] >= f_n_lim)
-
- #color = matplotlib.cm.spectral(data[:,6]/f_n_max)
- for i in I[0]:
-
- x1 = data[i,0]
- #y1 = data[i,1]
- z1 = data[i,2]
- x2 = data[i,3]
- #y2 = data[i,4]
- z2 = data[i,5]
- f_n = data[i,6]
-
- lw_max = 1.0
- if f_n >= f_n_max:
- lw = lw_max
else:
- lw = (f_n - f_n_lim)/(f_n_max - f_n_lim)*lw_max
-
- #print lw
- axfc2.plot([x1,x2], [z1,z2], '-k', linewidth=lw)
- #axfc2.plot([x1,x2], [z1,z2], '-', linewidth=lw, color=forcemagnitudes[sc])
- #axfc2.plot([x1,x2], [z1,z2], '-', linewidth=lw, color=color)
-
- axfc2.spines['right'].set_visible(False)
- axfc2.spines['left'].set_visible(False)
- # Only show ticks on the left and bottom spines
- axfc2.xaxis.set_ticks_position('none')
- axfc2.yaxis.set_ticks_position('none')
- axfc2.set_xticklabels([])
- axfc2.set_yticklabels([])
- axfc2.set_xlim([numpy.min(data[I[0],0]), numpy.max(data[I[0],0])])
- axfc2.set_ylim([numpy.min(data[I[0],2]), numpy.max(data[I[0],2])])
- axfc2.set_aspect('equal')
-
-
-
-
- # Scatter plot 3
- sc=2
- #Lx=.50; Ly=.40;
- Lx=.65; Ly=.40;
- #xytext=(Lx+.5*dx, Ly+.5*dy)
- #xytext=(Lx+1.75*dx, Ly+0.01)
- xytext=(Lx+1.15*dx, 0.2)
- xy=(ts[sc]*scalingfactor, Ns[sc]/1000.)
- #print xytext
- #print xy
- ax1.annotate('',
- xytext=xytext, textcoords='axes fraction',
- xy=xy, xycoords='data',
- arrowprops=dict(arrowstyle='->', connectionstyle='arc3'))
-
- axsc2 = fig.add_axes([Lx, Ly, dx, dy], polar=True)
-
- cs = axsc2.scatter(strikelists[sc], 90. - diplists[sc], marker='o',
- c=forcemagnitudes[sc],
- s=forcemagnitudes[sc]/f_n_max*5.,
- alpha=alphas[sc],
- edgecolors='none',
- vmin=f_n_maxs[sc]*lower_limit,
- vmax=f_n_maxs[sc]*upper_limit,
- cmap=matplotlib.cm.get_cmap('afmhot_r'))#,
- #norm=matplotlib.colors.LogNorm())
- # tick locations
- #thetaticks = numpy.arange(0,360,90)
- # set ticklabels location at 1.3 times the axes' radius
- #ax.set_thetagrids(thetaticks, frac=1.3)
- axsc2.set_xticklabels([])
- axsc2.set_yticklabels([])
-
- axsc2.set_title('\\textbf{3. Slip}', fontsize=7)
- if upper_limit < 1.0:
- cbar = plt.colorbar(cs, extend='max', fraction=0.035, pad=0.04)
- else:
- cbar = plt.colorbar(cs, fraction=0.045, pad=0.04)
- cbar.set_label('$\\boldsymbol{f}_\\text{n}$ [N]')
- #cbar.set_label('Contact force [N]')
- cbar.locator = matplotlib.ticker.MaxNLocator(nbins=4)
- cbar.update_ticks()
-
- # plot defined max compressive stress from tau/N ratio
- axsc2.scatter(0., numpy.degrees(numpy.arctan(taus[sc]/Ns[sc])),
- marker='+', c='none', edgecolor='red', s=100)
- axsc2.scatter(0., numpy.degrees(numpy.arctan(taus[sc]/Ns[sc])),
- marker='o', c='none', edgecolor='red', s=100)
- '''
- ax.sc2scatter(0., # actual stress
- numpy.degrees(numpy.arctan(
- self.shearStress('effective')/
- self.currentNormalStress('effective'))),
- marker='+', color='blue', s=300)
- '''
-
- axsc2.set_rmax(90)
- axsc2.set_rticks([])
- #axsc2.grid(False)
-
- # force chain plot
-
- axfc2 = fig.add_axes([Lx-0.007, Ly-0.7*dy, dx, dy*0.7])
- #axfc2 = fig.add_axes([Lx+dx+0.05, Ly+0.03, dx, dy*0.7])
-
- data = datalists[sc]
-
- # find the max. value of the normal force
- #f_n_max = numpy.amax(data[:,6])
-
- # specify the lower limit of force chains to do statistics on
- f_n_lim = lower_limit * f_n_max
-
- # find the indexes of these contacts
- I = numpy.nonzero(data[:,6] >= f_n_lim)
-
- #color = matplotlib.cm.spectral(data[:,6]/f_n_max)
- for i in I[0]:
-
- x1 = data[i,0]
- #y1 = data[i,1]
- z1 = data[i,2]
- x2 = data[i,3]
- #y2 = data[i,4]
- z2 = data[i,5]
- f_n = data[i,6]
-
- lw_max = 1.0
- if f_n >= f_n_max:
- lw = lw_max
- else:
- lw = (f_n - f_n_lim)/(f_n_max - f_n_lim)*lw_max
-
- #print lw
- axfc2.plot([x1,x2], [z1,z2], '-k', linewidth=lw)
- #axfc2.plot([x1,x2], [z1,z2], '-', linewidth=lw, color=forcemagnitudes[sc])
- #axfc2.plot([x1,x2], [z1,z2], '-', linewidth=lw, color=color)
-
- axfc2.spines['right'].set_visible(False)
- axfc2.spines['left'].set_visible(False)
- # Only show ticks on the left and bottom spines
- axfc2.xaxis.set_ticks_position('none')
- axfc2.yaxis.set_ticks_position('none')
- axfc2.set_xticklabels([])
- axfc2.set_yticklabels([])
- axfc2.set_xlim([numpy.min(data[:,0]), numpy.max(data[:,0])])
- axfc2.set_ylim([numpy.min(data[:,2]), numpy.max(data[:,2])])
- axfc2.set_aspect('equal')
-
-
-
-
+ axhistload1 = fig.add_axes([Lx[sc], Ly[sc], dx, dy])
+ axhistload1.hist(datalists[sc][:,6], alpha=0.75, facecolor='gray', log=True)
+ #plt.yscale('log', nonposy='clip')
+ axhistload1.set_xlim([0,50.])
+ axhistload1.set_xlabel('Contact load [N]')
+ axhistload1.set_ylabel('Count')
+
+ # force chain plot
+
+ axfc1 = fig.add_axes([fc_x[sc], fc_y[sc], fc_dx[sc], fc_dy[sc]])
+
+ data = datalists[sc]
+
+ # find the max. value of the normal force
+ #f_n_max = numpy.amax(data[:,6])
+
+ # specify the lower limit of force chains to do statistics on
+ f_n_lim = lower_limit * f_n_max
+
+ # find the indexes of these contacts
+ I = numpy.nonzero(data[:,6] >= f_n_lim)
+
+ #color = matplotlib.cm.spectral(data[:,6]/f_n_max)
+ for i in I[0]:
+
+ x1 = data[i,0]
+ #y1 = data[i,1]
+ z1 = data[i,2]
+ x2 = data[i,3]
+ #y2 = data[i,4]
+ z2 = data[i,5]
+ f_n = data[i,6]
+
+ lw_max = 1.0
+ if f_n >= f_n_max:
+ lw = lw_max
+ else:
+ lw = (f_n - f_n_lim)/(f_n_max - f_n_lim)*lw_max
+
+ #print lw
+ axfc1.plot([x1,x2], [z1,z2], '-k', linewidth=lw)
+ #axfc1.plot([x1,x2], [z1,z2], '-', linewidth=lw, color=forcemagnitudes[sc])
+ #axfc1.plot([x1,x2], [z1,z2], '-', linewidth=lw, color=color)
+
+ axfc1.spines['right'].set_visible(False)
+ axfc1.spines['left'].set_visible(False)
+ # Only show ticks on the left and bottom spines
+ axfc1.xaxis.set_ticks_position('none')
+ axfc1.yaxis.set_ticks_position('none')
+ axfc1.set_xticklabels([])
+ axfc1.set_yticklabels([])
+ axfc1.set_xlim([numpy.min(data[I[0],0]), numpy.max(data[I[0],0])])
+ axfc1.set_ylim([numpy.min(data[I[0],2]), numpy.max(data[I[0],2])])
+ axfc1.set_aspect('equal')
(DIR) diff --git a/python/sphere.py b/python/sphere.py
t@@ -5372,6 +5372,7 @@ class sim:
#hist, bins = numpy.histogram(datadata[:,6], bins=10)
n, bins, patches = plt.hist(data[:,6], alpha=0.75, facecolor='gray')
#plt.xlabel('$\\boldsymbol{f}_\text{n}$ [N]')
+ plt.yscale('log', nonposy='clip')
plt.xlabel('Contact load [N]')
plt.ylabel('Count $N$')
plt.grid(True)