%%time
length = 10000.0
dx = 100.0
imax = int(length/dx) + 1
dt = 10.0
totalTime = 2000.1*3600.0
outTimeStep = 200.0*3600.0
RunUpTime = 3.0 * 3600.0
hini = 1.0
manning = 0.03
ib = 1.0/200.0
ib2 = 1.0/200.0
outputfilename = '1D_case3.json'
# grain diameter classification
screenclass = np.array( [ 2.0**1, 2.0**2, 2.0**3, 2.0**4, 2.0**5, 2.0**6, 2.0**7], 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
dmax = screenclass[-1]
P = (screenclass/dmax)**1
dratioStandard1 = P[1:] - P[:-1]
dratioStandard1 = np.full_like(dsize, 1/len(dsize), dtype=float)
dratioStandard = np.full( (imax, len(dsize) ), dratioStandard1, dtype=float)
dratioStandard[0,:] = 0.0
dratioStandard[0,-1] = 1.0
# 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]
# WLdown = 2.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[-1]**2 )**0.3 * B[-1]
def Qup(time):
return Q[0]
S1Dbed.bedvariation(
dx,dt,manning,totalTime,outTimeStep,RunUpTime
,dsize ,dratioStandard ,dratio ,hExlayer ,A ,Q ,B ,zb ,dAb ,Qup ,Adown
,outputfilename, screenclass
)