In [3]:

    

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np


def avalanche(x, z0, dzdx_max=1., max_iter=1000, method=1):
    
    # TODO:
    # - 2D
    # - varying grid size
    
    z = z0.copy()

    for n in range(max_iter):
        
        dzs = np.zeros(z.shape)
        dzi = np.zeros(z.shape)

        if method == 1:
            dzdx = np.zeros(z.shape)

            dz = np.diff(z)
            dx = np.diff(x)
            dzdx[:-1] = dz / dx

            ixd = abs(dzdx) > abs(dzdx_max)
            ixe = np.roll(ixd, 1)

            dzi[ixd] = np.sign(dzdx[ixd]) * (abs(dzdx[ixd]) - abs(dzdx_max)) / 2.
            dzs[ixe] -= dzi[ixd]
            dzs[ixd] += dzi[ixd]

        elif method == 2:
            for i in range(len(x)-1):
                dz = z[i+1] - z[i]
                dx = x[i+1] - x[i]
                dzdx = dz / dx

                if abs(dzdx) > abs(dzdx_max):

                    dzi = np.sign(dzdx) * (abs(dzdx) - abs(dzdx_max)) / 2.
                    dzs[i+1] -= dzi
                    dzs[i] += dzi
        
        else:
            ValueError('Unknown method: %d' % method)
        
        z += dzs
        
        yield z

        if max(abs(dzs)) < 1e-3:
            break
            
            
def plot_profiles(x, z0, z, dzdx_max=1.):
    
    e = int(np.ceil(max(z0.max()/dzdx_max, x.max())))
    
    fig, axs = plt.subplots(figsize=(10,4))
    for x0 in range(-e, e, 10):
        axs.plot(x0+np.arange(e), dzdx_max * np.arange(e), ':k', alpha=.5)
        axs.plot(x0+np.arange(e), -dzdx_max * np.arange(e) + dzdx_max * e - 1, ':k', alpha=.5)

    axs.plot(x, z0, label='initial profile (%0.2f $\mathdefault{m^3/m}$)' % np.trapz(z0, x))
    axs.plot(x, z, label='final profile (%0.2f $\mathdefault{m^3/m}$)' % np.trapz(z, x))

    axs.set_xlim((0, x.max()))
    axs.set_ylim((0, z0.max()))
    axs.set_xlabel('horizontal position [m]')
    axs.set_ylabel('bed level [m]')
    axs.legend(loc='upper right')

    return fig, axs


def plot_avalanching(x, z0, dzdx_max=1.):
    profiles = []
    for z in avalanche(x, z0, dzdx_max=dzdx_max):
        profiles.append(z)
    fig, axs = plot_profiles(x, z0, profiles[-1], dzdx_max=dzdx_max)
    axs.set_title('iterations: %d' % len(profiles))
    for z in profiles:
        axs.plot(x, z, '-k', alpha=.1, zorder=-1)
        
    return fig, axs


    

In [4]:

    

x = np.arange(101)
dzdx_max = 1.

z0 = np.zeros(x.shape)

plot_avalanching(x, z0, dzdx_max=dzdx_max)

z0 = np.zeros(x.shape)
z0[41:51] = np.linspace(0, 5, 10)
z0[60:50:-1] = np.linspace(0, 5, 10)

plot_avalanching(x, z0, dzdx_max=dzdx_max)

z0 = np.zeros(x.shape)
z0[41:51] = np.linspace(0, 15, 10)
z0[60:50:-1] = np.linspace(0, 15, 10)

plot_avalanching(x, z0, dzdx_max=dzdx_max)

z0 = np.zeros(x.shape)
z0[50] = 2

plot_avalanching(x, z0, dzdx_max=dzdx_max)

z0 = np.zeros(x.shape)
z0[50] = 200

plot_avalanching(x, z0, dzdx_max=dzdx_max)

z0 = np.zeros(x.shape)
z0[10:20] = np.linspace(0, 30, 10)
z0[25:35] = np.linspace(0, 20, 10)
z0[40:50] = np.linspace(0, 10, 10)
z0[55:65] = np.linspace(0, 5, 10)
z0[90:80:-1] = np.linspace(0, 30, 10)
z0[75:65:-1] = np.linspace(0, 20, 10)
z0[60:50:-1] = np.linspace(0, 10, 10)
z0[45:35:-1] = np.linspace(0, 5, 10)

plot_avalanching(x, z0, dzdx_max=dzdx_max)
plot_avalanching(x, z0, dzdx_max=5)
plot_avalanching(x, z0, dzdx_max=.5)
plot_avalanching(x, z0, dzdx_max=.1)
plot_avalanching(x, z0, dzdx_max=.01)


    

    




/Users/hoonhout/Python/VirtualEnvs/aeolis/lib/python3.6/site-packages/matplotlib/axes/_base.py:3193: UserWarning: Attempting to set identical bottom==top results
in singular transformations; automatically expanding.
bottom=0, top=0.0
  'bottom=%s, top=%s') % (bottom, top))






    
Out[4]:





(<matplotlib.figure.Figure at 0x105b91320>,
 <matplotlib.axes._subplots.AxesSubplot at 0x105b91b00>)






    













    













    













    













    













    













    













    













    













    

In [ ]: