commit 7fc6a939b43d9f7333a72bf298b6359747f1b5c7
parent 9c1dfa711717600fb0bc43f6d6503294c5cddd9b
Author: Anders Damsgaard Christensen <adc@geo.au.dk>
Date: Tue, 31 Jan 2017 21:56:11 -0800
transition from effective channel pressure to water pressure
Diffstat:
| M | 1d-test.py | | | 67 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++---------- |
1 file changed, 57 insertions(+), 10 deletions(-)
diff --git a/1d-test.py b/1d-test.py
@@ -21,11 +21,12 @@ 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]
+dt = 60.*60. # 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
+#tol_N_c = 1e-3 # Tolerance criteria for the normalized max. residual for N_c
+tol_P_c = 1e-3 # Tolerance criteria for the normalized max. residual for P_c
max_iter = 1e4 # Maximum number of solver iterations before failure
output_convergence = True
@@ -36,6 +37,11 @@ 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]
+# Water source term [m/s]
+#m_dot = 7.93e-11
+m_dot = 5.79e-9
+#m_dot = 4.5e-8
+
# 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
@@ -75,10 +81,6 @@ 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]
@@ -87,6 +89,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]
+P_c = numpy.zeros_like(S) # Water 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]
@@ -165,6 +168,7 @@ def flux_and_pressure_solver(S):
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:
@@ -192,6 +196,41 @@ def flux_and_pressure_solver(S):
'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:
@@ -199,14 +238,16 @@ def flux_and_pressure_solver(S):
'Iterative solution not found for Q')
it_Q += 1
- return Q, N_c
+ #return Q, N_c
+ return Q, P_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_c/1000., N_c/1000., '--r', label='$N_c$')
+ plt.plot(s_c/1000., P_c/1000., '--r', label='$P_c$')
#plt.plot(s, b, ':k', label='$b$')
#plt.plot(s, b, ':k', label='$b$')
plt.legend()
@@ -224,8 +265,10 @@ 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 effective pressure in channels [Pa]
+# 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)
@@ -263,9 +306,13 @@ while t < t_end:
#import pdb; pdb.set_trace()
# Find new fluxes and effective pressures
- Q, N_c = flux_and_pressure_solver(S)
+ #Q, N_c = flux_and_pressure_solver(S)
+ Q, P_c = flux_and_pressure_solver(S)
# TODO: Q is zig zag
+ # Find new effective pressure in channel segments
+ N_c = rho_i*g*H_c - P_c
+
#import ipdb; ipdb.set_trace()
plot_state(step)
