commit e1e23e617f46e30a99b7ef75ff8010c3813ab1b2
parent d2b05a1797b69a2d3cad693f89c9b9560e825ea2
Author: Anders Damsgaard Christensen <adc@geo.au.dk>
Date: Wed, 1 Feb 2017 13:53:51 -0800
show time info to stdout
Diffstat:
| M | 1d-test.py | | | 64 | +++++++++++++++++++++++++++++++++++++++------------------------- |
1 file changed, 39 insertions(+), 25 deletions(-)
diff --git a/1d-test.py b/1d-test.py
@@ -14,33 +14,29 @@
import numpy
import matplotlib.pyplot as plt
-
+import sys
## Model parameters
-Ns = 10 # Number of nodes [-]
+#Ns = 10 # Number of nodes [-]
+Ns = 25 # 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. # 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_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
+t_end = 24.*60.*60.*7. # 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
+output_convergence = False
# 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]
+rho_s = 2600. # 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-5
-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
@@ -64,7 +60,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-5
+S_min = 1e-2
@@ -75,7 +71,8 @@ ds = s[1:] - s[:-1]
# Ice thickness and bed topography
#H = 6.*(numpy.sqrt(Ls - s + 5e3) - numpy.sqrt(5e3)) + 1.0
-H = 0.6*(numpy.sqrt(Ls - s + 5e3) - numpy.sqrt(5e3)) + 1.0
+#H = 0.6*(numpy.sqrt(Ls - s + 5e3) - numpy.sqrt(5e3)) + 1.0
+H = 1.*(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]
@@ -223,7 +220,7 @@ def pressure_solver(psi, F, Q, S):
#import ipdb; ipdb.set_trace()
return P_c
-def plot_state(step):
+def plot_state(step, time):
# Plot parameters along profile
plt.plot(s_c/1000., S, '-k', label='$S$')
plt.plot(s_c/1000., S_max, '--k', label='$S_{max}$')
@@ -234,6 +231,7 @@ def plot_state(step):
#plt.plot(s, b, ':k', label='$b$')
plt.legend()
plt.xlabel('Distance from terminus [km]')
+ plt.title('Day: {:.2}'.format(time/(60.*60.*24.)))
plt.tight_layout()
if step == -1:
plt.savefig('chan-0.init.pdf')
@@ -241,6 +239,17 @@ def plot_state(step):
plt.savefig('chan-' + str(step) + '.pdf')
plt.clf()
+def find_new_timestep(ds, Q, S):
+ # Determine the timestep using the Courant-Friedrichs-Lewy condition
+ return numpy.minimum(60.*60.*24., numpy.min(numpy.abs(ds/(Q*S))))
+
+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.flush()
+
s_c = avg_midpoint(s) # Channel section midpoint coordinates [m]
## Initialization
@@ -261,13 +270,18 @@ 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)
+plot_state(-1, 0.0)
#import ipdb; ipdb.set_trace()
+
## Time loop
-t = 0.; step = 0
-while t < t_end:
+time = 0.; step = 0
+while time < t_end:
+
+ dt = find_new_timestep(ds, Q, S)
+
+ print_status_to_stdout(time, dt)
# Find average shear stress from water flux for each channel segment
tau = channel_shear_stress(Q, S)
@@ -299,12 +313,12 @@ while t < t_end:
#import ipdb; ipdb.set_trace()
- plot_state(step)
-
- #find_new_timestep(
+ plot_state(step, time)
# Update time
- t += dt
+ time += dt
step += 1
#print(step)
#break
+
+print('')
