commit 709930aac7cbe7585ab6f956aed5f370b41d7166
parent 464ba42ad1e569c9c86c0b29e272ad9555b32feb
Author: Anders Damsgaard Christensen <adc@geo.au.dk>
Date: Thu, 2 Feb 2017 21:47:16 -0800
working example
Diffstat:
| M | 1d-channel.py | | | 60 | +++++++++++++++++++++++++++++++++--------------------------- |
1 file changed, 33 insertions(+), 27 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.*365 # Total simulation time [s]
+t_end = 24.*60.*60.*30. # 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
@@ -39,13 +39,15 @@ theta = 30. # Angle of internal friction in sediment [deg]
# Water source term [m/s]
# m_dot = 7.93e-11
-m_dot = 4.5e-7
+# m_dot = 4.5e-7
+m_dot = 4.5e-6
# m_dot = 5.79e-5
+# m_dot = 1.8/(1000.*365.*24.*60.*60.) # Whillan's melt rate from Joughin 2004
# Walder and Fowler 1994 sediment transport parameters
-K_e = 0.1 # Erosion constant [-], disabled when 0.0
+K_e = 6.0 # Erosion constant [-], disabled when 0.0
# K_d = 6.0 # Deposition constant [-], disabled when 0.0
-K_d = 0.1 # Deposition constant [-], disabled when 0.0
+K_d = K_e # 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
@@ -74,9 +76,10 @@ 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 = 6.*(numpy.sqrt(Ls - s + 5e3) - numpy.sqrt(5e3)) + 1.0 # glacier
+slope = 0.1 # Surface slope [%]
+H = 1000. + -slope/100.*s
+
b = numpy.zeros_like(H)
N = H*0.1*rho_i*g # Initial effective stress [Pa]
@@ -200,10 +203,9 @@ def channel_growth_rate(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.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)
+ S_max = numpy.maximum((c_1*numpy.maximum(N, 0.)/1000. + c_2) *
+ numpy.tan(numpy.deg2rad(theta)), S_min)
+ S = numpy.maximum(numpy.minimum(S + dSdt*dt, S_max), S_min)
W = S/numpy.tan(numpy.deg2rad(theta)) # Assume no channel floor wedge
return S, W, S_max
@@ -276,39 +278,43 @@ def pressure_solver(psi, F, Q, S):
return P_c
-def plot_state(step, time):
+def plot_state(step, time, S_, S_max_, title=True):
# Plot parameters along profile
fig = plt.gcf()
- fig.set_size_inches(3.3, 3.3)
+ fig.set_size_inches(3.3*1.1, 3.3*1.1)
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_Pa.plot(s_c/1000., P_c/1000., '--r', label='$P_c$')
+ ax_Pa.plot(s/1000., N/1000., '--r', label='$N$')
ax_m3s = ax_Pa.twinx() # axis with m3/s as y-axis unit
ax_m3s.plot(s_c/1000., Q, '-b', label='$Q$')
- plt.title('Day: {:.3}'.format(time/(60.*60.*24.)))
+ if title:
+ plt.title('Day: {:.3}'.format(time/(60.*60.*24.)))
ax_Pa.legend(loc=2)
- ax_m3s.legend(loc=1)
+ ax_m3s.legend(loc=3)
ax_Pa.set_ylabel('[kPa]')
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_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_, '--', color='#666666', 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 = ax_m2.twinx()
ax_ms.plot(s_c/1000., e_dot, '--r', label='$\dot{e}$')
ax_ms.plot(s_c/1000., d_dot, ':b', label='$\dot{d}$')
- ax_m.legend(loc=2)
- ax_ms.legend(loc=1)
- ax_m.set_xlabel('$s$ [km]')
- ax_m.set_ylabel('[m]')
+ ax_m2.legend(loc=2)
+ ax_ms.legend(loc=3)
+ ax_m2.set_xlabel('$s$ [km]')
+ ax_m2.set_ylabel('[m$^2$]')
ax_ms.set_ylabel('[m/s]')
+ ax_Pa.set_xlim([s.min()/1000., s.max()/1000.])
+
plt.setp(ax_Pa.get_xticklabels(), visible=False)
plt.tight_layout()
if step == -1:
@@ -352,7 +358,7 @@ 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)
+plot_state(-1, 0.0, S, S_max)
# # Time loop
@@ -394,8 +400,8 @@ while time <= t_end:
# Find new effective pressure in channel segments
N_c = rho_i*g*H_c - P_c
- if step % 10 == 0:
- plot_state(step, time)
+ if step + 1 % 10 == 0:
+ plot_state(step, time, S, S_max)
# import ipdb; ipdb.set_trace()
# if step > 0: break
