setup.grid.1D.R
##==============================================================================
## Creation of a one-dimensional finite difference grid
##==============================================================================
setup.grid.1D <- function(x.up=0, x.down=NULL, L=NULL,
N=NULL, dx.1 =NULL, p.dx.1 = rep(1,length(L)),
max.dx.1 = L, dx.N =NULL, p.dx.N = rep(1,length(L)),
max.dx.N = L) {
## Check on the input
for (i in 1:length(L)) {
if (is.null(x.down[i]) && is.null(L[i]))
stop (paste("Error in setup.grid.1D! Grid zone",i,": either the length (L) or the end point (x.down) should be specified"))
}
## If the interval lengths are given, create the end points of each zone
if (is.null(x.down[1])) {
if (any(L<0)) stop (paste("Error in setup.grid.1D! L[",which(L<0),"] < 0",sep=""))
x.down <- x.up + cumsum(L)
}
## If the interval lengths are given, create the end points of each zone
if (is.null(L[1])) {
L <- diff(c(x.up,x.down))
if (any(L<0)) stop (paste("Error in setup.grid.1D! L[",which(L<0),"] < 0",sep=""))
}
## Check wether all info is available in each zone
for (i in 1:length(L)) {
if (is.null(N[i]) && is.null(dx.1[i]) && is.null(dx.N[i]))
stop (paste("Error in setup.grid.1D! Both N[",i,"], dx.1[",i,"] and dx.N[",i,"] are NULL"))
}
## Calculation of the grid cell sizes
dx <- vector()
## Power law function that controls the increase in grid size. The
## root of this function determines the power factor that corresponds to a
## desired number of cells N
p1 <- 1.00001
p2 <- 3
f.root <- function(p,dx,N,L) dx*(p^(N)-1)/(p-1) - L
## Loop over all grid zones
for (i in 1:length(L)) {
## Option 1: only N[i] is specified, all grid cells are equal
if (!(is.null(N[i])) && (is.null(dx.1[i])) && (is.null(dx.N[i]))) {
if (N[i] < 1) stop (paste("Error in setup.grid.1D! N[",i,"] < 1",sep=""))
A.dx <- rep(L[i]/N[i],N[i])
}
## Option 2: dx.1[i] specified, N[i] = NULL and dx.N[i] = NULL
if ((is.null(N[i])) && !(is.null(dx.1[i])) && (is.null(dx.N[i]))) {
if (dx.1[i]<=0)
stop (paste("Error in setup.grid.1D! dx.1[",i,"] <= 0",sep=""))
if (dx.1[i] > L[i])
stop (paste("Error in setup.grid.1D! dx.1[",i,"] > L[",i,"]",sep=""))
# use gradual increase of grid cell size at upstream interface
A.dx <- vector()
pos <- vector()
A.dx[1] <- dx.1[i]
pos[1] <- dx.1[i]
j <- 1
while (pos[j] < L[i]) {
j <- j + 1
A.dx[j] <- min(max.dx.1[i],A.dx[j-1]*p.dx.1[i])
pos[j] <- pos[j-1] + A.dx[j]
}
# rescaling to fit the interval length
A.dx <- A.dx*(L[i]/pos[length(pos)])
}
## Option 3: dx.N[i] specified, N[i] = NULL and dx.1[i] = NULL
if ((is.null(N[i])) && (is.null(dx.1[i])) && !(is.null(dx.N[i]))) {
if (dx.N[i]<=0)
stop (paste("Error in setup.grid.1D! dx.N[",i,"] <= 0",sep=""))
if (dx.N[i] > L[i])
stop (paste("Error in setup.grid.1D! dx.N[",i,"] > L[",i,"]",sep=""))
# use gradual increase of grid cell size at downstream interface
A.dx <- vector()
pos <- vector()
A.dx[1] <- dx.N[i]
pos[1] <- dx.N[i]
j <- 1
while (pos[j] < L[i]) {
j <- j + 1
A.dx[j] <- min(max.dx.N[i],A.dx[j-1]*p.dx.N[i])
pos[j] <- pos[j-1] + A.dx[j]
}
# rescaling to fit the interval length
A.dx <- A.dx*(L[i]/pos[length(pos)])
# reversing the A.dx vector
A.dx <- A.dx[length(A.dx):1]
}
## Option 4: N[i] and dx.1[i] are specified, dx.N[i] = NULL
if (!(is.null(N[i])) && !(is.null(dx.1[i])) && (is.null(dx.N[i]))) {
if (N[i] < 1) stop (paste("Error in setup.grid.1D! N[",i,"] < 1",sep=""))
if (dx.1[i]<=0)
stop (paste("Error in setup.grid.1D! dx.1[",i,"] <= 0",sep=""))
if (dx.1[i] > L[i])
stop (paste("Error in setup.grid.1D! dx.1[",i,"] > L[",i,"]",sep=""))
# estimate power in power law
if (f.root(p=p1,dx.1[i],N[i],L[i])*f.root(p=p2,dx.1[i],N[i],L[i]) >= 0)
{
A.dx <- rep(L[i]/N[i],times=N[i])
} else {
p.estim <- uniroot(f=f.root,dx=dx.1[i],N=N[i],L=L[i],lower=p1,upper=p2)$root
# use gradual increase of grid cell size at upper interface
A.dx <- vector()
pos <- vector()
A.dx[1] <- dx.1[i]
pos[1] <- dx.1[i]
j <- 1
while (j < N[i]) {
j <- j + 1
A.dx[j] <- min(max.dx.1[i],A.dx[j-1]*p.estim)
pos[j] <- pos[j-1] + A.dx[j]
}
# rescaling to fit the interval length
A.dx <- A.dx*(L[i]/pos[length(pos)])
}
}
## Option 5: N[i] and dx.N[i] are specified, dx.1[i] = NULL
if (!(is.null(N[i])) && (is.null(dx.1[i])) && !(is.null(dx.N[i]))) {
if (N[i] < 1) stop (paste("Error in setup.grid.1D! N[",i,"] < 1",sep=""))
if (dx.N[i]<=0)
stop (paste("Error in setup.grid.1D! dx.N[",i,"] <= 0",sep=""))
if (dx.N[i] > L[i])
stop (paste("Error in setup.grid.1D! dx.N[",i,"] > L[",i,"]",sep=""))
if (f.root(p=p1,dx.N[i],N[i],L[i])*f.root(p=p2,dx.N[i],N[i],L[i]) >= 0)
{
A.dx <- rep(L[i]/N[i],times=N[i])
} else {
# estimate power in power law
p.estim <- uniroot(f=f.root,dx=dx.N[i],N=N[i],L=L[i],lower=p1,upper=p2)$root
# use gradual increase of grid cell size at upper interface
A.dx <- vector()
pos <- vector()
A.dx[1] <- dx.N[i]
pos[1] <- dx.N[i]
j <- 1
while (j < N[i]) {
j <- j + 1
A.dx[j] <- min(max.dx.N[i],A.dx[j-1]*p.estim)
pos[j] <- pos[j-1] + A.dx[j]
}
# rescaling to fit the interval length
A.dx <- A.dx*(L[i]/pos[length(pos)])
}
# reversing the A.dx vector
A.dx <- A.dx[length(A.dx):1]
}
## Option 6: dx.1[i] and dx.N[i] are specified, N[i] can be NULL
if (!(is.null(dx.1[i])) && !(is.null(dx.N[i]))) {
if (dx.1[i]<=0)
stop (paste("Error in setup.grid.1D! dx.1[",i,"] <= 0",sep=""))
if (dx.1[i] > L[i]/2)
stop (paste("Error in setup.grid.1D! dx.1[",i,"] > L[",i,"]",sep=""))
if (dx.N[i]<=0)
stop (paste("Error in setup.grid.1D! dx.N[",i,"] <= 0",sep=""))
if (dx.N[i] > L[i]/2)
stop (paste("Error in setup.grid.1D! dx.N[",i,"] > L[",i,"]",sep=""))
L.A <- 0.5*L[i]
L.B <- 0.5*L[i]
# estimate power in power law
if (!is.null(N[i])) {
N.A <- as.integer(N[i]/2)
N.B <- N[i] - N.A
L.A <- (N.A/N[i])*L[i]
L.B <- L[i] - L.A
if (f.root(p=p1,dx.1[i],N.A,L.A)*f.root(p=p2,dx.1[i],N.A,L.A) >= 0)
{
p.dx.1[i] <- 1
} else {
p.dx.1[i] <- uniroot(f=f.root,dx=dx.1[i],N=N.A,L=L.A,lower=p1,upper=p2)$root
}
if (f.root(p=p1,dx.N[i],N.B,L.B)*f.root(p=p2,dx.N[i],N.B,L.B) >= 0)
{
p.dx.N[i] <- 1
} else {
p.dx.N[i] <- uniroot(f=f.root,dx=dx.N[i],N=N.B,L=L.B,lower=p1,upper=p2)$root
}
}
# use gradual increase of grid cell size at upsteam interface
A.dx <- vector()
pos <- vector()
A.dx[1] <- dx.1[i]
pos[1] <- dx.1[i]
j <- 1
while ((is.null(N[i])&&(pos[j] < L.A))||(!is.null(N[i])&&(j < N.A))) {
j <- j + 1
A.dx[j] <- min(max.dx.1[i],A.dx[j-1]*p.dx.1[i])
pos[j] <- pos[j-1] + A.dx[j]
}
A.dx <- A.dx*(L.A/pos[length(pos)])
# use gradual increase of grid cell size at downsteam interface
B.dx <- vector()
pos <- vector()
B.dx[1] <- dx.N[i]
pos[1] <- dx.N[i]
j <- 1
while ((is.null(N[i])&&(pos[j] < L.B))||(!is.null(N[i])&&(j < N.B))) {
j <- j + 1
B.dx[j] <- min(max.dx.N[i],B.dx[j-1]*p.dx.N[i])
pos[j] <- pos[j-1] + B.dx[j]
}
B.dx <- B.dx*(L.B/pos[length(pos)])
# assemble the distance vector
A.dx <- c(A.dx,B.dx[length(B.dx):1])
}
## calculation of present zone has finished, add distances to existing array
dx <- c(dx,A.dx)
}
## positions and distances
x.int <- x.up + diffinv(dx)
x.mid <- colMeans(rbind(x.int[-1],x.int[-(length(dx)+1)]))
dx.aux <- c(dx[1]/2,diff(x.mid),dx[length(dx)]/2)
## Packaging of results
Res <- list(x.up = x.up,
x.down = x.down[length(x.down)],
x.mid = x.mid, # position of centre of the grid cells, vector of length N
x.int = x.int, # position of the grid cell interfaces , vector of length N+1
dx = dx , # thickness of the grid cells , vector length N
dx.aux = dx.aux, # auxiliary vector with distances between centre of adjacent cells, first and last: half of cell size, vector of length N+1
N = length(dx)) # total number of grid cells
class(Res) <- "grid.1D"
return(Res)
}
##==============================================================================
## S3 method: Plotting of a one-dimensional finite difference grid
##==============================================================================
plot.grid.1D <- function(x,...) {
mf <- par(mfrow=c(2,1))
on.exit(par(mf))
plot(x$x.int,main="position of cells",ylab="x",xlab="index",...)
points(0.5+(1:x$N),x$x.mid,pch=16)
legend("topleft",c("x.int","x.mid"),pch=c(1,16))
plot(x$dx.aux,main="box thickness",ylab="dx",xlab="index",...)
points(0.5+(1:x$N),x$dx,pch=16)
legend("topleft",c("dx.aux","dx.mid"),pch=c(1,16))
}
