granular-channel-hydro

commit 464ba42ad1e569c9c86c0b29e272ad9555b32feb
parent 400da52ca6dea9a381f902ba653ca1c267ac286f
Author: Anders Damsgaard Christensen <adc@geo.au.dk>
Date:   Thu,  2 Feb 2017 16:46:46 -0800

experiment with different parameter values for erosion/deposition

Diffstat:
M
1d-channel.py
 | 
21
++++++++++++---------

1 file changed, 12 insertions(+), 9 deletions(-)
diff --git a/1d-channel.py b/1d-channel.py
@@ -23,7 +23,7 @@ import sys
 # # Model parameters
 Ns = 25                # Number of nodes [-]
 Ls = 100e3             # Model length [m]
-t_end = 24.*60.*60.*2  # Total simulation time [s]
+t_end = 24.*60.*60.*365  # 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
@@ -44,8 +44,9 @@ m_dot = 4.5e-7
 
 # 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
-alpha = 2e5  # Geometric correction factor (Carter et al 2017)
+# K_d = 6.0   # Deposition constant [-], disabled when 0.0
+K_d = 0.1   # Deposition constant [-], disabled when 0.0
+alpha = 1e5  # 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
 d15 = 1e-3  # Characteristic grain size [m]
@@ -85,7 +86,7 @@ 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.empty_like(S)   # Transient in channel cross-sect. area [m^2/s]
+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]
@@ -124,9 +125,10 @@ def channel_shear_stress(Q, S):
 
 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)**(3./2)
     # Carter et al 2017
-    return K_e*v_s/alpha*(tau - tau_c).clip(min=0.)/(g*(rho_s - rho_w)*d15)
+    return K_e*v_s/alpha*(tau - tau_c).clip(min=0.) / \
+        (g*(rho_s - rho_w)*d15)**(3./2.)
 
 
 def channel_deposition_rate_kernel(tau, c_bar, ix):
@@ -200,6 +202,7 @@ 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_max[:] = 1000.
     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
@@ -391,11 +394,11 @@ while time <= t_end:
     # Find new effective pressure in channel segments
     N_c = rho_i*g*H_c - P_c
 
-    plot_state(step, time)
+    if step % 10 == 0:
+        plot_state(step, time)
 
     # import ipdb; ipdb.set_trace()
-    #if step > 0:
-        #break
+    # if step > 0: break
 
     # Update time
     time += dt