tupdated 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 57d5ce5829d8daebbb03d7ce6e164b36fec3ff88
 (DIR) parent 099a12bc8d1dbbdebff7e24a174c82db6055f79b
 (HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
       Date:   Fri, 17 Oct 2014 10:09:29 +0200
       
       updated plots
       
       Diffstat:
         M python/shear-results-internals.py   |      11 ++++++++---
         M python/shear-results-strain.py      |      36 +++++++++++++++----------------
         M python/shear-results.py             |      19 +++++++++++++++++--
       
       3 files changed, 42 insertions(+), 24 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@@ -50,6 +50,7 @@ v_z_f = numpy.zeros((len(steps), sim.num[0], sim.num[1], sim.num[2]))
        
        # pressure - hydrostatic pressure
        dev_p = numpy.zeros((len(steps), sim.num[2]))
       +p     = numpy.zeros((len(steps), sim.num[2]))
        
        # mean per-particle values
        v_z_p_bar = numpy.zeros((len(steps), sim.num[2]))
       t@@ -85,7 +86,7 @@ for step_str in steps:
                    if step + substep > sim.status():
                        raise Exception(
                                'Simulation step %d not available (sim.status = %d).'
       -                        % (step, sim.status()))
       +                        % (step + substep, sim.status()))
        
                    sim.readstep(step + substep, verbose=False)
        
       t@@ -111,6 +112,9 @@ for step_str in steps:
                                + sim.p_f[0,0,-1])) \
                                /nsteps_avg
        
       +            p[s,:] += numpy.average(numpy.average(sim.p_f[:,:,:], axis=0),\
       +                    axis=0)/nsteps_avg
       +
                    v_z_f[s,:] += sim.v_f[:,:,:,2]/nsteps_avg
        
                    v_z_f_bar[s,:] += \
       t@@ -213,6 +217,7 @@ for s in numpy.arange(len(steps)):
        
            # 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
 (DIR) diff --git a/python/shear-results-strain.py b/python/shear-results-strain.py
       t@@ -14,10 +14,10 @@ import matplotlib.pyplot as plt
        from matplotlib.ticker import MaxNLocator
        
        sigma0 = 20000.0
       -cvals = ['dry', 1.0, 0.1]
       +cvals = ['dry', 1.0, 0.1, 0.01]
       +#cvals = ['dry', 1.0, 0.1]
        #cvals = ['dry', 1.0]
       -step = 800
       -nsteps_avg = 1 # no. of steps to average over
       +#step = 1999
        
        sim = sphere.sim('halfshear-sigma0=' + str(sigma0) + '-shear')
        sim.readfirst(verbose=False)
       t@@ -51,20 +51,13 @@ for c in cvals:
        
            if os.path.isfile('../output/' + sid + '.status.dat'):
        
       -        for substep in numpy.arange(nsteps_avg):
       +        sim.readlast(verbose=False)
        
       -            if step + substep > sim.status():
       -                raise Exception(
       -                        'Simulation step %d not available (sim.status = %d).'
       -                        % (step, sim.status()))
       +        zpos_p[s,:] = sim.x[:,2]
        
       -            sim.readstep(step + substep, verbose=False)
       +        xdisp[s,:] = sim.xyzsum[:,0]
        
       -            zpos_p[s,:] += sim.x[:,2]/nsteps_avg
       -
       -            xdisp[s,:] += sim.xyzsum[:,0]/nsteps_avg
       -
       -            #shear_strain[s] += sim.shearStrain()/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]):
       t@@ -74,6 +67,10 @@ for c in cvals:
                    if len(I) > 0:
                        xdisp_mean[s,iz] = numpy.mean(xdisp[s,I])
        
       +        # normalize distance
       +        max_dist = numpy.nanmax(xdisp_mean[s])
       +        xdisp_mean[s] /= max_dist
       +
            else:
                print(sid + ' not found')
            s += 1
       t@@ -85,8 +82,8 @@ fig = plt.figure(figsize=(8,6))
        
        ax = []
        #linetype = ['-', '--', '-.']
       -linetype = ['-', '-', '-']
       -color = ['b','g','r','c']
       +linetype = ['-', '-', '-', '-']
       +color = ['b','g','c','y']
        for s in numpy.arange(len(cvals)):
        
            ax.append(plt.subplot(111))
       t@@ -102,11 +99,12 @@ for s in numpy.arange(len(cvals)):
        
            #ax[0].plot(xdisp[s], zpos_p[s], ',', color = '#888888')
            #ax[0].plot(xdisp[s], zpos_p[s], ',', color=color[s], alpha=0.5)
       -    ax[0].plot(xdisp_mean[s], zpos_c[s], linetype[s], color=color[s],
       -            label=legend, linewidth=1)
       +    ax[0].plot(xdisp_mean[s], zpos_c[s], linetype[s],
       +            color=color[s], label=legend, linewidth=1)
        
            ax[0].set_ylabel('Vertical position $z$ [m]')
       -    ax[0].set_xlabel('$\\boldsymbol{x}^x_\\text{p}$ [m]')
       +    #ax[0].set_xlabel('$\\boldsymbol{x}^x_\\text{p}$ [m]')
       +    ax[0].set_xlabel('Normalized horizontal distance')
        
            #ax[s*4+0].get_xaxis().set_major_locator(MaxNLocator(nbins=5))
            #ax[s*4+1].get_xaxis().set_major_locator(MaxNLocator(nbins=5))
 (DIR) diff --git a/python/shear-results.py b/python/shear-results.py
       t@@ -16,6 +16,7 @@ import matplotlib.pyplot as plt
        smoothed_results = False
        contact_forces = False
        pressures = False
       +zflow = True
        
        #sigma0_list = numpy.array([1.0e3, 2.0e3, 4.0e3, 10.0e3, 20.0e3, 40.0e3])
        #sigma0 = 10.0e3
       t@@ -104,6 +105,7 @@ p_mean = [[], [], [], []]
        p_max = [[], [], [], []]
        f_n_mean = [[], [], [], []]
        f_n_max  = [[], [], [], []]
       +v_f_z_mean  = [[], [], [], []]
        
        fluid=True
        
       t@@ -163,6 +165,7 @@ 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@@ -176,6 +179,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:
       +                    v_f_z_mean[c][i] = numpy.mean(sim.v_f[:,:,:,2])
       +
            else:
                print(sid + ' not found')
        
       t@@ -195,8 +201,11 @@ fig.subplots_adjust(hspace=0.0)
        #plt.subplot(3,1,1)
        #plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
        
       -ax1 = plt.subplot(211)
       -ax2 = plt.subplot(212, sharex=ax1)
       +#ax1 = plt.subplot(211)
       +#ax2 = plt.subplot(212, sharex=ax1)
       +ax1 = plt.subplot(311)
       +ax2 = plt.subplot(312, sharex=ax1)
       +ax3 = plt.subplot(313, sharex=ax1)
        #ax3 = plt.subplot(413, sharex=ax1)
        #ax4 = plt.subplot(414, sharex=ax1)
        alpha = 0.5
       t@@ -222,6 +231,11 @@ for c in numpy.arange(1,len(cvals)+1):
            ax2.plot(shear_strain[c][1:], dilation[c][1:], \
                    label='$c$ = %.2f' % (cvals[c-1]), linewidth=1)
        
       +    if zflow:
       +        ax3.plot(shear_strain[c][1:], v_f_z_mean[c][1:],
       +            label='$c$ = %.2f' % (cvals[c-1]), linewidth=1)
       +
       +
            '''
            alpha = 0.5
            ax3.plot(shear_strain[c][1:], p_max[c][1:], '-' + color[c], alpha=alpha)
       t@@ -243,6 +257,7 @@ for c in numpy.arange(1,len(cvals)+1):
        
        ax1.set_ylabel('Shear friction $\\tau/\\sigma\'$ [-]')
        ax2.set_ylabel('Dilation $\\Delta h/(2r)$ [-]')
       +ax3.set_ylabel('$\\boldsymbol{v}_\\text{f}^z h$ [ms$^{-1}$]')
        #ax3.set_ylabel('Fluid pressure $p_\\text{f}$ [kPa]')
        #ax4.set_ylabel('Particle contact force $||\\boldsymbol{f}_\\text{p}||$ [N]')