GopherProxy

	tadd script to plot shear zone internals - sphere - GPU-based 3D discrete eleme…
	git clone git://src.adamsgaard.dk/sphere
	Log
	Files
	Refs
	LICENSE
	---
	commit c3f43ac112b3dbefc930e2c82e437b6def40446b
	parent e4530ec4d724ea43f1e7ec7695c426c0021247db
	Author: Anders Damsgaard <[email protected]>
	Date: Tue, 9 Sep 2014 14:42:21 +0200

	add script to plot shear zone internals

	Diffstat:
	A python/shear-results-forces.py \| 107 +++++++++++++++++++++++++++++…
	M python/shear-results.py \| 2 +-

	2 files changed, 108 insertions(+), 1 deletion(-)
	---
	diff --git a/python/shear-results-forces.py b/python/shear-results-forces.py
	t@@ -0,0 +1,107 @@
	+#!/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
	+
	+#steps = [5, 10, 100]
	+steps = [3]
	+nsteps_avg = 1 # no. of steps to average over
	+
	+sigma0 = 10.0e3
	+c_grad_p = 1.0
	+c_phi = 1.0
	+
	+sid = 'shear-sigma0=' + str(sigma0) + '-c_phi=' + \
	+ str(c_phi) + '-c_grad_p=' + str(c_grad_p) + '-hi_mu-lo_visc'
	+sim = sphere.sim(sid, fluid=True)
	+sim.readfirst(verbose=False)
	+
	+# 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] = idz + 0.5dz
	+
	+# particle x displacements
	+xdisp = numpy.zeros((len(steps), sim.np))
	+
	+# particle-fluid force per particle
	+f_pf = numpy.zeros_like(xdisp)
	+
	+# pressure - hydrostatic pressure
	+dev_p = numpy.zeros((len(steps), sim.num[2]))
	+
	+s = 0
	+for step in steps:
	+
	+ if os.path.isfile('../output/' + sid + '.status.dat'):
	+
	+ if step > sim.status():
	+ raise Exception(
	+ 'Simulation step %d not available (sim.status = %d).'
	+ % (step, sim.status()))
	+
	+ for substep in numpy.arange(nsteps_avg):
	+
	+ sim.readstep(step + substep, verbose=False)
	+
	+ zpos_p[s,:] += sim.x[:,2]/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
	+
	+ dev_p[s,:] += \
	+ numpy.average(numpy.average(sim.p_f, axis=0), axis=0)\
	+ /nsteps_avg
	+
	+ else:
	+ print(sid + ' not found')
	+ s += 1
	+
	+fig = plt.figure(figsize=(8,4*(len(steps))))
	+
	+for s in numpy.arange(len(steps)):
	+ ax1 = plt.subplot((s+1)*100 + 31)
	+ ax2 = plt.subplot((s+1)*100 + 32, sharey=ax1)
	+ ax3 = plt.subplot((s+1)*100 + 33, sharey=ax1)
	+
	+ ax1.plot(xdisp[s], zpos_p[s], '.')
	+ ax2.plot(f_pf[s], zpos_p[s], '.')
	+ ax3.plot(dev_p[s], zpos_c[s], '.')
	+ #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]))
	+ #plt.loglog(dpdz[c], K[c], 'o-', label='$c$ = %.2f' % (cvals[c]))
	+
	+ ax1.set_ylabel('Vertical position $z$ [m]')
	+ ax1.set_xlabel('Horizontal particle displacement [m]')
	+ ax2.set_xlabel('Average fluid pressure $\\bar{p_f}$ [Pa]')
	+ ax3.set_xlabel('Average fluid-particle interaction force '
	+ + '$\|\|\\bar{\\boldsymbol{f_{pf}}}\|\|$ [N]')
	+ plt.setp(ax2.get_yticklabels(), visible=False)
	+ plt.setp(ax3.get_yticklabels(), visible=False)
	+ #ax1.grid()
	+ #ax2.grid()
	+ #ax1.legend(loc='lower right', prop={'size':18})
	+ #ax2.legend(loc='lower right', prop={'size':18})
	+
	+plt.tight_layout()
	+filename = 'shear-10kPa-forces.pdf'
	+plt.savefig(filename)
	+shutil.copyfile(filename, '/home/adc/articles/own/2-org/' + filename)
	+print(filename)
	diff --git a/python/shear-results.py b/python/shear-results.py
	t@@ -104,7 +104,7 @@ ax3.set_xlabel('Shear strain $\\gamma$ [-]')

	ax1.set_ylabel('Shear friction $\\tau/\\sigma\'$ [-]')
	ax2.set_ylabel('Dilation $\\Delta h/(2r)$ [-]')
	-ax3.set_ylabel('Fluid pressure $p_\\text{f}$ [Pa]')
	+ax3.set_ylabel('Fluid pressure $p_f$ [Pa]')

	plt.setp(ax1.get_xticklabels(), visible=False)
	plt.setp(ax2.get_xticklabels(), visible=False)