1d-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