timprove internals plot - 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 fd4eb344b8af1b8de25c5a2c587c01f94fd9adb4
(DIR) parent 9d9e5ba4f59c0dbafc66116261bd5766e48f380f
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Thu, 16 Apr 2015 11:48:55 +0200
improve internals plot
Diffstat:
M python/halfshear-darcy-internals.py | 373 +++++++++++++++++--------------
M python/sphere.py | 9 +++++++--
2 files changed, 207 insertions(+), 175 deletions(-)
---
(DIR) diff --git a/python/halfshear-darcy-internals.py b/python/halfshear-darcy-internals.py
t@@ -13,213 +13,240 @@ from permeabilitycalculator import *
import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocator
+import seaborn as sns
+#sns.set(style='ticks', palette='Set2')
+#sns.set(style='ticks', palette='colorblind')
+#sns.set(style='ticks', palette='Set2')
+sns.set(style='ticks', palette='Blues')
+
#sigma0 = float(sys.argv[1])
sigma0 = 20000.0
#k_c = 3.5e-13
-k_c = float(sys.argv[1])
-
-if k_c == 3.5e-15:
- steps = [1232, 1332, 1433, 1534, 1635]
-elif k_c == 3.5e-13:
- steps = [100, 200, 300, 410, 515]
-else:
- steps = [10, 50, 100, 1000, 1999]
+#k_c = float(sys.argv[1])
+k_c_list = [3.5e-13, 3.5e-14, 3.5e-15]
+
+#if k_c == 3.5e-15:
+# steps = [1232, 1332, 1433, 1534, 1635]
+#elif k_c == 3.5e-13:
+# steps = [100, 200, 300, 410, 515]
+#else:
+# steps = [10, 50, 100, 1000, 1999]
nsteps_avg = 1 # no. of steps to average over
#nsteps_avg = 100 # no. of steps to average over
+steps = [1, 50, 100, 150, 200]
-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)
+for k_c in k_c_list:
-# particle z positions
-zpos_p = numpy.zeros((len(steps), sim.np))
+ 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)
-# 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 z positions
+ zpos_p = numpy.zeros((len(steps), sim.np))
-# particle x displacements
-xdisp = numpy.zeros((len(steps), sim.np))
+ # cell midpoint cell positions
+ zpos_c = numpy.zeros((len(steps), sim.num[2]))
+ dz = sim.L[2]/sim.num[2]
+ for i in numpy.arange(sim.num[2]):
+ zpos_c[:,i] = i*dz + 0.5*dz
-# particle z velocity
-v_z_p = numpy.zeros((len(steps), sim.np))
+ # particle x displacements
+ xdisp = 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]))
+ # particle z velocity
+ v_z_p = numpy.zeros((len(steps), sim.np))
-# pressure
-p = 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]))
-# 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]))
+ # pressure
+ p = numpy.zeros((len(steps), sim.num[2]))
-# particle-fluid force per particle
-f_pf = numpy.zeros_like(xdisp)
+ # 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]))
-# pressure - hydrostatic pressure
-#dev_p = numpy.zeros((len(steps), sim.num[2]))
+ # particle-fluid force per particle
+ f_pf = numpy.zeros_like(xdisp)
-# mean porosity
-phi_bar = numpy.zeros((len(steps), sim.num[2]))
+ # pressure - hydrostatic pressure
+ #dev_p = numpy.zeros((len(steps), sim.num[2]))
-# mean porosity change
-dphi_bar = numpy.zeros((len(steps), sim.num[2]))
+ # mean porosity
+ phi_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 porosity change
+ dphi_bar = numpy.zeros((len(steps), sim.num[2]))
-shear_strain_start = numpy.zeros(len(steps))
-shear_strain_end = numpy.zeros(len(steps))
+ # mean per-particle values
+ xdisp_mean = numpy.zeros((len(steps), sim.num[2]))
+ f_pf_mean = numpy.zeros((len(steps), sim.num[2]))
-#fig = plt.figure(figsize=(8,4*(len(steps))+1))
-fig = plt.figure(figsize=(8,4.5))
-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
+ shear_strain_start = numpy.zeros(len(steps))
+ shear_strain_end = numpy.zeros(len(steps))
-s = 0
-for step_str in steps:
+ #fig = plt.figure(figsize=(8,4*(len(steps))+1))
+ #fig = plt.figure(figsize=(8,4.5))
+ fig = plt.figure(figsize=(3.74*2,4.5))
+ 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
- step = int(step_str)
+ s = 0
+ for step_str in steps:
- if os.path.isfile('../output/' + sid + '.status.dat'):
+ step = int(step_str)
- for substep in numpy.arange(nsteps_avg):
+ if os.path.isfile('../output/' + sid + '.status.dat'):
- if step + substep > sim.status():
- raise Exception(
- 'Simulation step %d not available (sim.status = %d).'
- % (step + substep, sim.status()))
+ for substep in numpy.arange(nsteps_avg):
- sim.readstep(step + substep, verbose=False)
+ if step + substep > sim.status():
+ raise Exception(
+ 'Simulation step %d not available (sim.status = %d).'
+ % (step + substep, sim.status()))
- zpos_p[s,:] += sim.x[:,2]/nsteps_avg
+ sim.readstep(step + substep, verbose=False)
- xdisp[s,:] += sim.xyzsum[:,0]/nsteps_avg
+ zpos_p[s,:] += sim.x[:,2]/nsteps_avg
- '''
- for i in numpy.arange(sim.np):
- f_pf[s,i] += \
- sim.f_sum[i].dot(sim.f_sum[i])/nsteps_avg
- '''
- f_pf[s,:] += sim.f_p[:,2]
+ xdisp[s,:] += sim.xyzsum[:,0]/nsteps_avg
- dz = sim.L[2]/sim.num[2]
- wall0_iz = int(sim.w_x[0]/dz)
+ '''
+ 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]
- p[s,:] += numpy.average(numpy.average(sim.p_f[:,:,:], axis=0),\
- axis=0)/nsteps_avg
+ dz = sim.L[2]/sim.num[2]
+ wall0_iz = int(sim.w_x[0]/dz)
- sim.findPermeabilities()
- k[s,:] += sim.k[:,:,:]/nsteps_avg
+ p[s,:] += numpy.average(numpy.average(sim.p_f[:,:,:], axis=0),\
+ axis=0)/nsteps_avg
- k_bar[s,:] += \
- numpy.average(numpy.average(sim.k[:,:,:], axis=0), axis=0)\
- /nsteps_avg
+ sim.findPermeabilities()
+ k[s,:] += sim.k[:,:,:]/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,:] += \
+ numpy.average(numpy.average(sim.k[:,:,:], axis=0), axis=0)\
+ /nsteps_avg
- #ax[0].plot(xdisp[s], zpos_p[s], ',', color = '#888888')
- ax[0].plot(xdisp_mean[s], zpos_c[s], label='$\gamma$ = %.2f' %
- (shear_strain_start[s]))
+ 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])
- ax[1].semilogx(k_bar[s], zpos_c[s], label='$\gamma$ = %.2f' %
- (shear_strain_start[s]))
-
- ax[2].plot(p[s]/1000.0, zpos_c[s], label='$\gamma$ = %.2f' %
- (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]
-
- #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$ = %.2f' %
- (shear_strain_start[s]))
-
- else:
- print(sid + ' not found')
- s += 1
-
-
-
-max_z = numpy.max(zpos_p)
-ax[0].set_ylim([0, max_z])
-ax[0].set_xlim([0, 0.5])
-
-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-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()
-
-legend_alpha=0.5
-ax[0].legend(loc='lower center', prop={'size':12}, fancybox=True,
- framealpha=legend_alpha)
-
-#plt.subplots_adjust(wspace = .05) # doesn't work with tight_layout()
-plt.tight_layout()
-#plt.MaxNLocator(nbins=1) # doesn't work?
-ax[0].locator_params(nbins=3)
-ax[2].locator_params(nbins=3)
-ax[3].locator_params(nbins=3)
-
-filename = 'halfshear-darcy-internals-k_c=%.0e.pdf' % (k_c)
-plt.savefig(filename)
-shutil.copyfile(filename, '/home/adc/articles/own/2/graphics/' + filename)
-print(filename)
+ 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$ = %.2f' %
+ (shear_strain_start[s]))
+
+ ax[1].semilogx(k_bar[s], zpos_c[s], label='$\gamma$ = %.2f' %
+ (shear_strain_start[s]))
+
+ ax[2].plot(p[s]/1000.0, zpos_c[s], label='$\gamma$ = %.2f' %
+ (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]
+
+ #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$ = %.2f' %
+ (shear_strain_start[s]))
+
+ else:
+ print(sid + ' not found')
+ s += 1
+
+
+ max_z = numpy.max(zpos_p)
+ 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()
+ '''
+
+ 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)
+
+ legend_alpha=0.5
+ ax[0].legend(loc='lower center', prop={'size':12}, fancybox=True,
+ framealpha=legend_alpha)
+
+ #plt.subplots_adjust(wspace = .05) # doesn't work with tight_layout()
+ plt.tight_layout()
+ #plt.MaxNLocator(nbins=1) # doesn't work?
+ ax[0].locator_params(nbins=5)
+ ax[2].locator_params(nbins=5)
+ ax[3].locator_params(nbins=5)
+
+ sns.despine() # remove chartjunk
+
+ filename = 'halfshear-darcy-internals-k_c=%.0e.pdf' % (k_c)
+ plt.savefig(filename)
+ shutil.copyfile(filename, '/home/adc/articles/own/2/graphics/' + filename)
+ print(filename)
(DIR) diff --git a/python/sphere.py b/python/sphere.py
t@@ -6785,7 +6785,8 @@ class sim:
i_max = sb.status()
# use largest difference in p from 0 as +/- limit on colormap
#print i_min, i_max
- p_ext = numpy.max(numpy.abs(pres))
+ #p_ext = numpy.max(numpy.abs(pres))
+ p_ext = numpy.max(numpy.abs(pres[0:9,:])) # for article2
if sb.wmode[0] == 3:
x = t
t@@ -6803,7 +6804,8 @@ class sim:
else:
im1 = ax.pcolormesh(
x, zpos_c, pres,
- cmap=matplotlib.cm.get_cmap('bwr'),
+ #cmap=matplotlib.cm.get_cmap('bwr'),
+ cmap=matplotlib.cm.get_cmap('RdBu_r'),
#cmap=matplotlib.cm.get_cmap('coolwarm'),
vmin=-p_ext, vmax=p_ext,
rasterized=True)
t@@ -6824,6 +6826,9 @@ class sim:
if xlim:
ax.set_xlim([x[0], x[-1]])
+ # for article2
+ ax.set_ylim([zpos_c[0], zpos_c[9]])
+
cb = plt.colorbar(im1)
cb.set_label('$p_\\text{f}$ [kPa]')
cb.solids.set_rasterized(True)