commit 6afcc6591b23c355ae2cae50953afbba1de9e0d6
parent 4bb1e340b01ce37312e12942c05715d30529e2b0
Author: Anders Damsgaard <andersd@riseup.net>
Date: Mon, 15 May 2017 16:11:43 -0400
implement Meyer-Peter and Muller 1948 relationship for sediment transport
Diffstat:
| M | 1d-channel.py | | | 141 | ++++++++++++++++++++++++------------------------------------------------------- |
1 file changed, 43 insertions(+), 98 deletions(-)
diff --git a/1d-channel.py b/1d-channel.py
@@ -13,7 +13,7 @@
# or CUDA.
#
# License: Gnu Public License v3
-# Author: Anders Damsgaard, adamsgaard@ucsd.edu, https://adamsgaard.dk
+# Author: Anders Damsgaard, andersd@princeton.edu, https://adamsgaard.dk
import numpy
import matplotlib.pyplot as plt
@@ -22,9 +22,10 @@ import sys
# # Model parameters
Ns = 25 # Number of nodes [-]
-Ls = 10e3 # Model length [m]
+Ls = 1e3 # Model length [m]
total_days = 60. # Total simulation time [d]
t_end = 24.*60.*60.*total_days # Total simulation time [s]
+tol_S = 1e-3 # Tolerance criteria for the norm. max. residual for Q
tol_Q = 1e-3 # Tolerance criteria for the norm. max. residual for Q
tol_N_c = 1e-3 # Tolerance criteria for the norm. max. residual for N_c
max_iter = 1e2*Ns # Maximum number of solver iterations before failure
@@ -32,7 +33,9 @@ print_output_convergence = False # Display convergence in nested loops
print_output_convergence_main = True # Display convergence in main loop
safety = 0.01 # Safety factor ]0;1] for adaptive timestepping
plot_interval = 20 # Time steps between plots
+plot_during_iterations = False # Generate plots for intermediate results
speedup_factor = 1. # Speed up channel growth to reach steady state faster
+# relax = 0.05 # Relaxation parameter for effective pressure ]0;1]
# Physical parameters
rho_w = 1000. # Water density [kg/m^3]
@@ -40,15 +43,12 @@ 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]
-sand_fraction = 0.5 # Initial volumetric fraction of sand relative to gravel
-D_g = 5e-3 # Mean grain size in gravel fraction (> 2 mm) [m]
-D_s = 5e-4 # Mean grain size in sand fraction (<= 2 mm) [m]
-#D_g = 1
-#D_g = 0.1
+D = 1.15e-3 # Mean grain size [m], Lajeuness et al 2010, series 1
+tau_c = 0.016 # Critical Shields stress, Lajeunesse et al 2010, series 1
# Boundary conditions
P_terminus = 0. # Water pressure at terminus [Pa]
-m_dot = 1e-6 # Water source term [m/s]
+m_dot = numpy.linspace(0., 1e-5, Ns-1) # Water source term [m/s]
Q_upstream = 1e-5 # Water influx upstream (must be larger than 0) [m^3/s]
# Channel hydraulic properties
@@ -61,8 +61,6 @@ c_2 = 4.60 # [m]
# Minimum channel size [m^2], must be bigger than 0
S_min = 1e-2
-# S_min = 1e-1
-# S_min = 1.
# # Initialize model arrays
@@ -81,7 +79,7 @@ b = numpy.zeros_like(H)
N = H*0.1*rho_i*g # Initial effective stress [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 = numpy.ones(len(s) - 1)*0.1 # 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
@@ -92,10 +90,6 @@ P_c = numpy.zeros_like(S) # Water pressure in channel segments [Pa]
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
-Q_t = numpy.zeros_like(S) # Total sediment flux [m3/s]
-Q_s = numpy.zeros_like(S) # Sediment flux where D <= 2 mm [m3/s]
-Q_g = numpy.zeros_like(S) # Sediment flux where D > 2 mm [m3/s]
-f_s = numpy.ones_like(S)*sand_fraction # Initial sediment fraction of sand [-]
# # Helper functions
@@ -123,75 +117,23 @@ def channel_shear_stress(Q, S):
return 1./8.*friction_factor*rho_w*u_bar**2.
-def channel_sediment_flux_sand(tau, W, f_s, D_s):
- # Parker 1979, Wilcock 1997, 2001, Egholm 2013
+def channel_sediment_flux(tau, W):
+ # Meyer-Peter and Mueller 1948
# tau: Shear stress by water flow
# W: Channel width
- # f_s: Sand volume fraction
- # D_s: Mean sand fraction grain size
-
- # Piecewise linear functions for nondimensional critical shear stresses
- # dependent on sand fraction from Gasparini et al 1999 of Wilcock 1997
- # data.
- ref_shear_stress = numpy.ones_like(f_s)*0.04
- ref_shear_stress[numpy.nonzero(f_s <= 0.1)] = 0.88
- I = numpy.nonzero((0.1 < f_s) & (f_s <= 0.4))
- ref_shear_stress[I] = 0.88 - 2.8*(f_s[I] - 0.1)
-
- # Non-dimensionalize shear stress
- shields_stress = tau/((rho_s - rho_w)*g*D_s)
-
- # import ipdb; ipdb.set_trace()
- Q_c = 11.2*f_s*W/((rho_s - rho_w)/rho_w*g) \
- * (tau/rho_w)**1.5 \
- * numpy.maximum(0.0,
- (1.0 - 0.846*numpy.sqrt(ref_shear_stress/shields_stress))
- )**4.5
-
- return Q_c
-
-
-def channel_sediment_flux_gravel(tau, W, f_g, D_g):
- # Parker 1979, Wilcock 1997, 2001, Egholm 2013
- # tau: Shear stress by water flow
- # W: Channel width
- # f_g: Gravel volume fraction
- # D_g: Mean gravel fraction grain size
-
- # Piecewise linear functions for nondimensional critical shear stresses
- # dependent on sand fraction from Gasparini et al 1999 of Wilcock 1997
- # data.
- ref_shear_stress = numpy.ones_like(f_g)*0.01
- ref_shear_stress[numpy.nonzero(f_g <= 0.1)] = 0.04
- I = numpy.nonzero((0.1 < f_g) & (f_g <= 0.4))
- ref_shear_stress[I] = 0.04 - 0.1*(f_g[I] - 0.1)
# Non-dimensionalize shear stress
- shields_stress = tau/((rho_s - rho_w)*g*D_g)
-
- # From Wilcock 2001, eq. 3
- Q_g = 11.2*f_g*W/((rho_s - rho_w)/rho_w*g) \
- * (tau/rho_w)**1.5 \
- * numpy.maximum(0.0,
- (1.0 - 0.846*ref_shear_stress/shields_stress))**4.5
-
- # From Wilcock 2001, eq. 4
- I = numpy.nonzero(ref_shear_stress/shields_stress < 1.)
- Q_g[I] = f_g[I]*W[I]/((rho_s - rho_w)/rho_w*g) \
- * (tau[I]/rho_w)**1.5 \
- * 0.0025*(shields_stress[I]/ref_shear_stress[I])**14.2
+ shields_stress = tau/((rho_s - rho_w)*g*D)
- return Q_g
+ stress_excess = shields_stress - tau_c
+ stress_excess[stress_excess < 0.] = 0.
+ return 8.*stress_excess**(3./2.)*W \
+ * numpy.sqrt((rho_s - rho_w)/rho_w*g*D**3.)
-def channel_growth_rate(e_dot, d_dot, W):
+def channel_growth_rate_sedflux(Q_s, porosity, s_c):
# Damsgaard et al, in prep
- return (e_dot - d_dot)*W
-
-
-def channel_growth_rate_sedflux(Q_t, porosity, s_c):
- # Damsgaard et al, in prep
- return 1./porosity[1:] * gradient(Q_t, s_c)
+ return 1./porosity[1:] * gradient(Q_s, s_c)
def update_channel_size_with_limit(S, S_old, dSdt, dt, N_c):
@@ -213,13 +155,13 @@ def flux_solver(m_dot, ds):
# Iteratively find solution, do not settle for less iterations than the
# number of nodes
- while max_res > tol_Q or it < Ns:
+ while max_res > tol_Q:
Q_old = Q.copy()
# dQ/ds = m_dot -> Q_out = m*delta(s) + Q_in
# Upwind information propagation (upwind)
Q[0] = Q_upstream
- Q[1:] = m_dot*ds[1:] + Q[:-1]
+ Q[1:] = m_dot[1:]*ds[1:] + Q[:-1]
max_res = numpy.max(numpy.abs((Q - Q_old)/(Q + 1e-16)))
if print_output_convergence:
@@ -240,7 +182,7 @@ def pressure_solver(psi, f, Q, S):
it = 0
max_res = 1e9 # arbitrary large value
- while max_res > tol_N_c or it < Ns:
+ while max_res > tol_N_c:
N_c_old = N_c.copy()
@@ -250,7 +192,7 @@ def pressure_solver(psi, f, Q, S):
N_c[:-1] = N_c[1:] \
+ psi[:-1]*ds[:-1] \
- f[:-1]*rho_w*g*Q[:-1]*numpy.abs(Q[:-1]) \
- /(S[:-1]**(8./3.))*ds[:-1]
+ / (S[:-1]**(8./3.))*ds[:-1]
max_res = numpy.max(numpy.abs((N_c - N_c_old)/(N_c + 1e-16)))
@@ -263,6 +205,7 @@ def pressure_solver(psi, f, Q, S):
it += 1
return N_c
+ # return N_c_old*(1 - relax_N_c) + N_c*relax_N_c
def plot_state(step, time, S_, S_max_, title=True):
@@ -286,9 +229,7 @@ def plot_state(step, time, S_, S_max_, title=True):
ax_m3s.set_ylabel('[m$^3$/s]')
ax_m3s_sed = plt.subplot(3, 1, 2, sharex=ax_Pa)
- ax_m3s_sed.plot(s_c/1000., Q_g, ':', label='$Q_{gravel}$')
- ax_m3s_sed.plot(s_c/1000., Q_s, '-', label='$Q_{sand}$')
- ax_m3s_sed.plot(s_c/1000., Q_t, '--', label='$Q_{total}$')
+ ax_m3s_sed.plot(s_c/1000., Q_s, '-', label='$Q_{s}$')
ax_m3s_sed.set_ylabel('[m$^3$/s]')
ax_m3s_sed.legend(loc=2)
@@ -314,11 +255,14 @@ def plot_state(step, time, S_, S_max_, title=True):
else:
plt.savefig('chan-' + str(step) + '.pdf')
plt.clf()
+ plt.close()
-def find_new_timestep(ds, Q, S):
+def find_new_timestep(ds, Q, Q_s, S):
# Determine the timestep using the Courant-Friedrichs-Lewy condition
- dt = safety*numpy.minimum(60.*60.*24., numpy.min(numpy.abs(ds/(Q*S))))
+ dt = safety*numpy.minimum(60.*60.*24.,
+ numpy.min(numpy.abs(ds/(Q*S),
+ ds/(Q_s*S)+1e-16)))
if dt < 1.0:
raise Exception('Error: Time step less than 1 second at step '
@@ -334,6 +278,8 @@ def print_status_to_stdout(step, time, dt):
.format(time, time/(60.*60.*24.), dt))
sys.stdout.flush()
+
+# Initialize remaining values before entering time loop
s_c = avg_midpoint(s) # Channel section midpoint coordinates [m]
H_c = avg_midpoint(H)
@@ -351,7 +297,7 @@ step = 0
while time <= t_end:
# Determine time step length from water flux
- dt = find_new_timestep(ds, Q, S)
+ dt = find_new_timestep(ds, Q, Q_s, S)
# Display current simulation status
print_status_to_stdout(step, time, dt)
@@ -363,10 +309,8 @@ while time <= t_end:
max_res = 1e9
S_old = S.copy()
- # Iteratively find solution, do not settle for less iterations than the
- # number of nodes to make sure information has had a chance to pass through
- # the system
- while max_res > tol_Q or it < Ns:
+ # Iteratively find solution with the Jacobi relaxation method
+ while max_res > tol_S:
S_prev_it = S.copy()
@@ -377,23 +321,21 @@ while time <= t_end:
# Find average shear stress from water flux for each channel segment
tau = channel_shear_stress(Q, S)
- # Determine sediment fluxes for each size fraction
- f_g = 1./f_s # gravel volume fraction is reciprocal to sand
- Q_s = channel_sediment_flux_sand(tau, W, f_s, D_s)
- Q_g = channel_sediment_flux_gravel(tau, W, f_g, D_g)
- Q_t = Q_s + Q_g
+ # Determine sediment flux
+ Q_s = channel_sediment_flux(tau, W)
# Determine change in channel size for each channel segment.
# Use backward differences and assume dS/dt=0 in first segment.
- #dSdt[1:] = channel_growth_rate_sedflux(Q_t, porosity, s_c)
- #dSdt *= speedup_factor
+ dSdt[1:] = channel_growth_rate_sedflux(Q_s, porosity, s_c)
+ # dSdt *= speedup_factor * relax
# Update channel cross-sectional area and width according to growth
# rate and size limit for each channel segment
+ # S_prev = S.copy()
S, W, S_max, dSdt = \
update_channel_size_with_limit(S, S_old, dSdt, dt, N_c)
+ # S = S_prev*(1.0 - relax) + S*relax
- # Find hydraulic roughness
f = channel_hydraulic_roughness(manning, S, W, theta)
# Find new water pressures consistent with the flow law
@@ -402,6 +344,9 @@ while time <= t_end:
# Find new effective pressure in channel segments
P_c = rho_i*g*H_c - N_c
+ if plot_during_iterations:
+ plot_state(step + it/1e4, time, S, S_max)
+
# Find new maximum normalized residual value
max_res = numpy.max(numpy.abs((S - S_prev_it)/(S + 1e-16)))
if print_output_convergence_main:
