twater flux debugging in progress - 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 518b45d9299a64bda252848ed1d6ee2c8ba71aad
 (DIR) parent 7fc6a939b43d9f7333a72bf298b6359747f1b5c7
 (HTM) Author: Anders Damsgaard Christensen <adc@geo.au.dk>
       Date:   Tue, 31 Jan 2017 22:14:00 -0800
       
       water flux debugging in progress
       
       Diffstat:
         M 1d-test.py                          |      11 +++++++----
       
       1 file changed, 7 insertions(+), 4 deletions(-)
       ---
 (DIR) diff --git a/1d-test.py b/1d-test.py
       t@@ -74,7 +74,8 @@ 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)) + 1.0
       +#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
        b = numpy.zeros_like(H)
        
        N = H*0.1*rho_i*g            # Initial effective stress [Pa]
       t@@ -160,9 +161,11 @@ def flux_and_pressure_solver(S):
            while max_res_Q > tol_Q:
        
                Q_old = Q.copy()
       -        # dQ/ds = m_dot  ->  Q_out = m*delta(s) - Q_in
       +        # dQ/ds = m_dot  ->  Q_out = m*delta(s) + Q_in
                # Propagate information along drainage direction (upwind)
       -        Q[1:] = m_dot*ds[1:] - Q[:-1]
       +        #Q[1:] = m_dot*ds[1:] - Q[:-1]
       +        Q[0] = 0.
       +        Q[1:] = m_dot*ds[1:] + Q[:-1]
                max_res_Q = numpy.max(numpy.abs((Q - Q_old)/(Q + 1e-16)))
        
                if output_convergence:
       t@@ -281,7 +284,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