timportant fix for water velocity and shear stress - 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
---
(DIR) commit 8119c0ae2e91d175f51c319d32a10abc108ebcd3
(DIR) parent cdd2e02befacebb825dd8718658f44c28fe27e82
(HTM) Author: Anders Damsgaard Christensen <adc@geo.au.dk>
Date: Wed, 1 Feb 2017 20:52:27 -0800
important fix for water velocity and shear stress
Diffstat:
A 1d-channel.py | 340 +++++++++++++++++++++++++++++++
D 1d-test.py | 322 -------------------------------
2 files changed, 340 insertions(+), 322 deletions(-)
---
(DIR) diff --git a/1d-channel.py b/1d-channel.py
t@@ -0,0 +1,340 @@
+#!/usr/bin/env python
+
+## ABOUT THIS FILE
+# The following script uses basic Python and Numpy functionality to solve the
+# coupled systems of equations describing subglacial channel development in
+# soft beds as presented in `Damsgaard et al. "Sediment plasticity controls
+# channelization of subglacial meltwater in soft beds"`, submitted to Journal
+# of Glaciology.
+#
+# High performance is not the goal for this implementation, which is instead
+# intended as a heavily annotated example on the solution procedure without
+# relying on solver libraries, suitable for low-level languages like C, Fortran
+# or CUDA.
+#
+# License: Gnu Public License v3
+# Author: Anders Damsgaard, adamsgaard@ucsd.edu, https://adamsgaard.dk
+
+import numpy
+import matplotlib.pyplot as plt
+import sys
+
+
+## Model parameters
+Ns = 25 # Number of nodes [-]
+Ls = 100e3 # Model length [m]
+t_end = 24.*60.*60.*2 # Total simulation time [s]
+tol_Q = 1e-3 # Tolerance criteria for the normalized max. residual for Q
+tol_P_c = 1e-3 # Tolerance criteria for the normalized max. residual for P_c
+max_iter = 1e2*Ns # Maximum number of solver iterations before failure
+output_convergence = False # Display convergence statistics during run
+
+# Physical parameters
+rho_w = 1000. # Water density [kg/m^3]
+rho_i = 910. # Ice density [kg/m^3]
+rho_s = 2600. # Sediment density [kg/m^3]
+g = 9.8 # Gravitational acceleration [m/s^2]
+theta = 30. # Angle of internal friction in sediment [deg]
+
+# Water source term [m/s]
+#m_dot = 7.93e-11
+#m_dot = 4.5e-8
+m_dot = 5.79e-5
+
+# Walder and Fowler 1994 sediment transport parameters
+K_e = 0.1 # Erosion constant [-], disabled when 0.0
+K_d = 6.0 # Deposition constant [-], disabled when 0.0
+#D50 = 1e-3 # Median grain size [m]
+#tau_c = 0.5*g*(rho_s - rho_i)*D50 # Critical shear stress for transport
+d15 = 1e-3 # Characteristic grain size [m]
+tau_c = 0.025*d15*g*(rho_s - rho_i) # Critical shear stress (Carter 2016)
+#tau_c = 0.
+mu_w = 1.787e-3 # Water viscosity [Pa*s]
+froude = 0.1 # Friction factor [-]
+v_s = d15**2.*g*2.*(rho_s - rho_i)/(9.*mu_w) # Settling velocity (Carter 2016)
+
+# Hewitt 2011 channel flux parameters
+manning = 0.1 # Manning roughness coefficient [m^{-1/3} s]
+F = rho_w*g*manning*(2.*(numpy.pi + 2)**2./numpy.pi)**(2./3.)
+
+# Channel growth-limit parameters
+c_1 = -0.118 # [m/kPa]
+c_2 = 4.60 # [m]
+
+# Minimum channel size [m^2], must be bigger than 0
+S_min = 1e-1
+
+
+## Initialize model arrays
+# Node positions, terminus at Ls
+s = numpy.linspace(0., Ls, Ns)
+ds = s[1:] - s[:-1]
+
+# Ice thickness and bed topography
+H = 6.*(numpy.sqrt(Ls - s + 5e3) - numpy.sqrt(5e3)) + 1.0 # max: 1.5 km
+#H = 1.*(numpy.sqrt(Ls - s + 5e3) - numpy.sqrt(5e3)) + 1.0 # max: 255 m
+#H = 0.6*(numpy.sqrt(Ls - s + 5e3) - numpy.sqrt(5e3)) + 1.0
+b = numpy.zeros_like(H)
+
+N = H*0.1*rho_i*g # Initial effective stress [Pa]
+p_w = rho_i*g*H - N # Initial guess of water pressure [Pa]
+hydro_pot = rho_w*g*b + p_w # Initial guess of hydraulic potential [Pa]
+
+# Initialize arrays for channel segments between nodes
+S = numpy.ones(len(s) - 1)*S_min # Cross-sectional area of channel segments[m^2]
+S_max = numpy.zeros_like(S) # Max. channel size [m^2]
+dSdt = numpy.empty_like(S) # Transient in channel cross-sectional area [m^2/s]
+W = S/numpy.tan(numpy.deg2rad(theta)) # Assuming no channel floor wedge
+Q = numpy.zeros_like(S) # Water flux in channel segments [m^3/s]
+Q_s = numpy.zeros_like(S) # Sediment flux in channel segments [m^3/s]
+N_c = numpy.zeros_like(S) # Effective pressure in channel segments [Pa]
+P_c = numpy.zeros_like(S) # Water pressure in channel segments [Pa]
+e_dot = numpy.zeros_like(S) # Sediment erosion rate in channel segments [m/s]
+d_dot = numpy.zeros_like(S) # Sediment deposition rate in chan. segments [m/s]
+c_bar = numpy.zeros_like(S) # Vertically integrated sediment content [m]
+tau = numpy.zeros_like(S) # Avg. shear stress from current [Pa]
+porosity = numpy.ones_like(S)*0.3 # Sediment porosity [-]
+res = numpy.zeros_like(S) # Solution residual during solver iterations
+
+
+## Helper functions
+def gradient(arr, arr_x):
+ # Central difference gradient of an array ``arr`` with node positions at
+ # ``arr_x``.
+ return (arr[:-1] - arr[1:])/(arr_x[:-1] - arr_x[1:])
+
+def avg_midpoint(arr):
+ # Averaged value of neighboring array elements
+ return (arr[:-1] + arr[1:])/2.
+
+def channel_water_flux(S, hydro_pot_grad):
+ # Hewitt 2011
+ return numpy.sqrt(1./F*S**(8./3.)*-hydro_pot_grad)
+
+def channel_shear_stress(Q, S):
+ # Weertman 1972, Walder and Fowler 1994
+ u_bar = Q/S
+ return 1./8.*froude*rho_w*u_bar**2.
+
+def channel_erosion_rate(tau):
+ # Parker 1979, Walder and Fowler 1994
+ return K_e*v_s*(tau - tau_c).clip(0.)/(g*(rho_s - rho_w)*d15)
+
+def channel_deposition_rate_kernel(tau, c_bar, ix):
+ # Parker 1979, Walder and Fowler 1994
+ return K_d*v_s*c_bar[ix]*(g*(rho_s - rho_w)*d15/tau[ix])**0.5
+
+def channel_deposition_rate(tau, c_bar, d_dot, Ns):
+ # Parker 1979, Walder and Fowler 1994
+ # Find deposition rate from upstream to downstream, margin at is=0
+
+ # No sediment deposition at upstream end
+ c_bar[0] = 0.
+ d_dot[0] = 0.
+ for ix in numpy.arange(1, Ns - 1):
+
+ # Net erosion in upstream cell
+ c_bar[ix] += numpy.maximum((e_dot[ix - 1] - d_dot[ix - 1])*dt, 0.)
+
+ d_dot[ix] = channel_deposition_rate_kernel(tau, c_bar, ix)
+
+ return d_dot, c_bar
+
+def channel_growth_rate(e_dot, d_dot, porosity, W):
+ # Damsgaard et al, in prep
+ return (e_dot - d_dot)/porosity*W
+
+def update_channel_size_with_limit(S, dSdt, dt, N):
+ # Damsgaard et al, in prep
+ S_max = ((c_1*N/1000. + c_2)*\
+ numpy.tan(numpy.deg2rad(theta))).clip(min=S_min)
+ S = numpy.minimum(S + dSdt*dt, S_max).clip(min=S_min)
+ W = S/numpy.tan(numpy.deg2rad(theta)) # Assume no channel floor wedge
+ return S, W, S_max
+
+def flux_solver(m_dot, ds):
+ # Iteratively find new fluxes
+ it = 0
+ max_res = 1e9 # arbitrary large value
+
+ # Iteratively find solution, do not settle for less iterations than the
+ # number of nodes
+ while max_res > tol_Q or it < Ns:
+
+ Q_old = Q.copy()
+ # dQ/ds = m_dot -> Q_out = m*delta(s) + Q_in
+ # Upwind information propagation (upwind)
+ Q[0] = 1e-2 # Ng 2000
+ Q[1:] = m_dot*ds[1:] + Q[:-1]
+ max_res = numpy.max(numpy.abs((Q - Q_old)/(Q + 1e-16)))
+
+ if output_convergence:
+ print('it = {}: max_res = {}'.format(it, max_res))
+
+ #import ipdb; ipdb.set_trace()
+ if it >= max_iter:
+ raise Exception('t = {}, step = {}:'.format(time, step) +
+ 'Iterative solution not found for Q')
+ it += 1
+
+ return Q
+
+def suspended_sediment_flux(c_bar, Q, S):
+ # Find the fluvial sediment flux through the system
+ # Q_s = c_bar * u * S, where u = Q/S
+ return c_bar*Q
+
+def pressure_solver(psi, F, Q, S):
+ # Iteratively find new water pressures
+ # dP_c/ds = psi - FQ^2/S^{8/3}
+
+ it = 0
+ max_res = 1e9 # arbitrary large value
+ while max_res > tol_P_c or it < Ns*40:
+
+ P_c_old = P_c.copy()
+
+ # Upwind finite differences
+ P_c[:-1] = -psi[:-1]*ds[:-1] \
+ + F*Q[:-1]**2./(S[:-1]**(8./3.))*ds[:-1] \
+ + P_c[1:] # Upstream
+
+ # Dirichlet BC (fixed pressure) at terminus
+ P_c[-1] = 0.
+
+ # von Neumann BC (no gradient = no flux) at s=0
+ P_c[0] = P_c[1]
+
+ max_res = numpy.max(numpy.abs((P_c - P_c_old)/(P_c + 1e-16)))
+
+ if output_convergence:
+ print('it = {}: max_res = {}'.format(it, max_res))
+
+ if it >= max_iter:
+ raise Exception('t = {}, step = {}:'.format(time, step) +
+ 'Iterative solution not found for P_c')
+ it += 1
+
+ return P_c
+
+def plot_state(step, time):
+ # Plot parameters along profile
+ fig = plt.gcf()
+ fig.set_size_inches(3.3, 3.3)
+
+ ax_Pa = plt.subplot(2, 1, 1) # axis with Pascals as y-axis unit
+ ax_Pa.plot(s_c/1000., P_c/1000., '--r', label='$P_c$')
+
+ ax_m3s = ax_Pa.twinx() # axis with m3/s as y-axis unit
+ ax_m3s.plot(s_c/1000., Q, '-b', label='$Q$')
+
+ plt.title('Day: {:.3}'.format(time/(60.*60.*24.)))
+ ax_Pa.legend(loc=2)
+ ax_m3s.legend(loc=1)
+ ax_Pa.set_ylabel('[kPa]')
+ ax_m3s.set_ylabel('[m$^3$/s]')
+
+ ax_m = plt.subplot(2, 1, 2, sharex=ax_Pa)
+ #ax_m.plot(s_c/1000., S, '-k', label='$S$')
+ #ax_m.plot(s_c/1000., S_max, '--k', label='$S_{max}$')
+ ax_m.semilogy(s_c/1000., S, '-k', label='$S$')
+ ax_m.semilogy(s_c/1000., S_max, '--k', label='$S_{max}$')
+
+ ax_ms = ax_m.twinx()
+ ax_ms.plot(s_c/1000., e_dot, '--r', label='$\dot{e}$')
+ ax_ms.plot(s_c/1000., d_dot, ':b', label='$\dot{d}$')
+
+ ax_m.legend(loc=2)
+ ax_ms.legend(loc=1)
+ ax_m.set_xlabel('$s$ [km]')
+ ax_m.set_ylabel('[m]')
+ ax_ms.set_ylabel('[m/s]')
+
+ plt.setp(ax_Pa.get_xticklabels(), visible=False)
+ plt.tight_layout()
+ if step == -1:
+ plt.savefig('chan-0.init.pdf')
+ else:
+ plt.savefig('chan-' + str(step) + '.pdf')
+ plt.clf()
+
+def find_new_timestep(ds, Q, S):
+ # Determine the timestep using the Courant-Friedrichs-Lewy condition
+ safety = 0.2
+ dt = safety*numpy.minimum(60.*60.*24., numpy.min(numpy.abs(ds/(Q*S))))
+
+ if dt < 1.0:
+ raise Exception('Error: Time step less than 1 second at step '
+ + '{}, time '.format(step)
+ + '{:.3} s/{:.3} d'.format(time, time/(60.*60.*24.)))
+
+ return dt
+
+def print_status_to_stdout(time, dt):
+ sys.stdout.write('\rt = {:.2} s or {:.4} d, dt = {:.2} s '\
+ .format(time, time/(60.*60.*24.), dt))
+ sys.stdout.flush()
+
+s_c = avg_midpoint(s) # Channel section midpoint coordinates [m]
+
+# Find gradient in hydraulic potential between the nodes
+hydro_pot_grad = gradient(hydro_pot, s)
+
+# Find field values at the middle of channel segments
+N_c = avg_midpoint(N)
+H_c = avg_midpoint(N)
+
+# Find fluxes in channel segments [m^3/s]
+Q = channel_water_flux(S, hydro_pot_grad)
+
+# Water-pressure gradient from geometry [Pa/m]
+psi = -rho_i*g*gradient(H, s) - (rho_w - rho_i)*g*gradient(b, s)
+
+# Prepare figure object for plotting during the simulation
+fig = plt.figure('channel')
+plot_state(-1, 0.0)
+
+
+## Time loop
+time = 0.; step = 0
+while time <= t_end:
+
+ dt = find_new_timestep(ds, Q, S)
+
+ print_status_to_stdout(time, dt)
+
+ # Find average shear stress from water flux for each channel segment
+ tau = channel_shear_stress(Q, S)
+
+ # Find sediment erosion and deposition rates for each channel segment
+ e_dot = channel_erosion_rate(tau)
+ d_dot, c_bar = channel_deposition_rate(tau, c_bar, d_dot, Ns)
+
+ # Determine change in channel size for each channel segment
+ dSdt = channel_growth_rate(e_dot, d_dot, porosity, W)
+
+ # Update channel cross-sectional area and width according to growth rate
+ # and size limit for each channel segment
+ S, W, S_max = update_channel_size_with_limit(S, dSdt, dt, N_c)
+
+ # Find new water fluxes consistent with mass conservation and local
+ # meltwater production (m_dot)
+ Q = flux_solver(m_dot, ds)
+
+ # Find the corresponding sediment flux
+ #Q_b = bedload_sediment_flux(
+ Q_s = suspended_sediment_flux(c_bar, Q, S)
+
+ # Find new water pressures consistent with the flow law
+ P_c = pressure_solver(psi, F, Q, S)
+
+ # Find new effective pressure in channel segments
+ N_c = rho_i*g*H_c - P_c
+
+ plot_state(step, time)
+
+ if step > 0:
+ break
+ # Update time
+ time += dt
+ step += 1
(DIR) diff --git a/1d-test.py b/1d-test.py
t@@ -1,322 +0,0 @@
-#!/usr/bin/env python
-
-## ABOUT THIS FILE
-# The following script uses basic Python and Numpy functionality to solve the
-# coupled systems of equations describing subglacial channel development in
-# soft beds as presented in `Damsgaard et al. "Sediment plasticity controls
-# channelization of subglacial meltwater in soft beds"`, submitted to Journal
-# of Glaciology.
-# High performance is not the goal for this implementation, which is instead
-# intended as a heavily annotated example on the solution procedure without
-# relying on solver libraries, suitable for low-level languages like C, Fortran
-# or CUDA.
-
-
-import numpy
-import matplotlib.pyplot as plt
-import sys
-
-## Model parameters
-Ns = 25 # Number of nodes [-]
-Ls = 100e3 # Model length [m]
-t_end = 24.*60.*60.*2 # Total simulation time [s]
-tol_Q = 1e-3 # Tolerance criteria for the normalized max. residual for Q
-tol_P_c = 1e-3 # Tolerance criteria for the normalized max. residual for P_c
-max_iter = 1e2*Ns # Maximum number of solver iterations before failure
-output_convergence = False
-
-# Physical parameters
-rho_w = 1000. # Water density [kg/m^3]
-rho_i = 910. # Ice density [kg/m^3]
-rho_s = 2600. # Sediment density [kg/m^3]
-g = 9.8 # Gravitational acceleration [m/s^2]
-theta = 30. # Angle of internal friction in sediment [deg]
-
-# Water source term [m/s]
-#m_dot = 7.93e-11
-m_dot = 5.79e-5
-#m_dot = 4.5e-8
-
-# Walder and Fowler 1994 sediment transport parameters
-K_e = 0.1 # Erosion constant [-], disabled when 0.0
-K_d = 6.0 # Deposition constant [-], disabled when 0.0
-#D50 = 1e-3 # Median grain size [m]
-#tau_c = 0.5*g*(rho_s - rho_i)*D50 # Critical shear stress for transport
-d15 = 1e-3 # Characteristic grain size [m]
-tau_c = 0.025*d15*g*(rho_s - rho_i) # Critical shear stress (Carter 2016)
-#tau_c = 0.
-mu_w = 1.787e-3 # Water viscosity [Pa*s]
-froude = 0.1 # Friction factor [-]
-v_s = d15**2.*g*2.*(rho_s - rho_i)/(9.*mu_w) # Settling velocity (Carter 2016)
-
-# Hewitt 2011 channel flux parameters
-manning = 0.1 # Manning roughness coefficient [m^{-1/3} s]
-F = rho_w*g*manning*(2.*(numpy.pi + 2)**2./numpy.pi)**(2./3.)
-
-# Channel growth-limit parameters
-c_1 = -0.118 # [m/kPa]
-c_2 = 4.60 # [m]
-
-# Minimum channel size [m^2], must be bigger than 0
-S_min = 1e-1
-
-
-## Initialize model arrays
-# Node positions, terminus at Ls
-s = numpy.linspace(0., Ls, Ns)
-ds = s[1:] - s[:-1]
-
-# Ice thickness and bed topography
-H = 6.*(numpy.sqrt(Ls - s + 5e3) - numpy.sqrt(5e3)) + 1.0 # max: 1.5 km
-#H = 1.*(numpy.sqrt(Ls - s + 5e3) - numpy.sqrt(5e3)) + 1.0 # max: 255 m
-#H = 0.6*(numpy.sqrt(Ls - s + 5e3) - numpy.sqrt(5e3)) + 1.0
-b = numpy.zeros_like(H)
-
-N = H*0.1*rho_i*g # Initial effective stress [Pa]
-p_w = rho_i*g*H - N # Initial guess of water pressure [Pa]
-hydro_pot = rho_w*g*b + p_w # Initial guess of hydraulic potential [Pa]
-
-# Initialize arrays for channel segments between nodes
-S = numpy.ones(len(s) - 1)*S_min # Cross-sectional area of channel segments[m^2]
-S_max = numpy.zeros_like(S) # Max. channel size [m^2]
-dSdt = numpy.empty_like(S) # Transient in channel cross-sectional area [m^2/s]
-W = S/numpy.tan(numpy.deg2rad(theta)) # Assuming no channel floor wedge
-Q = numpy.zeros_like(S) # Water flux in channel segments [m^3/s]
-Q_s = numpy.zeros_like(S) # Sediment flux in channel segments [m^3/s]
-N_c = numpy.zeros_like(S) # Effective pressure in channel segments [Pa]
-P_c = numpy.zeros_like(S) # Water pressure in channel segments [Pa]
-e_dot = numpy.zeros_like(S) # Sediment erosion rate in channel segments [m/s]
-d_dot = numpy.zeros_like(S) # Sediment deposition rate in chan. segments [m/s]
-c_bar = numpy.zeros_like(S) # Vertically integrated sediment content [m]
-tau = numpy.zeros_like(S) # Avg. shear stress from current [Pa]
-porosity = numpy.ones_like(S)*0.3 # Sediment porosity [-]
-res = numpy.zeros_like(S) # Solution residual during solver iterations
-
-
-## Helper functions
-def gradient(arr, arr_x):
- # Central difference gradient of an array ``arr`` with node positions at
- # ``arr_x``.
- return (arr[:-1] - arr[1:])/(arr_x[:-1] - arr_x[1:])
-
-def avg_midpoint(arr):
- # Averaged value of neighboring array elements
- return (arr[:-1] + arr[1:])/2.
-
-def channel_water_flux(S, hydro_pot_grad):
- # Hewitt 2011
- return numpy.sqrt(1./F*S**(8./3.)*-hydro_pot_grad)
-
-def channel_shear_stress(Q, S):
- # Weertman 1972, Walder and Fowler 1994
- u_bar = Q*S
- return 1./8.*froude*rho_w*u_bar**2.
-
-def channel_erosion_rate(tau):
- # Parker 1979, Walder and Fowler 1994
- return K_e*v_s*(tau - tau_c).clip(0.)/(g*(rho_s - rho_w)*d15)
-
-def channel_deposition_rate_kernel(tau, c_bar, ix):
- # Parker 1979, Walder and Fowler 1994
- return K_d*v_s*c_bar[ix]*(g*(rho_s - rho_w)*d15/tau[ix])**0.5
-
-def channel_deposition_rate(tau, c_bar, d_dot, Ns):
- # Parker 1979, Walder and Fowler 1994
- # Find deposition rate from upstream to downstream, margin at is=0
- #for ix in numpy.arange(Ns-2, -1, -1):
- for ix in numpy.arange(Ns - 1):
- if ix == 0: # No sediment deposition at upstream end
- c_bar[ix] = 0.
- d_dot[ix] = 0.
- else:
- c_bar[ix] = (e_dot[ix - 1] - d_dot[ix - 1])*dt # Net erosion upstr.
- d_dot[ix] = channel_deposition_rate_kernel(tau, c_bar, ix)
- return d_dot, c_bar
-
-def channel_growth_rate(e_dot, d_dot, porosity, W):
- # Damsgaard et al, in prep
- return (e_dot - d_dot)/porosity*W
-
-def update_channel_size_with_limit(S, dSdt, dt, N):
- # Damsgaard et al, in prep
- S_max = ((c_1*N/1000. + c_2)*\
- numpy.tan(numpy.deg2rad(theta))).clip(min=S_min)
- S = numpy.minimum(S + dSdt*dt, S_max).clip(min=S_min)
- W = S/numpy.tan(numpy.deg2rad(theta)) # Assume no channel floor wedge
- return S, W, S_max
-
-def flux_solver(m_dot, ds):
- # Iteratively find new fluxes
- it = 0
- max_res = 1e9 # arbitrary large value
-
- # Iteratively find solution, do not settle for less iterations than the
- # number of nodes
- while max_res > tol_Q or it < Ns:
-
- Q_old = Q.copy()
- # dQ/ds = m_dot -> Q_out = m*delta(s) + Q_in
- # Upwind information propagation (upwind)
- Q[0] = 1e-2 # Ng 2000
- Q[1:] = m_dot*ds[1:] + Q[:-1]
- max_res = numpy.max(numpy.abs((Q - Q_old)/(Q + 1e-16)))
-
- if output_convergence:
- print('it = {}: max_res = {}'.format(it, max_res))
-
- #import ipdb; ipdb.set_trace()
- if it >= max_iter:
- raise Exception('t = {}, step = {}:'.format(time, step) +
- 'Iterative solution not found for Q')
- it += 1
-
- return Q
-
-def pressure_solver(psi, F, Q, S):
- # Iteratively find new water pressures
- # dP_c/ds = psi - FQ^2/S^{8/3}
-
- it = 0
- max_res = 1e9 # arbitrary large value
- while max_res > tol_P_c or it < Ns*40:
-
- P_c_old = P_c.copy()
-
- # Upwind finite differences
- P_c[:-1] = -psi[:-1]*ds[:-1] \
- + F*Q[:-1]**2./(S[:-1]**(8./3.))*ds[:-1] \
- + P_c[1:] # Upstream
-
- # Dirichlet BC (fixed pressure) at terminus
- P_c[-1] = 0.
-
- # von Neumann BC (no gradient = no flux) at s=0
- P_c[0] = P_c[1]
-
- max_res = numpy.max(numpy.abs((P_c - P_c_old)/(P_c + 1e-16)))
-
- if output_convergence:
- print('it = {}: max_res = {}'.format(it, max_res))
-
- if it >= max_iter:
- raise Exception('t = {}, step = {}:'.format(time, step) +
- 'Iterative solution not found for P_c')
- it += 1
-
- return P_c
-
-def plot_state(step, time):
- # Plot parameters along profile
- fig = plt.gcf()
- fig.set_size_inches(3.3, 3.3)
-
- ax_Pa = plt.subplot(2, 1, 1) # axis with Pascals as y-axis unit
- ax_Pa.plot(s_c/1000., P_c/1000., '--r', label='$P_c$')
-
- ax_m3s = ax_Pa.twinx() # axis with m3/s as y-axis unit
- ax_m3s.plot(s_c/1000., Q, '-b', label='$Q$')
-
- plt.title('Day: {:.3}'.format(time/(60.*60.*24.)))
- ax_Pa.legend(loc=2)
- ax_m3s.legend(loc=1)
- ax_Pa.set_ylabel('[kPa]')
- ax_m3s.set_ylabel('[m$^3$/s]')
-
- ax_m = plt.subplot(2, 1, 2, sharex=ax_Pa)
- #ax_m.plot(s_c/1000., S, '-k', label='$S$')
- #ax_m.plot(s_c/1000., S_max, '--k', label='$S_{max}$')
- ax_m.semilogy(s_c/1000., S, '-k', label='$S$')
- ax_m.semilogy(s_c/1000., S_max, '--k', label='$S_{max}$')
-
- ax_ms = ax_m.twinx()
- ax_ms.plot(s_c/1000., e_dot, '--r', label='$\dot{e}$')
- ax_ms.plot(s_c/1000., d_dot, ':b', label='$\dot{d}$')
-
- ax_m.legend(loc=2)
- ax_ms.legend(loc=1)
- ax_m.set_xlabel('$s$ [km]')
- ax_m.set_ylabel('[m]')
- ax_ms.set_ylabel('[m/s]')
-
- plt.setp(ax_Pa.get_xticklabels(), visible=False)
- plt.tight_layout()
- if step == -1:
- plt.savefig('chan-0.init.pdf')
- else:
- plt.savefig('chan-' + str(step) + '.pdf')
- plt.clf()
-
-def find_new_timestep(ds, Q, S):
- # Determine the timestep using the Courant-Friedrichs-Lewy condition
- safety = 0.2
- dt = safety*numpy.minimum(60.*60.*24., numpy.min(numpy.abs(ds/(Q*S))))
-
- if dt < 1.0:
- raise Exception('Error: Time step less than 1 second at step '
- + '{}, time '.format(step)
- + '{:.3} s/{:.3} d'.format(time, time/(60.*60.*24.)))
-
- return dt
-
-def print_status_to_stdout(time, dt):
- sys.stdout.write('\rt = {:.2} s or {:.4} d, dt = {:.2} s '\
- .format(time, time/(60.*60.*24.), dt))
- sys.stdout.flush()
-
-s_c = avg_midpoint(s) # Channel section midpoint coordinates [m]
-
-# Find gradient in hydraulic potential between the nodes
-hydro_pot_grad = gradient(hydro_pot, s)
-
-# Find field values at the middle of channel segments
-N_c = avg_midpoint(N)
-H_c = avg_midpoint(N)
-
-# Find fluxes in channel segments [m^3/s]
-Q = channel_water_flux(S, hydro_pot_grad)
-
-# Water-pressure gradient from geometry [Pa/m]
-psi = -rho_i*g*gradient(H, s) - (rho_w - rho_i)*g*gradient(b, s)
-
-# Prepare figure object for plotting during the simulation
-fig = plt.figure('channel')
-plot_state(-1, 0.0)
-
-
-## Time loop
-time = 0.; step = 0
-while time <= t_end:
-
- dt = find_new_timestep(ds, Q, S)
-
- print_status_to_stdout(time, dt)
-
- # Find average shear stress from water flux for each channel segment
- tau = channel_shear_stress(Q, S)
-
- # Find sediment erosion and deposition rates for each channel segment
- e_dot = channel_erosion_rate(tau)
- d_dot, c_bar = channel_deposition_rate(tau, c_bar, d_dot, Ns)
-
- # Determine change in channel size for each channel segment
- dSdt = channel_growth_rate(e_dot, d_dot, porosity, W)
-
- # Update channel cross-sectional area and width according to growth rate
- # and size limit for each channel segment
- S, W, S_max = update_channel_size_with_limit(S, dSdt, dt, N_c)
-
- # Find new water fluxes consistent with mass conservation and local
- # meltwater production (m_dot)
- Q = flux_solver(m_dot, ds)
-
- # Find new water pressures consistent with the flow law
- P_c = pressure_solver(psi, F, Q, S)
-
- # Find new effective pressure in channel segments
- N_c = rho_i*g*H_c - P_c
-
- plot_state(step, time)
-
- # Update time
- time += dt
- step += 1