tuse upwind differences for pressure and a larger S_min - 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 b927f7b0194424fbefffb922391b0d999459702a
(DIR) parent 3a8080534dbb60ccc47435d15ed6ce9fd909fc59
(HTM) Author: Anders Damsgaard Christensen <adc@geo.au.dk>
Date: Wed, 1 Feb 2017 16:05:32 -0800
use upwind differences for pressure and a larger S_min
Diffstat:
M 1d-test.py | 89 +++++++++++++++----------------
1 file changed, 42 insertions(+), 47 deletions(-)
---
(DIR) diff --git a/1d-test.py b/1d-test.py
t@@ -19,7 +19,7 @@ import sys
## Model parameters
Ns = 25 # Number of nodes [-]
Ls = 100e3 # Model length [m]
-t_end = 24.*60.*60.*1.9 # Total simulation time [s]
+t_end = 24.*60.*60.*2 # Total simulation time [s]
tol_Q = 1e-3 # Tolerance criteria for the normalized max. residual for Q
tol_P_c = 1e-3 # Tolerance criteria for the normalized max. residual for P_c
max_iter = 1e2*Ns # Maximum number of solver iterations before failure
t@@ -44,8 +44,8 @@ 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.
+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)
t@@ -59,7 +59,7 @@ 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-1
t@@ -151,69 +151,61 @@ def update_channel_size_with_limit(S, dSdt, dt, N):
def flux_solver(m_dot, ds):
# Iteratively find new fluxes
- it_Q = 0
- max_res_Q = 1e9 # arbitrary large value
+ it = 0
+ max_res = 1e9 # arbitrary large value
# Iteratively find solution, do not settle for less iterations than the
# number of nodes
- while max_res_Q > tol_Q and it_Q < Ns:
+ while max_res > tol_Q or it < Ns:
Q_old = Q.copy()
# dQ/ds = m_dot -> Q_out = m*delta(s) + Q_in
# Upwind information propagation (upwind)
- #Q[1:] = m_dot*ds[1:] - Q[:-1]
- #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)))
+ max_res = 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))
+ print('it = {}: max_res = {}'.format(it, max_res))
#import ipdb; ipdb.set_trace()
- if it_Q >= max_iter:
- raise Exception('t = {}, step = {}:'.format(t, step) +
+ if it >= max_iter:
+ raise Exception('t = {}, step = {}:'.format(time, step) +
'Iterative solution not found for Q')
- it_Q += 1
+ it += 1
return Q
def pressure_solver(psi, F, Q, S):
# Iteratively find new water pressures
+ # dP_c/ds = psi - FQ^2/S^{8/3}
- 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:
+ it = 0
+ max_res = 1e9 # arbitrary large value
+ while max_res > tol_P_c or it < Ns*40:
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
+
+ # Upwind finite differences
+ P_c[:-1] = -psi[:-1]*ds[:-1] \
+ + F*Q[:-1]**2./(S[:-1]**(8./3.))*ds[:-1] \
+ + P_c[1:] # Upstream
+
+ # Dirichlet BC (fixed pressure) at terminus
P_c[-1] = 0.
- max_res_P_c = numpy.max(numpy.abs((P_c - P_c_old)/(P_c + 1e-16)))
+ # von Neumann BC (no gradient = no flux) at s=0
+ P_c[0] = P_c[1]
+
+ max_res = 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))
+ print('it = {}: max_res = {}'.format(it, max_res))
- if it_P_c >= max_iter:
- raise Exception('t = {}, step = {}:'.format(t, step) +
+ if it >= max_iter:
+ raise Exception('t = {}, step = {}:'.format(time, step) +
'Iterative solution not found for P_c')
- it_P_c += 1
+ it += 1
#import ipdb; ipdb.set_trace()
#import ipdb; ipdb.set_trace()
t@@ -250,7 +242,6 @@ def plot_state(step, time):
ax_m.set_ylabel('[m]')
ax_ms.set_ylabel('[m/s]')
-
plt.setp(ax_Pa.get_xticklabels(), visible=False)
plt.tight_layout()
if step == -1:
t@@ -261,14 +252,18 @@ def plot_state(step, time):
def find_new_timestep(ds, Q, S):
# Determine the timestep using the Courant-Friedrichs-Lewy condition
- safety = 0.8
- return safety*numpy.minimum(60.*60.*24., numpy.min(numpy.abs(ds/(Q*S))))
+ safety = 0.2
+ 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 '
+ + '{:.3} s/{:.3} d'.format(time, time/(60.*60.*24.)))
+
+ return dt
def print_status_to_stdout(time, dt):
- sys.stdout.write('\rt = {:.2} s or {:.4} d, dt = {:.2} s'.format(\
- time,
- time/(60.*60.*24.),
- dt))
+ sys.stdout.write('\rt = {:.2} s or {:.4} d, dt = {:.2} s '\
+ .format(time, time/(60.*60.*24.), dt))
sys.stdout.flush()
s_c = avg_midpoint(s) # Channel section midpoint coordinates [m]
t@@ -298,7 +293,7 @@ plot_state(-1, 0.0)
## Time loop
time = 0.; step = 0
-while time < t_end:
+while time <= t_end:
dt = find_new_timestep(ds, Q, S)