tdo not solve flux and pressure in nested loops, fix ds - 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 d2b05a1797b69a2d3cad693f89c9b9560e825ea2
(DIR) parent 518b45d9299a64bda252848ed1d6ee2c8ba71aad
(HTM) Author: Anders Damsgaard Christensen <adc@geo.au.dk>
Date: Wed, 1 Feb 2017 12:29:37 -0800
do not solve flux and pressure in nested loops, fix ds
Diffstat:
M 1d-test.py | 145 ++++++++++++++-----------------
1 file changed, 64 insertions(+), 81 deletions(-)
---
(DIR) diff --git a/1d-test.py b/1d-test.py
t@@ -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 = 5.79e-9
-#m_dot = 4.5e-8
+#m_dot = 5.79e-5
+m_dot = 4.5e-8
# Walder and Fowler 1994 sediment transport parameters
K_e = 0.1 # Erosion constant [-], disabled when 0.0
t@@ -64,14 +64,14 @@ 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
+S_min = 1e-5
## Initialize model arrays
# Node positions, terminus at Ls
s = numpy.linspace(0., Ls, Ns)
-ds = s[:-1] - s[1:]
+ds = s[1:] - s[:-1]
# Ice thickness and bed topography
#H = 6.*(numpy.sqrt(Ls - s + 5e3) - numpy.sqrt(5e3)) + 1.0
t@@ -153,96 +153,75 @@ def update_channel_size_with_limit(S, dSdt, dt, N):
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
-
+def flux_solver(m_dot, ds):
+ # Iteratively find new fluxes
it_Q = 0
max_res_Q = 1e9 # arbitrary large value
- while max_res_Q > tol_Q:
+
+ # Iteratively find solution, do not settle for less iterations than the
+ # number of nodes
+ while max_res_Q > tol_Q and it_Q < Ns:
Q_old = Q.copy()
# dQ/ds = m_dot -> Q_out = m*delta(s) + Q_in
- # Propagate information along drainage direction (upwind)
+ # Upwind information propagation (upwind)
#Q[1:] = m_dot*ds[1:] - Q[:-1]
- Q[0] = 0.
+ #Q[0] = 0.
+ Q[0] = 1e-2 # Ng 2000
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()
- '''
- it_P_c = 0
- max_res_P_c = 1e9 # arbitrary large value
- while max_res_P_c > tol_P_c:
-
- P_c_old = P_c.copy()
- # dP_c/ds = psi - FQ^2/S^{8/3}
- #if it_P_c % 2 == 0: # Alternate direction between iterations
- #P_c[1:] = psi[1:]*ds[1:] \
- #- F*Q[1:]**2./(S[1:]**(8./3.))*ds[1:] \
- #+ P_c[:-1] # Downstream
- #else:
-
- #P_c[:-1] = -psi[:-1]*ds[:-1] \
- #+ F*Q[:-1]**2./(S[:-1]**(8./3.))*ds[:-1] \
- #+ P_c[1:] # Upstream
-
- P_c[:-1] = -avg_midpoint(psi)*ds[:-1] \
- + F*avg_midpoint(Q)**2./(S[:-1]**(8./3.))*ds[:-1] \
- + P_c[1:] # Upstream
-
- # Dirichlet BC at terminus
- P_c[-1] = 0.
-
- max_res_P_c = numpy.max(numpy.abs((P_c - P_c_old)/(P_c + 1e-16)))
-
- if output_convergence:
- print('it_P_c = {}: max_res_P_c = {}'.format(it_P_c,
- max_res_P_c))
-
- if it_P_c >= max_iter:
- raise Exception('t = {}, step = {}:'.format(t, step) +
- 'Iterative solution not found for P_c')
- it_P_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
- return Q, P_c
+ return Q
+
+def pressure_solver(psi, F, Q, S):
+ # Iteratively find new water pressures
+
+ it_P_c = 0
+ max_res_P_c = 1e9 # arbitrary large value
+ while max_res_P_c > tol_P_c and it_P_c < Ns:
+
+ P_c_old = P_c.copy()
+ # dP_c/ds = psi - FQ^2/S^{8/3}
+ #if it_P_c % 2 == 0: # Alternate direction between iterations
+ #P_c[1:] = psi[1:]*ds[1:] \
+ #- F*Q[1:]**2./(S[1:]**(8./3.))*ds[1:] \
+ #+ P_c[:-1] # Downstream
+ #else:
+
+ #P_c[:-1] = -psi[:-1]*ds[:-1] \
+ #+ F*Q[:-1]**2./(S[:-1]**(8./3.))*ds[:-1] \
+ #+ P_c[1:] # Upstream
+
+ P_c[:-1] = -avg_midpoint(psi)*ds[:-1] \
+ + F*avg_midpoint(Q)**2./(S[:-1]**(8./3.))*ds[:-1] \
+ + P_c[1:] # Upstream
+
+ # Dirichlet BC at terminus
+ P_c[-1] = 0.
+
+ max_res_P_c = numpy.max(numpy.abs((P_c - P_c_old)/(P_c + 1e-16)))
+
+ if output_convergence:
+ print('it_P_c = {}: max_res_P_c = {}'.format(it_P_c,
+ max_res_P_c))
+
+ if it_P_c >= max_iter:
+ raise Exception('t = {}, step = {}:'.format(t, step) +
+ 'Iterative solution not found for P_c')
+ it_P_c += 1
+ #import ipdb; ipdb.set_trace()
+
+ #import ipdb; ipdb.set_trace()
+ return P_c
def plot_state(step):
# Plot parameters along profile
t@@ -284,7 +263,7 @@ psi = -rho_i*g*gradient(H, s) - (rho_w - rho_i)*g*gradient(b, s)
fig = plt.figure('channel', figsize=(3.3, 2.))
plot_state(-1)
-import ipdb; ipdb.set_trace()
+#import ipdb; ipdb.set_trace()
## Time loop
t = 0.; step = 0
t@@ -307,11 +286,13 @@ while t < t_end:
# 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)
+
#import pdb; pdb.set_trace()
- # Find new fluxes and effective pressures
- #Q, N_c = flux_and_pressure_solver(S)
- Q, P_c = flux_and_pressure_solver(S)
- # TODO: Q is zig zag
+ # 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
t@@ -320,6 +301,8 @@ while t < t_end:
plot_state(step)
+ #find_new_timestep(
+
# Update time
t += dt
step += 1