tadd new plot - sphere - GPU-based 3D discrete element method algorithm with op… | |
git clone git://src.adamsgaard.dk/sphere | |
Log | |
Files | |
Refs | |
LICENSE | |
--- | |
commit f747e43ef3abf99c7365e9eb6451f7bbd8828304 | |
parent 3e077e59134905ebc4b377df485aeb0f422682a9 | |
Author: Anders Damsgaard <[email protected]> | |
Date: Mon, 6 Oct 2014 14:53:53 +0200 | |
add new plot | |
Diffstat: | |
M python/shear-results-pressures.py | 4 ++-- | |
A python/shear-results-strain.py | 148 +++++++++++++++++++++++++++++… | |
2 files changed, 150 insertions(+), 2 deletions(-) | |
--- | |
diff --git a/python/shear-results-pressures.py b/python/shear-results-pressures… | |
t@@ -47,9 +47,8 @@ for i in numpy.arange(sim.status()): | |
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[… | |
pres_static[z,i] = \ | |
- (wall0_iz*dz - zpos_c[z])*sim.rho_f*numpy.abs(sim.g[2])\ | |
+ (wall0_iz*dz - zpos_c[z] + 0.5*dz)*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@@ -77,6 +76,7 @@ ax1 = plt.subplot(311) | |
# cmap=cmap, norm=norm) | |
im1 = ax1.pcolormesh(shear_strain, zpos_c, dev_pres/1000.0, vmin=min_p, | |
vmax=max_p, rasterized=True) | |
+#im1 = ax1.pcolormesh(shear_strain, zpos_c, dev_pres/1000.0, rasterized=True) | |
#ax1.set_xlim([0, shear_strain[-1]]) | |
#ax1.set_ylim([zpos_c[0], sim.w_x[0]]) | |
ax1.set_xlabel('Shear strain $\\gamma$ [-]') | |
diff --git a/python/shear-results-strain.py b/python/shear-results-strain.py | |
t@@ -0,0 +1,148 @@ | |
+#!/usr/bin/env python | |
+import matplotlib | |
+matplotlib.use('Agg') | |
+matplotlib.rcParams.update({'font.size': 18, 'font.family': 'serif'}) | |
+matplotlib.rc('text', usetex=True) | |
+matplotlib.rcParams['text.latex.preamble']=[r"\usepackage{amsmath}"] | |
+import shutil | |
+ | |
+import os | |
+import numpy | |
+import sphere | |
+from permeabilitycalculator import * | |
+import matplotlib.pyplot as plt | |
+from matplotlib.ticker import MaxNLocator | |
+ | |
+#cvals = ['dry', 1.0, 0.1] | |
+sigma0 = 20000.0 | |
+cvals = ['dry', 1.0] | |
+step = 1000 | |
+nsteps_avg = 1 # no. of steps to average over | |
+ | |
+sim = sphere.sim('halfshear-sigma0=' + str(sigma0) + '-shear') | |
+sim.readfirst(verbose=False) | |
+ | |
+ | |
+# particle z positions | |
+zpos_p = numpy.zeros((len(cvals), sim.np)) | |
+ | |
+# cell midpoint cell positions | |
+zpos_c = numpy.zeros((len(cvals), 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(cvals), sim.np)) | |
+ | |
+xdisp_mean = numpy.zeros((len(cvals), sim.num[2])) | |
+ | |
+s = 0 | |
+for c in cvals: | |
+ | |
+ if c == 'dry': | |
+ fluid = False | |
+ sid = 'halfshear-sigma0=' + str(sigma0) + '-shear' | |
+ else: | |
+ fluid = True | |
+ sid = 'halfshear-sigma0=' + str(sigma0) + '-c=' + str(c) + '-shear' | |
+ | |
+ sim = sphere.sim(sid, fluid=fluid) | |
+ | |
+ if os.path.isfile('../output/' + sid + '.status.dat'): | |
+ | |
+ for substep in numpy.arange(nsteps_avg): | |
+ | |
+ if step + substep > sim.status(): | |
+ raise Exception( | |
+ 'Simulation step %d not available (sim.status = %d).' | |
+ % (step, sim.status())) | |
+ | |
+ sim.readstep(step + substep, verbose=False) | |
+ | |
+ zpos_p[s,:] += sim.x[:,2]/nsteps_avg | |
+ | |
+ xdisp[s,:] += sim.xyzsum[:,0]/nsteps_avg | |
+ | |
+ #shear_strain[s] += sim.shearStrain()/nsteps_avg | |
+ | |
+ # 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]) | |
+ | |
+ 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,6)) | |
+ | |
+ax = [] | |
+linetype = ['-', '--', '-.'] | |
+for s in numpy.arange(len(cvals)): | |
+ | |
+ 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' | |
+ else: | |
+ legend = 'c = ' + str(cvals[s]) | |
+ | |
+ ax[0].plot(xdisp[s], zpos_p[s], ',', color = '#888888') | |
+ ax[0].plot(xdisp_mean[s], zpos_c[s], linetype[s], color='k', label = legen… | |
+ linewidth=2) | |
+ | |
+ ax[0].set_ylabel('Vertical position $z$ [m]') | |
+ ax[0].set_xlabel('$\\boldsymbol{x}^x_\\text{p}$ [m]') | |
+ | |
+ #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 fractio… | |
+ #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}) | |
+ | |
+legend_alpha=0.5 | |
+ax[0].legend(loc='best', prop={'size':18}, fancybox=True, framealpha=legend_al… | |
+ax[0].grid() | |
+plt.tight_layout() | |
+plt.subplots_adjust(wspace = .05) | |
+plt.MaxNLocator(nbins=4) | |
+ | |
+filename = 'shear-' + str(int(sigma0/1000.0)) + 'kPa-strain.pdf' | |
+plt.savefig(filename) | |
+shutil.copyfile(filename, '/home/adc/articles/own/2-org/' + filename) | |
+print(filename) | |
+ | |
+ |