commit 633970e66a0eb2d60b96eab40e8a9f891201ecf0
parent 6afcc6591b23c355ae2cae50953afbba1de9e0d6
Author: Anders Damsgaard <andersd@riseup.net>
Date: Mon, 15 May 2017 16:27:58 -0400
fix figure and tweak parameters
Diffstat:
1 file changed, 12 insertions(+), 17 deletions(-)
diff --git a/1d-channel.py b/1d-channel.py
@@ -23,19 +23,17 @@ import sys
# # Model parameters
Ns = 25 # Number of nodes [-]
Ls = 1e3 # Model length [m]
-total_days = 60. # Total simulation time [d]
+total_days = 2. # 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
+tol_S = 1e-2 # Tolerance criteria for the norm. max. residual for S
+tol_Q = 1e-2 # Tolerance criteria for the norm. max. residual for Q
+tol_N_c = 1e-2 # 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 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]
@@ -48,8 +46,8 @@ 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 = 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]
+m_dot = numpy.linspace(1e-6, 1e-5, Ns-1) # Water source term [m/s]
+Q_upstream = 1e-3 # Water influx upstream (must be larger than 0) [m^3/s]
# Channel hydraulic properties
manning = 0.1 # Manning roughness coefficient [m^{-1/3} s]
@@ -70,8 +68,6 @@ ds = s[1:] - s[:-1]
# Ice thickness [m]
H = 6.*(numpy.sqrt(Ls - s + 5e3) - numpy.sqrt(5e3)) + 1.0
-# slope = 0.1 # Surface slope [%]
-# H = 1000. + -slope/100.*s
# Bed topography [m]
b = numpy.zeros_like(H)
@@ -205,10 +201,9 @@ 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):
+def plot_state(step, time, S_, S_max_, title=False):
# Plot parameters along profile
fig = plt.gcf()
fig.set_size_inches(3.3*1.1, 3.3*1.1*1.5)
@@ -224,11 +219,14 @@ def plot_state(step, time, S_, S_max_, title=True):
if title:
plt.title('Day: {:.3}'.format(time/(60.*60.*24.)))
ax_Pa.legend(loc=2)
- ax_m3s.legend(loc=4)
+ ax_m3s.legend(loc=1)
ax_Pa.set_ylabel('[MPa]')
ax_m3s.set_ylabel('[m$^3$/s]')
- ax_m3s_sed = plt.subplot(3, 1, 2, sharex=ax_Pa)
+ # ax_m3s_sed = plt.subplot(3, 1, 2, sharex=ax_Pa)
+ ax_m3s_sed_blank = plt.subplot(3, 1, 2, sharex=ax_Pa)
+ ax_m3s_sed_blank.get_yaxis().set_visible(False)
+ ax_m3s_sed = ax_m3s_sed_blank.twinx()
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)
@@ -327,14 +325,11 @@ while time <= t_end:
# 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_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
f = channel_hydraulic_roughness(manning, S, W, theta)
