tremove legacy sediment model, solve for N_c instead of P_c - 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 d0ec40ca4dde58229f1c40d4fc0d3789a3d679a2
(DIR) parent 736f50c0a19c6d0ab65821b49e60a0fc010e0f84
(HTM) Author: Anders Damsgaard <andersd@riseup.net>
Date: Wed, 8 Mar 2017 17:23:29 -0800
remove legacy sediment model, solve for N_c instead of P_c
Diffstat:
M 1d-channel.py | 102 ++++++--------------------------
1 file changed, 17 insertions(+), 85 deletions(-)
---
(DIR) diff --git a/1d-channel.py b/1d-channel.py
t@@ -27,7 +27,7 @@ Ls = 1e3 # Model length [m]
total_days = 60. # Total simulation time [d]
t_end = 24.*60.*60.*total_days # 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 norm. max. residual for P_c
+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
print_output_convergence = False # Display convergence statistics during run
safety = 0.01 # Safety factor ]0;1] for adaptive timestepping
t@@ -44,8 +44,8 @@ sand_fraction = 0.5 # Initial volumetric fraction of sand relative to gravel
D_g = 1. # Mean grain size in gravel fraction (> 2 mm)
D_s = 0.01 # Mean grain size in sand fraction (<= 2 mm)
-# Water source term [m/s]
-m_dot = 2.5e-8
+Q_terminus = 0.01/2. # Desired water flux at terminus [m^3/s]
+m_dot = Q_terminus/Ls # Water source term [m/s]
mu_w = 1.787e-3 # Water viscosity [Pa*s]
friction_factor = 0.1 # Darcy-Weisbach friction factor [-]
t@@ -88,7 +88,6 @@ 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 [-]
t@@ -126,70 +125,6 @@ def channel_shear_stress(Q, S):
return 1./8.*friction_factor*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)**(3./2)
-
- # Carter et al 2017
- # return K_e*v_s/alpha*(tau - tau_c).clip(min=0.) / \
- # (g*(rho_s - rho_w)*d15)**(3./2.)
-
- # Ng 2000
- return 0.092*(tau/(2.*(rho_s - rho_w)*g*d15))**(3./2.)
-
-
-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
-
- # Carter et al. 2017
- return K_d*v_s/alpha*c_bar[ix]*(g*(rho_s - rho_w)*d15/tau[ix])**0.5
-
-
-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*friction_factor))*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.
- 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*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_sediment_flux_sand(tau, W, f_s, D_s):
# Parker 1979, Wilcock 1997, 2001, Egholm 2013
# tau: Shear stress by water flow
t@@ -305,40 +240,40 @@ def flux_solver(m_dot, ds):
def pressure_solver(psi, f, Q, S):
# Iteratively find new water pressures
- # dP_c/ds = psi - f*rho_w*g*Q^2/S^{8/3} (Kingslake and Ng 2013)
+ # dN_c/ds = f*rho_w*g*Q^2/S^{8/3} - psi (Kingslake and Ng 2013)
it = 0
max_res = 1e9 # arbitrary large value
- while max_res > tol_P_c or it < Ns:
+ while max_res > tol_N_c or it < Ns:
- P_c_old = P_c.copy()
+ N_c_old = N_c.copy()
# P_downstream = P_upstream + dP
- # P_c[1:] = P_c[:-1] \
+ # N_c[1:] = N_c[:-1] \
# + psi[:-1]*ds[:-1] \
# - f[:-1]*rho_w*g*Q[:-1]**2./(S[:-1]**(8./3.))*ds[:-1] \
# Dirichlet BC (fixed pressure) at terminus
- P_c[-1] = 0.
+ N_c[-1] = 0.
# P_upstream = P_downstream - dP
- P_c[:-1] = P_c[1:] \
- - psi[:-1]*ds[:-1] \
- + f[:-1]*rho_w*g*Q[:-1]**2./(S[:-1]**(8./3.))*ds[:-1]
+ N_c[:-1] = N_c[1:] \
+ + psi[:-1]*ds[:-1] \
+ - f[:-1]*rho_w*g*Q[:-1]**2./(S[:-1]**(8./3.))*ds[:-1]
# + psi[:-1]*ds[:-1] \
# - f[:-1]*rho_w*g*Q[:-1]**2./(S[:-1]**(8./3.))*ds[:-1]
- max_res = numpy.max(numpy.abs((P_c - P_c_old)/(P_c + 1e-16)))
+ max_res = numpy.max(numpy.abs((N_c - N_c_old)/(N_c + 1e-16)))
if print_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')
+ 'Iterative solution not found for N_c')
it += 1
- return P_c
+ return N_c
def plot_state(step, time, S_, S_max_, title=True):
t@@ -351,7 +286,7 @@ def plot_state(step, time, S_, S_max_, title=True):
# ax_Pa.plot(s/1000., N/1000., '--r', label='$N$')
ax_Pa.plot(s_c/1000., N_c/1e6, '-k', label='$N$')
ax_Pa.plot(s_c/1000., H_c*rho_i*g/1e6, '--r', label='$P_i$')
- ax_Pa.plot(s_c/1000., P_c/1e6, ':y', label='$P_c$')
+ #ax_Pa.plot(s_c/1000., P_c/1e6, ':y', 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$')
t@@ -486,11 +421,8 @@ while time <= t_end:
# Find hydraulic roughness
f = channel_hydraulic_roughness(manning, S, W, theta)
- # 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
+ # Find new effective pressures consistent with the flow law
+ N_c = pressure_solver(psi, f, Q, S)
# Find new maximum normalized residual value
max_res = numpy.max(numpy.abs((S - S_prev_it)/(S + 1e-16)))