import S1DBedloadSolver as S1Dbed
import numpy as np
%%time
length = 10000.0
dx = 100.0
imax = int(length/dx) + 1
dt = 10.0
totalTime = 5000.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/1000.0
outputfilename = '1D_case2.json'
# grain diameter classification
screenclass = np.array( [4.0], dtype=float )/1000
dsize = np.array( [4.0/1000.0], 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[i] = zb[i-1] + ib*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]
S1Dbed.bedvariation(
dx,dt,manning,totalTime,outTimeStep,RunUpTime
,dsize ,dratioStandard ,dratio ,hExlayer ,A ,Q ,B ,zb ,dAb ,Qup ,Adown
,outputfilename, screenclass
)
0 second 720000 second 1440000 second 2160000 second 2880000 second 3600000 second 4320000 second 5040000 second 5760000 second 6480000 second 7200000 second 7920000 second 8640000 second 9360000 second 10080000 second 10800000 second 11520000 second 12240000 second 12960000 second 13680000 second 14400000 second 15120000 second 15840000 second 16560000 second 17280000 second 18000000 second Wall time: 2min 42s
import xarray as xr
import json
import numpy as np
outputfile = '1D_case2.json'
data = json.load( open(outputfile, 'r') )
cond = data['condition']
d = data['output']
time = [ p['time']/3600 for p in d ]
x = [ p['distance'] for p in d[0]['profile'] ]
screenclass = np.array( cond['screenclass'] )*1000
zb = np.array( cond['elevation'] )
dep = []
dzb = []
for dr in d :
A = np.array( [ dd['A'] for dd in dr['profile'] ] )
dAb = np.array( [ dd['dAb'] for dd in dr['profile'] ] )
dep.append( zb + A + dAb)
dzb.append( zb + dAb)
dep = np.array(dep)
dzb = np.array(dzb)
ds = xr.Dataset({'H': (['t','x'], dep), 'zb': (['t','x'], dzb) }, coords={'x': x , 't': time})
out = ds.to_netcdf('case2.nc')
del out, ds
import xarray as xr
import hvplot.xarray
dsin =xr.open_dataset('case2.nc')
dini = dsin.isel(t=0)
dini['H_ini'] = dini['H']
dini['zb_ini'] = dini['zb']
dini = dini.drop(['H','zb'])
fini = dini.hvplot()
f = dsin.hvplot(groupby='t',x='x')
f * fini
del dsin, f, fini