granular-channel-hydro

commit 25bddd0fd670ac412295cf87cdd4ed1759875e12
parent 63f99a6d5a5fdb264a19d29281a8b8a373358632
Author: Anders Damsgaard <andersd@riseup.net>
Date:   Wed,  8 Mar 2017 20:52:39 -0800

change geometry and upstream flux BC

Diffstat:
M
1d-channel.py
 | 
12
+++++++-----

1 file changed, 7 insertions(+), 5 deletions(-)
diff --git a/1d-channel.py b/1d-channel.py
@@ -23,7 +23,7 @@ import sys
 # # Model parameters
 Ns = 25             # Number of nodes [-]
 # Ls = 100e3        # Model length [m]
-Ls = 1e3            # Model length [m]
+Ls = 10e3           # 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]
@@ -47,7 +47,9 @@ D_s = 0.01     # Mean grain size in sand fraction (<= 2 mm)
 
 # Boundary conditions
 P_terminus = 0.      # Water pressure at terminus [Pa]
-m_dot = 1.0e-5  # Water source term [m/s]
+m_dot = 2.0e-7  # Water source term [m/s]
+Q_upstream = 1e-5  # Water influx, upstream end (must be larger than 0) [m^3/s]
+
 
 # Channel hydraulic properties
 manning = 0.1  # Manning roughness coefficient [m^{-1/3} s]
@@ -69,7 +71,7 @@ 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)) + 10.0  # glacier
+H = 6.*(numpy.sqrt(Ls - s + 5e3) - numpy.sqrt(5e3)) + 1.0  # glacier
 # slope = 0.1  # Surface slope [%]
 # H = 1000. + -slope/100.*s
 
@@ -219,7 +221,7 @@ def flux_solver(m_dot, ds):
         Q_old = Q.copy()
         # dQ/ds = m_dot  ->  Q_out = m*delta(s) + Q_in
         # Upwind information propagation (upwind)
-        Q[0] = 1e-2  # Ng 2000
+        Q[0] = Q_upstream
         Q[1:] = m_dot*ds[1:] + Q[:-1]
         max_res = numpy.max(numpy.abs((Q - Q_old)/(Q + 1e-16)))
 
@@ -312,7 +314,7 @@ def plot_state(step, time, S_, S_max_, title=True):
     ax_m2s = ax_m2.twinx()
     # ax_ms.plot(s_c/1000., e_dot, '--r', label='$\dot{e}$')
     # ax_ms.plot(s_c/1000., d_dot, ':b', label='$\dot{d}$')
-    ax_m2s.plot(s_c/1000., dSdt, ':', label='$\partial S/\partial t$')
+    ax_m2s.plot(s_c/1000., dSdt, ':', label='$dS/dt$')
 
     ax_m2.legend(loc=2)
     ax_m2s.legend(loc=3)