# License
https://raw.githubusercontent.com/computational-sediment-hyd/1DBedVariationMixed/master/LICENSE

# import module


```python
import numpy as np
import pandas as pd
import numba
import sys
import os
import json
import glob as glob
```


```python
screenclass = np.array( [ 2.0**0, 2.0**1, 2.0**2, 2.0**3, 2.0**4, 2.0**5, 2.0**6] )
dsize = np.array( [ 10**(0.5*(np.log10(screenclass[i]) + np.log10(screenclass[i+1]))) for i in range(len(screenclass)-1) ] )

h=1.0
ib = 1.0/200.0
hini = 1.0
n = 0.03
v = ib**0.5 * h**(2/3) / n
fr = v/np.sqrt(9.8*h)
print('Froude Number:' + str(fr))

columns = ['diameter[mm]' ]
df = pd.DataFrame(dsize, columns=columns)
df['tau_star'] = h*ib/dsize/1.65*1000
df
```

    Froude Number:0.7529232524210427





<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>diameter[mm]</th>
      <th>tau_star</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>1.414214</td>
      <td>2.142748</td>
    </tr>
    <tr>
      <th>1</th>
      <td>2.828427</td>
      <td>1.071374</td>
    </tr>
    <tr>
      <th>2</th>
      <td>5.656854</td>
      <td>0.535687</td>
    </tr>
    <tr>
      <th>3</th>
      <td>11.313708</td>
      <td>0.267843</td>
    </tr>
    <tr>
      <th>4</th>
      <td>22.627417</td>
      <td>0.133922</td>
    </tr>
    <tr>
      <th>5</th>
      <td>45.254834</td>
      <td>0.066961</td>
    </tr>
  </tbody>
</table>
</div>



# Governing Equation of River flow 

\begin{align}
    \frac{\partial A}{\partial t}+\frac{\partial Q}{\partial x} = 0 \\
\end{align}

\begin{align}
    \frac{\partial Q}{\partial t}+\frac{\partial}{\partial x}\left(\frac{ Q^2}{ A}\right) + gA\frac{\partial H}{\partial x}+gAi_e = 0
\end{align}


```python
@numba.jit(nopython=True, parallel=False)
def flowmodel(A, Q, Adown, Qup, dAb, zb, B, dx, dt, g, manning):
    imax = len(A)
    Anew, Qnew = np.zeros(imax), np.zeros(imax)
    
# continuous equation
    for i in numba.prange(1, imax-1) : Anew[i] = A[i] - dt * ( Q[i] - Q[i-1] ) / dx
        
    Anew[0], Anew[-1] = Anew[1], Adown
    
# moumentum equation
    for i in numba.prange(1,imax-1): 
        ip, ic, im = (i+1, i, i-1) 
        Cr1 = 0.5*( Q[ic]/A[ic] + Q[ip]/A[ip] )*dt/dx
        Cr2 = 0.5*( Q[ic]/A[ic] + Q[im]/A[im] )*dt/dx
        dHdx1 = ( Anew[ip]/B[ip] + zb[ip] + dAb[ip]/B[ip] - Anew[ic]/B[ic] - zb[ic] - dAb[ic]/B[ic] ) / dx
        dHdx2 = ( Anew[ic]/B[ic] + zb[ic] + dAb[ic]/B[ic] - Anew[im]/B[im] - zb[im] - dAb[im]/B[im] ) / dx
        dHdx = (1.0 - Cr1) * dHdx1 + Cr2 * dHdx2
        
        Qnew[ic] = Q[ic] - dt * ( Q[ic]**2/A[ic] - Q[im]**2/A[im] ) / dx  \
                   - dt * g * Anew[ic] * dHdx \
                   - dt * g * A[ic] * manning**2 * Q[ic]**2 / B[ic]**2 / ( A[ic]/B[ic] )**(10.0/3.0)
                
    Qnew[0], Qnew[-1] = Qup, Qnew[-2]
            
    return Anew, Qnew 
```

# Governing Equation of River Bed : Hirano model


\begin{align}
    (1-\lambda)\frac{\partial A_{b}}{\partial t}+\frac{\partial }{\partial x} \sum_{i=1}^n ( Q_{bi}P_i)= 0 \\
\end{align}

\begin{align}
    \frac{\partial P_i}{\partial t} &= - \frac{1}{(1-\lambda)E_dB} \frac{\partial (Q_{bi}P_i)}{\partial x} - \frac{P_{si}}{E_d B}\frac{\partial A_b}{\partial t} \\
    &= - \frac{1}{E_d B}\left(\frac{\partial A_{bi}}{\partial t} + P_{si}\frac{\partial A_b}{\partial t}\right) 
\end{align}

\begin{align}
&P_{si} = 
\left\{ \begin{array}{ll}
    P_i \mbox{ in exchange layer} & \left(\dfrac{\partial A_b}{\partial t} > 0 \right) \\
    P_i \mbox{ under exchange layer} & \left(\dfrac{\partial A_b}{\partial t} < 0 \right) \\
\end{array} \right. \\
\end{align}


$Q_{bi}$ : Ashida-Michiue Eq. and  Egiazaroff Eq.

$E_d$ : thickness of exchange layer

\begin{align}
    \frac{\tau_{*c i}}{\tau_{*cm}} &= 0.85 \frac{D_m}{D_i} \qquad & \left(\frac{D_i}{D_m} < 0.4 \right) \\
    \dfrac{\tau_{*ci}}{\tau_{*cm}} &= \left( \dfrac{ \displaystyle \log_e 19 }{ \displaystyle \log_e \left(19\dfrac{D_i}{D_m} \right) } \right)^2 \qquad & \left(\frac{D_i}{D_m} \geq 0.4 \right)
\end{align}


```python
@numba.jit(nopython=True, parallel=False)
def sedimentmodel(dAb, dratio, A, Q, B, dx, dt, g, manning, dsize, dratioStandard, hExlayer, Qbup):
    rhosw = 1.65 # grain specific gravity in water
    porosity = 0.4
    tscAve = 0.05 # critical tractive force of average grain size
    
    dAbnew = np.zeros_like( dAb )
    drationew = np.zeros_like( dratio )
    Qbsub = np.zeros_like( dratio )
    dAbsub = np.zeros_like( dratio )
    
    imax, lmax = len(dratio), len(dratio[0])
    
    for i in numba.prange(imax):
        Ap, Qp, dr, Bp = A[i], Q[i], dratio[i], B[i]
        dAve = np.sum(dsize * dr)
        us = np.sqrt(g * Ap/Bp) * Qp/Ap * manning / (Ap/Bp)**(2/3)
        tsAve = us**2.0/rhosw/g/dAve
        use = Qp/Ap / ( 6.0 + 2.5 * np.log( Ap/Bp/dAve/( 1.0+2.0*tsAve) ) )
        Kc=1.0
        
        for l in range(lmax):
            dri, dsi = dr[l], dsize[l]
            ts  = us**2.0 /rhosw/g/dsi
            tse = use**2.0/rhosw/g/dsi
            
# Egiazaroff Eq. 
            x = dsi/dAve
            tscbytscm = (np.log(19)/np.log(19*x))**2 if x > 0.4 else 0.85/x
            tsc = Kc * tscbytscm * tscAve 
       
            if ts > tsc :
# Ashida-Michiue Eq. 
                Qbsub[i][l] = np.sqrt(rhosw * g * dsi**3.0) \
                        * 17.0 * tse**1.5 * ( 1.0 - tsc / ts) * ( 1.0 - np.sqrt(tsc / ts) ) \
                        * dri * B[i]
            else:
                Qbsub[i][l] = 0.0
                
    for i in numba.prange(imax):
        Qbs = Qbsub[i]
        if i == 0:
            if np.min(Qbup) < 0.0 : # equilibrium sand supply
                dAbsub[i][:] = 0.0
            else:
                for l in range(lmax):
                    dAbsub[i][l] = - dt / (1.0 - porosity) * ( Qbsub[i][l]-Qbup[l] )/dx
                    if np.abs(dAbsub[i][l]) > hExlayer :
                        print(i) ; print('dzsub-error')
        else:
            for l in range(lmax):
                dAbsub[i][l] = - dt / (1.0 - porosity) * ( Qbsub[i][l] - Qbsub[i-1][l] )/dx
                if np.abs(dAbsub[i][l]) > hExlayer : 
                    print(i) ; print('dzsub-error')
                
    for i in numba.prange(imax):
        Qbs = Qbsub[i]
# update dAb 
        dAball = np.sum( dAbsub[i] )
        if np.abs(dAball) > hExlayer:
            print(i) ; print('dz-error')
            
        dAbnew[i] = dAb[i] + dAball
        
# update dratio
        dratioIn = dratioStandard[i][:] if dAball < 0.0 else dratio[i][:]
            
        for l in range(lmax):
            drationew[i][l] = dratio[i][l] + ( dAbsub[i][l] - dAball * dratioIn[l] ) /hExlayer/B[i]
            if drationew[i][l] < 0.0 : drationew[i][l] = 0.0
            
# correct ration so that sum of ration become 100% 
        sumdr = np.sum( drationew[i] )
        drationew[i] /= sumdr
        
    return dAbnew, drationew, Qbsub
```

# main


```python
def bedvariation(
dx,dt,manning,totalTime,outTimeStep,RunUpTime
,dsize ,dratioStandard ,dratio ,hExlayer ,A ,Q ,B ,zb ,dAb ,Qup ,Adown 
,outputfilename
):
    
    g = 9.8
    Qbup = np.full( ( len(dsize) ), -9999.0 ) # when equilibrium sand supply, set qbup to minus value
        
# run-up calculation
    for i in range(int(RunUpTime/dt)):
        ib = ( ( dAb[-2]/B[-2] + zb[-2] ) - ( dAb[-1]/B[-1] + zb[-1] ) )/dx
        Adownp = Adown(0.0, Q[-1], dAb[-1]/B[-1], ib)
        Qupp = Qup(0.0)
        A, Q = flowmodel(A, Q, Adownp, Qupp, dAb, zb, B, dx, dt, g, manning)
        
# main calculation
    for i in range( int(totalTime/dt) ):
        
# cal bed variation
        dAb, dratio, Qbsub = sedimentmodel(dAb, dratio, A, Q, B, dx, dt, g, manning, dsize, dratioStandard, hExlayer, Qbup)
        
# cal water profile
        ib = ( ( dAb[-2]/B[-2] + zb[-2] ) - ( dAb[-1]/B[-1] + zb[-1] ) )/dx
        Adownp = Adown(i*dt, Q[-1], dAb[-1]/B[-1], ib)
        Qupp = Qup(i*dt)
        A, Q = flowmodel(A, Q, Adownp, Qupp, dAb, zb, B, dx, dt, g, manning)
        
# output 
        if( int(i*dt) % int(outTimeStep) ) == 0 :
            print( str( int(i*dt) ) + ' second') 
            profile = []
            for ii, (ap, qp, z, r, q) in enumerate(zip(A, Q, dAb, dratio, Qbsub)) : 
                profile.append({'distance' : int(ii*dx), 'A':ap, 'Q':qp, 'dAb':z, 'ratio':list(r), 'Qb':list(q)})
            json.dump( {'time':int(i*dt), 'profile' : profile}, open('%010d' % int(i*dt) + 'sec.json', 'w') )
            
# join json files
    d = [ json.load( open(f, 'r') ) for f in glob.glob('*sec.json') ]
# delete json files
    for f in glob.glob('*sec.json') : os.remove(f)
        
    d.sort(key=lambda x: x['time'])
    
    cond = {'width':list(B), 'elevation':list(zb), 'screenclass':list(screenclass), 'dsize':list(dsize) }
    json.dump( {'condition':cond, 'output':d}, open(outputfilename, 'w') )
```

# Sample
## case1


```python
%%time
length = 10000.0
dx = 100.0
imax = int(length/dx) + 1
dt = 10.0
totalTime = 1000.1*3600.0
outTimeStep = 100.0*3600.0
RunUpTime = 3.0 * 3600.0
hini = 1.0
manning  = 0.03
ib = 1.0/200.0
ib2 = 1.0/1000.0
outputfilename = 'case01.json'

# grain diameter classification
screenclass = np.array( [ 2.0**0, 2.0**1, 2.0**2, 2.0**3, 2.0**4, 2.0**5, 2.0**6], dtype=float )/1000
dsize = np.array( [ 10**(0.5*(np.log10(screenclass[i]) + np.log10(screenclass[i+1]))) for i in range(len(screenclass)-1) ], dtype=float )

# percentage of grain size under exchange layer
dratioStandard1 = np.full_like(dsize, 1.0/len(dsize), dtype=float)
dratioStandard = np.full( (imax, len(dsize) ), dratioStandard1, dtype=float)

# initial percentage of grain size in exchange layer
dratio = np.copy(dratioStandard)

# thickness of exchange layer 
hExlayer = dsize[-1]

# Initial & Boundary condition
B = np.full(imax, 1.0, dtype=float)
A = hini*B
Q = ib**0.5*(hini)**(5.0/3.0)/manning*B #normal flow
zb = np.zeros(imax)
for i in range(1,imax):
    zb[i] = zb[i-1] + ib2*dx if i < 50 else zb[i-1] + ib*dx
zb = zb[::-1]

dAb = np.zeros(imax)

Qup = Q[0]

def Adown(time, Q, dzb, ib):
    return ( manning**2*Q**2/ib/B[-1]**2 )**0.3 * B[-1]

def Qup(time):
    return Q[0]


bedvariation(
dx,dt,manning,totalTime,outTimeStep,RunUpTime
,dsize ,dratioStandard ,dratio ,hExlayer ,A ,Q ,B ,zb ,dAb ,Qup ,Adown
,outputfilename
    )
```

    0 second
    360000 second
    720000 second
    1080000 second
    1440000 second
    1800000 second
    2160000 second
    2520000 second
    2880000 second
    3240000 second
    3600000 second
    Wall time: 37.2 s


## case2


```python
%%time
length = 10000.0
dx = 100.0
imax = int(length/dx) + 1
dt = 10.0
totalTime = 3000.1*3600.0
outTimeStep = 300.0*3600.0
RunUpTime = 3.0 * 3600.0
hini = 1.0
manning  = 0.03
ib = 1.0/200.0
ib2 = 1.0/1000.0
outputfilename = 'case02.json'

# grain diameter classification
screenclass = np.array( [ 2.0**0, 2.0**1, 2.0**2, 2.0**3, 2.0**4, 2.0**5, 2.0**6], dtype=float )/1000
dsize = np.array( [ 10**(0.5*(np.log10(screenclass[i]) + np.log10(screenclass[i+1]))) for i in range(len(screenclass)-1) ], dtype=float )

# percentage of grain size under exchange layer
dratioStandard1 = np.full_like(dsize, 1.0/len(dsize), dtype=float)
dratioStandard2 = np.zeros(len(dsize))
dratioStandard2[-1] = 1.0

dratioStandard = np.full( (imax, len(dsize) ), dratioStandard1, dtype=float)
dratioStandard[:50] = dratioStandard2

# initial percentage of grain size in exchange layer
#     dratio = np.full( (imax, len(dsize) ), dratioStandard1, dtype=float)
dratio = np.copy(dratioStandard)

# thickness of exchange layer 
hExlayer = dsize[-1]

# Initial & Boundary condition
B = np.full(imax, 1.0, dtype=float)
A = hini*B
Q = ib**0.5*(hini)**(5.0/3.0)/manning*B #normal flow
zb = np.zeros(imax)
for i in range(1,imax):
    zb[i] = zb[i-1] + ib2*dx if i < 50 else zb[i-1] + ib*dx
zb = zb[::-1]

dAb = np.zeros(imax)

Qup = Q[0]

def Adown(time, Q, dzb, ib):
    return ( manning**2*Q**2/ib/B[-1]**2 )**0.3 * B[-1]

def Qup(time):
    return Q[0]

bedvariation(
dx,dt,manning,totalTime,outTimeStep,RunUpTime
,dsize ,dratioStandard ,dratio ,hExlayer ,A ,Q ,B ,zb ,dAb ,Qup ,Adown
,outputfilename
    )
```

    0 second
    1080000 second
    2160000 second
    3240000 second
    4320000 second
    5400000 second
    6480000 second
    7560000 second
    8640000 second
    9720000 second
    10800000 second
    Wall time: 2min 31s


## case3


```python
%%time
length = 10000.0
dx = 100.0
imax = int(length/dx) + 1
dt = 10.0
totalTime = 3000.1*3600.0
outTimeStep = 300.0*3600.0
RunUpTime = 3.0 * 3600.0
hini = 1.0
manning  = 0.03
ib = 1.0/1000.0
outputfilename = 'case03.json'

# grain diameter classification
screenclass = np.array( [ 2.0**0, 2.0**1, 2.0**2, 2.0**3, 2.0**4, 2.0**5, 2.0**6], dtype=float )/1000
dsize = np.array( [ 10**(0.5*(np.log10(screenclass[i]) + np.log10(screenclass[i+1]))) for i in range(len(screenclass)-1) ], dtype=float )

# percentage of grain size under exchange layer
dratioStandard1 = np.full_like(dsize, 1.0/len(dsize), dtype=float)
dratioStandard = np.full( (imax, len(dsize) ), dratioStandard1, dtype=float)

# initial percentage of grain size in exchange layer
dratio = np.copy(dratioStandard)

# thickness of exchange layer 
hExlayer = dsize[-1]

# Initial & Boundary condition
B = np.full(imax, 100.0, dtype=float)
for i in range(40, 50):
    B[i] = B[i-1] - 5.0
for i in range(50, 60):
    B[i] = B[i-1] + 5.0
    
A = hini*B
Q = ib**0.5*(hini)**(5.0/3.0)/manning*B #normal flow
zb = np.zeros(imax)
for i in range(1,imax):
    zb[i] = zb[i-1] + ib*dx
zb = zb[::-1]

dAb = np.zeros(imax)

Qup = Q[0]

def Adown(time, Q, dzb, ib):
    return ( manning**2*Q**2/ib/B[-1]**2 )**0.3 * B[-1]

def Qup(time):
    return Q[0]

bedvariation(
dx,dt,manning,totalTime,outTimeStep,RunUpTime
,dsize ,dratioStandard ,dratio ,hExlayer ,A ,Q ,B ,zb ,dAb ,Qup ,Adown
,outputfilename
    )
```

    0 second
    1080000 second
    2160000 second
    3240000 second
    4320000 second
    5400000 second
    6480000 second
    7560000 second
    8640000 second
    9720000 second
    10800000 second
    Wall time: 1min 53s


## case4


```python
%%time
length = 10000.0
dx = 100.0
imax = int(length/dx) + 1
dt = 10.0
totalTime = 5000.1*3600.0
outTimeStep = 500.0*3600.0
RunUpTime = 3.0 * 3600.0
hini = 1.0
manning  = 0.03
ib = 1.0/500.0
outputfilename = 'case04.json'

# grain diameter classification
screenclass = np.array( [ 2.0**0, 2.0**1, 2.0**2, 2.0**3, 2.0**4, 2.0**5, 2.0**6], dtype=float )/1000
dsize = np.array( [ 10**(0.5*(np.log10(screenclass[i]) + np.log10(screenclass[i+1]))) for i in range(len(screenclass)-1) ], dtype=float )

# percentage of grain size under exchange layer
dratioStandard1 = np.full_like(dsize, 1.0/len(dsize), dtype=float)
dratioStandard = np.full( (imax, len(dsize) ), dratioStandard1, dtype=float)

# initial percentage of grain size in exchange layer
dratio = np.copy(dratioStandard)

# thickness of exchange layer 
hExlayer = dsize[-1]

# Initial & Boundary condition
B = np.full(imax, 1.0, dtype=float)
A = hini*B
Q = ib**0.5*(hini)**(5.0/3.0)/manning*B #normal flow
zb = np.zeros(imax)
for i in range(1,imax):
    zb[i] = zb[i-1] + ib*dx
zb = zb[::-1]

WLdown = 5.0
for i, (Ap, zbp) in enumerate(zip(A, zb)):
    if( Ap/B[i] + zbp < WLdown) : A[i] = (WLdown - zbp)*B[i]

dAb = np.zeros(imax)

Qup = Q[0]

def Adown(time, Q, dzb, ib):
    return (WLdown - (dzb + zb[-1]))*B[-1]
#     return ( manning**2*Q**2/ib/B[0]**2 )**0.3 * B[0]

def Qup(time):
    return Q[0]

bedvariation(
dx,dt,manning,totalTime,outTimeStep,RunUpTime
,dsize ,dratioStandard ,dratio ,hExlayer ,A ,Q ,B ,zb ,dAb ,Qup ,Adown
,outputfilename
    )
```

    0 second
    1800000 second
    3600000 second
    5400000 second
    7200000 second
    9000000 second
    10800000 second
    12600000 second
    14400000 second
    16200000 second
    18000000 second
    Wall time: 3min 46s

