vignettes/publication-reproduction.Rmd
publication-reproduction.RmdThis document provides the code needed to reproduce all results and figures from
Michael C. Koohafkan and Bassam A. Younis (2015). Open-Channel Computation with R. The R Journal, 7(2), 249-262. URL https://journal.r-project.org/archive/2015-2/koohafkan-younis.pdf
The .Rnw file used to generate the full publication is hosted on the package development repository.
## ---- load-libs -------------------------------------------------------
library(knitr)
opts_chunk$set(fig.width = 6, fig.height = 4, fig.path = "fig/",
warning = FALSE, message = FALSE, error = FALSE, results = "asis",
out.width = '\\textwidth', cache = TRUE, cache.path = "cache/")
library(rivr)
library(ggplot2)
library(scales)
library(RColorBrewer)
library(dplyr)
library(xtable)
## ---- startparams -----------------------------------------------------
plotopts = list(theme_bw(), xlab(expression(Distance~from~control~section~~(ft))))
g = 32.2
Cm = 1.486
slope = 0.001
mannings = 0.045
flow = 250
width = 100
sideslope = 0
# calculate control depth as 1ft above the normal depth
depth.m1 = round(1 + normal_depth(slope, mannings, flow, 2, Cm,
width, sideslope), 2)
depth.m2 = round(1.1*critical_depth(flow, 1, g, width, sideslope), 2)
rivdist = 3000
## ---- restest ---------------------------------------------------------
# Test sensitivity of step size
stepsizes = c(500, 100, 50, 10, 1)
resolution.test.m1 = list()
resolution.test.m2 = list()
for(r in stepsizes){
resolution.test.m1[[paste0('dx=', r)]] = compute_profile(slope, mannings,
flow, depth.m1, Cm, g, width, sideslope, stepdist = r,
totaldist = rivdist)
resolution.test.m2[[paste0('dx=', r)]] = compute_profile(slope, mannings,
flow, depth.m2, Cm, g, width, sideslope, stepdist = r,
totaldist = rivdist)
}
resolution.plot = NULL
for(lbl in names(resolution.test.m1)){
f1 = resolution.test.m1[[lbl]]
f2 = resolution.test.m2[[lbl]]
f1['type'] = 'M1'
f2['type'] = 'M2'
f = rbind(f1, f2)
f['run'] = lbl
f['res'] = as.numeric(substr(lbl, 4, 10L))
resolution.plot = rbind(resolution.plot, f)
}
resolution.plot['res'] = factor(resolution.plot$res, levels = rev(sort(stepsizes)))
resolution.plot['run'] = factor(resolution.plot$run,
levels = names(resolution.test.m1))
resolution.plot['type'] = factor(resolution.plot$type, levels = c('M1', 'M2'))
stepsize = 10
## ----resfig ----------------------------------------------------------
ggplot(resolution.plot, aes(x = x, y = y + z)) + plotopts +
geom_line(aes(color = res, linetype = res), size = 1) +
facet_wrap(~type) + ylab(expression(River~stage~~(ft))) +
scale_linetype_manual("resolution (ft)", labels = levels(resolution.plot$res),
values = c("solid", "longdash", "dashed", "dotdash", "dotted")) +
scale_color_manual("resolution (ft)", values = brewer.pal(7, "YlGnBu")[3:7],
labels = levels(resolution.plot$res))
## ----roughness -------------------------------------------------------
roughness.test.m1 = list()
roughness.test.m2 = list()
roughnesses = 0.045*seq(0.5, 1.5, length = 5)
for(n in roughnesses){
roughness.test.m1[[paste0('n=',n)]] = compute_profile(slope, n, flow,
depth.m1, Cm, g, width, sideslope, stepdist=stepsize, totaldist=rivdist)
roughness.test.m2[[paste0('n=',n)]] = compute_profile(slope, n, flow,
depth.m2, Cm, g, width, sideslope, stepdist=stepsize, totaldist=rivdist)
}
roughness.plot = NULL
for(lbl in names(roughness.test.m1)){
f1 = roughness.test.m1[[lbl]]
f2 = roughness.test.m2[[lbl]]
f1['type'] = 'M1'
f2['type'] = 'M2'
f = rbind(f1, f2)
f['run'] = lbl
f['mannings'] = as.numeric(substr(lbl, 4, 10L))
roughness.plot = rbind(roughness.plot, f)
}
roughness.plot['run'] = factor(roughness.plot$run,
levels = names(roughness.test.m1))
roughness.plot['mannings'] = factor(roughness.plot$mannings,
levels = sort(roughnesses))
roughness.plot['type'] = factor(roughness.plot$type, levels=c('M1', 'M2'))
## ---- roughplot -------------------------------------------------------
ggplot(roughness.plot, aes(x = x, y = y + z)) +
geom_line(aes(linetype = mannings, color = mannings), size = 1) + plotopts +
facet_wrap(~ type) + ylab(expression(River~stage~~(ft))) +
scale_linetype_manual("Bed roughness", labels = levels(roughness.plot$mannings),
values = c("solid", "longdash", "dashed", "dotdash", "dotted")) +
scale_color_manual("Bed roughness", labels = levels(roughness.plot$mannings),
values = brewer.pal(6, "Oranges")[2:6])
## ---- relroughplot ----------------------------------------------------
roughness.rel.m1 = list()
roughness.rel.m2 = list()
for(n in roughnesses){
thisyn = normal_depth(slope, n, flow, 2, Cm, width, sideslope)
roughness.rel.m1[[paste0('n=',n)]] = compute_profile(slope, n, flow,
1.25*thisyn, Cm, g, width, sideslope, stepdist=stepsize, totaldist=rivdist)
roughness.rel.m2[[paste0('n=',n)]] = compute_profile(slope, n, flow,
0.75*thisyn, Cm, g, width, sideslope, stepdist=stepsize, totaldist=rivdist)
roughness.rel.m1[[paste0('n=',n)]]['pd.yn'] =
(roughness.rel.m1[[paste0('n=',n)]]$y - thisyn)/thisyn
roughness.rel.m2[[paste0('n=',n)]]['pd.yn'] =
(roughness.rel.m2[[paste0('n=',n)]]$y - thisyn)/thisyn
}
roughness.rel.plot = NULL
for(lbl in names(roughness.rel.m1)){
f1 = roughness.rel.m1[[lbl]]
f2 = roughness.rel.m2[[lbl]]
f1['type'] = 'M1'
f2['type'] = 'M2'
f = rbind(f1, f2)
f['run'] = lbl
f['mannings'] = as.numeric(substr(lbl, 4, 10L))
roughness.rel.plot = rbind(roughness.rel.plot, f)
}
roughness.rel.plot['mannings'] = factor(roughness.plot$mannings,
levels = sort(roughnesses))
roughness.rel.plot['run'] = factor(roughness.rel.plot$run,
levels = names(roughness.rel.m1))
roughness.rel.plot['type'] = factor(roughness.rel.plot$type,
levels = c('M1', 'M2'))
ggplot(roughness.rel.plot, aes(x = x, y = pd.yn)) +
geom_line(aes(linetype = mannings, color = mannings), size = 1) +
plotopts + facet_wrap(~type, scales = 'free_y') +
scale_y_continuous(expression(Percent~difference~from~normal~depth), labels = percent) +
scale_linetype_manual("Bed roughness", labels = levels(roughness.rel.plot$mannings),
values = c("solid", "longdash", "dashed", "dotdash", "dotted")) +
scale_color_manual("Bed roughness", labels = levels(roughness.rel.plot$mannings),
values = brewer.pal(6, "Oranges")[2:6])
## ---- loadup-kwm ------------------------------------------------------
oldscipen = options('scipen')
options(scipen = 1000)
slope = 0.001
extent = 150000
mannings = 0.045
B = 100
SS = 0
Cm = 1.486
g = 32.2
iflow = 250
# keep Courant number at 0.7 to balance temporal and spatial resolution
idepth = normal_depth(slope, mannings, 250, 10, 1.485, 100, SS)
iarea = channel_geom(idepth, B, SS)["A"]
cn = 0.7*iarea/iflow
# define upstream boundary condition assuming a timestep in seconds
bcfunc = function(x)
ifelse(x < 9000, 250 + (750/pi)*(1 - cos(pi*x/(60*75))), 250)
bctime = 76000
xnodes = extent/c(25, 50, 125, 250, 500, 1000, 2000, 5000) + 1
myxres = extent/(xnodes - 1)
mytres = cn*myxres
plotopts = list(theme_bw())
# run the model
modtime = list()
modelresults = list()
for(i in seq(length(xnodes))){
numnodes = xnodes[i]
xstep = myxres[i]
tstep = mytres[i]
bc = bcfunc(seq(0, bctime, by=tstep))
mp = c(1, as.integer(50000/xstep + 1), as.integer(100000/xstep + 1), numnodes)
mt = as.integer(round(seq(1, length(bc), length.out=10)))[seq(7)]
if(xstep == min(myxres)){
mt = as.integer(round(seq(1, length(bc), length.out=125)))
mtslice = mt[round(seq(1, length(mt), length.out=7))]
}
lbl = paste0('dx=',xstep)
modtime[[lbl]] <- system.time(modelresults[[lbl]] <-
route_wave(slope, mannings, Cm, g, B, SS, iflow, bc,
timestep=tstep, spacestep=xstep, numnodes=numnodes,
monitor.nodes=mp, monitor.times=mt, engine="Kinematic"))
modelresults[[lbl]]['deltax'] = xstep
modelresults[[lbl]]['deltat'] = tstep
modelresults[[lbl]]['computationtime'] = modtime[[lbl]][[3]]
modelresults[[lbl]]['label'] = lbl
}
allresults = do.call(rbind.data.frame, modelresults)
row.names(allresults) = NULL
names(allresults)[which(names(allresults)=='flow')] = "Q"
names(allresults)[which(names(allresults)=='distance')] = "x"
names(allresults)[which(names(allresults)=='time')] = "t"
## ---- floodhydrograph -------------------------------------------------
floodhydro = data.frame(time = seq(15000), flow=bcfunc(seq(15000)))
ggplot(floodhydro, aes(x=time/60, y=flow)) + geom_line(color="darkblue",
size = 1) + plotopts +
xlab(expression(time~~(minutes))) + ylab(expression(flow~~(ft^3~~s^-1)))
## ---- routeplot-kinematic ---------------------------------------------
# still plot
surfaces = group_by(filter(allresults, monitor.type=='timestep',
deltax==min(deltax)))
ggplot(filter(surfaces, step %in% mtslice[2:5]),
aes(x=x, y=Q, linetype=factor(round(t/60)), color=factor(round(t/60)))) +
geom_line(size = 1) + plotopts +
scale_color_manual("time\n(minutes)", values = brewer.pal(7, "YlGnBu")[3:6]) +
scale_linetype_manual("time\n(minutes)", values = c("solid",
"dashed", "dotdash", "dotted")) +
xlab(expression(distance~~downstream~~(ft))) +
ylab(expression(flow~~(ft^3~~s^-1)))
# animation
routeylims = c(min(surfaces$Q), max(surfaces$Q))
routexlims = c(min(surfaces$x), max(surfaces$x))
thisdt = unique(surfaces$deltat)
for(n in unique(surfaces$step))
print(
ggplot(filter(surfaces, step==n), aes(x=x, y=Q)) +
geom_line(color = "darkblue", size = 1) + plotopts +
scale_linetype_discrete("time\n(minutes)") +
scale_x_continuous(expression(distance~~downstream~~(ft)),
limits=routexlims) +
scale_y_continuous(expression(flow~~(ft^3~~s^-1)), limits=routeylims) +
ggtitle(paste(format(round((n - 1)*thisdt/60), width = 4), "minutes"))
)
## ---- xsectionplot-kinematic ------------------------------------------
xsections = filter(allresults, monitor.type == "node",
x %in% c(0,50000))[c("x","t","Q", "deltax")]
xsections['deltax'] = factor(paste(xsections$deltax, "feet"),
levels=paste(sort(unique(xsections$deltax)), "feet"))
ggplot(xsections, aes(x=t/60, y=Q, linetype=factor(format(x, big.mark=",",
trim=TRUE), levels = format(sort(unique(x)), big.mark=",", trim=TRUE)),
color = deltax)) +
geom_path(size = 1) + facet_grid(deltax~.) + plotopts +
scale_y_continuous(expression(flow~~(ft^3~~s^-1)), breaks=c(300, 500,700)) +
scale_linetype_manual("downstream\ndistance (ft)", values = c("solid",
"dashed", "dotted")) +
scale_x_continuous(expression(time~~(min)), limits=c(0,600)) +
scale_color_manual(values = c("#7570B3", "#8DA0CB", "#66C2A5", "#A6D854",
"#E5C494", "#E6AB02", "#FC8D62", "#E78AC3"), guide = FALSE)
## ---- peakstable-kinematic --------------------------------------------
peaks = as.data.frame(summarize(group_by(filter(allresults, x==50000,
monitor.type=='node'), computationtime, deltax, deltat, monitor.type),
peak.flow=max(Q), time.to.peak=t[which(Q==peak.flow)]))
real.peak = max(floodhydro$flow)
peaks['flow.percent.error'] = (peaks$peak.flow - real.peak)/real.peak
tbld = as.data.frame(peaks[c("deltax", "deltat", "flow.percent.error",
"computationtime")])
tbld["deltax"] = round(tbld$deltax)
tbld["deltat"] = round(tbld$deltat, 2)
tbld["flow.percent.error"] = round(100*tbld$flow.percent.error, 2)
names(tbld) = c("$\\Delta x$ (ft)", "$\\Delta t$ (s)", "\\% error (peak flow)",
"cost (s)")
ptable = xtable(tbld[order(tbld[,1]),], digits = c(0,0,2,2,2))
print(ptable, include.rownames = FALSE, sanitize.colnames.function = identity,
sanitize.text.function = function(x){x}, floating=FALSE, hline.after=NULL,
add.to.row=list(pos=list(-1,0, nrow(ptable)),
command=c('\\toprule ', '\\midrule ', '\\bottomrule ')),
format.args = list(big.mark = ","))
## ---- loadup-dwm ------------------------------------------------------
oldscipen = options('scipen')
options(scipen = 1000)
slope = 0.001
extent = 150000
mannings = 0.045
B = 100
SS = 0
g = 32.2
Cm = 1.486
iflow = 250
# define upstream boundary condition assuming a timestep in seconds
idepth = normal_depth(slope, mannings, 250, 10, Cm, 100, SS)
iarea = channel_geom(idepth, B, SS)[["A"]]
# keep Courant number at 0.06 to balance temporal and spatial resolution
cn = 0.06*iarea/250
bcfunc = function(x)
ifelse(x < 9000, 250 + (750/pi)*(1 - cos(pi*x/(60*75))), 250)
bctime = 76000
xnodes = extent/c(50, 125, 250, 500, 1000) + 1
myxres = extent/(xnodes - 1)
mytres = cn*myxres
plotopts = list(theme_bw())
# run the model
modtime = list()
modelresults = list()
for(i in seq(length(xnodes))){
numnodes = xnodes[i]
xstep = myxres[i]
tstep = mytres[i]
bc = bcfunc(seq(0, bctime, by=tstep))
dc = rep(-1, length(bc))
mp = c(1, as.integer(c(1000, 50000, 100000, 149000)/xstep + 1), numnodes)
mt = as.integer(round(seq(1, length(bc), length.out=10)))[seq(7)]
if(xstep == min(myxres)){
mt = as.integer(round(seq(1, length(bc), length.out=125)))
mtslice = mt[round(seq(1, length(mt), length.out=7))]
}
lbl = paste0('dx=',xstep)
modtime[[lbl]] <- system.time(modelresults[[lbl]] <-
route_wave(slope, mannings, Cm, g, B, SS, iflow, bc, dc, timestep=tstep,
spacestep=xstep, numnodes=numnodes, monitor.nodes=mp, monitor.times=mt,
engine="Dynamic", boundary.type="QQ"))
modelresults[[lbl]]['deltax'] = xstep
modelresults[[lbl]]['deltat'] = tstep
modelresults[[lbl]]['computationtime'] = modtime[[lbl]][[3]]
modelresults[[lbl]]['label'] = lbl
}
# combine results
allresults = do.call(rbind.data.frame, modelresults)
row.names(allresults) = NULL
names(allresults)[which(names(allresults)=='flow')] = "Q"
names(allresults)[which(names(allresults)=='distance')] = "x"
names(allresults)[which(names(allresults)=='time')] = "t"
allresults['t'] = (allresults$step - 1)*allresults$deltat
allresults['x'] = (allresults$node - 1)*allresults$deltax
## ---- routeplot-characteristic ----------------------------------------
# still plot
surfaces = filter(allresults, monitor.type=='timestep',
deltax==min(deltax))
ggplot(filter(surfaces, step %in% mtslice[2:5]),
aes(x=x, y=Q, linetype=factor(round(t/60)), color=factor(round(t/60)))) +
geom_line(size = 1) + plotopts +
scale_color_manual("time\n(minutes)", values = brewer.pal(7, "YlGnBu")[3:6]) +
scale_linetype_manual("time\n(minutes)", values = c("solid",
"dashed", "dotdash", "dotted")) +
xlab(expression(distance~~downstream~~(ft))) +
ylab(expression(flow~~(ft^3~~s^-1)))
# animation
routeylims = c(min(surfaces$Q), max(surfaces$Q))
routexlims = c(min(surfaces$x), max(surfaces$x))
thisdt = unique(surfaces$deltat)
for(n in unique(surfaces$step))
print(
ggplot(filter(surfaces, step==n), aes(x=x, y=Q)) +
geom_line(color = "darkblue", size = 1) + plotopts +
scale_linetype_discrete("time\n(minutes)") +
scale_x_continuous(expression(distance~~downstream~~(ft)),
limits=routexlims) +
scale_y_continuous(expression(flow~~(ft^3~~s^-1)), limits = routeylims) +
ggtitle(paste(format(round((n - 1)*thisdt/60), width = 4), "minutes"))
)
## ---- xsectionplot-characteristic -------------------------------------
xsections = filter(allresults, monitor.type == "node",
x %in% c(0, 50000, 150000))[c("x","t","Q", "deltax")]
xsections['deltax'] = factor(paste(xsections$deltax, "feet"),
levels=paste(sort(unique(xsections$deltax)), "feet"))
data(waterolympics)
realdat = NULL
for(x in unique(xsections$deltax)){
thisdat = waterolympics[waterolympics$t > 200*60,]
thisdat['deltax'] = x
realdat = rbind(realdat, thisdat)
}
realdat['deltax'] = factor(realdat$deltax)
ggplot(xsections, aes(x=t/60, y=Q, linetype=factor(format(x, big.mark=",",
trim=TRUE), levels = format(sort(unique(x)), big.mark=",", trim=TRUE)),
color = deltax)) +
geom_path(size = 1) + geom_point(data = realdat[realdat$t > 200*60,], color = "black") +
facet_grid(deltax~.) +scale_y_continuous(expression(flow~~(ft^3~~s^-1)),
breaks=c(300, 500,700)) +
scale_linetype_manual("downstream\ndistance (ft)", values = c("solid",
"dashed", "dotted")) + plotopts +
scale_x_continuous(expression(time~~(min)), limits=c(0,1250)) +
scale_color_manual(values = c("#8DA0CB", "#66C2A5", "#A6D854", "#E5C494",
"#E6AB02"), guide = FALSE)
## ---- peakstable-characteristic ---------------------------------------
peaks = as.data.frame(summarize(group_by(filter(allresults, x==50000,
monitor.type=='node'), computationtime, deltax, deltat, monitor.type),
peak.flow=max(Q), time.to.peak=t[which(Q==peak.flow)]))
real.peak.flow = max(waterolympics$Q)
real.time.to.peak = waterolympics$t[which(waterolympics$Q==real.peak.flow)][1]
peaks['flow.percent.error'] =
(peaks$peak.flow - real.peak.flow)/real.peak.flow
peaks['time.percent.error'] =
(peaks$time.to.peak - real.time.to.peak)/real.time.to.peak
tbld = as.data.frame(peaks[c("deltax", "deltat", "flow.percent.error",
"time.percent.error", "computationtime")])
tbld["flow.percent.error"] = round(100*tbld$flow.percent.error, 2)
tbld["time.percent.error"] = round(100*tbld$time.percent.error, 2)
names(tbld) = c("$\\Delta x$ (ft)", "$\\Delta t$ (s)", "\\% error (peak flow)",
"\\% error (time to peak)", "cost (s)")
ptable = xtable(tbld[order(tbld[,1]),], digits = c(0,0,2,2,2,2))
print(ptable, include.rownames=FALSE, sanitize.colnames.function=identity,
sanitize.text.function = function(x){x}, floating=FALSE, hline.after=NULL,
add.to.row=list(pos=list(-1,0, nrow(ptable)),
command=c('\\toprule ', '\\midrule ', '\\bottomrule ')),
format.args = list(big.mark = ","))
## ---- loadup-boundaries -----------------------------------------------
oldscipen = options('scipen')
options(scipen = 1000)
plotopts = list(theme_bw())
slope = 0.00008
extent = 5000
mannings = 0.013
B = 6.1
SS = 1.5
g = 9.81
Cm = 1
ic = 126
id = 5.79
ia = channel_geom(id, B, SS)[["A"]]
CN = 0.9
dx = 10
dt = round(dx*CN/(ic/ia + sqrt(id*g)), 2)
numnodes = extent/dx + 1
bctime = 2000
bc = rep(id, round(bctime/dt) + 1)
dc = rep(0, length(bc))
dt = round(bctime/(length(bc) - 1), 2)
CN = dt*(ic/ia + sqrt(id*g))/dx
mp = c(1, as.integer(c(1500, 2500, 3000, 5000)/dx + 1))
mt = as.integer(round(seq(0, length(bc)-1, by=25))) + 1L
mtslice = c(1, as.integer(c(500, 1000, 1500, 2000)/dt + 1))
dclose.lax = route_wave(slope, mannings, Cm, g, B, SS, ic, bc, dc,
timestep=dt, spacestep=dx, numnodes=numnodes, monitor.nodes=mp,
monitor.times=mt, engine="Dynamic", scheme="Lax", boundary.type="yQ")
dclose.mac = route_wave(slope, mannings, Cm, g, B, SS, ic, bc, dc,
timestep=dt, spacestep=dx, numnodes=numnodes, monitor.nodes=mp,
monitor.times=mt, engine="Dynamic", scheme="MacCormack", boundary.type="yQ")
dclose.lax['scheme'] = "Lax diffusive"
dclose.mac['scheme'] = "MacCormack predictor-corrector"
dclose = rbind(dclose.lax, dclose.mac)
row.names(dclose) = NULL
dclose["CN"] = dt*(dclose$flow/dclose$area + sqrt(dclose$depth*g))/dx
## ---- through-time ----------------------------------------------------
dclose.times = filter(dclose, monitor.type=="timestep")
ggplot(filter(dclose.times, step %in% mtslice),
aes(x=distance, y=depth, linetype = factor(round(time)),
color = factor(round(time)))) + geom_path(size = 1) +
scale_y_continuous(expression(depth~~(ft))) + facet_wrap(~scheme) +
scale_x_continuous(expression(distance~~downstream~~(ft))) + plotopts +
scale_color_manual("time (s)", values = brewer.pal(7, "YlGnBu")[3:7]) +
scale_linetype_manual("time (s)", values = c("solid", "longdash", "dashed",
"dotdash", "dotted"))
# animation
routeylims = c(min(dclose.times$depth), max(dclose.times$depth))
routexlims = c(min(dclose.times$distance), max(dclose.times$distance))
for(n in unique(dclose.times$step))
print(
ggplot(filter(dclose.times, step==n), aes(x=distance, y=depth)) +
geom_line(color = "darkblue", size = 1) + plotopts +
scale_linetype_discrete("time\n(minutes)") +
scale_x_continuous(expression(distance~~downstream~~(ft)),
limits=routexlims) +
scale_y_continuous(expression(depth~~(ft)), limits=routeylims) +
ggtitle(paste(format(round((n - 1)*dt), width = 4), "seconds")) +
facet_wrap(~scheme)
)