t1d-channel.py - granular-channel-hydro - subglacial hydrology model for sedimentary channels
(HTM) git clone git://src.adamsgaard.dk/granular-channel-hydro
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) 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 solve the
5 # coupled systems of equations describing subglacial channel development in
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 instead
11 # intended as a heavily annotated example on the solution procedure without
12 # relying on solver libraries, suitable for low-level languages like C, Fortran
13 # or CUDA.
14 #
15 # License: GNU Public License v3+ (see LICENSE.md)
16 # Author: Anders Damsgaard, andersd@princeton.edu, 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 for S
29 tol_Q = 1e-2 # Tolerance criteria for the norm. max. residual for Q
30 tol_P_c = 1e-2 # Tolerance criteria for the norm. max. residual for P_c
31 max_iter = 1e2*Ns # Maximum number of solver iterations before failure
32 print_output_convergence = False # Display convergence in nested loops
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 results
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 1
45 tau_c = 0.016 # Critical Shields stress, Lajeunesse et al 2010, series 1
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^3/s]
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 [m^2]
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^2/s]
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 [Pa]
85 P_w = numpy.zeros_like(S) # Effective pressure in channel segments [Pa]
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 positions at
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 channel
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 wedge
143 dSdt = S - S_old # adjust dSdt for clipping due to channel size limits
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 the
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.*24.)))
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 large
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 local
326 # meltwater production (m_dot)
327 Q = flux_solver(m_dot, ds)
328
329 # Find average shear stress from water flux for each channel segment
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 growth
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 water
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