tsupraglacial-plots.py - 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
---
tsupraglacial-plots.py (5145B)
---
1 #!/usr/bin/env python
2
3 import sphere
4 import numpy as np
5 import matplotlib
6 matplotlib.rcParams.update({'font.size': 12})
7 import matplotlib.pyplot as plt
8 import subprocess
9
10 dpdz_values = [0.0, -50.0, -100.0]
11 slope_angle_values = [5.0, 10.0, 15.0, 20.0]
12
13 scatter_color = '#666666'
14 scatter_alpha = 0.6
15 dp_dz_vals = np.empty(len(dpdz_values)*len(slope_angle_values))
16 slope_angles = np.empty_like(dp_dz_vals)
17 fluxes = np.empty_like(dp_dz_vals)
18
19 '''
20 j = 0
21 for dpdz in dpdz_values:
22 outfiles = ''
23 for slope_angle in slope_angle_values:
24
25 sim = sphere.sim("supraglacial-slope{}-dpdz{}".format(slope_angle, dpdz),
26 fluid=True)
27 print('### ' + sim.id())
28 print('Last output file: ' + str(sim.status()))
29 sim.readTime(4.99)
30
31 fig = plt.figure()
32
33 title = 'slope = ' + str(slope_angle) + '$^\circ$, ' + \
34 '$dp/dz$ = -' + str(dpdz) + ' Pa/m'
35
36 z = np.zeros(20)
37 v_x_space_avg = np.empty_like(z)
38 xsum_space_avg = np.empty_like(z)
39 dz = np.max(sim.x[:,2])/len(z)
40 for i in range(len(z)):
41 z[i] = i*dz + 0.5*dz
42 I = np.nonzero((sim.x[:,2] >= i*dz) & (sim.x[:,2] < (i+1)*dz))
43 v_x_space_avg[i] = np.mean(sim.vel[I,0])
44 xsum_space_avg[i] = np.mean(sim.xyzsum[I,0])
45
46 plt.plot(sim.vel[:,0], sim.x[:,2], '.',
47 color=scatter_color, alpha=scatter_alpha)
48 plt.plot(v_x_space_avg, z, '+-k')
49 plt.title(title)
50 plt.xlabel('Horizontal particle velocity $v_x$ [m/s]')
51 plt.ylabel('Vertical position $z$ [m]')
52 plt.savefig(sim.id() + '-vel.pdf')
53 plt.savefig(sim.id() + '-vel.png')
54 plt.clf()
55
56 plt.plot(sim.xyzsum[:,0]/sim.time_current, sim.x[:,2], '.',
57 color=scatter_color, alpha=scatter_alpha)
58 plt.plot(xsum_space_avg/sim.time_current, z, '+-k')
59 plt.title(title)
60 plt.xlabel('Average horizontal particle velocity $\\bar{v}_x$ [m/s]')
61 plt.ylabel('Vertical position $z$ [m]')
62 plt.savefig(sim.id() + '-avg_vel.pdf')
63 plt.savefig(sim.id() + '-avg_vel.png')
64 plt.clf()
65
66 plt.close()
67 outfiles += sim.id() + '-avg_vel.png '
68
69 dp_dz_vals[j] = dpdz
70 slope_angles[j] = slope_angle
71 fluxes[j] = np.trapz(xsum_space_avg/sim.time_current, dx=dz)
72 j += 1
73
74 subprocess.call('montage ' + outfiles +
75 '-geometry +0+0 ' +
76 'supraglacial-avg_vel-dpdz_' + str(dpdz) + '.png ',
77 shell=True)
78
79 for dpdz in dpdz_values:
80 I = np.nonzero(dp_dz_vals == dpdz)
81 plt.semilogy(slope_angles[I[0]], fluxes[I[0]], '+-', label='$dp/dz$ = ' + str(dpdz) + ' Pa/m')
82
83 plt.legend()
84 plt.xlabel('Slope [$^\circ$]')
85 plt.ylabel('Specific sediment flux [m$^2$/s]')
86 plt.savefig('supraglacial_flux.png')
87 plt.savefig('supraglacial_flux.pdf')
88 '''
89
90 # time series
91 for dpdz in dpdz_values:
92 for slope_angle in slope_angle_values:
93 sim = sphere.sim("supraglacial-slope{}-dpdz{}".format(slope_angle, dpdz), fluid=True)
94 print('### ' + sim.id())
95 sim.readlast()
96
97 N_time = min(100, sim.status())
98 timesteps = np.linspace(sim.time_file_dt[0], sim.time_current[0] - sim.time_file_dt[0], N_time)
99 porosity = np.empty(N_time)
100 velocity = np.empty(N_time)
101 displacement = np.empty(N_time)
102 flux = np.empty(N_time)
103 z = np.zeros(20)
104 v_x_space_avg = np.empty_like(z)
105 xsum_space_avg = np.empty_like(z)
106
107 for it in np.arange(N_time):
108 sim.readTime(timesteps[it])
109
110 dz = np.max(sim.x[:,2])/len(z)
111 for i in range(len(z)):
112 z[i] = i*dz + 0.5*dz
113 I = np.nonzero((sim.x[:,2] >= i*dz) & (sim.x[:,2] < (i+1)*dz))
114 xsum_space_avg[i] = np.mean(sim.xyzsum[I,0])
115
116 porosity[it] = np.mean(sim.phi)
117 velocity[it] = np.mean(sim.vel[:,0])
118 displacement[it] = np.mean(sim.xyzsum[:,0])
119 flux[it] = np.trapz(xsum_space_avg/sim.time_current, dx=dz)
120
121 fig = plt.figure(figsize=(6,8))
122 fig.suptitle('slope = ' + str(slope_angle) + '$^\circ$, ' + \
123 '$dp/dz$ = -' + str(dpdz) + ' Pa/m',
124 horizontalalignment='left')
125
126 ax1 = plt.subplot(3,1,1)
127 plt.plot(timesteps, porosity, '-')
128 plt.ylabel('Porosity [-]')
129 plt.setp(ax1.get_xticklabels(), visible=False)
130
131 #ax2 = plt.subplot(3,1,2)
132 #plt.semilogy(timesteps, velocity, '-')
133 #plt.ylabel('Avg. velocity [m/s]')
134 #plt.setp(ax2.get_xticklabels(), visible=False)
135
136 ax3 = plt.subplot(3,1,2)
137 plt.semilogy(timesteps, displacement, '-')
138 plt.ylabel('Avg. displacement [m]')
139 plt.setp(ax3.get_xticklabels(), visible=False)
140
141 ax1 = plt.subplot(3,1,3)
142 plt.semilogy(timesteps, flux, '-')
143 plt.ylabel('Cumulative flux [m$^2$/s]')
144 plt.xlabel('Time [s]')
145
146 plt.tight_layout()
147 plt.savefig(sim.id() + '-timeseries.png')
148 plt.savefig(sim.id() + '-timeseries.pdf')
149 plt.clf()
150 plt.close()