tvarious fixes, advection instability present - granular-channel-hydro - subglacial hydrology model for sedimentary channels
(HTM) git clone git://src.adamsgaard.dk/granular-channel-hydro
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) LICENSE
---
(DIR) commit cdd2e02befacebb825dd8718658f44c28fe27e82
(DIR) parent b927f7b0194424fbefffb922391b0d999459702a
(HTM) Author: Anders Damsgaard Christensen <adc@geo.au.dk>
Date: Wed, 1 Feb 2017 16:41:32 -0800
various fixes, advection instability present
Diffstat:
M 1d-test.py | 46 ++++++++++---------------------
1 file changed, 14 insertions(+), 32 deletions(-)
---
(DIR) diff --git a/1d-test.py b/1d-test.py
t@@ -39,8 +39,7 @@ m_dot = 5.79e-5
# Walder and Fowler 1994 sediment transport parameters
K_e = 0.1 # Erosion constant [-], disabled when 0.0
-#K_d = 6. # Deposition constant [-], disabled when 0.0
-K_d = 0.0 # Deposition constant [-], disabled when 0.0
+K_d = 6.0 # Deposition constant [-], disabled when 0.0
#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]
t@@ -62,21 +61,20 @@ c_2 = 4.60 # [m]
S_min = 1e-1
-
## Initialize model arrays
# Node positions, terminus at Ls
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
+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 = 1.*(numpy.sqrt(Ls - s + 5e3) - numpy.sqrt(5e3)) + 1.0
b = numpy.zeros_like(H)
N = H*0.1*rho_i*g # Initial effective stress [Pa]
-p_w = rho_i*g*H - N # Water pressure [Pa], here at floatation
-hydro_pot = rho_w*g*b + p_w # Hydraulic potential [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]
# Initialize arrays for channel segments between nodes
S = numpy.ones(len(s) - 1)*S_min # Cross-sectional area of channel segments[m^2]
t@@ -131,12 +129,10 @@ def channel_deposition_rate(tau, c_bar, d_dot, Ns):
c_bar[ix] = 0.
d_dot[ix] = 0.
else:
- c_bar[ix] = (e_dot[ix - 1] - d_dot[ix - 1])*dt
+ c_bar[ix] = (e_dot[ix - 1] - d_dot[ix - 1])*dt # Net erosion upstr.
d_dot[ix] = channel_deposition_rate_kernel(tau, c_bar, ix)
-
return d_dot, c_bar
-
def channel_growth_rate(e_dot, d_dot, porosity, W):
# Damsgaard et al, in prep
return (e_dot - d_dot)/porosity*W
t@@ -206,9 +202,7 @@ def pressure_solver(psi, F, Q, S):
raise Exception('t = {}, step = {}:'.format(time, step) +
'Iterative solution not found for P_c')
it += 1
- #import ipdb; ipdb.set_trace()
- #import ipdb; ipdb.set_trace()
return P_c
def plot_state(step, time):
t@@ -219,7 +213,7 @@ def plot_state(step, time):
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_m3s = ax_Pa.twinx() # axis with meters as y-axis unit
+ 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.)))
t@@ -229,8 +223,10 @@ def plot_state(step, time):
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_m.plot(s_c/1000., S, '-k', label='$S$')
+ #ax_m.plot(s_c/1000., S_max, '--k', 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.plot(s_c/1000., e_dot, '--r', label='$\dot{e}$')
t@@ -256,7 +252,8 @@ def find_new_timestep(ds, Q, S):
dt = safety*numpy.minimum(60.*60.*24., numpy.min(numpy.abs(ds/(Q*S))))
if dt < 1.0:
- raise Exception('Error: Time step less than 1 second at time '
+ raise Exception('Error: Time step less than 1 second at step '
+ + '{}, time '.format(step)
+ '{:.3} s/{:.3} d'.format(time, time/(60.*60.*24.)))
return dt
t@@ -268,28 +265,23 @@ def print_status_to_stdout(time, dt):
s_c = avg_midpoint(s) # Channel section midpoint coordinates [m]
-## Initialization
# Find gradient in hydraulic potential between the nodes
hydro_pot_grad = gradient(hydro_pot, s)
# Find field values at the middle of channel segments
N_c = avg_midpoint(N)
-#P_c = avg_midpoint(P)
H_c = avg_midpoint(N)
# Find fluxes in channel segments [m^3/s]
Q = channel_water_flux(S, hydro_pot_grad)
# Water-pressure gradient from geometry [Pa/m]
-#psi = -rho_i*g*gradient(H, s) - (rho_w - rho_i)*g*gradient(b, s)
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)
-#import ipdb; ipdb.set_trace()
-
## Time loop
time = 0.; step = 0
t@@ -302,12 +294,9 @@ while time <= t_end:
# Find average shear stress from water flux for each channel segment
tau = channel_shear_stress(Q, S)
- # Find erosion rates for each channel segment
+ # Find sediment erosion and deposition rates for each channel segment
e_dot = channel_erosion_rate(tau)
- # TODO: erosion law smooth for now with tau_c = 0.
d_dot, c_bar = channel_deposition_rate(tau, c_bar, d_dot, Ns)
- # TODO: d_dot and c_bar values unreasonably high
- # Deposition disabled for now with K_d = 0.
# Determine change in channel size for each channel segment
dSdt = channel_growth_rate(e_dot, d_dot, porosity, W)
t@@ -320,21 +309,14 @@ while time <= t_end:
# meltwater production (m_dot)
Q = flux_solver(m_dot, ds)
- #import pdb; pdb.set_trace()
# 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
- #import ipdb; ipdb.set_trace()
-
plot_state(step, time)
# Update time
time += dt
step += 1
- #print(step)
- #break
-
-print('')