Introduction
Introduction Statistics Contact Development Disclaimer Help
t1d-channel.py - granular-channel-hydro - subglacial hydrology model for sedime…
git clone git://src.adamsgaard.dk/granular-channel-hydro
Log
Files
Refs
README
LICENSE
---
t1d-channel.py (12863B)
---
1 #!/usr/bin/env python
2
3 # # ABOUT THIS FILE
4 # The following script uses basic Python and Numpy functionality to solv…
5 # coupled systems of equations describing subglacial channel development…
6 # soft beds as presented in `Damsgaard et al. "Sediment behavior controls
7 # equilibrium width of subglacial channels"`, accepted for publicaiton in
8 # Journal of Glaciology.
9 #
10 # High performance is not the goal for this implementation, which is ins…
11 # intended as a heavily annotated example on the solution procedure with…
12 # relying on solver libraries, suitable for low-level languages like C, …
13 # or CUDA.
14 #
15 # License: GNU Public License v3+ (see LICENSE.md)
16 # Author: Anders Damsgaard, [email protected], https://adamsgaard.dk
17
18 import numpy
19 import matplotlib.pyplot as plt
20 import sys
21
22
23 # # Model parameters
24 Ns = 25 # Number of nodes [-]
25 Ls = 1e3 # Model length [m]
26 total_days = 2. # Total simulation time [d]
27 t_end = 24.*60.*60.*total_days # Total simulation time [s]
28 tol_S = 1e-2 # Tolerance criteria for the norm. max. residual fo…
29 tol_Q = 1e-2 # Tolerance criteria for the norm. max. residual fo…
30 tol_P_c = 1e-2 # Tolerance criteria for the norm. max. residual fo…
31 max_iter = 1e2*Ns # Maximum number of solver iterations before failure
32 print_output_convergence = False # Display convergence in nested lo…
33 print_output_convergence_main = True # Display convergence in main loop
34 safety = 0.01 # Safety factor ]0;1] for adaptive timestepping
35 plot_interval = 20 # Time steps between plots
36 plot_during_iterations = False # Generate plots for intermediate …
37
38 # Physical parameters
39 rho_w = 1000. # Water density [kg/m^3]
40 rho_i = 910. # Ice density [kg/m^3]
41 rho_s = 2600. # Sediment density [kg/m^3]
42 g = 9.8 # Gravitational acceleration [m/s^2]
43 theta = 30. # Angle of internal friction in sediment [deg]
44 D = 1.15e-3 # Mean grain size [m], Lajeuness et al 2010, series…
45 tau_c = 0.016 # Critical Shields stress, Lajeunesse et al 2010, s…
46
47 # Boundary conditions
48 P_terminus = 0. # Water pressure at terminus [Pa]
49 m_dot = numpy.linspace(1e-6, 1e-5, Ns-1) # Water source term [m/s]
50 Q_upstream = 1e-3 # Water influx upstream (must be larger than 0) [m^…
51
52 # Channel hydraulic properties
53 manning = 0.1 # Manning roughness coefficient [m^{-1/3} s]
54 friction_factor = 0.1 # Darcy-Weisbach friction factor [-]
55
56 # Channel growth-limit parameters
57 c_1 = -0.118 # [m/kPa]
58 c_2 = 4.60 # [m]
59
60 # Minimum channel size [m^2], must be bigger than 0
61 S_min = 1e-2
62
63
64 # # Initialize model arrays
65 # Node positions, terminus at Ls
66 s = numpy.linspace(0., Ls, Ns)
67 ds = s[1:] - s[:-1]
68
69 # Ice thickness [m]
70 H = 6.*(numpy.sqrt(Ls - s + 5e3) - numpy.sqrt(5e3)) + 1.0
71
72 # Bed topography [m]
73 b = numpy.zeros_like(H)
74
75 N = H*0.1*rho_i*g # Initial effective stress [Pa]
76
77 # Initialize arrays for channel segments between nodes
78 S = numpy.ones(len(s) - 1)*0.1 # Cross-sect. area of channel segments …
79 S_max = numpy.zeros_like(S) # Max. channel size [m^2]
80 dSdt = numpy.zeros_like(S) # Transient in channel cross-sect. area [m^…
81 W = S/numpy.tan(numpy.deg2rad(theta)) # Assuming no channel floor wedge
82 Q = numpy.zeros_like(S) # Water flux in channel segments [m^3/s]
83 Q_s = numpy.zeros_like(S) # Sediment flux in channel segments [m^3/s]
84 P_c = numpy.zeros_like(S) # Effective pressure in channel segments [P…
85 P_w = numpy.zeros_like(S) # Effective pressure in channel segments [P…
86 tau = numpy.zeros_like(S) # Avg. shear stress from current [Pa]
87 porosity = numpy.ones_like(S)*0.3 # Sediment porosity [-]
88 res = numpy.zeros_like(S) # Solution residual during solver iterations
89
90
91 # # Helper functions
92 def gradient(arr, arr_x):
93 # Central difference gradient of an array ``arr`` with node position…
94 # ``arr_x``.
95 return (arr[1:] - arr[:-1])/(arr_x[1:] - arr_x[:-1])
96
97
98 def avg_midpoint(arr):
99 # Averaged value of neighboring array elements
100 return (arr[:-1] + arr[1:])/2.
101
102
103 def channel_hydraulic_roughness(manning, S, W, theta):
104 # Determine hydraulic roughness assuming that the channel has no cha…
105 # floor wedge.
106 l = W + W/numpy.cos(numpy.deg2rad(theta)) # wetted perimeter
107 return manning**2.*(l**2./S)**(2./3.)
108
109
110 def channel_shear_stress(Q, S):
111 # Determine mean wall shear stress from Darcy-Weisbach friction loss
112 u_bar = Q/S
113 return 1./8.*friction_factor*rho_w*u_bar**2.
114
115
116 def channel_sediment_flux(tau, W):
117 # Meyer-Peter and Mueller 1948
118 # tau: Shear stress by water flow
119 # W: Channel width
120
121 # Non-dimensionalize shear stress
122 shields_stress = tau/((rho_s - rho_w)*g*D)
123
124 stress_excess = shields_stress - tau_c
125 stress_excess[stress_excess < 0.] = 0.
126 return 8.*stress_excess**(3./2.)*W \
127 * numpy.sqrt((rho_s - rho_w)/rho_w*g*D**3.)
128
129
130 def channel_growth_rate_sedflux(Q_s, porosity, s_c):
131 # Damsgaard et al, in prep
132 return 1./porosity[1:] * gradient(Q_s, s_c)
133
134
135 def update_channel_size_with_limit(S, S_old, dSdt, dt, P_c):
136 # Damsgaard et al, in prep
137 S_max = numpy.maximum(
138 numpy.maximum(
139 1./4.*(c_1*numpy.maximum(P_c, 0.)/1000. + c_2), 0.)**2. *
140 numpy.tan(numpy.deg2rad(theta)), S_min)
141 S = numpy.maximum(numpy.minimum(S_old + dSdt*dt, S_max), S_min)
142 W = S/numpy.tan(numpy.deg2rad(theta)) # Assume no channel floor wed…
143 dSdt = S - S_old # adjust dSdt for clipping due to channel size li…
144 return S, W, S_max, dSdt
145
146
147 def flux_solver(m_dot, ds):
148 # Iteratively find new water fluxes
149 it = 0
150 max_res = 1e9 # arbitrary large value
151
152 # Iteratively find solution, do not settle for less iterations than …
153 # number of nodes
154 while max_res > tol_Q:
155
156 Q_old = Q.copy()
157 # dQ/ds = m_dot -> Q_out = m*delta(s) + Q_in
158 # Upwind information propagation (upwind)
159 Q[0] = Q_upstream
160 Q[1:] = m_dot[1:]*ds[1:] + Q[:-1]
161 max_res = numpy.max(numpy.abs((Q - Q_old)/(Q + 1e-16)))
162
163 if print_output_convergence:
164 print('it = {}: max_res = {}'.format(it, max_res))
165
166 # import ipdb; ipdb.set_trace()
167 if it >= max_iter:
168 raise Exception('t = {}, step = {}: '.format(time, step) +
169 'Iterative solution not found for Q')
170 it += 1
171
172 return Q
173
174
175 def pressure_solver(psi, f, Q, S):
176 # Iteratively find new water pressures
177 # dP_c/ds = f*rho_w*g*Q^2/S^{8/3} - psi (Kingslake and Ng 2013)
178
179 it = 0
180 max_res = 1e9 # arbitrary large value
181 while max_res > tol_P_c:
182
183 P_c_old = P_c.copy()
184
185 # Dirichlet BC (fixed pressure) at terminus
186 P_c[-1] = rho_i*g*H_c[-1] - P_terminus
187
188 P_c[:-1] = P_c[1:] \
189 + psi[:-1]*ds[:-1] \
190 - f[:-1]*rho_w*g*Q[:-1]*numpy.abs(Q[:-1]) \
191 / (S[:-1]**(8./3.))*ds[:-1]
192
193 max_res = numpy.max(numpy.abs((P_c - P_c_old)/(P_c + 1e-16)))
194
195 if print_output_convergence:
196 print('it = {}: max_res = {}'.format(it, max_res))
197
198 if it >= max_iter:
199 raise Exception('t = {}, step = {}:'.format(time, step) +
200 'Iterative solution not found for P_c')
201 it += 1
202
203 return P_c
204
205
206 def plot_state(step, time, S_, S_max_, title=False):
207 # Plot parameters along profile
208 fig = plt.gcf()
209 fig.set_size_inches(3.3*1.1, 3.3*1.1*1.5)
210
211 ax_Pa = plt.subplot(3, 1, 1) # axis with Pascals as y-axis unit
212 ax_Pa.plot(s_c/1000., P_c/1e6, '-k', label='$P_c$')
213 ax_Pa.plot(s_c/1000., H_c*rho_i*g/1e6, '--r', label='$P_i$')
214 ax_Pa.plot(s_c/1000., P_w/1e6, ':y', label='$P_w$')
215
216 ax_m3s = ax_Pa.twinx() # axis with m3/s as y-axis unit
217 ax_m3s.plot(s_c/1000., Q, '.-b', label='$Q$')
218
219 if title:
220 plt.title('Day: {:.3}'.format(time/(60.*60.*24.)))
221 ax_Pa.legend(loc=2)
222 ax_m3s.legend(loc=1)
223 ax_Pa.set_ylabel('[MPa]')
224 ax_m3s.set_ylabel('[m$^3$/s]')
225
226 # ax_m3s_sed = plt.subplot(3, 1, 2, sharex=ax_Pa)
227 ax_m3s_sed_blank = plt.subplot(3, 1, 2, sharex=ax_Pa)
228 ax_m3s_sed_blank.get_yaxis().set_visible(False)
229 ax_m3s_sed = ax_m3s_sed_blank.twinx()
230 #ax_m3s_sed.plot(s_c/1000., Q_s, '-', label='$Q_{s}$')
231 ax_m3s_sed.plot(s_c/1000., Q_s*1000., '-', label='$Q_{s}$')
232 ax_m3s_sed.legend(loc=2)
233
234 ax_m2 = plt.subplot(3, 1, 3, sharex=ax_Pa)
235 ax_m2.plot(s_c/1000., S_, '-k', label='$S$')
236 ax_m2.plot(s_c/1000., S_max_, '--', color='#666666', label='$S_{max}…
237
238 ax_m2s = ax_m2.twinx()
239 #ax_m2s.plot(s_c/1000., dSdt, ':', label='$dS/dt$')
240 ax_m2s.plot(s_c/1000., dSdt*1000., ':', label='$dS/dt$')
241
242 ax_m2.legend(loc=2)
243 ax_m2s.legend(loc=3)
244 ax_m2.set_xlabel('$s$ [km]')
245 ax_m2.set_ylabel('[m$^2$]')
246 #ax_m3s_sed.set_ylabel('[m$^3$/s]')
247 #ax_m2s.set_ylabel('[m$^2$/s]')
248 ax_m3s_sed.set_ylabel('[mm$^3$/s]')
249 ax_m2s.set_ylabel('[mm$^2$/s]')
250
251 # use scientific notation for m3s and m2s axes
252 #ax_m3s_sed.get_yaxis().set_major_formatter(plt.LogFormatter(10,
253 #labelOnlyBase=False))
254 #ax_m2s.get_yaxis().set_major_formatter(plt.LogFormatter(10,
255 #labelOnlyBase=False))
256
257 ax_Pa.set_xlim([s.min()/1000., s.max()/1000.])
258
259 plt.setp(ax_Pa.get_xticklabels(), visible=False)
260 plt.tight_layout()
261 if step == -1:
262 plt.savefig('chan-0.init.pdf')
263 else:
264 plt.savefig('chan-' + str(step) + '.pdf')
265 plt.clf()
266 plt.close()
267
268
269 def find_new_timestep(ds, Q, Q_s, S):
270 # Determine the timestep using the Courant-Friedrichs-Lewy condition
271 dt = safety*numpy.minimum(60.*60.*24.,
272 numpy.min(numpy.abs(ds/(Q*S),
273 ds/(Q_s*S)+1e-16)))
274
275 if dt < 1.0:
276 raise Exception('Error: Time step less than 1 second at step '
277 + '{}, time '.format(step)
278 + '{:.3} s/{:.3} d'.format(time, time/(60.*60.*2…
279
280 return dt
281
282
283 def print_status_to_stdout(step, time, dt):
284 sys.stdout.write('\rstep = {}, '.format(step) +
285 't = {:.2} s or {:.4} d, dt = {:.2} s '
286 .format(time, time/(60.*60.*24.), dt))
287 sys.stdout.flush()
288
289
290 # Initialize remaining values before entering time loop
291 s_c = avg_midpoint(s) # Channel section midpoint coordinates [m]
292 H_c = avg_midpoint(H)
293
294 # Water-pressure gradient from geometry [Pa/m]
295 psi = -rho_i*g*gradient(H, s) - (rho_w - rho_i)*g*gradient(b, s)
296
297 # Prepare figure object for plotting during the simulation
298 fig = plt.figure('channel')
299 plot_state(-1, 0.0, S, S_max)
300
301
302 # # Time loop
303 time = 0.
304 step = 0
305 while time <= t_end:
306
307 # Determine time step length from water flux
308 dt = find_new_timestep(ds, Q, Q_s, S)
309
310 # Display current simulation status
311 print_status_to_stdout(step, time, dt)
312
313 it = 0
314
315 # Initialize the maximum normalized residual for S to an arbitrary l…
316 # value
317 max_res = 1e9
318
319 S_old = S.copy()
320 # Iteratively find solution with the Jacobi relaxation method
321 while max_res > tol_S:
322
323 S_prev_it = S.copy()
324
325 # Find new water fluxes consistent with mass conservation and lo…
326 # meltwater production (m_dot)
327 Q = flux_solver(m_dot, ds)
328
329 # Find average shear stress from water flux for each channel seg…
330 tau = channel_shear_stress(Q, S)
331
332 # Determine sediment flux
333 Q_s = channel_sediment_flux(tau, W)
334
335 # Determine change in channel size for each channel segment.
336 # Use backward differences and assume dS/dt=0 in first segment.
337 dSdt[1:] = channel_growth_rate_sedflux(Q_s, porosity, s_c)
338
339 # Update channel cross-sectional area and width according to gro…
340 # rate and size limit for each channel segment
341 S, W, S_max, dSdt = \
342 update_channel_size_with_limit(S, S_old, dSdt, dt, P_c)
343
344 f = channel_hydraulic_roughness(manning, S, W, theta)
345
346 # Find new effective pressures consistent with the flow law and …
347 # pressures in channel segments
348 P_c = pressure_solver(psi, f, Q, S)
349 P_w = rho_i*g*H_c - P_c
350
351 if plot_during_iterations:
352 plot_state(step + it/1e4, time, S, S_max)
353
354 # Find new maximum normalized residual value
355 max_res = numpy.max(numpy.abs((S - S_prev_it)/(S + 1e-16)))
356 if print_output_convergence_main:
357 print('it = {}: max_res = {}'.format(it, max_res))
358
359 # import ipdb; ipdb.set_trace()
360 if it >= max_iter:
361 raise Exception('t = {}, step = {}: '.format(time, step) +
362 'Iterative solution not found')
363 it += 1
364
365 # Generate an output figure for every n time steps
366 if step % plot_interval == 0:
367 plot_state(step, time, S, S_max)
368
369 # Update time
370 time += dt
371 step += 1
You are viewing proxied material from mx1.adamsgaard.dk. The copyright of proxied material belongs to its original authors. Any comments or complaints in relation to proxied material should be directed to the original authors of the content concerned. Please see the disclaimer for more details.