commit ed3e36d656e883a7527ac866dad7e6d1ebffcf57
parent 25bddd0fd670ac412295cf87cdd4ed1759875e12
Author: Anders Damsgaard <andersd@riseup.net>
Date: Mon, 13 Mar 2017 15:59:48 -0700
solve for N_c
Diffstat:
| M | 1d-channel.py | | | 125 | +++++++++++++++++++++++++++++++++---------------------------------------------- |
1 file changed, 52 insertions(+), 73 deletions(-)
diff --git a/1d-channel.py b/1d-channel.py
@@ -21,43 +21,43 @@ import sys
# # Model parameters
-Ns = 25 # Number of nodes [-]
-# Ls = 100e3 # Model length [m]
-Ls = 10e3 # Model length [m]
-# Ls = 1e3 # Model length [m]
-total_days = 60. # Total simulation time [d]
+Ns = 25 # Number of nodes [-]
+Ls = 10e3 # 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
-max_iter = 1e2*Ns # Maximum number of solver iterations before failure
-print_output_convergence = False # Display convergence statistics during run
+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
+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_interval = 1 # Time steps between plots
+plot_interval = 20 # Time steps between plots
+speedup_factor = 1. # Speed up channel growth to reach steady state faster
# Physical parameters
-rho_w = 1000. # Water density [kg/m^3]
-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]
+rho_w = 1000. # Water density [kg/m^3]
+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 = 1. # Mean grain size in gravel fraction (> 2 mm)
-D_s = 0.01 # Mean grain size in sand fraction (<= 2 mm)
+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
# Boundary conditions
P_terminus = 0. # Water pressure at terminus [Pa]
-m_dot = 2.0e-7 # Water source term [m/s]
-Q_upstream = 1e-5 # Water influx, upstream end (must be larger than 0) [m^3/s]
-
+m_dot = 1e-6 # Water source term [m/s]
+Q_upstream = 1e-5 # Water influx upstream (must be larger than 0) [m^3/s]
# Channel hydraulic properties
-manning = 0.1 # Manning roughness coefficient [m^{-1/3} s]
+manning = 0.1 # Manning roughness coefficient [m^{-1/3} s]
friction_factor = 0.1 # Darcy-Weisbach friction factor [-]
# Channel growth-limit parameters
-c_1 = -0.118 # [m/kPa]
-c_2 = 4.60 # [m]
+c_1 = -0.118 # [m/kPa]
+c_2 = 4.60 # [m]
# Minimum channel size [m^2], must be bigger than 0
S_min = 1e-2
@@ -70,16 +70,15 @@ S_min = 1e-2
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 # glacier
+# 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)
-N = H*0.1*rho_i*g # Initial effective stress [Pa]
-#p_w = rho_i*g*H - N # Initial guess of water pressure [Pa]
-#hydro_pot = rho_w*g*b + p_w # Initial guess of hydraulic potential [Pa]
+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]
@@ -149,8 +148,6 @@ def channel_sediment_flux_sand(tau, W, f_s, D_s):
(1.0 - 0.846*numpy.sqrt(ref_shear_stress/shields_stress))
)**4.5
- # The above relation gives 'nan' values for low values of tau
-
return Q_c
@@ -239,40 +236,33 @@ 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()
-
- # P_downstream = P_upstream + 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_old = N_c.copy()
# Dirichlet BC (fixed pressure) at terminus
- P_c[-1] = P_terminus
+ N_c[-1] = rho_i*g*H_c[-1] - P_terminus
- # 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]
- # + 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]*numpy.abs(Q[:-1]) \
+ /(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):
@@ -281,8 +271,6 @@ def plot_state(step, time, S_, S_max_, title=True):
fig.set_size_inches(3.3*1.1, 3.3*1.1*1.5)
ax_Pa = plt.subplot(3, 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/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$')
@@ -294,7 +282,6 @@ def plot_state(step, time, S_, S_max_, title=True):
plt.title('Day: {:.3}'.format(time/(60.*60.*24.)))
ax_Pa.legend(loc=2)
ax_m3s.legend(loc=4)
- # ax_Pa.set_ylabel('[kPa]')
ax_Pa.set_ylabel('[MPa]')
ax_m3s.set_ylabel('[m$^3$/s]')
@@ -308,12 +295,8 @@ def plot_state(step, time, S_, S_max_, title=True):
ax_m2 = plt.subplot(3, 1, 3, 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_m2s = 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_m2s.plot(s_c/1000., dSdt, ':', label='$dS/dt$')
ax_m2.legend(loc=2)
@@ -345,17 +328,15 @@ def find_new_timestep(ds, Q, S):
return dt
-def print_status_to_stdout(time, dt):
- sys.stdout.write('\rt = {:.2} s or {:.4} d, dt = {:.2} s '
+def print_status_to_stdout(step, time, dt):
+ sys.stdout.write('\rstep = {}, '.format(step) +
+ 't = {:.2} s or {:.4} d, dt = {:.2} s '
.format(time, time/(60.*60.*24.), dt))
sys.stdout.flush()
s_c = avg_midpoint(s) # Channel section midpoint coordinates [m]
H_c = avg_midpoint(H)
-# Find water flux
-Q = flux_solver(m_dot, ds)
-
# Water-pressure gradient from geometry [Pa/m]
psi = -rho_i*g*gradient(H, s) - (rho_w - rho_i)*g*gradient(b, s)
@@ -373,7 +354,7 @@ while time <= t_end:
dt = find_new_timestep(ds, Q, S)
# Display current simulation status
- print_status_to_stdout(time, dt)
+ print_status_to_stdout(step, time, dt)
it = 0
@@ -389,6 +370,10 @@ while time <= t_end:
S_prev_it = S.copy()
+ # Find new water fluxes consistent with mass conservation and local
+ # meltwater production (m_dot)
+ Q = flux_solver(m_dot, ds)
+
# Find average shear stress from water flux for each channel segment
tau = channel_shear_stress(Q, S)
@@ -400,29 +385,26 @@ 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_t, porosity, s_c)
+ #dSdt[1:] = channel_growth_rate_sedflux(Q_t, porosity, s_c)
+ #dSdt *= speedup_factor
# Update channel cross-sectional area and width according to growth
# rate and size limit for each channel segment
S, W, S_max, dSdt = \
update_channel_size_with_limit(S, S_old, dSdt, dt, N_c)
- # Find new water fluxes consistent with mass conservation and local
- # meltwater production (m_dot)
- Q = flux_solver(m_dot, ds)
-
# 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)
+ N_c = pressure_solver(psi, f, Q, S)
# Find new effective pressure in channel segments
- N_c = rho_i*g*H_c - P_c
+ P_c = rho_i*g*H_c - N_c
# Find new maximum normalized residual value
max_res = numpy.max(numpy.abs((S - S_prev_it)/(S + 1e-16)))
- if print_output_convergence:
+ if print_output_convergence_main:
print('it = {}: max_res = {}'.format(it, max_res))
# import ipdb; ipdb.set_trace()
@@ -435,9 +417,6 @@ while time <= t_end:
if step % plot_interval == 0:
plot_state(step, time, S, S_max)
- # import ipdb; ipdb.set_trace()
- # if step > 0: break
-
# Update time
time += dt
step += 1
