Skip to content

Commit

Permalink
pre-release of v1.0: implementation of S3 class rivr and methods for
Browse files Browse the repository at this point in the history
`plot`, `summary`, `print`, `head`, `tail`. Also updated namespace to
import from base packages `graphics` and `utils`.
  • Loading branch information
mkoohafkan committed Aug 14, 2015
1 parent f9699b2 commit 09fce29
Show file tree
Hide file tree
Showing 5 changed files with 189 additions and 48 deletions.
4 changes: 2 additions & 2 deletions DESCRIPTION
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
Package: rivr
Type: Package
Title: Steady and Unsteady Open-Channel Flow Computation
Version: 0.9.3
Version: 1.0.0.9000
Date: 2015-04-01
Author: Michael C Koohafkan [aut, cre]
Maintainer: Michael C Koohafkan <michael.koohafkan@gmail.com>
Expand All @@ -12,7 +12,7 @@ Description: A tool for undergraduate courses in open-channel hydraulics.
URL: https://github.com/mkoohafkan/rivr
BugReports: https://github.com/mkoohafkan/rivr/issues
License: GPL (>= 3)
Depends: R (>= 3.1.2)
Depends: R (>= 3.1.2), utils, graphics
Imports: Rcpp (>= 0.11.3)
Suggests: dplyr, ggplot2, knitr
LinkingTo: Rcpp
Expand Down
11 changes: 11 additions & 0 deletions NAMESPACE
Original file line number Diff line number Diff line change
@@ -1,5 +1,10 @@
# Generated by roxygen2 (4.1.1): do not edit by hand

S3method(head,rivr)
S3method(plot,rivr)
S3method(print,rivr)
S3method(summary,rivr)
S3method(tail,rivr)
export(channel_geom)
export(compute_profile)
export(conveyance)
Expand All @@ -8,4 +13,10 @@ export(froude)
export(normal_depth)
export(route_wave)
importFrom(Rcpp,evalCpp)
importFrom(graphics,legend)
importFrom(graphics,par)
importFrom(graphics,plot)
importFrom(utils,head)
importFrom(utils,str)
importFrom(utils,tail)
useDynLib(rivr)
91 changes: 51 additions & 40 deletions R/graduallyvariedflow_v2.r
Original file line number Diff line number Diff line change
Expand Up @@ -12,47 +12,50 @@
#' @importFrom Rcpp evalCpp
NULL

check_profile = function(So, n, Q, g, Cm, B, SS, y0, stepdist){
get_profile = function(So, n, Q, g, Cm, B, SS, y0){
yc = critical_depth(Q, y0, g, B, SS)
yn = normal_depth(So, n, Q, y0, Cm, B, SS)
if(So < 0){ # adverse slope
if(y0 > yc){
message("A2 profile specified. Computing upstream profile")
return(-abs(stepdist))
} else
stop("A3 profile cannot be computed (rapidly-varied flow)")
} else if (So == 0){ # horizontal slope
if(y0 > yc){
message("H2 profile specified. Computing upstream profile")
return(-abs(stepdist))
} else
stop("H3 profile cannot be computed (rapidly-varied flow)")
} else if(yn > yc){ # Mild slope
if(y0 > yn)
message("M1 profile specified. Computing upstream profile")
else if(y0 > yc)
message("M2 profile specified. Computing upstream profile")
if(So < 0) # adverse slope
if (y0 > yc)
return("A2")
else
stop("M3 profile cannot be computed (rapidly-varied flow)")
return(-abs(stepdist))
} else if(yc > yn){ # steep slope
return("A3")
else if (So == 0) # horizontal slope
if(y0 > yc)
stop("S1 profile cannot be computed (rapidly-varied flow)")
else if(yn > y0)
message("S3 profile specified. Computing downstream profile")
return("H2")
else
return("H3")
else if (yn > yc) # Mild slope
if(y0 > yn)
return("M1")
else if (y0 > yc)
return("M2")
else
return("M3")
else if (yc > yn) # steep slope
if (y0 > yc)
return("S1")
else if (yn > y0)
return("S3")
else
message("S2 profile specified. Computing downstream profile")
return(abs(stepdist))
} else { # critical profile
if(y0 > yc){
message("C1 profile specified. Computing upstream profile")
return(-abs(stepdist))
}
else {
message("C3 profile specified. Computing downstream profile")
return(abs(stepdist))
}
}
return("S2")
else # critical profile
if (y0 > yc)
return("C1")
else
return("C3")
}

check_profile = function(p){
if(any(p == c("S1", "M3", "H3", "A3"))){
stop(p, " profile cannot be computed (rapidly-varied flow)")
} else if(any(p == c("A2", "H2", "M1", "M2", "C1"))){
message(p, " profile specified. Computing upstream profile")
function(x) -abs(x)
} else {
message(p, " profile specified. Computing downstream profile")
function(x) abs(x)
}
}

#' @title Gradually-varied flow profiles
Expand Down Expand Up @@ -97,10 +100,18 @@ check_profile = function(So, n, Q, g, Cm, B, SS, y0, stepdist){
#' compute_profile(0.001, 0.045, 250, 0.64, 1.486, 32.2, 100, 0, stepdist=50, totaldist=3000)
#' @export
compute_profile = function(So, n, Q, y0, Cm, g, B, SS, z0=0, x0=0, stepdist, totaldist){
stepsize = stepdist
# determine profile type
stepsize = check_profile(So, n, Q, g, Cm, B, SS, y0, stepdist)
proftype = get_profile(So, n, Q, g, Cm, B, SS, y0)
stepsize = check_profile(proftype)(stepdist)
res = as.data.frame(loop_step(So, n, Q, Cm, g, y0, B, SS, z0, x0, stepsize, totaldist))
names(res) = c("x", "z", "y", "v", "A", "Sf", "E", "Fr")
return(res)
names(res) = c("x", "z", "y", "v", "A", "Sf", "E", "Fr")
retobj = res
class(retobj) = c("rivr", "data.frame")
attr(retobj, "call") = match.call()
attr(retobj, "simtype") = "gvf"
attr(retobj, "proftype") = proftype
attr(retobj, "modspec") = list(delta.x = stepdist, channel.length = totaldist,
normal.depth = normal_depth(So, n, Q, y0, Cm, B, SS))
attr(retobj, "channel.geometry") = list(So = So, n = n, B = B, SS = SS)
return(retobj)
}
109 changes: 109 additions & 0 deletions R/rivr-methods.r
Original file line number Diff line number Diff line change
@@ -0,0 +1,109 @@
#' @export
summary.rivr = function(object, ...){
if(attr(object, "simtype") == "gvf")
summarize_gvf(object, ...)
else if(attr(object, "simtype") == "usf")
summarize_usf(object, ...)
else
stop("Attribute 'simtype' not recognized. Supported simulations are ",
"'usf' (unsteady flow) and 'gvf' (gradually-varied flow).")
}

summarize_gvf = function(x, ...){
is = unique(round(seq(1, nrow(x), length.out = 5)))
r = data.frame(x = x$x[is], wse = x$y[is] + x$z[is])
r["pfn"] = round((x$y[is]/attr(x, "modspec")$normal.depth - 1)*100, 2)
r["Fr"] = x$Fr[is]
colnames(r) = c("Distance from control section", "Water-surface elevation",
"Percent from normal", "Froude number")
print(r, rownames = FALSE)
invisible(r)
}

summarize_usf = function(x, ...){
r = list()
sx = x[x$monitor.type == "node",]
r[[1]] = data.frame(distance.downstream = unique(sx$distance),
peak.flow = tapply(sx$flow, sx$node, max),
time.to.peak = sx$time[tapply(sx$flow, sx$node, which.max)])
colnames(r[[1]]) = gsub("\\.", " ", names(r))
st = x[x$monitor.type == "timestep",]
r[[2]] = data.frame(time.since.start = unique(st$time),
peak.flow = tapply(st$flow, st$time, max),
distance.travelled = st$distance[tapply(st$flow, st$time, which.max)])
colnames(r[[2]]) = gsub("\\.", " ", names(r))
names(r) = paste("Summary by", c("timestep", "node"))
lapply(r, print, row.names = FALSE)
invisible(r)
}

#' @importFrom utils head
#' @export
head.rivr = function(x, ...){
class(x) = "data.frame"
NextMethod(x, ...)
}

#' @importFrom utils tail
#' @export
tail.rivr = function(x, ...){
class(x) = "data.frame"
NextMethod(x, ...)
}

#' @importFrom utils str
#' @export
print.rivr = function(x, ...){
cat("Simulation type:\n ")
if(attr(x, "simtype") == "gvf")
cat("Gradually-varied flow")
else
cat("Unsteady flow", paste0("(", attr(x, "engine"), ")"))
cat("\n\nCall:\n ", paste(deparse(attr(x, "call")), sep = "\n",
collapse = "\n"), "\n\n", sep = "")
cat("\nChannel geometry:\n")
str(attr(x, "channel.geometry"), comp.str = " ", no.list = TRUE,
give.head = FALSE)
cat("\nModel specification:\n")
str(attr(x, "modspec"), comp.str = " ", no.list = TRUE,
give.head = FALSE)
cat("\nData:\n")
str(unclass(x), give.attr = FALSE, comp.str = " ", no.list = TRUE)
}

#' @importFrom graphics plot
#' @importFrom graphics legend
#' @importFrom graphics par
#' @export
plot.rivr = function(x, ...){
if (attr(x, "simtype") == "gvf") {
with(x, plot(x, y + z, type = "l",
xlab = "Distance from control section", ylab = "Water-surface elevation",
main = "Water-surface elevation profile"))
} else {
par(mfrow = c(1, 2))
with(x[x$monitor.type == "node",], {
xvals <- tapply(time, distance, function(x) return(x))
yvals <- tapply(flow, distance, function(x) return(x))
plot(min(unlist(xvals)):max(unlist(xvals)),
ylim = (c(min(unlist(yvals)), max(unlist(yvals)))),type = "n",
xlab = "Time since start", ylab = "Flow",
main = "Flow at monitored nodes")
mapply(lines, xvals, yvals, lty = 1:length(unique(node)))
legend("topright", paste("x =", unique(distance)),
lty = 1:length(unique(node)), bty = "n")
})
with(x[x$monitor.type == "timestep",], {
xvals <- tapply(distance, time, function(x) return(x))
yvals <- tapply(flow, time, function(x) return(x))
plot(min(unlist(xvals)):max(unlist(xvals)),
ylim = (c(min(unlist(yvals)), max(unlist(yvals)))),type = "n",
xlab = "Distance downstream", ylab = "Flow",
main = "Flow at monitored timesteps")
mapply(lines, xvals, yvals, lty = 1:length(unique(step)))
legend("topright", paste("t =", unique(time)),
lty = 1:length(unique(step)), bty = "n")
})
}
par(mfrow = c(1, 1))
}
22 changes: 16 additions & 6 deletions R/route_wave_v4.r
Original file line number Diff line number Diff line change
Expand Up @@ -139,6 +139,7 @@ route_wave = function(So, n, Cm, g, B, SS,
if(engine == 'Kinematic'){
reslist = kinematic_wave(So, n, Cm, g, B, SS, numnodes, boundary.condition,
initial.condition, timestep, spacestep, monitor.nodes - 1L, monitor.times - 1L)
scheme = NULL
} else {
scheme = match.arg(scheme, c("MacCormack", "Lax"))
if(scheme == "MacCormack"){
Expand Down Expand Up @@ -172,15 +173,15 @@ route_wave = function(So, n, Cm, g, B, SS,
mpnames = c("mpdepth", "mpvelocity", "mparea")
mtnames= c("mtdepth", "mtvelocity", "mtarea")
addnames = c("depth", "velocity", "area")
for(n in mpnames){
thismat = reslist[[n]]
for(nm in mpnames){
thismat = reslist[[nm]]
colnames(thismat) = as.integer(monitor.nodes)
rownames(thismat) = as.integer(seq(nrow(thismat)))
thisdf = stackmat(thismat)
mpdf = cbind(mpdf, thisdf[, 3])
}
for(n in mtnames){
thismat = reslist[[n]]
for(nm in mtnames){
thismat = reslist[[nm]]
colnames(thismat) = as.integer(seq(ncol(thismat)))
rownames(thismat) = as.integer(monitor.times)
thisdf = stackmat(thismat)
Expand All @@ -194,6 +195,15 @@ route_wave = function(So, n, Cm, g, B, SS,
allres["distance"] = (allres$node - 1)*spacestep
allres["time"] = (allres$step - 1)*timestep
retnames = c("step", "node", "time", "distance", "flow",
addnames, "monitor.type")
return(allres[retnames])
addnames, "monitor.type")
retobj = allres[retnames]
class(retobj) = c("rivr", "data.frame")
attr(retobj, "call") = match.call()
attr(retobj, "simtype") = "usf"
attr(retobj, "modspec") = list(dx = spacestep, dt = timestep,
nodes.monitored = monitor.nodes,
steps.monitored = monitor.times)
attr(retobj, "engine") = paste(c(engine, scheme), collapse = ":")
attr(retobj, "channel.geometry") = list(So = So, n = n, B = B, SS = SS)
return(retobj)
}

0 comments on commit 09fce29

Please sign in to comment.