commit 9c1dfa711717600fb0bc43f6d6503294c5cddd9b
parent ea147780409ef9eacaa3e3cf29f1f10f4a8b8817
Author: Anders Damsgaard Christensen <adc@geo.au.dk>
Date: Mon, 30 Jan 2017 22:36:55 -0800
add simple 1d example, WIP
Diffstat:
2 files changed, 295 insertions(+), 0 deletions(-)
diff --git a/1d-test.py b/1d-test.py
@@ -0,0 +1,277 @@
+#!/usr/bin/env python
+
+## ABOUT THIS FILE
+# The following script uses basic Python and Numpy functionality to solve the
+# coupled systems of equations describing subglacial channel development in
+# soft beds as presented in `Damsgaard et al. "Sediment plasticity controls
+# channelization of subglacial meltwater in soft beds"`, submitted to Journal
+# of Glaciology.
+# High performance is not the goal for this implementation, which is instead
+# intended as a heavily annotated example on the solution procedure without
+# relying on solver libraries, suitable for low-level languages like C, Fortran
+# or CUDA.
+
+
+import numpy
+import matplotlib.pyplot as plt
+
+
+## Model parameters
+Ns = 10 # Number of nodes [-]
+Ls = 100e3 # Model length [m]
+#dt = 24.*60.*60. # Time step length [s]
+#dt = 1. # Time step length [s]
+dt = 60.*60.*24 # Time step length [s]
+#t_end = 24.*60.*60.*7. # Total simulation time [s]
+t_end = dt*4
+tol_Q = 1e-3 # Tolerance criteria for the normalized max. residual for Q
+tol_N_c = 1e-3 # Tolerance criteria for the normalized max. residual for N_c
+max_iter = 1e4 # Maximum number of solver iterations before failure
+output_convergence = True
+
+# Physical parameters
+rho_w = 1000. # Water density [kg/m^3]
+rho_i = 910. # Ice density [kg/m^3]
+rho_s = 2700. # Sediment density [kg/m^3]
+g = 9.8 # Gravitational acceleration [m/s^2]
+theta = 30. # Angle of internal friction in sediment [deg]
+
+# 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
+#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]
+#tau_c = 0.025*d15*g*(rho_s - rho_i) # Critical shear stress (Carter 2016)
+tau_c = 0.
+mu_w = 1.787e-3 # Water viscosity [Pa*s]
+froude = 0.1 # Friction factor [-]
+v_s = d15**2.*g*2.*(rho_s - rho_i)/(9.*mu_w) # Settling velocity (Carter 2016)
+
+# Hewitt 2011 channel flux parameters
+manning = 0.1 # Manning roughness coefficient [m^{-1/3} s]
+F = rho_w*g*manning*(2.*(numpy.pi + 2)**2./numpy.pi)**(2./3.)
+
+# Channel growth-limit parameters
+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
+
+
+
+## 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
+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]
+
+#m_dot = 7.93e-11 # Water source term [m/s]
+m_dot = 5.79e-9 # Water source term [m/s]
+#m_dot = 4.5e-8 # Water source term [m/s]
+
+# Initialize arrays for channel segments between nodes
+S = numpy.ones(len(s) - 1)*S_min # Cross-sectional area of channel segments[m^2]
+S_max = numpy.zeros_like(S) # Max. channel size [m^2]
+dSdt = numpy.empty_like(S) # Transient in channel cross-sectional area [m^2/s]
+W = S/numpy.tan(numpy.deg2rad(theta)) # Assuming no channel floor wedge
+Q = numpy.zeros_like(S) # Water flux in channel segments [m^3/s]
+Q_s = numpy.zeros_like(S) # Sediment flux in channel segments [m^3/s]
+N_c = numpy.zeros_like(S) # Effective pressure in channel segments [Pa]
+e_dot = numpy.zeros_like(S) # Sediment erosion rate in channel segments [m/s]
+d_dot = numpy.zeros_like(S) # Sediment deposition rate in chan. segments [m/s]
+c_bar = numpy.zeros_like(S) # Vertically integrated sediment content [m]
+tau = numpy.zeros_like(S) # Avg. shear stress from current [Pa]
+porosity = numpy.ones_like(S)*0.3 # Sediment porosity [-]
+res = numpy.zeros_like(S) # Solution residual during solver iterations
+
+
+## Helper functions
+def gradient(arr, arr_x):
+ # Central difference gradient of an array ``arr`` with node positions at
+ # ``arr_x``.
+ return (arr[:-1] - arr[1:])/(arr_x[:-1] - arr_x[1:])
+
+def avg_midpoint(arr):
+ # Averaged value of neighboring array elements
+ return (arr[:-1] + arr[1:])/2.
+
+def channel_water_flux(S, hydro_pot_grad):
+ # Hewitt 2011
+ return numpy.sqrt(1./F*S**(8./3.)*-hydro_pot_grad)
+
+def channel_shear_stress(Q, S):
+ # Weertman 1972, Walder and Fowler 1994
+ u_bar = Q*S
+ return 1./8.*froude*rho_w*u_bar**2.
+
+def channel_erosion_rate(tau):
+ # Parker 1979, Walder and Fowler 1994
+ return K_e*v_s*(tau - tau_c).clip(0.)/(g*(rho_s - rho_w)*d15)
+
+def channel_deposition_rate_kernel(tau, c_bar, ix):
+ # Parker 1979, Walder and Fowler 1994
+ return K_d*v_s*c_bar[ix]*(g*(rho_s - rho_w)*d15/tau[ix])**0.5
+
+def channel_deposition_rate(tau, c_bar, d_dot, Ns):
+ # Parker 1979, Walder and Fowler 1994
+ # Find deposition rate from upstream to downstream, margin at is=0
+ #for ix in numpy.arange(Ns-2, -1, -1):
+ for ix in numpy.arange(Ns - 1):
+ if ix == 0: # No sediment deposition at upstream end
+ c_bar[ix] = 0.
+ d_dot[ix] = 0.
+ else:
+ c_bar[ix] = (e_dot[ix - 1] - d_dot[ix - 1])*dt
+ 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
+
+def update_channel_size_with_limit(S, dSdt, dt, N):
+ # Damsgaard et al, in prep
+ S_max = ((c_1*N/1000. + c_2)*\
+ numpy.tan(numpy.deg2rad(theta))).clip(min=S_min)
+ S = numpy.minimum(S + dSdt*dt, S_max).clip(min=S_min)
+ W = S/numpy.tan(numpy.deg2rad(theta)) # Assume no channel floor wedge
+ return S, W, S_max
+
+def flux_and_pressure_solver(S):
+ # Iteratively find new fluxes and effective pressures in nested loops
+
+ it_Q = 0
+ max_res_Q = 1e9 # arbitrary large value
+ while max_res_Q > tol_Q:
+
+ Q_old = Q.copy()
+ # dQ/ds = m_dot -> Q_out = m*delta(s) - Q_in
+ # Propagate information along drainage direction (upwind)
+ Q[1:] = m_dot*ds[1:] - Q[:-1]
+ max_res_Q = numpy.max(numpy.abs((Q - Q_old)/(Q + 1e-16)))
+
+ if output_convergence:
+ print('it_Q = {}: max_res_Q = {}'.format(it_Q, max_res_Q))
+
+ it_N_c = 0
+ max_res_N_c = 1e9 # arbitrary large value
+ while max_res_N_c > tol_N_c:
+
+ N_c_old = N_c.copy()
+ # dN_c/ds = FQ^2/S^{8/3} - psi ->
+ #if it_N_c % 2 == 0: # Alternate direction between iterations
+ #N_c[1:] = F*Q[1:]**2./(S[1:]**(8./3.))*ds[1:] \
+ #- psi[1:]*ds[1:] + N_c[:-1] # Downstream
+ #else:
+ N_c[:-1] = -F*Q[:-1]**2./(S[:-1]**(8./3.))*ds[:-1] \
+ + psi[:-1]*ds[:-1] + N_c[1:] # Upstream
+
+ # Dirichlet BC at terminus
+ N_c[-1] = 0.
+
+ max_res_N_c = numpy.max(numpy.abs((N_c - N_c_old)/(N_c + 1e-16)))
+
+ if output_convergence:
+ print('it_N_c = {}: max_res_N_c = {}'.format(it_N_c,
+ max_res_N_c))
+
+ if it_N_c >= max_iter:
+ raise Exception('t = {}, step = {}:'.format(t, step) +
+ 'Iterative solution not found for N_c')
+ it_N_c += 1
+ #import ipdb; ipdb.set_trace()
+
+ #import ipdb; ipdb.set_trace()
+ if it_Q >= max_iter:
+ raise Exception('t = {}, step = {}:'.format(t, step) +
+ 'Iterative solution not found for Q')
+ it_Q += 1
+
+ return Q, N_c
+
+def plot_state(step):
+ # Plot parameters along profile
+ plt.plot(s_c/1000., S, '-k', label='$S$')
+ plt.plot(s_c/1000., S_max, '--k', label='$S_{max}$')
+ plt.plot(s_c/1000., Q, '-b', label='$Q$')
+ plt.plot(s_c/1000., N_c/1000., '--r', label='$N_c$')
+ #plt.plot(s, b, ':k', label='$b$')
+ #plt.plot(s, b, ':k', label='$b$')
+ plt.legend()
+ plt.xlabel('Distance from terminus [km]')
+ plt.tight_layout()
+ if step == -1:
+ plt.savefig('chan-0.init.pdf')
+ else:
+ plt.savefig('chan-' + str(step) + '.pdf')
+ plt.clf()
+
+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 effective pressure in channels [Pa]
+N_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', figsize=(3.3, 2.))
+plot_state(-1)
+
+#import ipdb; ipdb.set_trace()
+
+## Time loop
+t = 0.; step = 0
+while t < 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
+ 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)
+
+ # 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)
+
+ #import pdb; pdb.set_trace()
+ # Find new fluxes and effective pressures
+ Q, N_c = flux_and_pressure_solver(S)
+ # TODO: Q is zig zag
+
+ #import ipdb; ipdb.set_trace()
+
+ plot_state(step)
+
+ # Update time
+ t += dt
+ step += 1
+ #print(step)
+ #break
diff --git a/granular_channel_drainage/model.py b/granular_channel_drainage/model.py
@@ -14,6 +14,24 @@ class model:
'''
self.name = name
+ def genreateRegularGrid(self, Lx, Ly, Nx, Ny):
+ '''
+ Generate a uniform, regular and orthogonal grid using Landlab.
+
+ :param Lx: A tuple containing the length along x of the model
+ domain.
+ :type Lx: float
+ :param Ly: A tuple containing the length along y of the model
+ domain.
+ :type Ly: float
+ :param Nx: The number of random model nodes along ``x`` in the model.
+ :type Nx: int
+ :param Ny: The number of random model nodes along ``y`` in the model.
+ :type Ny: int
+ '''
+ self.grid_type = 'Regular'
+ self.grid = landlab.grid.RasterModelGrid(shape=(Nx, Ny), spacing=Lx/Nx)
+
def generateVoronoiDelaunayGrid(self, Lx, Ly, Nx, Ny,
structure='pseudorandom',
distribution='uniform'):
