tscale pressure gradient force with c_grad_p. improved plots - 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 3e077e59134905ebc4b377df485aeb0f422682a9
(DIR) parent fd79cc4d343e7c453c80d070111e68dacbb3f1d2
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Mon, 6 Oct 2014 13:22:05 +0200
scale pressure gradient force with c_grad_p. improved plots
Diffstat:
M python/permeability-results.py | 1 +
M python/shear-results-forces.py | 20 ++++++++++++++------
M python/shear-results-pressures.py | 43 +++++++++++++++++--------------
M python/shear-results.py | 6 +++---
M src/device.cu | 1 +
M src/navierstokes.cuh | 3 ++-
6 files changed, 45 insertions(+), 29 deletions(-)
---
(DIR) diff --git a/python/permeability-results.py b/python/permeability-results.py
t@@ -145,6 +145,7 @@ ax2.legend(loc='best', prop={'size':18}, fancybox=True, framealpha=0.5)
ax3.legend(loc='best', prop={'size':18}, fancybox=True, framealpha=0.5)
plt.tight_layout()
+plt.subplots_adjust(hspace = .12)
filename = 'permeability-dpdz-vs-K-vs-c.pdf'
#print(os.getcwd() + '/' + filename)
plt.savefig(filename)
(DIR) diff --git a/python/shear-results-forces.py b/python/shear-results-forces.py
t@@ -144,20 +144,28 @@ for s in numpy.arange(len(steps)):
#ax[s*4+1].plot([0.0, 0.0], [0.0, sim.L[2]], '--', color='k')
#ax[s*4+2].plot(dev_p[s]/1000.0, zpos_c[s], 'k')
- ax[s*4+2].plot(phi_bar[s,1:], zpos_c[s,1:], '-k', linewidth=3)
+ #ax[s*4+2].plot(phi_bar[s,1:], zpos_c[s,1:], '-k', linewidth=3)
+ ax[s*4+2].plot(phi_bar[s,1:], zpos_c[s,1:], '-k')
#phicolor = '#888888'
#ax[s*4+3].plot(phi_bar[s], zpos_c[s], '-', color = phicolor)
#for tl in ax[s*4+3].get_xticklabels():
#tl.set_color(phicolor)
- ax[s*4+3].plot(dphi_bar[s,1:], zpos_c[s,1:], '-k', linewidth=3)
- ax[s*4+3].plot(dphi_bar[s,1:], zpos_c[s,1:], '-w', linewidth=2)
+ ax[s*4+3].plot(dphi_bar[s,1:], zpos_c[s,1:], '--k')
+ #ax[s*4+3].plot(dphi_bar[s,1:], zpos_c[s,1:], '-k', linewidth=3)
+ #ax[s*4+3].plot(dphi_bar[s,1:], zpos_c[s,1:], '-w', linewidth=2)
max_z = numpy.max(zpos_p)
ax[s*4+0].set_ylim([0, max_z])
- ax[s*4+1].set_xlim([0.15, 0.46])
+
+ #ax[s*4+1].set_xlim([0.15, 0.46]) # f_pf
+ ax[s*4+1].set_xlim([0.235, 0.409]) # f_pf
ax[s*4+1].set_ylim([0, max_z])
+
ax[s*4+2].set_ylim([0, max_z])
+ ax[s*4+2].set_xlim([0.33, 0.6]) # phi
+ ax[s*4+3].set_xlim([-0.09, 0.024]) # dphi/dt
+
#plt.plot(dpdz[c], K[c], 'o-', label='$c$ = %.2f' % (cvals[c]))
#plt.semilogx(dpdz[c], K[c], 'o-', label='$c$ = %.2f' % (cvals[c]))
#plt.semilogy(dpdz[c], K[c], 'o-', label='$c$ = %.2f' % (cvals[c]))
t@@ -168,8 +176,8 @@ for s in numpy.arange(len(steps)):
ax[s*4+1].set_xlabel('$\\boldsymbol{f}^z_\\text{pf}$ [N]')
#ax[s*4+2].set_xlabel('$\\bar{p_\\text{f}}$ [kPa]')
#ax[s*4+3].set_xlabel('$\\bar{\\phi}$ [-]', color=phicolor)
- ax[s*4+2].set_xlabel('$\\bar{\\phi}$ [-]')
- ax[s*4+3].set_xlabel('$\\delta \\bar{\\phi}/\\delta t$ [-]')
+ ax[s*4+2].set_xlabel('$\\bar{\\phi}$ [-] (solid)')
+ ax[s*4+3].set_xlabel('$\\delta \\bar{\\phi}/\\delta t$ [-] (dashed)')
plt.setp(ax[s*4+1].get_yticklabels(), visible=False)
plt.setp(ax[s*4+2].get_yticklabels(), visible=False)
(DIR) diff --git a/python/shear-results-pressures.py b/python/shear-results-pressures.py
t@@ -35,7 +35,7 @@ for i in numpy.arange(sim.num[2]):
shear_strain = numpy.zeros(sim.status())
dev_pres = numpy.zeros((sim.num[2], sim.status()))
-pres_static = numpy.ones_like(dev_pres)*100.0e3
+pres_static = numpy.ones_like(dev_pres)*600.0e3
pres = numpy.zeros_like(dev_pres)
for i in numpy.arange(sim.status()):
t@@ -44,10 +44,12 @@ for i in numpy.arange(sim.status()):
pres[:,i] = numpy.average(numpy.average(sim.p_f, axis=0), axis=0)
- wall0_iz = int(sim.w_x[0]/(sim.L[2]/sim.num[2]))
+ dz = sim.L[2]/sim.num[2]
+ wall0_iz = int(sim.w_x[0]/dz)
for z in numpy.arange(0, wall0_iz+1):
+ #(wall0_iz*dz - zpos_c[z] + 0.5*dz)*sim.rho_f*numpy.abs(sim.g[2])\
pres_static[z,i] = \
- (sim.w_x[0] - zpos_c[z])*sim.rho_f*numpy.abs(sim.g[2])\
+ (wall0_iz*dz - zpos_c[z])*sim.rho_f*numpy.abs(sim.g[2])\
+ sim.p_f[0,0,-1]
#pres_static[z,i] = zpos_c[z]
#pres_static[z,i] = z
t@@ -56,8 +58,9 @@ for i in numpy.arange(sim.status()):
dev_pres = pres - pres_static
-fig = plt.figure(figsize=(8,6))
-#fig = plt.figure(figsize=(8,15))
+#fig = plt.figure(figsize=(8,6))
+#fig = plt.figure(figsize=(8,12))
+fig = plt.figure(figsize=(8,15))
min_p = numpy.min(dev_pres)/1000.0
#max_p = numpy.min(dev_pres)
t@@ -67,14 +70,15 @@ max_p = numpy.abs(min_p)
#bounds = [min_p, 0, max_p]
#norm = matplotlib.colors.BoundaryNorm(bounds, cmap.N)
-#ax1 = plt.subplot(311)
-ax1 = plt.subplot(111)
+ax1 = plt.subplot(311)
+#ax1 = plt.subplot(111)
+#ax1 = plt.subplot(211)
#im1 = ax1.pcolormesh(shear_strain, zpos_c, dev_pres/1000.0, rasterized=True,
# cmap=cmap, norm=norm)
im1 = ax1.pcolormesh(shear_strain, zpos_c, dev_pres/1000.0, vmin=min_p,
vmax=max_p, rasterized=True)
-ax1.set_xlim([0, shear_strain[-1]])
-ax1.set_ylim([zpos_c[0], sim.w_x[0]])
+#ax1.set_xlim([0, shear_strain[-1]])
+#ax1.set_ylim([zpos_c[0], sim.w_x[0]])
ax1.set_xlabel('Shear strain $\\gamma$ [-]')
ax1.set_ylabel('Vertical position $z$ [m]')
cb1 = plt.colorbar(im1)
t@@ -83,33 +87,34 @@ cb1.set_label('$p_\\text{f} - p^\\text{hyd}_\\text{f}$ [kPa]')
cb1.solids.set_rasterized(True)
# annotate plot
-ax1.text(0.02, 0.15, 'compressive',
- bbox={'facecolor':'white', 'alpha':0.5, 'pad':10})
+#ax1.text(0.02, 0.15, 'compressive',
+ #bbox={'facecolor':'white', 'alpha':0.5, 'pad':10})
-ax1.text(0.12, 0.25, 'dilative',
- bbox={'facecolor':'white', 'alpha':0.5, 'pad':10})
+#ax1.text(0.12, 0.25, 'dilative',
+ #bbox={'facecolor':'white', 'alpha':0.5, 'pad':10})
-'''
+#'''
ax2 = plt.subplot(312)
im2 = ax2.pcolormesh(shear_strain, zpos_c, pres/1000.0, rasterized=True)
-ax2.set_xlim([0, shear_strain[-1]])
-ax2.set_ylim([zpos_c[0], sim.w_x[0]])
+#ax2.set_xlim([0, shear_strain[-1]])
+#ax2.set_ylim([zpos_c[0], sim.w_x[0]])
ax2.set_xlabel('Shear strain $\\gamma$ [-]')
ax2.set_ylabel('Vertical position $z$ [m]')
cb2 = plt.colorbar(im2)
cb2.set_label('Pressure $p_\\text{f}$ [kPa]')
cb2.solids.set_rasterized(True)
+#'''
ax3 = plt.subplot(313)
im3 = ax3.pcolormesh(shear_strain, zpos_c, pres_static/1000.0, rasterized=True)
-ax3.set_xlim([0, shear_strain[-1]])
-ax3.set_ylim([zpos_c[0], sim.w_x[0]])
+#ax3.set_xlim([0, shear_strain[-1]])
+#ax3.set_ylim([zpos_c[0], sim.w_x[0]])
ax3.set_xlabel('Shear strain $\\gamma$ [-]')
ax3.set_ylabel('Vertical position $z$ [m]')
cb3 = plt.colorbar(im3)
cb3.set_label('Static Pressure $p_\\text{f}$ [kPa]')
cb3.solids.set_rasterized(True)
-'''
+#'''
#plt.MaxNLocator(nbins=4)
(DIR) diff --git a/python/shear-results.py b/python/shear-results.py
t@@ -16,8 +16,8 @@ import matplotlib.pyplot as plt
#sigma0_list = numpy.array([1.0e3, 2.0e3, 4.0e3, 10.0e3, 20.0e3, 40.0e3])
#sigma0 = 10.0e3
sigma0 = float(sys.argv[1])
-#cvals = [1.0, 0.1]
-cvals = [1.0]
+cvals = [1.0, 0.1]
+#cvals = [1.0]
shear_strain = [[], [], []]
friction = [[], [], []]
t@@ -169,7 +169,7 @@ ax4.legend(loc='best', prop={'size':18}, fancybox=True,
framealpha=legend_alpha)
plt.tight_layout()
-plt.subplots_adjust(hspace=0.0)
+plt.subplots_adjust(hspace=0.05)
filename = 'shear-' + str(int(sigma0/1000.0)) + 'kPa-stress-dilation.pdf'
#print(os.getcwd() + '/' + filename)
plt.savefig(filename)
(DIR) diff --git a/src/device.cu b/src/device.cu
t@@ -1125,6 +1125,7 @@ __host__ void DEM::startTime()
dev_ns_div_tau_x,
dev_ns_div_tau_y,
dev_ns_div_tau_z,
+ ns.c_grad_p,
dev_ns_f_pf,
dev_force,
dev_ns_f_d,
(DIR) diff --git a/src/navierstokes.cuh b/src/navierstokes.cuh
t@@ -3132,6 +3132,7 @@ __global__ void findInteractionForce(
const Float* __restrict__ dev_ns_div_tau_x,// in
const Float* __restrict__ dev_ns_div_tau_y,// in
const Float* __restrict__ dev_ns_div_tau_z,// in
+ const Float c_grad_p, // in
Float3* __restrict__ dev_ns_f_pf, // out
Float4* __restrict__ dev_force, // out
Float4* __restrict__ dev_ns_f_d, // out
t@@ -3211,7 +3212,7 @@ __global__ void findInteractionForce(
// Pressure gradient force
const Float3 f_p =
- -1.0*gradient(dev_ns_p, i_x, i_y, i_z, dx, dy, dz)*V_p;
+ -c_grad_p*gradient(dev_ns_p, i_x, i_y, i_z, dx, dy, dz)*V_p;
// Viscous force
const Float3 f_v = -1.0*div_tau*V_p;