tCalculate flux for supraglacial creep experiments - 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 497a00115fb29356886ae81684f62e6f80d021eb
 (DIR) parent 4250795968a1883585822f0eae9d1c2060cf91f9
 (HTM) Author: Anders Damsgaard <anders@adamsgaard.dk>
       Date:   Thu, 21 Nov 2019 11:19:01 +0100
       
       Calculate flux for supraglacial creep experiments
       
       Diffstat:
         M python/supraglacial-plots.py        |      35 ++++++++++++++++++++++++-------
       
       1 file changed, 27 insertions(+), 8 deletions(-)
       ---
 (DIR) diff --git a/python/supraglacial-plots.py b/python/supraglacial-plots.py
       t@@ -1,7 +1,7 @@
        #!/usr/bin/env python
        
        import sphere
       -import numpy
       +import numpy as np
        import matplotlib
        matplotlib.rcParams.update({'font.size': 12})
        import matplotlib.pyplot as plt
       t@@ -12,7 +12,11 @@ slope_angle_values = [5.0, 10.0, 15.0, 20.0]
        
        scatter_color = '#666666'
        scatter_alpha = 0.6
       +dp_dz_vals = np.empty(len(dpdz_values)*len(slope_angle_values))
       +slope_angles = np.empty_like(dp_dz_vals)
       +fluxes = np.empty_like(dp_dz_vals)
        
       +j = 0
        for dpdz in dpdz_values:
            outfiles = ''
            for slope_angle in slope_angle_values:
       t@@ -28,15 +32,15 @@ for dpdz in dpdz_values:
                title = 'slope = ' + str(slope_angle) + '$^\circ$, ' + \
                        '$dp/dz$ = -' + str(dpdz) + ' Pa/m'
        
       -        z = numpy.zeros(20)
       -        v_x_space_avg = numpy.empty_like(z)
       -        xsum_space_avg = numpy.empty_like(z)
       -        dz = numpy.max(sim.x[:,2])/len(z)
       +        z = np.zeros(20)
       +        v_x_space_avg = np.empty_like(z)
       +        xsum_space_avg = np.empty_like(z)
       +        dz = np.max(sim.x[:,2])/len(z)
                for i in range(len(z)):
                    z[i] = i*dz + 0.5*dz
       -            I = numpy.nonzero((sim.x[:,2] >= i*dz) & (sim.x[:,2] < (i+1)*dz))
       -            v_x_space_avg[i] = numpy.mean(sim.vel[I,0])
       -            xsum_space_avg[i] = numpy.mean(sim.xyzsum[I,0])
       +            I = np.nonzero((sim.x[:,2] >= i*dz) & (sim.x[:,2] < (i+1)*dz))
       +            v_x_space_avg[i] = np.mean(sim.vel[I,0])
       +            xsum_space_avg[i] = np.mean(sim.xyzsum[I,0])
        
                plt.plot(sim.vel[:,0], sim.x[:,2], '.',
                         color=scatter_color, alpha=scatter_alpha)
       t@@ -61,7 +65,22 @@ for dpdz in dpdz_values:
                plt.close()
                outfiles += sim.id() + '-avg_vel.png '
        
       +        dp_dz_vals[j] = dpdz
       +        slope_angles[j] = slope_angle
       +        fluxes[j] = np.trapz(xsum_space_avg/sim.time_current, dx=dz)
       +        j += 1
       +
            subprocess.call('montage ' + outfiles +
                            '-geometry +0+0 ' +
                            'supraglacial-avg_vel-dpdz_' + str(dpdz) + '.png ',
                             shell=True)
       +
       +for dpdz in dpdz_values:
       +    I = np.nonzero(dp_dz_vals == dpdz)
       +    plt.semilogy(slope_angles[I[0]], fluxes[I[0]], '+-', label='$dp/dz$ = ' + str(dpdz) + ' Pa/m')
       +
       +plt.legend()
       +plt.xlabel('Slope [$^\circ$]')
       +plt.ylabel('Specific sediment flux [m/s]')
       +plt.savefig('supraglacial_flux.png')
       +plt.savefig('supraglacial_flux.pdf')