commit 97a5cad4554cc24b354750d75d98d7e661bab5ca
parent 709930aac7cbe7585ab6f956aed5f370b41d7166
Author: Anders Damsgaard <andersd@riseup.net>
Date: Tue, 14 Feb 2017 13:16:36 -0800
use an iterative solution for channel dynamics
Diffstat:
| M | 1d-channel.py | | | 84 | +++++++++++++++++++++++++++++++++++++++++++++++++++---------------------------- |
1 file changed, 54 insertions(+), 30 deletions(-)
diff --git a/1d-channel.py b/1d-channel.py
@@ -39,8 +39,8 @@ 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-6
+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
@@ -198,14 +198,14 @@ def channel_deposition_rate(tau, c_bar, d_dot, Ns):
def channel_growth_rate(e_dot, d_dot, porosity, W):
# Damsgaard et al, in prep
- return (e_dot - d_dot)/porosity*W
+ return (e_dot - d_dot)*W
-def update_channel_size_with_limit(S, dSdt, dt, N):
+def update_channel_size_with_limit(S, S_old, dSdt, dt, N):
# Damsgaard et al, in prep
- S_max = numpy.maximum((c_1*numpy.maximum(N, 0.)/1000. + c_2) *
+ S_max = numpy.maximum((c_1*numpy.maximum(N, 0.)/1000. + c_2)**2./4. *
numpy.tan(numpy.deg2rad(theta)), S_min)
- S = numpy.maximum(numpy.minimum(S + dSdt*dt, S_max), S_min)
+ S = numpy.maximum(numpy.minimum(S_old + dSdt*dt, S_max), S_min)
W = S/numpy.tan(numpy.deg2rad(theta)) # Assume no channel floor wedge
return S, W, S_max
@@ -285,16 +285,19 @@ def plot_state(step, time, S_, S_max_, title=True):
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/1000., N/1000., '--r', label='$N$')
+ #ax_Pa.plot(s/1000., N/1000., '--r', label='$N$')
+ ax_Pa.plot(s/1000., H*rho_i*g/1e6, '--r', label='$P_i$')
+ ax_Pa.plot(s_c/1000., P_c/1e6, ':b', 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$')
+ ax_m3s.plot(s_c/1000., Q, '-k', label='$Q$')
if title:
plt.title('Day: {:.3}'.format(time/(60.*60.*24.)))
ax_Pa.legend(loc=2)
ax_m3s.legend(loc=3)
- ax_Pa.set_ylabel('[kPa]')
+ #ax_Pa.set_ylabel('[kPa]')
+ ax_Pa.set_ylabel('[MPa]')
ax_m3s.set_ylabel('[m$^3$/s]')
ax_m2 = plt.subplot(2, 1, 2, sharex=ax_Pa)
@@ -372,35 +375,56 @@ while time <= t_end:
print_status_to_stdout(time, dt)
- # Find average shear stress from water flux for each channel segment
- tau = channel_shear_stress(Q, S)
+ it = 0
+ max_res = 1e9 # arbitrary large value
+
+ S_old = S.copy()
+ # Iteratively find solution, do not settle for less iterations than the
+ # number of nodes
+ while max_res > tol_Q or it < Ns:
+
+ S_prev_it = S.copy()
+
+ # Find average shear stress from water flux for each channel segment
+ tau = channel_shear_stress(Q, S)
+
+ # Find sediment erosion and deposition rates for each channel segment
+ e_dot = channel_erosion_rate(tau)
+ d_dot, c_bar = channel_deposition_rate(tau, c_bar, d_dot, Ns)
- # Find sediment erosion and deposition rates for each channel segment
- e_dot = channel_erosion_rate(tau)
- d_dot, c_bar = channel_deposition_rate(tau, c_bar, d_dot, Ns)
+ # Determine change in channel size for each channel segment
+ dSdt = channel_growth_rate(e_dot, d_dot, porosity, W)
- # Determine change in channel size for each channel segment
- dSdt = channel_growth_rate(e_dot, d_dot, porosity, W)
+ # Update channel cross-sectional area and width according to growth rate
+ # and size limit for each channel segment
+ S, W, S_max = update_channel_size_with_limit(S, S_old, dSdt, dt, N_c)
- # Update channel cross-sectional area and width according to growth rate
- # and size limit for each channel segment
- S, W, S_max = update_channel_size_with_limit(S, 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 new water fluxes consistent with mass conservation and local
- # meltwater production (m_dot)
- Q = flux_solver(m_dot, ds)
+ # Find the corresponding sediment flux
+ # Q_b = bedload_sediment_flux(
+ Q_s = suspended_sediment_flux(c_bar, Q, S)
- # Find the corresponding sediment flux
- # Q_b = bedload_sediment_flux(
- Q_s = suspended_sediment_flux(c_bar, Q, S)
+ # Find new water pressures consistent with the flow law
+ P_c = pressure_solver(psi, F, Q, S)
- # 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 pressure in channel segments
- N_c = rho_i*g*H_c - P_c
+ # Find new maximum normalized residual value
+ max_res = numpy.max(numpy.abs((S - S_prev_it)/(S + 1e-16)))
+ if output_convergence:
+ print('it = {}: max_res = {}'.format(it, max_res))
+
+ # import ipdb; ipdb.set_trace()
+ if it >= max_iter:
+ raise Exception('t = {}, step = {}:'.format(time, step) +
+ 'Iterative solution not found for Q')
+ it += 1
- if step + 1 % 10 == 0:
+ if step % 10 == 0:
plot_state(step, time, S, S_max)
# import ipdb; ipdb.set_trace()
