tshear-results-strain.py - 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
---
tshear-results-strain.py (4510B)
---
1 #!/usr/bin/env python
2 import matplotlib
3 matplotlib.use('Agg')
4 matplotlib.rcParams.update({'font.size': 18, 'font.family': 'serif'})
5 matplotlib.rc('text', usetex=True)
6 matplotlib.rcParams['text.latex.preamble']=[r"\usepackage{amsmath}"]
7 import shutil
8
9 import os
10 import numpy
11 import sphere
12 from permeabilitycalculator import *
13 import matplotlib.pyplot as plt
14 from matplotlib.ticker import MaxNLocator
15
16 sigma0 = 20000.0
17 #cvals = ['dry', 1.0, 0.1, 0.01]
18 cvals = ['dry', 1.0, 0.1]
19 #cvals = ['dry', 1.0]
20 #step = 1999
21
22 sim = sphere.sim('halfshear-sigma0=' + str(sigma0) + '-shear')
23 sim.readfirst(verbose=False)
24
25
26 # particle z positions
27 zpos_p = numpy.zeros((len(cvals), sim.np))
28
29 # cell midpoint cell positions
30 zpos_c = numpy.zeros((len(cvals), sim.num[2]))
31 dz = sim.L[2]/sim.num[2]
32 for i in numpy.arange(sim.num[2]):
33 zpos_c[:,i] = i*dz + 0.5*dz
34
35 # particle x displacements
36 xdisp = numpy.zeros((len(cvals), sim.np))
37
38 xdisp_mean = numpy.zeros((len(cvals), sim.num[2]))
39
40 s = 0
41 for c in cvals:
42
43 if c == 'dry':
44 fluid = False
45 sid = 'halfshear-sigma0=' + str(sigma0) + '-shear'
46 else:
47 fluid = True
48 sid = 'halfshear-sigma0=' + str(sigma0) + '-c=' + str(c) + '-shear'
49
50 sim = sphere.sim(sid, fluid=fluid)
51
52 if os.path.isfile('../output/' + sid + '.status.dat'):
53
54 sim.readlast(verbose=False)
55
56 zpos_p[s,:] = sim.x[:,2]
57
58 xdisp[s,:] = sim.xyzsum[:,0]
59
60 #shear_strain[s] += sim.shearStrain()/nsteps_avg
61
62 # calculate mean values of xdisp and f_pf
63 for iz in numpy.arange(sim.num[2]):
64 z_bot = iz*dz
65 z_top = (iz+1)*dz
66 I = numpy.nonzero((zpos_p[s,:] >= z_bot) & (zpos_p[s,:] < z_top))
67 if len(I) > 0:
68 xdisp_mean[s,iz] = numpy.mean(xdisp[s,I])
69
70 # normalize distance
71 max_dist = numpy.nanmax(xdisp_mean[s])
72 xdisp_mean[s] /= max_dist
73
74 else:
75 print(sid + ' not found')
76 s += 1
77
78
79 #fig = plt.figure(figsize=(8,4*(len(steps))+1))
80 #fig = plt.figure(figsize=(8,5*(len(steps))+1))
81 fig = plt.figure(figsize=(8,6))
82
83 ax = []
84 #linetype = ['-', '--', '-.']
85 linetype = ['-', '-', '-', '-']
86 color = ['b','g','c','y']
87 for s in numpy.arange(len(cvals)):
88
89 ax.append(plt.subplot(111))
90 #ax.append(plt.subplot(len(steps)*100 + 31 + s*3))
91 #ax.append(plt.subplot(len(steps)*100 + 32 + s*3, sharey=ax[s*4+0]))
92 #ax.append(plt.subplot(len(steps)*100 + 33 + s*3, sharey=ax[s*4+0]))
93 #ax.append(ax[s*4+2].twiny())
94
95 if cvals[s] == 'dry':
96 legend = 'dry'
97 else:
98 legend = 'wet, c = ' + str(cvals[s])
99
100 #ax[0].plot(xdisp[s], zpos_p[s], ',', color = '#888888')
101 #ax[0].plot(xdisp[s], zpos_p[s], ',', color=color[s], alpha=0.5)
102 ax[0].plot(xdisp_mean[s], zpos_c[s], linetype[s],
103 color=color[s], label=legend, linewidth=1)
104
105 ax[0].set_ylabel('Vertical position $z$ [m]')
106 #ax[0].set_xlabel('$\\boldsymbol{x}^x_\\text{p}$ [m]')
107 ax[0].set_xlabel('Normalized horizontal distance')
108
109 #ax[s*4+0].get_xaxis().set_major_locator(MaxNLocator(nbins=5))
110 #ax[s*4+1].get_xaxis().set_major_locator(MaxNLocator(nbins=5))
111 #ax[s*4+2].get_xaxis().set_major_locator(MaxNLocator(nbins=5))
112
113 #plt.setp(ax[s*4+0].xaxis.get_majorticklabels(), rotation=90)
114 #plt.setp(ax[s*4+1].xaxis.get_majorticklabels(), rotation=90)
115 #plt.setp(ax[s*4+2].xaxis.get_majorticklabels(), rotation=90)
116 #plt.setp(ax[s*4+3].xaxis.get_majorticklabels(), rotation=90)
117
118 #if s == 0:
119 #y = 0.95
120 #if s == 1:
121 #y = 0.55
122
123 #strain_str = 'Shear strain $\\gamma = %.3f$' % (shear_strain[s])
124 #fig.text(0.1, y, strain_str, horizontalalignment='left', fontsize=22)
125 #ax[s*4+0].annotate(strain_str, xytext=(0,1.1), textcoords='figure fraction',
126 #horizontalalignment='left', fontsize=22)
127 #plt.text(0.05, 1.06, strain_str, horizontalalignment='left', fontsize=22,
128 #transform=ax[s*4+0].transAxes)
129 #ax[s*4+0].set_title(strain_str)
130
131 #ax[s*4+0].grid()
132 #ax[s*4+1].grid()
133 #ax[s*4+2].grid()
134 #ax1.legend(loc='lower right', prop={'size':18})
135 #ax2.legend(loc='lower right', prop={'size':18})
136
137 legend_alpha=0.5
138 ax[0].legend(loc='lower right', prop={'size':18}, fancybox=True, framealpha=legend_alpha)
139 ax[0].grid()
140 plt.tight_layout()
141 plt.subplots_adjust(wspace = .05)
142 plt.MaxNLocator(nbins=4)
143
144 filename = 'shear-' + str(int(sigma0/1000.0)) + 'kPa-strain.pdf'
145 plt.savefig(filename)
146 shutil.copyfile(filename, '/home/adc/articles/own/2-org/' + filename)
147 print(filename)
148
149