commit 400da52ca6dea9a381f902ba653ca1c267ac286f
parent 26bc00c665916c00d8fc195170f6df2062d77332
Author: Anders Damsgaard Christensen <adc@geo.au.dk>
Date: Thu, 2 Feb 2017 16:22:16 -0800
transport sediment for more than one cell
Diffstat:
| A | 1d-channel-flux.py | | | 292 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| M | 1d-channel.py | | | 153 | ++++++++++++++++++++++++++++++++++++++++++++----------------------------------- |
2 files changed, 377 insertions(+), 68 deletions(-)
diff --git a/1d-channel-flux.py b/1d-channel-flux.py
@@ -0,0 +1,292 @@
+#!/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.*120. # 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
+safety = 0.1 # Safety factor ]0;1] for adaptive timestepping
+
+# 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
+
+# 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
+S_min = 1.5e-2
+
+
+# # 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-sect. area of channel segments [m^2]
+S_max = numpy.zeros_like(S) # Max. channel size [m^2]
+dSdt = numpy.zeros_like(S) # Transient in channel cross-sect. 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]
+dQ_s_ds = numpy.empty_like(S) # Transient in channel cross-sect. area [m^2/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]
+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 update_channel_size_with_limit(S, dSdt, dt, N):
+ # Damsgaard et al, in prep
+ S_max = ((c_1*N.clip(min=0.)/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 sediment_flux(Q):
+ #return Q**(3./2.)
+ return Q/2.
+
+def sediment_flux_divergence(Q_s, ds):
+ # Damsgaard et al, in prep
+ return (Q_s[1:] - Q_s[:-1])/ds[1:]
+
+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$')
+ ax_m3s.plot(s_c/1000., Q_s, ':b', label='$Q_s$')
+
+ 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_m2 = plt.subplot(2, 1, 2, sharex=ax_Pa)
+ ax_m2.plot(s_c/1000., S, '-k', label='$S$')
+ ax_m2.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_m2s = ax_m2.twinx()
+ ax_m2s.plot(s_c/1000., dSdt, ':b', label='$dS/dt$')
+
+ ax_m2.legend(loc=2)
+ ax_m2s.legend(loc=1)
+ ax_m2.set_xlabel('$s$ [km]')
+ ax_m2.set_ylabel('[m$^2$]')
+ ax_m2s.set_ylabel('[m$^2$/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
+ 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:
+
+ #print('@ @ @ step ' + str(step))
+
+ dt = find_new_timestep(ds, Q, S)
+
+ print_status_to_stdout(time, dt)
+
+ # Find the sediment flux
+ Q_s = sediment_flux(Q)
+
+ # Find sediment flux divergence which determines channel growth, no growth
+ # in first segment
+ dSdt[1:] = sediment_flux_divergence(Q_s, ds)
+
+ # 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)
+
+ #import ipdb; ipdb.set_trace()
+ #if step > 0:
+ #break
+
+ # Update time
+ time += dt
+ step += 1
diff --git a/1d-channel.py b/1d-channel.py
@@ -1,6 +1,6 @@
#!/usr/bin/env python
-## ABOUT THIS FILE
+# # 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
@@ -20,38 +20,37 @@ 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
+# # 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
safety = 0.1 # Safety factor ]0;1] for adaptive timestepping
# 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]
+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 = 7.93e-11
m_dot = 4.5e-7
-#m_dot = 5.79e-5
+# 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
-K_d = 1e-1 # Deposition constant [-], disabled when 0.0
+K_d = 6.0 # Deposition constant [-], disabled when 0.0
alpha = 2e5 # Geometric correction factor (Carter et al 2017)
-#D50 = 1e-3 # Median grain size [m]
-#tau_c = 0.5*g*(rho_s - rho_i)*D50 # Critical shear stress for transport
+# 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 2017)
-#tau_c = 0.
+# 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 2017)
@@ -61,22 +60,22 @@ 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]
+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
+# # 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
+# 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]
@@ -84,79 +83,88 @@ 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]
+S = numpy.ones(len(s) - 1)*S_min # Cross-sect. 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-sect. 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 concentration [-]
-tau = numpy.zeros_like(S) # Avg. shear stress from current [Pa]
-porosity = numpy.ones_like(S)*0.3 # Sediment porosity [-]
+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 concentration [-]
+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
+# # 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(min=0.)/(g*(rho_s - rho_w)*d15)
+ # return K_e*v_s*(tau - tau_c).clip(min=0.)/(g*(rho_s - rho_w)*d15)
# Carter et al 2017
return K_e*v_s/alpha*(tau - tau_c).clip(min=0.)/(g*(rho_s - rho_w)*d15)
+
def channel_deposition_rate_kernel(tau, c_bar, ix):
# Parker 1979, Walder and Fowler 1994
- #result = K_d*v_s*c_bar[ix]*(g*(rho_s - rho_w)*d15/tau[ix])**0.5
+ # result = K_d*v_s*c_bar[ix]*(g*(rho_s - rho_w)*d15/tau[ix])**0.5
# Carter et al. 2017
result = K_d*v_s/alpha*c_bar[ix]*(g*(rho_s - rho_w)*d15/tau[ix])**0.5
+ '''
print('tau[{}] = {}'.format(ix, tau[ix]))
print('c_bar[{}] = {}'.format(ix, c_bar[ix]))
print('e_dot[{}] = {}'.format(ix, e_dot[ix]))
print('d_dot[{}] = {}'.format(ix, result))
print('')
+ '''
return result
+
def channel_deposition_rate_kernel_ng(c_bar, ix):
# Ng 2000
h = W[ix]/2.*numpy.tan(numpy.deg2rad(theta))
- epsilon = numpy.sqrt((psi[ix] - (P_c[ix] - P_c[ix - 1])/ds[ix])\
- /(rho_w*froude))*h**(3./2.)
+ epsilon = numpy.sqrt((psi[ix] - (P_c[ix] - P_c[ix - 1])/ds[ix])
+ / (rho_w*froude))*h**(3./2.)
return v_s/epsilon*c_bar[ix]
+
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
-
+ '''
print("\n## Before loop:")
print(c_bar)
print(d_dot)
print('')
-
+ '''
# No sediment deposition at upstream end
c_bar[0] = 0.
@@ -164,36 +172,39 @@ def channel_deposition_rate(tau, c_bar, d_dot, Ns):
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*ds[ix-1], 0.)
- c_bar[ix] = numpy.maximum(
- W[ix - 1]*ds[ix - 1]*rho_s/rho_w*
- (e_dot[ix - 1] - d_dot[ix - 1])/Q[ix - 1]
- , 0.)
-
- #d_dot[ix] = channel_deposition_rate_kernel(tau, c_bar, ix)
- d_dot[ix] = channel_deposition_rate_kernel_ng(c_bar, ix)
+ # c_bar[ix] = numpy.maximum((e_dot[ix-1]-d_dot[ix-1])*dt*ds[ix-1], 0.)
+ c_bar[ix] = c_bar[ix - 1] + \
+ numpy.maximum(
+ W[ix - 1]*ds[ix - 1]*rho_s/rho_w *
+ (e_dot[ix - 1] - d_dot[ix - 1])/Q[ix - 1], 0.)
+ d_dot[ix] = channel_deposition_rate_kernel(tau, c_bar, ix)
+ # d_dot[ix] = channel_deposition_rate_kernel_ng(c_bar, ix)
+ '''
print("\n## After loop:")
print(c_bar)
print(d_dot)
print('')
-
+ '''
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)*\
+ S_max = ((c_1*N.clip(min=0.)/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
@@ -213,7 +224,7 @@ def flux_solver(m_dot, ds):
if output_convergence:
print('it = {}: max_res = {}'.format(it, max_res))
- #import ipdb; ipdb.set_trace()
+ # import ipdb; ipdb.set_trace()
if it >= max_iter:
raise Exception('t = {}, step = {}:'.format(time, step) +
'Iterative solution not found for Q')
@@ -221,11 +232,13 @@ def flux_solver(m_dot, ds):
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}
@@ -259,6 +272,7 @@ def pressure_solver(psi, F, Q, S):
return P_c
+
def plot_state(step, time):
# Plot parameters along profile
fig = plt.gcf()
@@ -277,10 +291,10 @@ def plot_state(step, time):
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_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}$')
@@ -300,6 +314,7 @@ def plot_state(step, time):
plt.savefig('chan-' + str(step) + '.pdf')
plt.clf()
+
def find_new_timestep(ds, Q, S):
# Determine the timestep using the Courant-Friedrichs-Lewy condition
dt = safety*numpy.minimum(60.*60.*24., numpy.min(numpy.abs(ds/(Q*S))))
@@ -311,12 +326,13 @@ def find_new_timestep(ds, Q, S):
return dt
+
def print_status_to_stdout(time, dt):
- sys.stdout.write('\rt = {:.2} s or {:.4} d, dt = {:.2} s '\
+ 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]
+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)
@@ -336,11 +352,12 @@ fig = plt.figure('channel')
plot_state(-1, 0.0)
-## Time loop
-time = 0.; step = 0
+# # Time loop
+time = 0.
+step = 0
while time <= t_end:
- #print('@ @ @ step ' + str(step))
+ # print('@ @ @ step ' + str(step))
dt = find_new_timestep(ds, Q, S)
@@ -365,7 +382,7 @@ while time <= t_end:
Q = flux_solver(m_dot, ds)
# Find the corresponding sediment flux
- #Q_b = bedload_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
@@ -376,9 +393,9 @@ while time <= t_end:
plot_state(step, time)
- #import ipdb; ipdb.set_trace()
- if step > 0:
- break
+ # import ipdb; ipdb.set_trace()
+ #if step > 0:
+ #break
# Update time
time += dt
