tfinish figure 3 - 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 e5de47d2e2884179b7956ee59e5b26d4ccf9ebd3
(DIR) parent f724591d19e45dc7456e70c22ba7158d3534533e
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Fri, 29 May 2015 15:40:46 +0200
finish figure 3
Diffstat:
M python/halfshear-darcy-creep-dynam… | 494 +++++++++++++++++--------------
M python/sphere.py | 72 ++++++++++++++++----------------
2 files changed, 313 insertions(+), 253 deletions(-)
---
(DIR) diff --git a/python/halfshear-darcy-creep-dynamics.py b/python/halfshear-darcy-creep-dynamics.py
t@@ -36,18 +36,25 @@ plotForceChains = True
plotHistograms = False
plotCombinedHistogram = True
-sids = ['halfshear-darcy-sigma0=80000.0' +
- '-k_c=3.5e-13-mu=1.04e-07-ss=10000.0-A=70000.0-f=0.2']
-timescalings = [5.787e-5]
+sid = 'halfshear-darcy-sigma0=80000.0' +\
+ '-k_c=3.5e-13-mu=1.04e-07-ss=10000.0-A=70000.0-f=0.2'
+timescaling = 5.787e-5
n = 6 # wave number
period = 500 # 1 period = 500 steps
steps = numpy.arange(n*period + 0.25*period, n*period + 0.75*period,
dtype=numpy.int)
+''' # 6.5, 6.6, 6.7 d
plotsteps = numpy.array([
- n*period + 0.40*period,
n*period + 0.50*period,
- n*period + 0.60*period],
+ n*period + 0.60*period,
+ n*period + 0.65*period],
+ dtype=numpy.int) # slow creep, fast creep, slip
+ '''
+plotsteps = numpy.array([
+ n*period + 0.50*period,
+ n*period + 0.60*period,
+ n*period + 0.65*period],
dtype=numpy.int) # slow creep, fast creep, slip
# steps[-1] - 10], dtype=numpy.int) # slow creep, fast creep, slip
# plotsteps = numpy.array([1670])
t@@ -63,6 +70,7 @@ f_n_maxs = [[], [], []]
taus = [[], [], []]
ts = [[], [], []]
Ns = [[], [], []]
+vs = [[], [], []]
# f_min = 1.0
# f_max = 1.0e16
t@@ -74,6 +82,10 @@ f_n_max = 500 # for force chain plots
N = numpy.zeros_like(steps, dtype=numpy.float64)
t = numpy.zeros_like(steps, dtype=numpy.float64)
+v = numpy.zeros_like(steps, dtype=numpy.float64)
+
+vel_avg = numpy.zeros_like(N)
+angvel_avg = numpy.zeros_like(N)
# insert plot positions
Lx = [.17, .37, .65]
t@@ -89,225 +101,273 @@ fc_dx = [dx, dx, dx]
fc_dy = [dy*.7, dy*.7, dy*.7]
-s = 0
-for sid in sids:
-
- sim = sphere.sim(sid, fluid=True)
- # print(sim.status())
- t_DEM_to_t_real = timescalings[s]
+sim = sphere.sim(sid, fluid=True)
+# print(sim.status())
+t_DEM_to_t_real = timescaling
- i = 0
- i_scatter = 0
- for step in steps:
+i = 0
+i_scatter = 0
+for step in steps:
- sim.readstep(step, verbose=False)
- if i == 0:
- L = sim.L
+ sim.readstep(step, verbose=False)
+ if i == 0:
+ L = sim.L
- N[i] = sim.currentNormalStress('defined')
- t[i] = sim.currentTime()
+ N[i] = sim.currentNormalStress('defined')
+ t[i] = sim.currentTime()
+ vel_avg[i] = numpy.average(sim.vel[:, 0])
+ angvel_avg[i] = numpy.average(sim.angvel[:, 1])
+ v[i] = sim.shearVelocity()
- if plotContacts:
- outfolder = '../img_out/' + sim.id() + '/'
- if not os.path.exists(outfolder):
- os.makedirs(outfolder)
+ 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():
- datalists[i_scatter], strikelists[i_scatter], diplists[i_scatter],\
- forcemagnitudes[i_scatter], alphas[i_scatter], \
- f_n_maxs[i_scatter] = sim.plotContacts(return_data=True,
- alpha=0.9,
- lower_limit=lower_limit,
- upper_limit=upper_limit)
- # contactfigs.append(
- # sim.plotContacts(return_fig=True,
- # f_min=f_min,
- # f_max=f_max))
- # datalists.append(data)
- # strikelists.append(strikelist)
- # diplists.append(diplist)
- # forcemagnitudes.append(forcemagnitude)
- # alphas.append(alpha)
- # f_n_maxs.append(f_n_max)
- ts[i_scatter] = t[i]
- Ns[i_scatter] = N[i]
- # taus.append(sim.shearStress('defined'))
- taus[i_scatter] = sim.shearStress('defined')
-
- # contactidx.append(step)
- i_scatter += 1
-
- i += 1
- s += 1
-
- # PLOTTING ######################################################
-
- # Time in days
- scalingfactor = 1./t_DEM_to_t_real / (24.*3600.)
- t_scaled = t*scalingfactor
-
- # Normal stress plot
- fig = plt.figure(figsize=[3.5, 3.5])
-
- ax1 = plt.subplot(1, 1, 1)
-
- ax1.plot(t_scaled, N/1000., '-k', label='$N$', clip_on=False)
- ax1.plot([t_scaled[int(len(t_scaled)*0.60)], 10],
- [taus[0]/.3/1000., taus[0]/.3/1000.], '--', color='gray')
- ax1.set_xlabel('Time [d]')
- ax1.set_ylabel('Effective normal stress $N$ [kPa]')
-
- ax1.spines['right'].set_visible(False)
- ax1.spines['top'].set_visible(False)
- # Only show ticks on the left and bottom spines
- ax1.yaxis.set_ticks_position('left')
- ax1.xaxis.set_ticks_position('bottom')
-
- ax1.set_xlim([numpy.min(t_scaled), numpy.max(t_scaled)])
- # ax1.locator_params(axis='x', nbins=5)
- ax1.locator_params(axis='y', nbins=4)
-
- # Contact scatter plots
-
- # Scatter plot 1
- 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)
+ if (step == plotsteps).any():
+ datalists[i_scatter], strikelists[i_scatter], diplists[i_scatter],\
+ forcemagnitudes[i_scatter], alphas[i_scatter], \
+ f_n_maxs[i_scatter] = sim.plotContacts(return_data=True,
+ alpha=0.9,
+ lower_limit=lower_limit,
+ upper_limit=upper_limit)
+ # contactfigs.append(
+ # sim.plotContacts(return_fig=True,
+ # f_min=f_min,
+ # f_max=f_max))
+ # datalists.append(data)
+ # strikelists.append(strikelist)
+ # diplists.append(diplist)
+ # forcemagnitudes.append(forcemagnitude)
+ # alphas.append(alpha)
+ # f_n_maxs.append(f_n_max)
+ ts[i_scatter] = t[i]
+ Ns[i_scatter] = N[i]
+ vs[i_scatter] = sim.shearVelocity()
+ # taus.append(sim.shearStress('defined'))
+ taus[i_scatter] = sim.shearStress('defined')
+
+ # contactidx.append(step)
+ i_scatter += 1
+
+ i += 1
+
+# PLOTTING ######################################################
+
+# Time in days
+scalingfactor = 1./t_DEM_to_t_real / (24.*3600.)
+t_scaled = t*scalingfactor
+
+# Normal stress plot
+fig = plt.figure(figsize=[3.5, 3.5])
+
+# ax1 = plt.subplot(1, 1, 1)
+ax1 = plt.subplot2grid((2, 3), (0, 0), colspan=3)
+
+ax1.add_patch(matplotlib.patches.Rectangle((6.61, -100.), .2, 500., alpha=0.7,
+ facecolor='orange',
+ edgecolor='none'))
+# shade stick periods
+# collection = matplotlib.collections.BrokenBarHCollection.span_where(
+# t_scaled, ymin=-20.0, ymax=200.0,
+# where=numpy.isclose(v, 0.0),
+# facecolor='black', alpha=0.2,
+# linewidth=0)
+# ax1.add_collection(collection)
+
+lns0 = ax1.plot(t_scaled, N/1000., '-k', label='$N$', clip_on=False)
+# ax1.plot([t_scaled[int(len(t_scaled)*0.60)], 10],
+# [taus[0]/.3/1000., taus[0]/.3/1000.], '--', color='gray')
+ax1.set_xlabel('Time $t$ [d]')
+ax1.set_ylabel('Effective normal stress $N$ [kPa]')
+
+ax2 = ax1.twinx()
+lns1 = ax2.semilogy(t_scaled, numpy.abs(vel_avg)*timescaling, '-b',
+ label='$\\bar{\\boldsymbol{v}}_x$',
+ clip_on=False)
+lns2 = ax2.semilogy(t_scaled, numpy.abs(angvel_avg)*timescaling, '-r',
+ label='$\\bar{\\boldsymbol{\\omega}}_y$',
+ clip_on=False)
+ax2.set_ylabel('Linear $\\bar{\\boldsymbol{v}}_x$ [m/s] or '
+ + 'angular velocity $\\bar{\\boldsymbol{\\omega}}_y$ [rad/s]')
+
+ax1.set_ylim([-20, 150])
+ax1.text(6.42, -5., 'Creep', horizontalalignment='center')
+ax1.text(6.68, -5., 'Slip', horizontalalignment='center')
+
+ax1.spines['right'].set_visible(False)
+ax1.spines['top'].set_visible(False)
+# Only show ticks on the left and bottom spines
+ax1.yaxis.set_ticks_position('left')
+ax1.xaxis.set_ticks_position('bottom')
+
+lns = lns0+lns1+lns2
+labs = [l.get_label() for l in lns]
+ax2.legend(lns, labs, loc='upper center', ncol=3, bbox_to_anchor=(0.5, 1.12),
+ fancybox=True, framealpha=1.0)
+
+ax1.set_xlim([numpy.min(t_scaled), numpy.max(t_scaled)])
+# ax1.locator_params(axis='x', nbins=5)
+ax1.locator_params(axis='y', nbins=4)
+
+# Contact scatter plots
+
+# Scatter plot 1
+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)
-
- if plotHistograms:
- axhistload1 = fig.add_axes([Lx[sc], Ly[sc], dx, dy*.5])
- 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]')
-
- axdipload1 = fig.add_axes([Lx[sc], Ly[sc]-.8*dy, dx, dy*.5])
- axdipload1.hist(diplists[sc], alpha=0.75, facecolor='gray',
- log=False)
- axdipload1.set_xlabel('Contact angle')
- # axdipload1.set_ylabel('Count')
-
- # force chain plot
-
- if plotForceChains:
- 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 = 2.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=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')
-
- # fig.tight_layout()
- # plt.subplots_adjust(hspace=0.05)
- plt.subplots_adjust(right=0.82)
-
- filename = sid + '-creep-dynamics.' + outformat
- plt.savefig(filename)
- plt.close()
- shutil.copyfile(filename, '/home/adc/articles/own/3/graphics/' + filename)
- print(filename)
+ # 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)
+
+ if plotHistograms:
+ axhistload1 = fig.add_axes([Lx[sc], Ly[sc], dx, dy*.5])
+ 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]')
+
+ axdipload1 = fig.add_axes([Lx[sc], Ly[sc]-.8*dy, dx, dy*.5])
+ axdipload1.hist(diplists[sc], alpha=0.75, facecolor='gray',
+ log=False)
+ axdipload1.set_xlabel('Contact angle')
+ # axdipload1.set_ylabel('Count')
+
+ # force chain plot
+
+ if plotForceChains:
+ axfc1 = plt.subplot2grid((2, 3), (1, sc))
+ # 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.2
+ 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=color)
+
+ if sc == 0:
+ title = 'Slow creep'
+ elif sc == 1:
+ title = 'Fast creep'
+ elif sc == 2:
+ title = 'Slip'
+ else:
+ title = 'Unknown'
+ axfc1.set_title(title, fontsize=7)
+ if sc == 0:
+ axfc1.set_ylabel('Stress network', fontsize=7)
+ axfc1.set_xlabel('$t$ = {:.3} d\n'.format(ts[sc]*scalingfactor) +
+ '$v$ = {:.3} m/s'.format(vs[sc]*scalingfactor))
+
+ 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')
+
+fig.tight_layout()
+# plt.subplots_adjust(hspace=0.05)
+# plt.subplots_adjust(right=0.82)
+
+filename = sid + '-creep-dynamics.' + outformat
+plt.savefig(filename)
+plt.close()
+shutil.copyfile(filename, '/home/adc/articles/own/3/graphics/' + filename)
+print(filename)
if plotCombinedHistogram:
fig = plt.figure(figsize=[3.5, 3.5])
t@@ -319,8 +379,8 @@ if plotCombinedHistogram:
for sc in range(len(Lx)):
if f_max < numpy.max(datalists[sc][:, 6]):
f_max = numpy.max(datalists[sc][:, 6])
- f_max = numpy.ceil(f_max/10.)*10.
- print(f_max)
+ f_max = numpy.ceil(f_max/10.)*10.
+ print(f_max)
for sc in range(len(Lx)):
# axhistload1.hist(datalists[sc][:,6], alpha=0.75, facecolor='gray',
t@@ -336,8 +396,8 @@ if plotCombinedHistogram:
ax1.set_ylabel('Number of contacts')
ax1.legend(loc='upper right', fancybox=True, framealpha=1.0)
- filename = sid + '-creep-dynamics-hist.' + outformat
- plt.savefig(filename)
- plt.close()
- shutil.copyfile(filename, '/home/adc/articles/own/3/graphics/' + filename)
- print(filename)
+filename = sid + '-creep-dynamics-hist.' + outformat
+plt.savefig(filename)
+plt.close()
+shutil.copyfile(filename, '/home/adc/articles/own/3/graphics/' + filename)
+print(filename)
(DIR) diff --git a/python/sphere.py b/python/sphere.py
t@@ -1438,7 +1438,7 @@ class sim:
1.0, choose "scalar" as the "Scale Mode". Check the "Edit" checkbox, and
set the "Set Scale Factor" to 1.0. The field "Maximum Number of Points"
may be increased if the number of particles exceed the default value.
- Finally press "Apply", and the particles will appear in the main window.
+ Finally press "Apply", and the particles will appear in the main window.
The sphere resolution may be adjusted ("Theta resolution", "Phi
resolution") to increase the quality and the computational requirements
t@@ -1527,7 +1527,7 @@ class sim:
as the "Scale Mode". Check the "Edit" checkbox, and set the "Set Scale
Factor" to 1.0. The field "Maximum Number of Points" may be increased if
the number of particles exceed the default value. Finally press "Apply",
- and the particles will appear in the main window.
+ and the particles will appear in the main window.
The sphere resolution may be adjusted ("Theta resolution", "Phi
resolution") to increase the quality and the computational requirements
t@@ -1573,7 +1573,7 @@ class sim:
fh.write('\n')
fh.write(' </DataArray>\n')
fh.write(' </Points>\n')
-
+
### Data attributes
fh.write(' <PointData Scalars="Diameter [m]" Vectors="vector">\n')
t@@ -1610,7 +1610,7 @@ class sim:
if self.cfd_solver == 0: # Navier Stokes
# Fluid interaction force
- fh.write(' <DataArray type="Float32" '
+ fh.write(' <DataArray type="Float32" '
+ 'Name="Fluid force total [N]" '
+ 'NumberOfComponents="3" format="ascii">\n')
fh.write(' ')
t@@ -2244,7 +2244,7 @@ class sim:
radius_max = mean + variance
self.radius = numpy.random.uniform(radius_min, radius_max, self.np)
else:
- raise Exception('Particle size distribution type not understood ('
+ raise Exception('Particle size distribution type not understood ('
+ str(psd) + '). Valid values are \'uni\' or \'logn\'')
# Show radii as histogram
t@@ -2345,7 +2345,7 @@ class sim:
def wall0iz(self):
'''
Returns the cell index of wall 0 along z.
-
+
:returns: z cell index
:return type: int
'''
t@@ -2357,7 +2357,7 @@ class sim:
def normalBoundariesXY(self):
'''
Set the x and y boundary conditions to be static walls.
-
+
See also :func:`periodicBoundariesXY()` and
:func:`periodicBoundariesX()`
'''
t@@ -2366,7 +2366,7 @@ class sim:
def periodicBoundariesXY(self):
'''
Set the x and y boundary conditions to be periodic.
-
+
See also :func:`normalBoundariesXY()` and
:func:`periodicBoundariesX()`
'''
t@@ -2810,7 +2810,7 @@ class sim:
:param z_adjust: Increase the world and grid size by this amount to
allow for wall movement.
:type z_adjust: float
- :param
+ :param
'''
if idx == 0:
t@@ -2844,7 +2844,7 @@ class sim:
'''
Setup consolidation experiment. Specify the upper wall normal stress in
Pascal, default value is 10 kPa.
-
+
:param normal_stress: The normal stress to apply from the upper wall
:type normal_stress: float
'''
t@@ -3217,7 +3217,7 @@ class sim:
Initialize the fluid arrays and the fluid viscosity. The default value
of ``mu`` equals the dynamic viscosity of water at 25 degrees Celcius.
The value for water at 0 degrees Celcius is 17.87e-4 kg/(m*s).
-
+
:param mu: The fluid dynamic viscosity [kg/(m*s)]
:type mu: float
:param rho: The fluid density [kg/(m^3)]
t@@ -3487,7 +3487,7 @@ class sim:
:type gamma_n: float
:param gamma_t: Particle-particle contact tangential viscosity [Ns/m]
:type gamma_t: float
- :param gamma_r: Particle-particle contact rolling viscosity *Parameter
+ :param gamma_r: Particle-particle contact rolling viscosity *Parameter
not used*
:type gamma_r: float
:param gamma_wn: Wall-particle contact normal viscosity [Ns/m]
t@@ -3622,7 +3622,7 @@ class sim:
print('Warning: The system is critically dampened (ratio = '
+ str(damping_ratio) + ') in the normal component. '
+ '\nCritical damping = ' + str(critical_gamma) + '.')
-
+
def setDampingTangential(self, gamma, over_damping=False):
'''
Set the dampening coefficient (gamma) in the tangential direction of the
t@@ -3802,7 +3802,7 @@ class sim:
Calculates the current magnitude of the defined or effective top wall
normal stress.
- :param type: Find the 'defined' (default) or 'effective' normal stress
+ :param type: Find the 'defined' (default) or 'effective' normal stress
:type type: str
:returns: The current top wall normal stress in Pascal
t@@ -3850,7 +3850,7 @@ class sim:
:return type: float
'''
return self.rho[0]*self.volume(idx)
-
+
def smallestMass(self):
'''
Returns the mass of the leightest particle.
t@@ -3861,7 +3861,7 @@ class sim:
:return type: float
'''
return V_sphere(numpy.min(self.radius))
-
+
def largestMass(self):
'''
Returns the mass of the heaviest particle.
t@@ -3872,7 +3872,7 @@ class sim:
:return type: float
'''
return V_sphere(numpy.max(self.radius))
-
+
def momentOfInertia(self, idx):
'''
Returns the moment of inertia of a particle.
t@@ -3940,7 +3940,7 @@ class sim:
:return type: float
'''
return self.ev[idx]
-
+
def totalViscousEnergy(self):
'''
Returns the total viscous dissipated energy for all particles.
t@@ -3963,7 +3963,7 @@ class sim:
:return type: float
'''
return self.es[idx]
-
+
def totalFrictionalEnergy(self):
'''
Returns the total frictional dissipated energy for all particles.
t@@ -4050,7 +4050,7 @@ class sim:
def voidRatio(self):
'''
Calculates the current void ratio
-
+
:returns: The void ratio, in [0:1]
:return type: float
'''
t@@ -4253,7 +4253,7 @@ class sim:
:param workdir: The working directory during the calculations, if
`use_workdir=True`
:type workdir: str
-
+
'''
filename = self.sid + ".sh"
t@@ -4367,7 +4367,7 @@ class sim:
'''
Using the built-in ray tracer, render all output files that belong to
the simulation, determined by the simulation id (``sid``).
-
+
:param method: The color visualization method to use for the particles.
Possible values are: 'normal': color all particles with the same
color, 'pres': color by pressure, 'vel': color by translational
t@@ -4420,7 +4420,7 @@ class sim:
'''
Uses ffmpeg to combine images to animation. All images should be
rendered beforehand using :func:`render()`.
-
+
:param out_folder: The output folder for the video file
:type out_folder: str
:param video_format: The format of the output video
t@@ -4687,7 +4687,7 @@ class sim:
'''
Visualizes the force chains in the system from the magnitude of the
normal contact forces, and produces an image of them. Warning: Will
- segfault if no contacts are found.
+ segfault if no contacts are found.
:param lc: Lower cutoff of contact forces. Contacts below are not
visualized
t@@ -4848,7 +4848,7 @@ class sim:
'''
Returns the current simulation status by using the simulation id
(``sid``) as an identifier.
-
+
:returns: The number of the last output file written
:return type: int
'''
t@@ -4882,7 +4882,7 @@ class sim:
Plot the particle x-axis displacement against the original vertical
particle position. The plot is saved in the current directory with the
file name '<simulation id>-sheardisp.<graphics_format>'.
-
+
:param graphics_format: Save the plot in this format
:type graphics_format: str
'''
t@@ -5446,7 +5446,7 @@ class sim:
graphics_format,\
transparent=False)
-
+
plt.close()
t@@ -5739,7 +5739,7 @@ class sim:
self.H100 = H[-1]
self.H50 = (self.H0 + self.H100)/2.0
T50 = 0.197 # case I
-
+
# find the time where 50% of the consolidation (H50) has happened by
# linear interpolation. The values in H are expected to be
# monotonically decreasing. See Numerical Recipies p. 115
t@@ -5775,7 +5775,7 @@ class sim:
plt.savefig(self.sid + '-loadcurve.' + graphics_format)
plt.clf()
plt.close(fig)
-
+
def convergence(self):
'''
Read the convergence evolution in the CFD solver. The values are stored
t@@ -5789,7 +5789,7 @@ class sim:
def plotConvergence(self, graphics_format='png'):
'''
- Plot the convergence evolution in the CFD solver. The plot is saved
+ Plot the convergence evolution in the CFD solver. The plot is saved
in the output folder with the file name
'<simulation id>-conv.<graphics_format>'.
t@@ -6077,7 +6077,7 @@ class sim:
Calculates the sum of shear stress values measured on any moving
particles with a finite and fixed velocity.
- :param type: Find the 'defined' or 'effective' (default) shear stress
+ :param type: Find the 'defined' or 'effective' (default) shear stress
:type type: str
:returns: The shear stress in Pa
t@@ -6653,7 +6653,7 @@ class sim:
ax5 = plt.subplot(3, 1, 2, sharex=ax1)
ax5.semilogy(time[1:], self.v[1:], label='Shear velocity')
ax5.set_ylabel('Shear velocity [ms$^{-1}$]')
-
+
# shade stick periods
collection = \
matplotlib.collections.BrokenBarHCollection.span_where(
t@@ -7044,7 +7044,7 @@ class sim:
outfolder = '../img_out/')
# render images to movie
- subprocess.call('cd ../img_out/ && ' +
+ subprocess.call('cd ../img_out/ && ' +
'ffmpeg -sameq -i {}.%05d-contacts.png '.format(self.sid) +
'{}-contacts.mp4'.format(self.sid),
#'convert -quality 100 {}.*.png {}-contacts.avi'.format(
t@@ -7105,7 +7105,7 @@ def render(binary,
verbose=True):
'''
Render target binary using the ``sphere`` raytracer.
-
+
:param method: The color visualization method to use for the particles.
Possible values are: 'normal': color all particles with the same
color, 'pres': color by pressure, 'vel': color by translational
t@@ -7153,7 +7153,7 @@ def video(project,
'''
Uses ffmpeg to combine images to animation. All images should be
rendered beforehand using :func:`render()`.
-
+
:param project: The simulation id of the project to render
:type project: str
:param out_folder: The output folder for the video file
t@@ -7350,7 +7350,7 @@ def status(project):
'''
Check the status.dat file for the target project, and return the last output
file number.
-
+
:param project: The simulation id of the target project
:type project: str