tadd figures from simulations with N=80 kPa - 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 c9144d2a6bda44c762028cb0a92d972ebf774a2f
(DIR) parent 379813ad7fcb70b3647d4026d1764225873cc186
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Thu, 23 Apr 2015 22:30:28 +0200
add figures from simulations with N=80 kPa
Diffstat:
M python/halfshear-darcy-fluid-press… | 8 ++++++++
M python/halfshear-darcy-internals.py | 398 ++++++++++++++++---------------
M python/halfshear-darcy-strain.py | 334 ++++++++++++++++---------------
M python/halfshear-darcy-strength-di… | 652 ++++++++++++++++---------------
M python/halfshear-darcy-strength-di… | 657 ++++++++++++++++---------------
5 files changed, 1034 insertions(+), 1015 deletions(-)
---
(DIR) diff --git a/python/halfshear-darcy-fluid-pressures.py b/python/halfshear-darcy-fluid-pressures.py
t@@ -25,6 +25,14 @@ sids = [
'halfshear-darcy-sigma0=20000.0-k_c=3.5e-15-mu=1.797e-07-velfac=1.0-shear',
'halfshear-darcy-sigma0=20000.0-k_c=3.5e-15-mu=1.797e-08-velfac=1.0-shear'#,
+
+
+'halfshear-darcy-sigma0=80000.0-k_c=3.5e-15-mu=1.797e-06-velfac=1.0-shear',
+'halfshear-darcy-sigma0=80000.0-k_c=3.5e-14-mu=1.797e-06-velfac=1.0-shear',
+'halfshear-darcy-sigma0=80000.0-k_c=3.5e-13-mu=1.797e-06-velfac=1.0-shear',
+
+'halfshear-darcy-sigma0=80000.0-k_c=3.5e-15-mu=1.797e-07-velfac=1.0-shear',
+'halfshear-darcy-sigma0=80000.0-k_c=3.5e-15-mu=1.797e-08-velfac=1.0-shear'#,
]
for sid in sids:
(DIR) diff --git a/python/halfshear-darcy-internals.py b/python/halfshear-darcy-internals.py
t@@ -21,7 +21,7 @@ import seaborn as sns
sns.set(style='ticks', palette='Blues')
#sigma0 = float(sys.argv[1])
-sigma0 = 20000.0
+sigma0_list = [20000.0, 80000.0]
#k_c = 3.5e-13
#k_c = float(sys.argv[1])
k_c_list = [3.5e-13, 3.5e-14, 3.5e-15]
t@@ -37,228 +37,232 @@ nsteps_avg = 1 # no. of steps to average over
steps = [1, 51, 101, 152, 204]
-for k_c in k_c_list:
+for sigma0 in sigma0_list:
+ for k_c in k_c_list:
- sid = 'halfshear-darcy-sigma0=' + str(sigma0) + '-k_c=' + str(k_c) + \
- '-mu=1.797e-06-velfac=1.0-shear'
- sim = sphere.sim(sid, fluid=True)
- sim.readfirst(verbose=False)
+ sid = 'halfshear-darcy-sigma0=' + str(sigma0) + '-k_c=' + str(k_c) + \
+ '-mu=1.797e-06-velfac=1.0-shear'
+ sim = sphere.sim(sid, fluid=True)
+ sim.readfirst(verbose=False)
- # particle z positions
- zpos_p = numpy.zeros((len(steps), sim.np))
+ # 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
+ # 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 x displacements
+ xdisp = numpy.zeros((len(steps), sim.np))
- # particle z velocity
- v_z_p = numpy.zeros((len(steps), sim.np))
+ # particle z velocity
+ v_z_p = numpy.zeros((len(steps), sim.np))
- # fluid permeability
- k = numpy.zeros((len(steps), sim.num[0], sim.num[1], sim.num[2]))
- k_bar = numpy.zeros((len(steps), sim.num[2]))
+ # fluid permeability
+ k = numpy.zeros((len(steps), sim.num[0], sim.num[1], sim.num[2]))
+ k_bar = numpy.zeros((len(steps), sim.num[2]))
- # pressure
- p = numpy.zeros((len(steps), sim.num[2]))
+ # pressure
+ 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]))
+ # 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)
+ # particle-fluid force per particle
+ f_pf = numpy.zeros_like(xdisp)
- # pressure - hydrostatic pressure
- #dev_p = numpy.zeros((len(steps), sim.num[2]))
+ # 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
+ phi_bar = numpy.zeros((len(steps), sim.num[2]))
- # mean porosity change
- dphi_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]))
+ # 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_start = numpy.zeros(len(steps))
- shear_strain_end = numpy.zeros(len(steps))
+ shear_strain_start = numpy.zeros(len(steps))
+ shear_strain_end = numpy.zeros(len(steps))
- #fig = plt.figure(figsize=(8,4*(len(steps))+1))
- #fig = plt.figure(figsize=(8,4.5))
- fig = plt.figure(figsize=(3.74*2,3.00))
- ax = []
- n = 4
- ax.append(plt.subplot(1, n, 1)) # 0: xdisp
- ax.append(plt.subplot(1, n, 2, sharey=ax[0])) # 3: k
- ax.append(plt.subplot(1, n, 3, sharey=ax[0])) # 5: p_f
- ax.append(plt.subplot(1, n, 4, sharey=ax[0])) # 6: f_pf_z
+ #fig = plt.figure(figsize=(8,4*(len(steps))+1))
+ #fig = plt.figure(figsize=(8,4.5))
+ fig = plt.figure(figsize=(3.74*2,3.00))
+ ax = []
+ n = 4
+ ax.append(plt.subplot(1, n, 1)) # 0: xdisp
+ ax.append(plt.subplot(1, n, 2, sharey=ax[0])) # 3: k
+ ax.append(plt.subplot(1, n, 3, sharey=ax[0])) # 5: p_f
+ ax.append(plt.subplot(1, n, 4, sharey=ax[0])) # 6: f_pf_z
- s = 0
- for step_str in steps:
+ s = 0
+ for step_str in steps:
- step = int(step_str)
+ step = int(step_str)
- if os.path.isfile('../output/' + sid + '.status.dat'):
+ if os.path.isfile('../output/' + sid + '.status.dat'):
- for substep in numpy.arange(nsteps_avg):
+ for substep in numpy.arange(nsteps_avg):
- if step + substep > sim.status():
- raise Exception(
- 'Simulation step %d not available (sim.status = %d).'
- % (step + substep, sim.status()))
+ if step + substep > sim.status():
+ raise Exception(
+ 'Simulation step %d not available (sim.status = %d).'
+ % (step + substep, sim.status()))
- sim.readstep(step + substep, verbose=False)
+ sim.readstep(step + substep, verbose=False)
- zpos_p[s,:] += sim.x[:,2]/nsteps_avg
+ zpos_p[s,:] += sim.x[:,2]/nsteps_avg
- xdisp[s,:] += sim.xyzsum[:,0]/nsteps_avg
+ xdisp[s,:] += sim.xyzsum[:,0]/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_p[:,2]
+ '''
+ 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_p[:,2]
- dz = sim.L[2]/sim.num[2]
- wall0_iz = int(sim.w_x[0]/dz)
+ dz = sim.L[2]/sim.num[2]
+ wall0_iz = int(sim.w_x[0]/dz)
- p[s,:] += numpy.average(numpy.average(sim.p_f[:,:,:], axis=0),\
- axis=0)/nsteps_avg
+ p[s,:] += numpy.average(numpy.average(sim.p_f[:,:,:], axis=0),\
+ axis=0)/nsteps_avg
- sim.findPermeabilities()
- k[s,:] += sim.k[:,:,:]/nsteps_avg
+ sim.findPermeabilities()
+ k[s,:] += sim.k[:,:,:]/nsteps_avg
- k_bar[s,:] += \
- numpy.average(numpy.average(sim.k[:,:,:], axis=0), axis=0)\
- /nsteps_avg
-
- if substep == 0:
- shear_strain_start[s] = sim.shearStrain()
- else:
- shear_strain_end[s] = sim.shearStrain()
+ k_bar[s,:] += \
+ numpy.average(numpy.average(sim.k[:,:,:], axis=0), axis=0)\
+ /nsteps_avg
+
+ if substep == 0:
+ shear_strain_start[s] = sim.shearStrain()
+ else:
+ shear_strain_end[s] = sim.shearStrain()
- # 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])
- f_pf_mean[s,iz] = numpy.mean(f_pf[s,I])
-
- k_bar[s][0] = k_bar[s][1]
-
-
-
- #ax[0].plot(xdisp[s], zpos_p[s], ',', color = '#888888')
- ax[0].plot(xdisp_mean[s], zpos_c[s], label='$\gamma$ = %.3f' %
- (shear_strain_start[s]))
-
- ax[1].semilogx(k_bar[s], zpos_c[s], label='$\gamma$ = %.3f' %
- (shear_strain_start[s]))
-
- ax[2].plot(p[s]/1000.0, zpos_c[s], label='$\gamma$ = %.3f' %
- (shear_strain_start[s]))
-
- # 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]
- #f_pf_mean_nonzero = numpy.append(f_pf_mean_nonzero, 0.0)
- #zpos_c_nonzero = numpy.append(zpos_c_nonzero, zpos_c[s][zpos_c_nonzero.size])
-
- #ax[3].plot(f_pf_nonzero, zpos_p_nonzero, ',', alpha=0.5,
- #color='#888888')
- ax[3].plot(f_pf_mean_nonzero, zpos_c_nonzero, label='$\gamma$ = %.3f' %
- (shear_strain_start[s]))
-
- else:
- print(sid + ' not found')
- s += 1
-
-
- #max_z = numpy.max(zpos_p)
- max_z = 0.5
- #max_z = numpy.max(zpos_c[-2])
- ax[0].set_ylim([0, max_z])
- #ax[0].set_xlim([0, 0.5])
- ax[0].set_xlim([0, 0.05])
-
- if k_c == 3.5e-15:
- #ax[1].set_xlim([1e-14, 1e-12])
- ax[1].set_xlim([1e-16, 1e-14])
- elif k_c == 3.5e-14:
- ax[1].set_xlim([1e-15, 1e-13])
- elif k_c == 3.5e-13:
- #ax[1].set_xlim([1e-12, 1e-10])
- ax[1].set_xlim([1e-14, 1e-12])
-
- ax[0].set_ylabel('Vertical position $z$ [m]')
- ax[0].set_xlabel('$\\bar{\\boldsymbol{x}}^x_\\text{p}$ [m]')
- ax[1].set_xlabel('$\\bar{k}$ [m$^{2}$]')
- ax[2].set_xlabel('$\\bar{p_\\text{f}}$ [kPa]')
- ax[3].set_xlabel('$\\boldsymbol{f}^z_\\text{i}$ [N]')
-
- # align x labels
- #labely = -0.3
- #ax[0].xaxis.set_label_coords(0.5, labely)
- #ax[1].xaxis.set_label_coords(0.5, labely)
- #ax[2].xaxis.set_label_coords(0.5, labely)
- #ax[3].xaxis.set_label_coords(0.5, labely)
-
- plt.setp(ax[1].get_yticklabels(), visible=False)
- plt.setp(ax[2].get_yticklabels(), visible=False)
- plt.setp(ax[3].get_yticklabels(), visible=False)
-
- plt.setp(ax[0].xaxis.get_majorticklabels(), rotation=90)
- plt.setp(ax[1].xaxis.get_majorticklabels(), rotation=90)
- plt.setp(ax[2].xaxis.get_majorticklabels(), rotation=90)
- plt.setp(ax[3].xaxis.get_majorticklabels(), rotation=90)
-
- '''
- ax[0].grid()
- ax[1].grid()
- ax[2].grid()
- ax[3].grid()
- '''
-
- sns.despine() # remove chartjunk
-
- for i in range(4):
- # vertical grid lines
- ax[i].get_xaxis().grid(True, linestyle=':', linewidth=0.5)
- # horizontal grid lines
- ax[i].get_yaxis().grid(True, linestyle=':', linewidth=0.5)
-
- if i>0:
- ax[i].spines['left'].set_visible(False)
- ax[i].get_yaxis().set_ticks_position('none')
-
- legend_alpha=0.5
- #ax[0].legend(loc='lower center', prop={'size':12}, fancybox=True,
- #framealpha=legend_alpha)
- ax[0].legend(loc='lower right')
-
- #plt.subplots_adjust(wspace = .05) # doesn't work with tight_layout()
- #plt.MaxNLocator(nbins=1) # doesn't work?
- ax[0].locator_params(nbins=4)
- ax[2].locator_params(nbins=5)
- ax[3].locator_params(nbins=5)
-
- plt.tight_layout()
-
- filename = 'halfshear-darcy-internals-k_c=%.0e.pdf' % (k_c)
- plt.savefig(filename)
- #subprocess.call('pdfcrop ' + filename, shell=True)
- shutil.copyfile(filename, '/home/adc/articles/own/2/graphics/' + filename)
- print(filename)
+ # 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])
+ f_pf_mean[s,iz] = numpy.mean(f_pf[s,I])
+
+ k_bar[s][0] = k_bar[s][1]
+
+
+
+ #ax[0].plot(xdisp[s], zpos_p[s], ',', color = '#888888')
+ ax[0].plot(xdisp_mean[s], zpos_c[s], label='$\gamma$ = %.3f' %
+ (shear_strain_start[s]))
+
+ ax[1].semilogx(k_bar[s], zpos_c[s], label='$\gamma$ = %.3f' %
+ (shear_strain_start[s]))
+
+ ax[2].plot(p[s]/1000.0, zpos_c[s], label='$\gamma$ = %.3f' %
+ (shear_strain_start[s]))
+
+ # 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]
+ #f_pf_mean_nonzero = numpy.append(f_pf_mean_nonzero, 0.0)
+ #zpos_c_nonzero = numpy.append(zpos_c_nonzero, zpos_c[s][zpos_c_nonzero.size])
+
+ #ax[3].plot(f_pf_nonzero, zpos_p_nonzero, ',', alpha=0.5,
+ #color='#888888')
+ ax[3].plot(f_pf_mean_nonzero, zpos_c_nonzero, label='$\gamma$ = %.3f' %
+ (shear_strain_start[s]))
+
+ else:
+ print(sid + ' not found')
+ s += 1
+
+
+ #max_z = numpy.max(zpos_p)
+ max_z = 0.5
+ #max_z = numpy.max(zpos_c[-2])
+ ax[0].set_ylim([0, max_z])
+ #ax[0].set_xlim([0, 0.5])
+ ax[0].set_xlim([0, 0.05])
+
+ if k_c == 3.5e-15:
+ #ax[1].set_xlim([1e-14, 1e-12])
+ ax[1].set_xlim([1e-16, 1e-14])
+ elif k_c == 3.5e-14:
+ ax[1].set_xlim([1e-15, 1e-13])
+ elif k_c == 3.5e-13:
+ #ax[1].set_xlim([1e-12, 1e-10])
+ ax[1].set_xlim([1e-14, 1e-12])
+
+ ax[0].set_ylabel('Vertical position $z$ [m]')
+ ax[0].set_xlabel('$\\bar{\\boldsymbol{x}}^x_\\text{p}$ [m]')
+ ax[1].set_xlabel('$\\bar{k}$ [m$^{2}$]')
+ ax[2].set_xlabel('$\\bar{p_\\text{f}}$ [kPa]')
+ ax[3].set_xlabel('$\\boldsymbol{f}^z_\\text{i}$ [N]')
+
+ # align x labels
+ #labely = -0.3
+ #ax[0].xaxis.set_label_coords(0.5, labely)
+ #ax[1].xaxis.set_label_coords(0.5, labely)
+ #ax[2].xaxis.set_label_coords(0.5, labely)
+ #ax[3].xaxis.set_label_coords(0.5, labely)
+
+ plt.setp(ax[1].get_yticklabels(), visible=False)
+ plt.setp(ax[2].get_yticklabels(), visible=False)
+ plt.setp(ax[3].get_yticklabels(), visible=False)
+
+ plt.setp(ax[0].xaxis.get_majorticklabels(), rotation=90)
+ plt.setp(ax[1].xaxis.get_majorticklabels(), rotation=90)
+ plt.setp(ax[2].xaxis.get_majorticklabels(), rotation=90)
+ plt.setp(ax[3].xaxis.get_majorticklabels(), rotation=90)
+
+ '''
+ ax[0].grid()
+ ax[1].grid()
+ ax[2].grid()
+ ax[3].grid()
+ '''
+
+ sns.despine() # remove chartjunk
+
+ for i in range(4):
+ # vertical grid lines
+ ax[i].get_xaxis().grid(True, linestyle=':', linewidth=0.5)
+ # horizontal grid lines
+ ax[i].get_yaxis().grid(True, linestyle=':', linewidth=0.5)
+
+ if i>0:
+ ax[i].spines['left'].set_visible(False)
+ ax[i].get_yaxis().set_ticks_position('none')
+
+ legend_alpha=0.5
+ #ax[0].legend(loc='lower center', prop={'size':12}, fancybox=True,
+ #framealpha=legend_alpha)
+ ax[0].legend(loc='lower right')
+
+ #plt.subplots_adjust(wspace = .05) # doesn't work with tight_layout()
+ #plt.MaxNLocator(nbins=1) # doesn't work?
+ ax[0].locator_params(nbins=4)
+ ax[2].locator_params(nbins=5)
+ ax[3].locator_params(nbins=5)
+
+ plt.tight_layout()
+
+ filename = 'halfshear-darcy-internals-k_c=%.0e.pdf' % (k_c)
+ if sigma0 == 80000.0:
+ filename = 'halfshear-darcy-internals-k_c=%.0e-N80.pdf' % (k_c)
+
+ plt.savefig(filename)
+ #subprocess.call('pdfcrop ' + filename, shell=True)
+ shutil.copyfile(filename, '/home/adc/articles/own/2/graphics/' + filename)
+ print(filename)
(DIR) diff --git a/python/halfshear-darcy-strain.py b/python/halfshear-darcy-strain.py
t@@ -19,181 +19,183 @@ import seaborn as sns
sns.set(style='ticks', palette='Set2')
sns.despine() # remove chartjunk
-sigma0 = 20000.0
+sigma0_list = [20000.0, 80000.0]
#cvals = ['dry', 1.0, 0.1, 0.01]
#cvals = ['dry', 3.5e-13, 3.5e-15]
cvals = ['dry', 3.5e-13, 3.5e-14, 3.5e-15]
#cvals = ['dry', 1.0]
#step = 1999
-sim = sphere.sim('halfshear-sigma0=' + str(sigma0) + '-shear')
-sim.readfirst(verbose=False)
-
+for sigma0 in sigma0_list:
-# particle z positions
-zpos_p = [[], [], [], []]
-
-# cell midpoint cell positions
-zpos_c = [[], [], [], []]
+ sim = sphere.sim('halfshear-sigma0=' + str(sigma0) + '-shear')
+ sim.readfirst(verbose=False)
-# particle x displacements
-xdisp = [[], [], [], []]
-xdisp_mean = [[], [], [], []]
-s = 0
-for c in cvals:
+ # particle z positions
+ zpos_p = [[], [], [], []]
- if c == 'dry':
- fluid = False
- sid = 'halfshear-sigma0=' + str(sigma0) + '-shear'
- else:
- fluid = True
- sid = 'halfshear-darcy-sigma0=' + str(sigma0) + '-k_c=' + str(c) + \
- '-mu=1.797e-06-velfac=1.0-shear'
-
- sim = sphere.sim(sid, fluid=fluid)
-
- if os.path.isfile('../output/' + sid + '.status.dat'):
-
- sim.readlast(verbose=False)
-
- zpos_c[s] = numpy.zeros(sim.num[2]*2)
- dz = sim.L[2]/(sim.num[2]*2)
- for i in numpy.arange(sim.num[2]*2):
- zpos_c[s][i] = i*dz + 0.5*dz
-
- xdisp[s] = numpy.zeros(sim.np)
- xdisp_mean[s] = numpy.zeros(sim.num[2]*2)
-
-
- zpos_p[s][:] = sim.x[:,2]
-
- xdisp[s][:] = sim.xyzsum[:,0]
-
- #shear_strain[s] += sim.shearStrain()/nsteps_avg
-
- # calculate mean values of xdisp and f_pf
- for iz in numpy.arange(sim.num[2]*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])
-
- # normalize distance
- max_dist = numpy.nanmax(xdisp_mean[s])
- xdisp_mean[s] /= max_dist
-
- 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=(8/2,6/2))
-fig = plt.figure(figsize=(3.74,3.47)) # 3.14 inch = 80 mm, 3.74 = 95 mm
-#fig = plt.figure(figsize=(8,6))
-
-ax = []
-#linetype = ['-', '--', '-.']
-#linetype = ['-', '-', '-', '-']
-linetype = ['-', '--', '-.', ':']
-#color = ['b','g','c','y']
-#color = ['k','g','c','y']
-color = ['y','g','c','k']
-#color = ['c','m','y','k']
-for s in numpy.arange(len(cvals)):
-#for s in numpy.arange(len(cvals)-1, -1, -1):
-
- ax.append(plt.subplot(111))
- #ax.append(plt.subplot(len(steps)*100 + 31 + s*3))
- #ax.append(plt.subplot(len(steps)*100 + 32 + s*3, sharey=ax[s*4+0]))
- #ax.append(plt.subplot(len(steps)*100 + 33 + s*3, sharey=ax[s*4+0]))
- #ax.append(ax[s*4+2].twiny())
-
- if cvals[s] == 'dry':
- legend = 'dry'
- elif cvals[s] == 3.5e-13:
- legend = 'wet, high permeability'
- elif cvals[s] == 3.5e-14:
- legend = 'wet, interm. permeability'
- elif cvals[s] == 3.5e-15:
- legend = 'wet, low permeability'
- else:
- legend = 'wet, $k_c$ = ' + str(cvals[s]) + ' m$^2$'
-
- #ax[0].plot(xdisp[s], zpos_p[s], ',', color = '#888888')
- #ax[0].plot(xdisp[s], zpos_p[s], ',', color=color[s], alpha=0.5)
- ax[0].plot(xdisp_mean[s], zpos_c[s], linetype[s],
- label=legend)#,
- #color=color[s],
- #linewidth=2.0)
-
- ax[0].set_ylabel('Vertical position $z$ [m]')
- #ax[0].set_xlabel('$\\boldsymbol{x}^x_\\text{p}$ [m]')
- ax[0].set_xlabel('Normalized horizontal movement')
-
- #ax[s*4+0].get_xaxis().set_major_locator(MaxNLocator(nbins=5))
- #ax[s*4+1].get_xaxis().set_major_locator(MaxNLocator(nbins=5))
- #ax[s*4+2].get_xaxis().set_major_locator(MaxNLocator(nbins=5))
-
- #plt.setp(ax[s*4+0].xaxis.get_majorticklabels(), rotation=90)
- #plt.setp(ax[s*4+1].xaxis.get_majorticklabels(), rotation=90)
- #plt.setp(ax[s*4+2].xaxis.get_majorticklabels(), rotation=90)
- #plt.setp(ax[s*4+3].xaxis.get_majorticklabels(), rotation=90)
-
- #if s == 0:
- #y = 0.95
- #if s == 1:
- #y = 0.55
-
- #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 fraction',
- #horizontalalignment='left', fontsize=22)
- #plt.text(0.05, 1.06, strain_str, horizontalalignment='left', fontsize=22,
- #transform=ax[s*4+0].transAxes)
- #ax[s*4+0].set_title(strain_str)
-
- #ax[s*4+0].grid()
- #ax[s*4+1].grid()
- #ax[s*4+2].grid()
- #ax1.legend(loc='lower right', prop={'size':18})
- #ax2.legend(loc='lower right', prop={'size':18})
-
-# remove box at top and right
-ax[0].spines['top'].set_visible(False)
-ax[0].spines['right'].set_visible(False)
-# remove ticks at top and right
-ax[0].get_xaxis().tick_bottom()
-ax[0].get_yaxis().tick_left()
-ax[0].get_xaxis().grid(False) # horizontal grid lines
-#ax[0].get_yaxis().grid(True, linestyle='--', linewidth=0.5) # vertical grid lines
-ax[0].get_xaxis().grid(True, linestyle=':', linewidth=0.5) # vertical grid lines
-ax[0].get_yaxis().grid(True, linestyle=':', linewidth=0.5) # vertical grid lines
-
-# reverse legend order
-handles, labels = ax[0].get_legend_handles_labels()
-ax[0].legend(handles[::-1], labels[::-1], loc='best')
-
-#legend_alpha=0.5
-#ax[0].legend(loc='lower right', prop={'size':18}, fancybox=True, framealpha=legend_alpha)
-#ax[0].legend(loc='best', prop={'size':18}, fancybox=True, framealpha=legend_alpha)
-#ax[0].legend(loc='best')
-#ax[0].grid()
-#ax[0].set_xlim([-0.05, 1.01])
-ax[0].set_xlim([-0.05, 1.05])
-#ax[0].set_ylim([0.0, 0.47])
-ax[0].set_ylim([0.20, 0.47])
-plt.tight_layout()
-plt.subplots_adjust(wspace = .05)
-plt.MaxNLocator(nbins=4)
-
-filename = 'halfshear-darcy-strain.pdf'
-plt.savefig(filename)
-#shutil.copyfile(filename, '/Users/adc/articles/own/2/graphics/' + filename)
-shutil.copyfile(filename, '/home/adc/articles/own/2/graphics/' + filename)
-print(filename)
+ # cell midpoint cell positions
+ zpos_c = [[], [], [], []]
+ # particle x displacements
+ xdisp = [[], [], [], []]
+ xdisp_mean = [[], [], [], []]
+ s = 0
+ for c in cvals:
+
+ if c == 'dry':
+ fluid = False
+ sid = 'halfshear-sigma0=' + str(sigma0) + '-shear'
+ else:
+ fluid = True
+ sid = 'halfshear-darcy-sigma0=' + str(sigma0) + '-k_c=' + str(c) + \
+ '-mu=1.797e-06-velfac=1.0-shear'
+
+ sim = sphere.sim(sid, fluid=fluid)
+
+ if os.path.isfile('../output/' + sid + '.status.dat'):
+
+ sim.readlast(verbose=False)
+
+ zpos_c[s] = numpy.zeros(sim.num[2]*2)
+ dz = sim.L[2]/(sim.num[2]*2)
+ for i in numpy.arange(sim.num[2]*2):
+ zpos_c[s][i] = i*dz + 0.5*dz
+
+ xdisp[s] = numpy.zeros(sim.np)
+ xdisp_mean[s] = numpy.zeros(sim.num[2]*2)
+
+
+ zpos_p[s][:] = sim.x[:,2]
+
+ xdisp[s][:] = sim.xyzsum[:,0]
+
+ #shear_strain[s] += sim.shearStrain()/nsteps_avg
+
+ # calculate mean values of xdisp and f_pf
+ for iz in numpy.arange(sim.num[2]*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])
+
+ # normalize distance
+ max_dist = numpy.nanmax(xdisp_mean[s])
+ xdisp_mean[s] /= max_dist
+
+ 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=(8/2,6/2))
+ fig = plt.figure(figsize=(3.74,3.47)) # 3.14 inch = 80 mm, 3.74 = 95 mm
+ #fig = plt.figure(figsize=(8,6))
+
+ ax = []
+ #linetype = ['-', '--', '-.']
+ #linetype = ['-', '-', '-', '-']
+ linetype = ['-', '--', '-.', ':']
+ #color = ['b','g','c','y']
+ #color = ['k','g','c','y']
+ color = ['y','g','c','k']
+ #color = ['c','m','y','k']
+ for s in numpy.arange(len(cvals)):
+ #for s in numpy.arange(len(cvals)-1, -1, -1):
+
+ ax.append(plt.subplot(111))
+ #ax.append(plt.subplot(len(steps)*100 + 31 + s*3))
+ #ax.append(plt.subplot(len(steps)*100 + 32 + s*3, sharey=ax[s*4+0]))
+ #ax.append(plt.subplot(len(steps)*100 + 33 + s*3, sharey=ax[s*4+0]))
+ #ax.append(ax[s*4+2].twiny())
+
+ if cvals[s] == 'dry':
+ legend = 'dry'
+ elif cvals[s] == 3.5e-13:
+ legend = 'wet, high permeability'
+ elif cvals[s] == 3.5e-14:
+ legend = 'wet, interm. permeability'
+ elif cvals[s] == 3.5e-15:
+ legend = 'wet, low permeability'
+ else:
+ legend = 'wet, $k_c$ = ' + str(cvals[s]) + ' m$^2$'
+
+ #ax[0].plot(xdisp[s], zpos_p[s], ',', color = '#888888')
+ #ax[0].plot(xdisp[s], zpos_p[s], ',', color=color[s], alpha=0.5)
+ ax[0].plot(xdisp_mean[s], zpos_c[s], linetype[s],
+ label=legend)#,
+ #color=color[s],
+ #linewidth=2.0)
+
+ ax[0].set_ylabel('Vertical position $z$ [m]')
+ #ax[0].set_xlabel('$\\boldsymbol{x}^x_\\text{p}$ [m]')
+ ax[0].set_xlabel('Normalized horizontal movement')
+
+ #ax[s*4+0].get_xaxis().set_major_locator(MaxNLocator(nbins=5))
+ #ax[s*4+1].get_xaxis().set_major_locator(MaxNLocator(nbins=5))
+ #ax[s*4+2].get_xaxis().set_major_locator(MaxNLocator(nbins=5))
+
+ #plt.setp(ax[s*4+0].xaxis.get_majorticklabels(), rotation=90)
+ #plt.setp(ax[s*4+1].xaxis.get_majorticklabels(), rotation=90)
+ #plt.setp(ax[s*4+2].xaxis.get_majorticklabels(), rotation=90)
+ #plt.setp(ax[s*4+3].xaxis.get_majorticklabels(), rotation=90)
+
+ #if s == 0:
+ #y = 0.95
+ #if s == 1:
+ #y = 0.55
+
+ #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 fraction',
+ #horizontalalignment='left', fontsize=22)
+ #plt.text(0.05, 1.06, strain_str, horizontalalignment='left', fontsize=22,
+ #transform=ax[s*4+0].transAxes)
+ #ax[s*4+0].set_title(strain_str)
+
+ #ax[s*4+0].grid()
+ #ax[s*4+1].grid()
+ #ax[s*4+2].grid()
+ #ax1.legend(loc='lower right', prop={'size':18})
+ #ax2.legend(loc='lower right', prop={'size':18})
+
+ # remove box at top and right
+ ax[0].spines['top'].set_visible(False)
+ ax[0].spines['right'].set_visible(False)
+ # remove ticks at top and right
+ ax[0].get_xaxis().tick_bottom()
+ ax[0].get_yaxis().tick_left()
+ ax[0].get_xaxis().grid(False) # horizontal grid lines
+ #ax[0].get_yaxis().grid(True, linestyle='--', linewidth=0.5) # vertical grid lines
+ ax[0].get_xaxis().grid(True, linestyle=':', linewidth=0.5) # vertical grid lines
+ ax[0].get_yaxis().grid(True, linestyle=':', linewidth=0.5) # vertical grid lines
+
+ # reverse legend order
+ handles, labels = ax[0].get_legend_handles_labels()
+ ax[0].legend(handles[::-1], labels[::-1], loc='best')
+
+ #legend_alpha=0.5
+ #ax[0].legend(loc='lower right', prop={'size':18}, fancybox=True, framealpha=legend_alpha)
+ #ax[0].legend(loc='best', prop={'size':18}, fancybox=True, framealpha=legend_alpha)
+ #ax[0].legend(loc='best')
+ #ax[0].grid()
+ #ax[0].set_xlim([-0.05, 1.01])
+ ax[0].set_xlim([-0.05, 1.05])
+ #ax[0].set_ylim([0.0, 0.47])
+ ax[0].set_ylim([0.20, 0.47])
+ plt.tight_layout()
+ plt.subplots_adjust(wspace = .05)
+ plt.MaxNLocator(nbins=4)
+
+ filename = 'halfshear-darcy-strain.pdf'
+ if sigma0 == 80000.0:
+ filename = 'halfshear-darcy-strain-N80.pdf'
+ plt.savefig(filename)
+ #shutil.copyfile(filename, '/Users/adc/articles/own/2/graphics/' + filename)
+ shutil.copyfile(filename, '/home/adc/articles/own/2/graphics/' + filename)
+ print(filename)
(DIR) diff --git a/python/halfshear-darcy-strength-dilation-rate.py b/python/halfshear-darcy-strength-dilation-rate.py
t@@ -27,7 +27,7 @@ smooth_friction = True
smooth_window = 100
#sigma0_list = numpy.array([1.0e3, 2.0e3, 4.0e3, 10.0e3, 20.0e3, 40.0e3])
-sigma0 = 20000.0
+sigma0_list = [20000.0, 80000.0]
#k_c_vals = [3.5e-13, 3.5e-15]
k_c = 3.5e-15
#k_c = 3.5e-13
t@@ -96,354 +96,356 @@ def smooth(x, window_len=10, window='hanning'):
y = numpy.convolve(w/w.sum(), s, mode='same')
return y[window_len-1:-window_len+1]
+for sigma0 in sigma0_list:
+ shear_strain = [[], [], [], []]
+ friction = [[], [], [], []]
+ dilation = [[], [], [], []]
+ p_min = [[], [], [], []]
+ p_mean = [[], [], [], []]
+ p_max = [[], [], [], []]
+ f_n_mean = [[], [], [], []]
+ f_n_max = [[], [], [], []]
+ v_f_z_mean = [[], [], [], []]
+
+ fluid=True
+
+ # wet shear
+ for c in numpy.arange(0,len(mu_f_vals)):
+ #for c in numpy.arange(len(mu_f_vals)-1, -1, -1):
+ mu_f = mu_f_vals[c]
+
+ # halfshear-darcy-sigma0=20000.0-k_c=3.5e-13-mu=1.797e-06-velfac=1.0-shear
+ sid = 'halfshear-darcy-sigma0=' + str(sigma0) + '-k_c=' + str(k_c) + \
+ '-mu=' + str(mu_f) + '-velfac=' + str(velfac) + '-shear'
+ #sid = 'halfshear-sigma0=' + str(sigma0) + '-c_v=' + str(c_v) +\
+ #'-c_a=0.0-velfac=1.0-shear'
+ if os.path.isfile('../output/' + sid + '.status.dat'):
+
+ sim = sphere.sim(sid, fluid=fluid)
+ n = sim.status()
+ #n = 20
+ shear_strain[c] = numpy.zeros(n)
+ friction[c] = numpy.zeros_like(shear_strain[c])
+ dilation[c] = numpy.zeros_like(shear_strain[c])
+
+ '''
+ sim.readlast(verbose=False)
+ #sim.visualize('shear')
+ shear_strain[c] = sim.shear_strain
+ #shear_strain[c] = numpy.arange(sim.status()+1)
+ #friction[c] = sim.tau/1000.0#/sim.sigma_eff
+ friction[c] = sim.shearStress('effective')/sim.currentNormalStress('defined')
+ dilation[c] = sim.dilation
+ '''
+
+ # fluid pressures and particle forces
+ p_mean[c] = numpy.zeros_like(shear_strain[c])
+ p_min[c] = numpy.zeros_like(shear_strain[c])
+ p_max[c] = numpy.zeros_like(shear_strain[c])
+ f_n_mean[c] = numpy.zeros_like(shear_strain[c])
+ f_n_max[c] = numpy.zeros_like(shear_strain[c])
-shear_strain = [[], [], [], []]
-friction = [[], [], [], []]
-dilation = [[], [], [], []]
-p_min = [[], [], [], []]
-p_mean = [[], [], [], []]
-p_max = [[], [], [], []]
-f_n_mean = [[], [], [], []]
-f_n_max = [[], [], [], []]
-v_f_z_mean = [[], [], [], []]
-
-fluid=True
-
-# wet shear
-for c in numpy.arange(0,len(mu_f_vals)):
-#for c in numpy.arange(len(mu_f_vals)-1, -1, -1):
- mu_f = mu_f_vals[c]
-
- # halfshear-darcy-sigma0=20000.0-k_c=3.5e-13-mu=1.797e-06-velfac=1.0-shear
- sid = 'halfshear-darcy-sigma0=' + str(sigma0) + '-k_c=' + str(k_c) + \
- '-mu=' + str(mu_f) + '-velfac=' + str(velfac) + '-shear'
- #sid = 'halfshear-sigma0=' + str(sigma0) + '-c_v=' + str(c_v) +\
- #'-c_a=0.0-velfac=1.0-shear'
- if os.path.isfile('../output/' + sid + '.status.dat'):
-
- sim = sphere.sim(sid, fluid=fluid)
- n = sim.status()
- #n = 20
- shear_strain[c] = numpy.zeros(n)
- friction[c] = numpy.zeros_like(shear_strain[c])
- dilation[c] = numpy.zeros_like(shear_strain[c])
-
- '''
- sim.readlast(verbose=False)
- #sim.visualize('shear')
- shear_strain[c] = sim.shear_strain
- #shear_strain[c] = numpy.arange(sim.status()+1)
- #friction[c] = sim.tau/1000.0#/sim.sigma_eff
- friction[c] = sim.shearStress('effective')/sim.currentNormalStress('defined')
- dilation[c] = sim.dilation
- '''
-
- # fluid pressures and particle forces
- p_mean[c] = numpy.zeros_like(shear_strain[c])
- p_min[c] = numpy.zeros_like(shear_strain[c])
- p_max[c] = numpy.zeros_like(shear_strain[c])
- f_n_mean[c] = numpy.zeros_like(shear_strain[c])
- f_n_max[c] = numpy.zeros_like(shear_strain[c])
-
- for i in numpy.arange(n):
-
- sim.readstep(i, verbose=False)
-
- shear_strain[c][i] = sim.shearStrain()
- friction[c][i] = sim.shearStress('effective')/sim.currentNormalStress('defined')
- dilation[c][i] = sim.w_x[0]
-
- if pressures:
- iz_top = int(sim.w_x[0]/(sim.L[2]/sim.num[2]))-1
- p_mean[c][i] = numpy.mean(sim.p_f[:,:,0:iz_top])/1000
- p_min[c][i] = numpy.min(sim.p_f[:,:,0:iz_top])/1000
- p_max[c][i] = numpy.max(sim.p_f[:,:,0:iz_top])/1000
+ for i in numpy.arange(n):
- if contact_forces:
- sim.findNormalForces()
- f_n_mean[c][i] = numpy.mean(sim.f_n_magn)
- f_n_max[c][i] = numpy.max(sim.f_n_magn)
+ sim.readstep(i, verbose=False)
- if zflow:
- v_f_z_mean[c] = numpy.zeros_like(shear_strain[c])
- for i in numpy.arange(n):
- v_f_z_mean[c][i] = numpy.mean(sim.v_f[:,:,:,2])
+ shear_strain[c][i] = sim.shearStrain()
+ friction[c][i] = sim.shearStress('effective')/sim.currentNormalStress('defined')
+ dilation[c][i] = sim.w_x[0]
- dilation[c] =\
- (dilation[c] - dilation[c][0])/(numpy.mean(sim.radius)*2.0)
+ if pressures:
+ iz_top = int(sim.w_x[0]/(sim.L[2]/sim.num[2]))-1
+ p_mean[c][i] = numpy.mean(sim.p_f[:,:,0:iz_top])/1000
+ p_min[c][i] = numpy.min(sim.p_f[:,:,0:iz_top])/1000
+ p_max[c][i] = numpy.max(sim.p_f[:,:,0:iz_top])/1000
- else:
- print(sid + ' not found')
+ if contact_forces:
+ sim.findNormalForces()
+ f_n_mean[c][i] = numpy.mean(sim.f_n_magn)
+ f_n_max[c][i] = numpy.max(sim.f_n_magn)
- # produce VTK files
- #for sid in sids:
- #sim = sphere.sim(sid, fluid=True)
- #sim.writeVTKall()
+ if zflow:
+ v_f_z_mean[c] = numpy.zeros_like(shear_strain[c])
+ for i in numpy.arange(n):
+ v_f_z_mean[c][i] = numpy.mean(sim.v_f[:,:,:,2])
+ dilation[c] =\
+ (dilation[c] - dilation[c][0])/(numpy.mean(sim.radius)*2.0)
-if zflow or pressures:
- #fig = plt.figure(figsize=(8,10))
- #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))
+ else:
+ print(sid + ' not found')
-#plt.subplot(3,1,1)
-#plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+ # produce VTK files
+ #for sid in sids:
+ #sim = sphere.sim(sid, fluid=True)
+ #sim.writeVTKall()
-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)
+ #fig = plt.figure(figsize=(8,10))
+ #fig = plt.figure(figsize=(3.74, 2*3.74))
+ fig = plt.figure(figsize=(2*3.74, 2*3.74))
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):
- mu_f = mu_f_vals[c]
-
- if numpy.isclose(mu_f, 1.797e-6):
- label = 'ref. shear velocity'
- #elif numpy.isclose(mu_f, 1.204e-6):
- #label = 'ref. shear velocity$\\times$0.67'
- #elif numpy.isclose(mu_f, 1.797e-8):
- #label = 'ref. shear velocity$\\times$0.01'
- else:
- #label = '$\\mu_\\text{{f}}$ = {:.3e} Pa s'.format(mu_f)
- 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.3, color=color[c], clip_on=False)
- #alpha=0.2, color='gray')
- #alpha=alpha, color=color[c])
-
- # smoothed
- ax1.plot(shear_strain[c][1:-50], smooth(friction[c], smooth_window)[1:-50], \
- label=label, linewidth=1,
- alpha=alpha, color=color[c])
+ fig = plt.figure(figsize=(8,8)) # (w,h)
+ #fig = plt.figure(figsize=(8,12))
+ #fig = plt.figure(figsize=(8,16))
+
+ #plt.subplot(3,1,1)
+ #plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+
+ 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=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):
+ mu_f = mu_f_vals[c]
+
+ if numpy.isclose(mu_f, 1.797e-6):
+ label = 'ref. shear velocity'
+ #elif numpy.isclose(mu_f, 1.204e-6):
+ #label = 'ref. shear velocity$\\times$0.67'
+ #elif numpy.isclose(mu_f, 1.797e-8):
+ #label = 'ref. shear velocity$\\times$0.01'
+ else:
+ #label = '$\\mu_\\text{{f}}$ = {:.3e} Pa s'.format(mu_f)
+ 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.3, color=color[c], clip_on=False)
+ #alpha=0.2, color='gray')
+ #alpha=alpha, color=color[c])
+
+ # smoothed
+ ax1.plot(shear_strain[c][1:-50], smooth(friction[c], smooth_window)[1:-50], \
+ label=label, linewidth=1,
+ alpha=alpha, color=color[c])
+
+
+ ax2.plot(shear_strain[c], dilation[c], \
+ label=label, linewidth=1,
+ color=color[c])
+ if zflow:
+ ax3.plot(shear_strain[c], v_f_z_mean[c],
+ label=label, linewidth=1)
- ax2.plot(shear_strain[c], dilation[c], \
- label=label, linewidth=1,
- color=color[c])
-
- if zflow:
- ax3.plot(shear_strain[c], v_f_z_mean[c],
- label=label, linewidth=1)
+ if pressures:
+ #ax3.plot(shear_strain[c], p_max[c], '-', color=color[c], alpha=0.5)
- if pressures:
- #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], p_mean[c], '-', color=color[c], \
- label=label, linewidth=1)
+ #ax3.plot(shear_strain[c], p_min[c], '-', color=color[c], 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], edgecolor='None',
+ interpolate=True, alpha=0.3)
- ax3.fill_between(shear_strain[c], p_min[c], p_max[c],
- where=p_min[c]<=p_max[c], facecolor=color[c], edgecolor='None',
- interpolate=True, alpha=0.3)
-
- #ax4.plot(shear_strain[c][1:], f_n_mean[c][1:], '-' + color[c],
+ #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.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$ [-]')
+
+ 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('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('Fluid pressure $\\bar{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])
- #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)
- 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()
+ '''
- '''
- 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,
+ 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)
- #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])
- ax1.set_xlim([0.0, 1.0])
- ax1.set_ylim([0.0, 1.0])
- ax2.set_ylim([0.0, 1.0])
- ax3.set_ylim([-150., 150.])
-
- #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'
-filename = 'halfshear-darcy-rate.pdf'
-#print(os.getcwd() + '/' + filename)
-plt.savefig(filename)
-shutil.copyfile(filename, '/home/adc/articles/own/2/graphics/' + filename)
-plt.close()
-print(filename)
+ #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])
+ ax1.set_xlim([0.0, 1.0])
+ ax1.set_ylim([0.0, 1.0])
+ ax2.set_ylim([0.0, 1.0])
+ ax3.set_ylim([-150., 150.])
+
+ #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'
+ filename = 'halfshear-darcy-rate.pdf'
+ if sigma0 == 80000.0:
+ filename = 'halfshear-darcy-rate-N80.pdf'
+ #print(os.getcwd() + '/' + filename)
+ plt.savefig(filename)
+ shutil.copyfile(filename, '/home/adc/articles/own/2/graphics/' + filename)
+ plt.close()
+ print(filename)
(DIR) diff --git a/python/halfshear-darcy-strength-dilation.py b/python/halfshear-darcy-strength-dilation.py
t@@ -28,7 +28,7 @@ smooth_window = 100
#smooth_window = 200
#sigma0_list = numpy.array([1.0e3, 2.0e3, 4.0e3, 10.0e3, 20.0e3, 40.0e3])
-sigma0 = 20000.0
+sigma0_list = [20000.0, 80000.0]
#k_c_vals = [3.5e-13, 3.5e-15]
k_c = 3.5e-15
t@@ -96,365 +96,368 @@ def smooth(x, window_len=10, window='hanning'):
return y[window_len-1:-window_len+1]
-shear_strain = [[], [], [], []]
-friction = [[], [], [], []]
-dilation = [[], [], [], []]
-p_min = [[], [], [], []]
-p_mean = [[], [], [], []]
-p_max = [[], [], [], []]
-f_n_mean = [[], [], [], []]
-f_n_max = [[], [], [], []]
-v_f_z_mean = [[], [], [], []]
+for sigma0 in sigma0_list:
+ shear_strain = [[], [], [], []]
+ friction = [[], [], [], []]
+ dilation = [[], [], [], []]
+ p_min = [[], [], [], []]
+ p_mean = [[], [], [], []]
+ p_max = [[], [], [], []]
+ f_n_mean = [[], [], [], []]
+ f_n_max = [[], [], [], []]
+ v_f_z_mean = [[], [], [], []]
+
+ fluid=True
+
+ for c in numpy.arange(0,len(k_c_vals)):
+ k_c = k_c_vals[c]
+
+ if k_c == 'dry':
+ sid = 'halfshear-sigma0=' + str(sigma0) + '-shear'
+ fluid = False
+ else:
+ sid = 'halfshear-darcy-sigma0=' + str(sigma0) + '-k_c=' + str(k_c) + \
+ '-mu=' + str(mu_f) + '-velfac=' + str(velfac) + '-shear'
+ fluid = True
+ #sid = 'halfshear-sigma0=' + str(sigma0) + '-c_v=' + str(c_v) +\
+ #'-c_a=0.0-velfac=1.0-shear'
+ if os.path.isfile('../output/' + sid + '.status.dat'):
+
+ sim = sphere.sim(sid, fluid=fluid)
+ n = sim.status()
+ #n = 20
+ shear_strain[c] = numpy.zeros(n)
+ friction[c] = numpy.zeros_like(shear_strain[c])
+ dilation[c] = numpy.zeros_like(shear_strain[c])
+
+ # fluid pressures and particle forces
+ if fluid:
+ p_mean[c] = numpy.zeros_like(shear_strain[c])
+ p_min[c] = numpy.zeros_like(shear_strain[c])
+ p_max[c] = numpy.zeros_like(shear_strain[c])
+ if contact_forces:
+ f_n_mean[c] = numpy.zeros_like(shear_strain[c])
+ f_n_max[c] = numpy.zeros_like(shear_strain[c])
-fluid=True
+ for i in numpy.arange(n):
-for c in numpy.arange(0,len(k_c_vals)):
- k_c = k_c_vals[c]
+ sim.readstep(i, verbose=False)
- if k_c == 'dry':
- sid = 'halfshear-sigma0=' + str(sigma0) + '-shear'
- fluid = False
- else:
- sid = 'halfshear-darcy-sigma0=' + str(sigma0) + '-k_c=' + str(k_c) + \
- '-mu=' + str(mu_f) + '-velfac=' + str(velfac) + '-shear'
- fluid = True
- #sid = 'halfshear-sigma0=' + str(sigma0) + '-c_v=' + str(c_v) +\
- #'-c_a=0.0-velfac=1.0-shear'
- if os.path.isfile('../output/' + sid + '.status.dat'):
-
- sim = sphere.sim(sid, fluid=fluid)
- n = sim.status()
- #n = 20
- shear_strain[c] = numpy.zeros(n)
- friction[c] = numpy.zeros_like(shear_strain[c])
- dilation[c] = numpy.zeros_like(shear_strain[c])
-
- # fluid pressures and particle forces
- if fluid:
- p_mean[c] = numpy.zeros_like(shear_strain[c])
- p_min[c] = numpy.zeros_like(shear_strain[c])
- p_max[c] = numpy.zeros_like(shear_strain[c])
- if contact_forces:
- f_n_mean[c] = numpy.zeros_like(shear_strain[c])
- f_n_max[c] = numpy.zeros_like(shear_strain[c])
+ shear_strain[c][i] = sim.shearStrain()
+ friction[c][i] = sim.shearStress('effective')/sim.currentNormalStress('defined')
+ dilation[c][i] = sim.w_x[0]
- for i in numpy.arange(n):
+ if fluid and pressures:
+ iz_top = int(sim.w_x[0]/(sim.L[2]/sim.num[2]))-1
+ p_mean[c][i] = numpy.mean(sim.p_f[:,:,0:iz_top])/1000
+ p_min[c][i] = numpy.min(sim.p_f[:,:,0:iz_top])/1000
+ p_max[c][i] = numpy.max(sim.p_f[:,:,0:iz_top])/1000
- sim.readstep(i, verbose=False)
+ if contact_forces:
+ sim.findNormalForces()
+ f_n_mean[c][i] = numpy.mean(sim.f_n_magn)
+ f_n_max[c][i] = numpy.max(sim.f_n_magn)
- shear_strain[c][i] = sim.shearStrain()
- friction[c][i] = sim.shearStress('effective')/sim.currentNormalStress('defined')
- dilation[c][i] = sim.w_x[0]
+ if fluid and zflow:
+ v_f_z_mean[c] = numpy.zeros_like(shear_strain[c])
+ for i in numpy.arange(n):
+ v_f_z_mean[c][i] = numpy.mean(sim.v_f[:,:,:,2])
- if fluid and pressures:
- iz_top = int(sim.w_x[0]/(sim.L[2]/sim.num[2]))-1
- p_mean[c][i] = numpy.mean(sim.p_f[:,:,0:iz_top])/1000
- p_min[c][i] = numpy.min(sim.p_f[:,:,0:iz_top])/1000
- p_max[c][i] = numpy.max(sim.p_f[:,:,0:iz_top])/1000
+ dilation[c] =\
+ (dilation[c] - dilation[c][0])/(numpy.mean(sim.radius)*2.0)
- if contact_forces:
- sim.findNormalForces()
- f_n_mean[c][i] = numpy.mean(sim.f_n_magn)
- f_n_max[c][i] = numpy.max(sim.f_n_magn)
+ else:
+ print(sid + ' not found')
- if fluid and zflow:
- v_f_z_mean[c] = numpy.zeros_like(shear_strain[c])
- for i in numpy.arange(n):
- v_f_z_mean[c][i] = numpy.mean(sim.v_f[:,:,:,2])
+ # produce VTK files
+ #for sid in sids:
+ #sim = sphere.sim(sid, fluid=True)
+ #sim.writeVTKall()
- dilation[c] =\
- (dilation[c] - dilation[c][0])/(numpy.mean(sim.radius)*2.0)
+ if zflow or pressures:
+ #fig = plt.figure(figsize=(8,10))
+ #fig = plt.figure(figsize=(3.74, 2*3.74))
+ fig = plt.figure(figsize=(2*3.74, 2*3.74))
else:
- print(sid + ' not found')
+ fig = plt.figure(figsize=(8,8)) # (w,h)
+ #fig = plt.figure(figsize=(8,12))
+ #fig = plt.figure(figsize=(8,16))
+
+ #plt.subplot(3,1,1)
+ #plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+
+ for c in numpy.arange(0,len(k_c_vals)):
+
+ if zflow or pressures:
+ ax1 = plt.subplot(3, len(k_c_vals), 1+c)
+ ax2 = plt.subplot(3, len(k_c_vals), 5+c, sharex=ax1)
+ if c > 0:
+ ax3 = plt.subplot(3, len(k_c_vals), 9+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=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):
+ k_c = k_c_vals[c]
- # produce VTK files
- #for sid in sids:
- #sim = sphere.sim(sid, fluid=True)
- #sim.writeVTKall()
+ fluid = True
+ if k_c == 'dry':
+ label = 'dry'
+ fluid = False
+ elif numpy.isclose(k_c, 3.5e-13, atol=1.0e-16):
+ label = 'high permeability'
+ elif numpy.isclose(k_c, 3.5e-14, atol=1.0e-16):
+ label = 'interm. permeability'
+ elif numpy.isclose(k_c, 3.5e-15, atol=1.0e-16):
+ label = 'low permeability'
+ else:
+ label = str(k_c)
+
+ # unsmoothed
+ ax1.plot(shear_strain[c][1:], friction[c][1:], \
+ label=label, linewidth=1,
+ alpha=0.3, color=color[c], clip_on=False)
+ #alpha=0.2, color='gray', clip_on=False)
+ #alpha=alpha, color=color[c])
+
+ # smoothed
+ ax1.plot(shear_strain[c][1:-50], smooth(friction[c], smooth_window)[1:-50],\
+ label=label, linewidth=1,
+ alpha=alpha, color=color[c])
+
+
+ ax2.plot(shear_strain[c], dilation[c], \
+ label=label, linewidth=1,
+ color=color[c])
+ if zflow:
+ ax3.plot(shear_strain[c], v_f_z_mean[c],
+ label=label, linewidth=1)
-if zflow or pressures:
- #fig = plt.figure(figsize=(8,10))
- #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))
+ if fluid and pressures:
+ #ax3.plot(shear_strain[c], p_max[c], '-', color=color[c], alpha=0.5)
-#plt.subplot(3,1,1)
-#plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
+ ax3.plot(shear_strain[c], p_mean[c], '-', color=color[c], \
+ label=label, linewidth=1)
-for c in numpy.arange(0,len(k_c_vals)):
+ #ax3.plot(shear_strain[c], p_min[c], '-', color=color[c], alpha=0.5)
- if zflow or pressures:
- ax1 = plt.subplot(3, len(k_c_vals), 1+c)
- ax2 = plt.subplot(3, len(k_c_vals), 5+c, sharex=ax1)
- if c > 0:
- ax3 = plt.subplot(3, len(k_c_vals), 9+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=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):
- k_c = k_c_vals[c]
-
- fluid = True
- if k_c == 'dry':
- label = 'dry'
- fluid = False
- elif numpy.isclose(k_c, 3.5e-13, atol=1.0e-16):
- label = 'high permeability'
- elif numpy.isclose(k_c, 3.5e-14, atol=1.0e-16):
- label = 'interm. permeability'
- elif numpy.isclose(k_c, 3.5e-15, atol=1.0e-16):
- label = 'low permeability'
- else:
- label = str(k_c)
- # unsmoothed
- ax1.plot(shear_strain[c][1:], friction[c][1:], \
- label=label, linewidth=1,
- alpha=0.3, color=color[c], clip_on=False)
- #alpha=0.2, color='gray', clip_on=False)
- #alpha=alpha, color=color[c])
+ ax3.fill_between(shear_strain[c], p_min[c], p_max[c],
+ where=p_min[c]<=p_max[c], facecolor=color[c], edgecolor='None',
+ interpolate=True, alpha=0.3)
- # smoothed
- ax1.plot(shear_strain[c][1:-50], smooth(friction[c], smooth_window)[1:-50],\
- label=label, linewidth=1,
- alpha=alpha, color=color[c])
+ #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)
- ax2.plot(shear_strain[c], dilation[c], \
- label=label, linewidth=1,
- color=color[c])
- if zflow:
- ax3.plot(shear_strain[c], v_f_z_mean[c],
- label=label, linewidth=1)
+ #ax4.set_xlabel('Shear strain $\\gamma$ [-]')
+ if fluid and (zflow or pressures):
+ ax3.set_xlabel('Shear strain $\\gamma$ [-]')
+ else:
+ ax2.set_xlabel('Shear strain $\\gamma$ [-]')
- if fluid and pressures:
- #ax3.plot(shear_strain[c], p_max[c], '-', color=color[c], alpha=0.5)
+ 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)$ [-]')
- ax3.plot(shear_strain[c], p_mean[c], '-', color=color[c], \
- label=label, linewidth=1)
+ if c == 1:
+ if zflow:
+ ax3.set_ylabel('$\\boldsymbol{v}_\\text{f}^z h$ [ms$^{-1}$]')
+ if pressures:
+ ax3.set_ylabel('Fluid pressure $\\bar{p}_\\text{f}$ [kPa]')
+ #ax4.set_ylabel('Particle contact force $||\\boldsymbol{f}_\\text{p}||$ [N]')
- #ax3.plot(shear_strain[c], p_min[c], '-', color=color[c], alpha=0.5)
+ #ax1.set_xlim([200,300])
+ #ax3.set_ylim([595,608])
+ plt.setp(ax1.get_xticklabels(), visible=False)
+ if fluid and (zflow or pressures):
+ plt.setp(ax2.get_xticklabels(), visible=False)
+ #plt.setp(ax2.get_xticklabels(), visible=False)
+ #plt.setp(ax3.get_xticklabels(), visible=False)
- ax3.fill_between(shear_strain[c], p_min[c], p_max[c],
- where=p_min[c]<=p_max[c], facecolor=color[c], edgecolor='None',
- interpolate=True, alpha=0.3)
+ '''
+ ax1.grid()
+ ax2.grid()
+ if zflow or pressures:
+ ax3.grid()
+ #ax4.grid()
+ '''
- #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)
+ 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(True)
+ # 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().tick_bottom()
+
+ '''
+ # remove box at top and right
+ ax3.spines['top'].set_visible(False)
+ ax3.spines['left'].set_visible(False)
+ ax3.spines['bottom'].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')
+ plt.setp(ax3.get_xticklabels(), visible=False)
+ plt.setp(ax3.get_yticklabels(), visible=False)
+ '''
+
+ elif c == len(k_c_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_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)
+ if c == 1:
+ #ax3.get_yaxis().tick_left()
+ ax3.spines['left'].set_visible(True)
+ # vertical grid lines
+ ax1.get_xaxis().grid(True, linestyle=':', linewidth=0.5)
+ ax2.get_xaxis().grid(True, linestyle=':', linewidth=0.5)
+ if fluid:
+ ax3.get_xaxis().grid(True, linestyle=':', linewidth=0.5)
- #ax4.set_xlabel('Shear strain $\\gamma$ [-]')
- if fluid and (zflow or pressures):
- ax3.set_xlabel('Shear strain $\\gamma$ [-]')
- else:
- ax2.set_xlabel('Shear strain $\\gamma$ [-]')
- 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)$ [-]')
+ # horizontal grid lines
+ ax1.get_yaxis().grid(True, linestyle=':', linewidth=0.5)
+ ax2.get_yaxis().grid(True, linestyle=':', linewidth=0.5)
+ if fluid:
+ ax3.get_yaxis().grid(True, linestyle=':', linewidth=0.5)
- if c == 1:
- if zflow:
- ax3.set_ylabel('$\\boldsymbol{v}_\\text{f}^z h$ [ms$^{-1}$]')
- if pressures:
- ax3.set_ylabel('Fluid pressure $\\bar{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])
-
- plt.setp(ax1.get_xticklabels(), visible=False)
- if fluid and (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(True)
- # 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().tick_bottom()
+ 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)
- '''
- # remove box at top and right
- ax3.spines['top'].set_visible(False)
- ax3.spines['left'].set_visible(False)
- ax3.spines['bottom'].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')
- plt.setp(ax3.get_xticklabels(), visible=False)
- plt.setp(ax3.get_yticklabels(), visible=False)
- '''
+ #ax1.set_xlim([0.0, 0.09])
+ #ax2.set_xlim([0.0, 0.09])
+ #ax2.set_xlim([0.0, 0.2])
- elif c == len(k_c_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_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)
- if c == 1:
- #ax3.get_yaxis().tick_left()
- ax3.spines['left'].set_visible(True)
-
-
- # vertical grid lines
- ax1.get_xaxis().grid(True, linestyle=':', linewidth=0.5)
- ax2.get_xaxis().grid(True, linestyle=':', linewidth=0.5)
- if fluid:
- 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)
- if fluid:
- 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([0.0, 0.09])
- #ax2.set_xlim([0.0, 0.09])
- #ax2.set_xlim([0.0, 0.2])
-
- #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])
- if fluid:
- ax3.set_ylim([-150., 150.])
-
- #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'
-filename = 'halfshear-darcy.pdf'
-#print(os.getcwd() + '/' + filename)
-plt.savefig(filename)
-shutil.copyfile(filename, '/home/adc/articles/own/2/graphics/' + filename)
-plt.close()
-print(filename)
+ #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])
+ if fluid:
+ ax3.set_ylim([-150., 150.])
+
+ #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'
+ filename = 'halfshear-darcy.pdf'
+ if sigma0 == 80000.0:
+ filename = 'halfshear-darcy-N80.pdf'
+ #print(os.getcwd() + '/' + filename)
+ plt.savefig(filename)
+ shutil.copyfile(filename, '/home/adc/articles/own/2/graphics/' + filename)
+ plt.close()
+ print(filename)