timprove plots, remove stray tabs - 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 b87092360b0bc18e75a6850737beb40dc60373f1
(DIR) parent 57d5ce5829d8daebbb03d7ce6e164b36fec3ff88
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Mon, 20 Oct 2014 13:29:11 +0200
improve plots, remove stray tabs
Diffstat:
M python/shear-results-internals.py | 58 ++++++++++++++++++-------------
M python/shear-results.py | 13 ++++++++-----
M src/device.cu | 14 +++++++-------
3 files changed, 48 insertions(+), 37 deletions(-)
---
(DIR) diff --git a/python/shear-results-internals.py b/python/shear-results-internals.py
t@@ -16,8 +16,8 @@ from matplotlib.ticker import MaxNLocator
#steps = [5, 10, 100]
#steps = [5, 10]
steps = sys.argv[3:]
-nsteps_avg = 5 # no. of steps to average over
-#nsteps_avg = 100 # no. of steps to average over
+#nsteps_avg = 5 # no. of steps to average over
+nsteps_avg = 100 # no. of steps to average over
sigma0 = float(sys.argv[1])
#c_grad_p = 1.0
t@@ -72,7 +72,8 @@ dphi_bar = numpy.zeros((len(steps), sim.num[2]))
xdisp_mean = numpy.zeros((len(steps), sim.num[2]))
f_pf_mean = numpy.zeros((len(steps), sim.num[2]))
-shear_strain = numpy.zeros(len(steps))
+shear_strain_start = numpy.zeros(len(steps))
+shear_strain_end = numpy.zeros(len(steps))
s = 0
for step_str in steps:
t@@ -130,7 +131,10 @@ for step_str in steps:
/nsteps_avg/sim.time_dt
- shear_strain[s] += sim.shearStrain()/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]):
t@@ -151,6 +155,7 @@ for step_str in steps:
fig = plt.figure(figsize=(16,5*(len(steps))+1))
def color(c):
+ return 'black'
if c == 1.0:
return 'green'
elif c == 0.1:
t@@ -166,11 +171,11 @@ for s in numpy.arange(len(steps)):
#strain_str = 'Shear strain\n $\\gamma = %.3f$' % (shear_strain[s])
if s == 0:
- strain_str = 'Dilating state\n $\\gamma = %.2f$, $c = %.2f$' % \
- (shear_strain[s], c_grad_p)
+ strain_str = 'Dilating state\n$\\gamma = %.2f$ to $%.2f$\n$c = %.2f$' % \
+ (shear_strain_start[s], shear_strain_end[s], c_grad_p)
else:
- strain_str = 'Steady state\n $\\gamma = %.2f$, $c = %.2f$' % \
- (shear_strain[s], c_grad_p)
+ strain_str = 'Steady state\n$\\gamma = %.2f$ to $%.2f$\n$c = %.2f$' % \
+ (shear_strain_start[s], shear_strain_end[s], c_grad_p)
n = 7
if s == 0:
t@@ -195,7 +200,7 @@ for s in numpy.arange(len(steps)):
ax.append(plt.subplot(len(steps), n-1, s*(n-1)+6, sharey=ax[s*n+0],
sharex=ax[6])) # 6: v_z^f
- #ax[s*n+0].plot(xdisp[s], zpos_p[s], ',', color = '#888888')
+ ax[s*n+0].plot(xdisp[s], zpos_p[s], ',', color = '#888888')
ax[s*n+0].plot(xdisp_mean[s], zpos_c[s], color=color(c_grad_p))
#ax[s*4+2].plot(dev_p[s]/1000.0, zpos_c[s], 'k')
t@@ -211,20 +216,21 @@ for s in numpy.arange(len(steps)):
#ax[s*4+3].plot(dphi_bar[s,1:], zpos_c[s,1:], '-w', linewidth=2)
ax[s*n+3].plot(v_z_p[s]*100.0, zpos_p[s], ',', alpha=0.5,
- color=color(c_grad_p))
+ color='#888888')
+ #color=color(c_grad_p))
ax[s*n+3].plot(v_z_p_bar[s]*100.0, zpos_c[s], color=color(c_grad_p))
#ax[s*n+0].plot([0.0,0.0], [0.0, sim.L[2]], '--', color='k')
# hydrostatic pressure distribution
- ax[s*n+4].plot(dev_p[s]/1000.0, zpos_c[s], color=color(c_grad_p))
- #ax[s*n+4].plot(p[s]/1000.0, zpos_c[s], color=color(c_grad_p))
- #dz = sim.L[2]/sim.num[2]
- #wall0_iz = int(sim.w_x[0]/dz)
- #y_top = wall0_iz*dz + 0.5*dz
- #x_top = sim.p_f[0,0,-1]
- #y_bot = 0.0
- #x_bot = x_top + (wall0_iz*dz - zpos_c[s][0] + 0.5*dz)*sim.rho_f*numpy.abs(sim.g[2])
- #ax[s*n+4].plot([x_top/1000.0, x_bot/1000.0], [y_top, y_bot], '--', color='k')
+ #ax[s*n+4].plot(dev_p[s]/1000.0, zpos_c[s], color=color(c_grad_p))
+ ax[s*n+4].plot(p[s]/1000.0, zpos_c[s], color=color(c_grad_p))
+ dz = sim.L[2]/sim.num[2]
+ wall0_iz = int(sim.w_x[0]/dz)
+ y_top = wall0_iz*dz + 0.5*dz
+ x_top = sim.p_f[0,0,-1]
+ y_bot = 0.0
+ x_bot = x_top + (wall0_iz*dz - zpos_c[s][0] + 0.5*dz)*sim.rho_f*numpy.abs(sim.g[2])
+ ax[s*n+4].plot([x_top/1000.0, x_bot/1000.0], [y_top, y_bot], '--', color='k')
#ax[s*n+1].set_title(strain_str)
#ax[s*n+1].set_title(' ')
t@@ -237,7 +243,8 @@ for s in numpy.arange(len(steps)):
zpos_c_nonzero = zpos_c[s][I]
ax[s*n+5].plot(f_pf_nonzero, zpos_p_nonzero, ',', alpha=0.5,
- color=color(c_grad_p))
+ color='#888888')
+ #color=color(c_grad_p))
#ax[s*4+1].plot(f_pf_mean[s][1:-2], zpos_c[s][1:-2], color = 'k')
ax[s*n+5].plot(f_pf_mean_nonzero, zpos_c_nonzero, color=color(c_grad_p))
#ax[s*4+1].plot([0.0, 0.0], [0.0, sim.L[2]], '--', color='k')
t@@ -276,9 +283,10 @@ for s in numpy.arange(len(steps)):
#ax[s*4+1].set_xlim([0.15, 0.46]) # f_pf
#ax[s*n+1].set_xlim([0.235, 0.409]) # f_pf
- ax[s*n+1].set_xlim([0.33, 0.6]) # phi
+ ax[s*n+1].set_xlim([0.33, 0.6]) # phi
ax[s*n+2].set_xlim([-0.09, 0.035]) # dphi/dt
- ax[s*n+3].set_xlim([-1.50, 1.50]) # v_z_p
+ ax[s*n+3].set_xlim([-1.50, 1.50]) # v_z_p
+ ax[s*n+5].set_xlim([5.0, 8.0]) # f_z_pf
#plt.plot(dpdz[c], K[c], 'o-', label='$c$ = %.2f' % (cvals[c]))
t@@ -291,8 +299,8 @@ for s in numpy.arange(len(steps)):
ax[s*n+1].set_xlabel('$\\bar{\\phi}$ [-] (solid)')
ax[s*n+2].set_xlabel('$\\delta \\bar{\\phi}/\\delta t$ [-] (dashed)')
ax[s*n+3].set_xlabel('$\\boldsymbol{v}^z_\\text{p}$ [cms$^{-1}$]')
- #ax[s*n+4].set_xlabel('$\\bar{p_\\text{f}}$ [kPa]')
- ax[s*n+4].set_xlabel('$\\bar{p_\\text{f}} - p_\\text{hyd}$ [kPa]')
+ ax[s*n+4].set_xlabel('$\\bar{p_\\text{f}}$ [kPa]')
+ #ax[s*n+4].set_xlabel('$\\bar{p_\\text{f}} - p_\\text{hyd}$ [kPa]')
ax[s*n+5].set_xlabel('$\\boldsymbol{f}^z_\\text{pf}$ [N]')
ax[s*n+6].set_xlabel('$\\bar{\\boldsymbol{v}}^z_\\text{f}$ [cms$^{-1}$]')
t@@ -364,7 +372,7 @@ for s in numpy.arange(len(steps)):
#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.38, 1.15, strain_str, horizontalalignment='left', fontsize=22,
+ plt.text(-0.38, 1.10, strain_str, horizontalalignment='left', fontsize=22,
transform=ax[s*n+0].transAxes)
#plt.title(' ')
(DIR) diff --git a/python/shear-results.py b/python/shear-results.py
t@@ -165,7 +165,6 @@ for c in numpy.arange(1,len(cvals)+1):
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])
- v_f_z_mean[c] = numpy.zeros_like(shear_strain[c])
for i in numpy.arange(sim.status()):
if pressures:
sim.readstep(i, verbose=False)
t@@ -179,7 +178,9 @@ for c in numpy.arange(1,len(cvals)+1):
f_n_mean[c][i] = numpy.mean(sim.f_n_magn)
f_n_max[c][i] = numpy.max(sim.f_n_magn)
- if zflow:
+ if zflow:
+ v_f_z_mean[c] = numpy.zeros_like(shear_strain[c])
+ for i in numpy.arange(sim.status()):
v_f_z_mean[c][i] = numpy.mean(sim.v_f[:,:,:,2])
else:
t@@ -270,7 +271,8 @@ plt.setp(ax1.get_xticklabels(), visible=False)
ax1.grid()
ax2.grid()
-#ax3.grid()
+if zflow:
+ ax3.grid()
#ax4.grid()
legend_alpha=0.5
t@@ -278,8 +280,9 @@ ax1.legend(loc='lower right', prop={'size':18}, fancybox=True,
framealpha=legend_alpha)
ax2.legend(loc='lower right', prop={'size':18}, fancybox=True,
framealpha=legend_alpha)
-#ax3.legend(loc='lower right', prop={'size':18}, fancybox=True,
- #framealpha=legend_alpha)
+if zflow:
+ ax3.legend(loc='lower right', prop={'size':18}, fancybox=True,
+ framealpha=legend_alpha)
#ax4.legend(loc='best', prop={'size':18}, fancybox=True,
#framealpha=legend_alpha)
(DIR) diff --git a/src/device.cu b/src/device.cu
t@@ -6,8 +6,8 @@
#include <cuda.h>
#include <helper_math.h>
-#include "vector_arithmetic.h" // for arbitrary prec. vectors
-//#include <vector_functions.h> // for single prec. vectors
+#include "vector_arithmetic.h" // for arbitrary prec. vectors
+//#include <vector_functions.h> // for single prec. vectors
#include "thrust/device_ptr.h"
#include "thrust/sort.h"
t@@ -17,7 +17,7 @@
#include "constants.cuh"
#include "debug.h"
-#include "sorting.cuh"
+#include "sorting.cuh"
#include "contactmodels.cuh"
#include "cohesion.cuh"
#include "contactsearch.cuh"
t@@ -223,7 +223,7 @@ __host__ void DEM::checkConstantMemory()
// Compare values between global and constant memory
// structures on the device.
- int* equal = new int; // The values are equal = 0, if not = 1
+ int* equal = new int; // The values are equal = 0, if not = 1
*equal = 0;
int* dev_equal;
cudaMalloc((void**)&dev_equal, sizeof(int));
t@@ -548,15 +548,15 @@ __host__ void DEM::transferToGlobalDeviceMemory(int statusmsg)
//cudaMemcpy(dev_time, &time, sizeof(Time), cudaMemcpyHostToDevice);
// Kinematic particle values
- cudaMemcpy( dev_x, k.x,
+ cudaMemcpy( dev_x, k.x,
memSizeF4, cudaMemcpyHostToDevice);
- cudaMemcpy( dev_xyzsum, k.xyzsum,
+ cudaMemcpy( dev_xyzsum, k.xyzsum,
memSizeF4, cudaMemcpyHostToDevice);
cudaMemcpy( dev_vel, k.vel,
memSizeF4, cudaMemcpyHostToDevice);
cudaMemcpy( dev_vel0, k.vel,
memSizeF4, cudaMemcpyHostToDevice);
- cudaMemcpy( dev_acc, k.acc,
+ cudaMemcpy( dev_acc, k.acc,
memSizeF4, cudaMemcpyHostToDevice);
cudaMemcpy( dev_force, k.force,
memSizeF4, cudaMemcpyHostToDevice);