import S1DBedloadSolver as S1Dbed
import numpy as np
%%time
length = 1000.0
dx = 10.0
imax = int(length/dx) + 1
dt = 2.0
totalTime = 300.0 * 3600.0
outTimeStep = 5.0*3600.0
RunUpTime = 3.0 * 3600.0
hini = 1.0
manning = 0.03
ib = 1.0/1000.0
outputfilename = '1D_case1.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]
zb[50:60] -= 0.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 18000 second 36000 second 54000 second 72000 second 90000 second 108000 second 126000 second 144000 second 162000 second 180000 second 198000 second 216000 second 234000 second 252000 second 270000 second 288000 second 306000 second 324000 second 342000 second 360000 second 378000 second 396000 second 414000 second 432000 second 450000 second 468000 second 486000 second 504000 second 522000 second 540000 second 558000 second 576000 second 594000 second 612000 second 630000 second 648000 second 666000 second 684000 second 702000 second 720000 second 738000 second 756000 second 774000 second 792000 second 810000 second 828000 second 846000 second 864000 second 882000 second 900000 second 918000 second 936000 second 954000 second 972000 second 990000 second 1008000 second 1026000 second 1044000 second 1062000 second Wall time: 37.5 s
import xarray as xr
import json
import numpy as np
outputfile = '1D_case1.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('case1.nc')
del out, ds
import xarray as xr
import hvplot.xarray
dsin =xr.open_dataset('case1.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