commit 63f99a6d5a5fdb264a19d29281a8b8a373358632
parent d0ec40ca4dde58229f1c40d4fc0d3789a3d679a2
Author: Anders Damsgaard <andersd@riseup.net>
Date: Wed, 8 Mar 2017 19:42:12 -0800
solve for water pressure
Diffstat:
| M | 1d-channel.py | | | 66 | ++++++++++++++++++++++++++++++------------------------------------ |
1 file changed, 30 insertions(+), 36 deletions(-)
diff --git a/1d-channel.py b/1d-channel.py
@@ -24,10 +24,11 @@ import sys
Ns = 25 # Number of nodes [-]
# Ls = 100e3 # Model length [m]
Ls = 1e3 # Model length [m]
+# Ls = 1e3 # 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_N_c = 1e-3 # Tolerance criteria for the norm. max. residual for N_c
+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
safety = 0.01 # Safety factor ]0;1] for adaptive timestepping
@@ -44,15 +45,13 @@ 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)
-Q_terminus = 0.01/2. # Desired water flux at terminus [m^3/s]
-m_dot = Q_terminus/Ls # Water source term [m/s]
-
-mu_w = 1.787e-3 # Water viscosity [Pa*s]
-friction_factor = 0.1 # Darcy-Weisbach friction factor [-]
+# Boundary conditions
+P_terminus = 0. # Water pressure at terminus [Pa]
+m_dot = 1.0e-5 # Water source term [m/s]
# Channel hydraulic properties
manning = 0.1 # Manning roughness coefficient [m^{-1/3} s]
-#F = rho_w*g*manning*(2.*(numpy.pi + 2)**2./numpy.pi)**(2./3.)
+friction_factor = 0.1 # Darcy-Weisbach friction factor [-]
# Channel growth-limit parameters
c_1 = -0.118 # [m/kPa]
@@ -70,15 +69,15 @@ 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
+H = 6.*(numpy.sqrt(Ls - s + 5e3) - numpy.sqrt(5e3)) + 10.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]
-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]
+#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-sect. area of channel segments[m^2]
@@ -88,9 +87,7 @@ 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 concentration [-]
+P_c = numpy.zeros_like(S) # Water pressure in channel segments [Pa]
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
@@ -240,40 +237,40 @@ def flux_solver(m_dot, ds):
def pressure_solver(psi, f, Q, S):
# Iteratively find new water pressures
- # dN_c/ds = f*rho_w*g*Q^2/S^{8/3} - psi (Kingslake and Ng 2013)
+ # dP_c/ds = psi - f*rho_w*g*Q^2/S^{8/3} (Kingslake and Ng 2013)
it = 0
max_res = 1e9 # arbitrary large value
- while max_res > tol_N_c or it < Ns:
+ while max_res > tol_P_c or it < Ns:
- N_c_old = N_c.copy()
+ P_c_old = P_c.copy()
# P_downstream = P_upstream + dP
- # N_c[1:] = N_c[:-1] \
+ # P_c[1:] = P_c[:-1] \
# + psi[:-1]*ds[:-1] \
# - f[:-1]*rho_w*g*Q[:-1]**2./(S[:-1]**(8./3.))*ds[:-1] \
# Dirichlet BC (fixed pressure) at terminus
- N_c[-1] = 0.
+ P_c[-1] = P_terminus
# P_upstream = P_downstream - dP
- N_c[:-1] = N_c[1:] \
- + psi[:-1]*ds[:-1] \
- - f[:-1]*rho_w*g*Q[:-1]**2./(S[:-1]**(8./3.))*ds[:-1]
+ 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]
- max_res = numpy.max(numpy.abs((N_c - N_c_old)/(N_c + 1e-16)))
+ max_res = numpy.max(numpy.abs((P_c - P_c_old)/(P_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 N_c')
+ 'Iterative solution not found for P_c')
it += 1
- return N_c
+ return P_c
def plot_state(step, time, S_, S_max_, title=True):
@@ -286,7 +283,7 @@ def plot_state(step, time, S_, S_max_, title=True):
# 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$')
+ ax_Pa.plot(s_c/1000., P_c/1e6, ':y', 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$')
@@ -300,9 +297,9 @@ def plot_state(step, time, S_, S_max_, title=True):
ax_m3s.set_ylabel('[m$^3$/s]')
ax_m3s_sed = plt.subplot(3, 1, 2, sharex=ax_Pa)
- ax_m3s_sed.plot(s_c/1000., Q_g, ':', label='$Q_g$')
- ax_m3s_sed.plot(s_c/1000., Q_s, '-', label='$Q_s$')
- ax_m3s_sed.plot(s_c/1000., Q_t, '--', label='$Q_t$')
+ ax_m3s_sed.plot(s_c/1000., Q_g, ':', label='$Q_{gravel}$')
+ ax_m3s_sed.plot(s_c/1000., Q_s, '-', label='$Q_{sand}$')
+ ax_m3s_sed.plot(s_c/1000., Q_t, '--', label='$Q_{total}$')
ax_m3s_sed.set_ylabel('[m$^3$/s]')
ax_m3s_sed.legend(loc=2)
@@ -352,12 +349,6 @@ def print_status_to_stdout(time, dt):
sys.stdout.flush()
s_c = avg_midpoint(s) # Channel section midpoint coordinates [m]
-
-# 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)
H_c = avg_midpoint(H)
# Find water flux
@@ -421,8 +412,11 @@ while time <= t_end:
# Find hydraulic roughness
f = channel_hydraulic_roughness(manning, S, W, theta)
- # Find new effective pressures consistent with the flow law
- N_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 maximum normalized residual value
max_res = numpy.max(numpy.abs((S - S_prev_it)/(S + 1e-16)))
