tran.volume.2D.R
##==============================================================================
## Transport in a two-dimensional finite volume grid
##==============================================================================
tran.volume.2D <- function(C,
C.x.up = C[1,], C.x.down = C[nrow(C),],
C.y.up = C[,1], C.y.down = C[,ncol(C)],
C.z = C, masscons = TRUE,
F.x.up = NULL, F.x.down = NULL,
F.y.up = NULL, F.y.down = NULL,
Disp.grid = NULL, Disp.x = NULL, Disp.y = Disp.x,
flow.grid = NULL, flow.x = NULL, flow.y = NULL,
AFDW.grid = NULL, AFDW.x = 1, AFDW.y = AFDW.x,
V = NULL, full.check = FALSE, full.output = FALSE)
{
DD <- dim(C)
Nx <- nrow(C)
Ny <- ncol(C)
if (is.null(V))
stop("'V' should be specified ")
#==============================================================================
# infilling of grids with only x.int and y.int needed
#==============================================================================
gridInt <- function(G.x, G.y, Name) { # define a function first
# check if G.x and G.y is not NULL
if (is.null(G.x) | is.null(G.y))
stop( (paste("error: ",Name,"and (",Name,".x and", Name,".y) cannot be NULL at the same time", del="")))
G.grid <- list()
# infilling of x-matrix
if (is.matrix(G.x)) {
if (sum(abs(dim(G.x) - c(Nx+1,Ny)))!=0)
stop (paste("error: ",Name,".x matrix not of correct (Nx+1) Ny dimensions", del=""))
G.grid$x.int <- G.x
} else if (class(G.x)=="prop.1D") {
G.grid$x.int <- matrix(data=G.x$int,nrow=(Nx+1),ncol=Ny)
} else if (length(G.x) == 1) {
G.grid$x.int <- matrix(data=G.x,nrow=(Nx+1),ncol=Ny)
} else if (length(G.x) != Nx+1) {
stop (paste("error: ",Name,".x should be a vector of length 1 or Nx+1", del=""))
} else { # correct length
G.grid$x.int <- matrix(data=G.x,nrow=(Nx+1),ncol=Ny)
}
# infilling of y-matrix
if (is.matrix(G.y)) {
if (sum(abs(dim(G.y) - c(Nx,Ny+1)))!=0)
stop (paste("error: ",Name,".y matrix not of correct Nx(Ny+1)dimensions", del=""))
G.grid$y.int <- G.y
} else if (class(G.y)=="prop.1D") {
G.grid$y.int <- t(matrix(data=G.y$int,ncol=Nx,nrow=(Ny+1))) #Karline: changed 007/2011
} else if (length(G.y) == 1) {
G.grid$y.int <- matrix(data=G.y,nrow=Nx,ncol=(Ny+1))
} else if (length(G.y) != Ny+1) {
stop (paste("error: ",Name,".y should be a vector of length 1 or Ny+1", del=""))
} else { # correct length
G.grid$y.int <- matrix(data=G.y,nrow=Nx,ncol=(Ny+1))
}
G.grid
}
# Need this for AFDW , Disp and flow
if (is.null(AFDW.grid)) AFDW.grid <- gridInt(AFDW.x, AFDW.y, "AFDW")
if (is.null(Disp.grid)) Disp.grid <- gridInt(Disp.x, Disp.y, "Disp")
if (is.null(flow.grid)) flow.grid <- gridInt(flow.x, flow.y, "flow")
#==============================================================================
# INPUT CHECKS
#==============================================================================
if (full.check) {
## check dimensions of input concentrations
if (!is.null(C.x.up)) {
if (!((length(C.x.up)==1) || (length(C.x.up)==(Ny))))
stop("error: C.x.up should be a vector of length 1 or ncol(C)")
}
if (!is.null(C.x.down)) {
if (!((length(C.x.down)==1) || (length(C.x.down)==(Ny))))
stop("error: C.x.down should be a vector of length 1 or ncol(C)")
}
if (!is.null(C.y.up)) {
if (!((length(C.y.up)==1) || (length(C.y.up)==(Nx))))
stop("error: C.y.up should be a vector of length 1 or nrow(C)")
}
if (!is.null(C.y.down)) {
if (!((length(C.y.down)==1) || (length(C.y.down)==(Nx))))
stop("error: C.y.down should be a vector of length 1 or nrow(C)")
}
# check dimensions of input fluxes
if (!is.null(F.x.up)) {
if (!((length(F.x.up)==1) || (length(F.x.up)==(Ny))))
stop("error: F.x.up should be a vector of length 1 or ncol(C)")
}
if (!is.null(F.x.down)) {
if (!((length(F.x.down)==1) || (length(F.x.down)==(Ny))))
stop("error: F.x.down should be a vector of length 1 or ncol(C)")
}
if (!is.null(F.y.up)) {
if (!((length(F.y.up)==1) || (length(F.y.up)==(Nx))))
stop("error: F.y.up should be a vector of length 1 or nrow(C)")
}
if (!is.null(F.y.down)) {
if (!((length(F.y.down)==1) || (length(F.y.down)==(Nx))))
stop("error: F.y.down should be a vector of length 1 or nrow(C)")
}
## check input of volumes
if (any(V <= 0))
stop("error: the volumes should always be positive")
if (nrow(V) != Nx || ncol(V) != Ny)
stop("error: the dimension of 'V' should be = dimension of 'C'")
## If mass should be conserved...
if (masscons) {
if (nrow(C.z) != Nx || ncol(C.z) != Ny)
stop("error: the dimension of 'C.z' should be = dimension of 'C' if 'masscons' = TRUE")
}
## check input of AFDW.grid
gn <- names(AFDW.grid)
if (! "x.int" %in% gn)
stop("error: AFDW.grid should be a list that contains 'x.int', the AFDW values at the interfaces of the grid cells in x-direction")
if (! "y.int" %in% gn)
stop("error: AFDW.grid should be a list that contains 'y.int', the AFDW values at the interfaces of the grid cells in y-direction")
if (is.null(AFDW.grid$x.int))
stop("error: AFDW.grid$x.int should be a list with (numeric) values")
if (is.null(AFDW.grid$y.int))
stop("error: AFDW.grid$y.int should be a list with (numeric) values")
if (any (AFDW.grid$x.int < 0)||any (AFDW.grid$x.int > 1))
stop("error: the AFDW should range between 0 and 1")
if (any (AFDW.grid$y.int < 0)||any (AFDW.grid$y.int > 1))
stop("error: the AFDW should range between 0 and 1")
## check input of Disp.grid
gn <- names(Disp.grid)
if (! "x.int" %in% gn)
stop("error: Disp.grid should be a list that contains 'x.int', the Disp values at the interfaces of the grid cells in x-direction")
if (! "y.int" %in% gn)
stop("error: Disp.grid should be a list that contains 'y.int', the Disp values at the interfaces of the grid cells in y-direction")
if (is.null(Disp.grid$x.int))
stop("error: Disp.grid$x.int should be a list with (numeric) values")
if (is.null(Disp.grid$y.int))
stop("error: Disp.grid$y.int should be a list with (numeric) values")
if (any (Disp.grid$x.int < 0)||any (Disp.grid$y.int < 0))
stop("error: the diffusion coefficient, Disp, should always be positive")
## check input of flow.grid
gn <- names(flow.grid)
if (! "x.int" %in% gn)
stop("error: flow.grid should be a list that contains 'x.int', the flows at the interfaces of the grid cells in x-direction")
if (! "y.int" %in% gn)
stop("error: flow.grid should be a list that contains 'y.int', the flows at the interfaces of the grid cells in y-direction")
if (is.null(flow.grid$x.int))
stop("error: the flows flow.grid$x.int should be a list with (numeric) values")
if (is.null(flow.grid$y.int))
stop("error: the flows flow.grid$y.int should be a list with (numeric) values")
}
## FUNCTION BODY: CALCULATIONS
## Calculate diffusive part of the flow
x.Dif.flow <- as.matrix(-Disp.grid$x.int *
diff(rbind(C.x.up, C, C.x.down, deparse.level = 0)))
y.Dif.flow <- as.matrix(-Disp.grid$y.int *
t(diff(t(cbind(C.y.up,C,C.y.down,deparse.level = 0)))))
## Calculate advective part of the flow
x.Adv.flow <- 0
if (any(flow.grid$x.int >0) ) {
vv <- flow.grid$x.int
vv[vv < 0]<-0
x.Adv.flow <- x.Adv.flow + as.matrix(vv * (
AFDW.grid$x.int * rbind(C.x.up,C,deparse.level = 0)
+ (1-AFDW.grid$x.int) * rbind(C,C.x.down,deparse.level = 0)))
}
if (any (flow.grid$x.int < 0)) {
vv <- flow.grid$x.int
vv[vv>0]<-0
x.Adv.flow <- x.Adv.flow + as.matrix(vv * (
(1-AFDW.grid$x.int) * rbind(C.x.up,C,deparse.level = 0)
+ AFDW.grid$x.int * rbind(C,C.x.down,deparse.level = 0)))
}
y.Adv.flow <- 0
if (any(flow.grid$y.int >0) ) {
vv <- flow.grid$y.int
vv[vv<0]<-0
y.Adv.flow <- y.Adv.flow + as.matrix(vv * (
AFDW.grid$y.int * cbind(C.y.up,C,deparse.level = 0)
+ (1-AFDW.grid$y.int) * cbind(C,C.y.down,deparse.level = 0)))
}
if (any (flow.grid$y.int < 0)) {
vv <- flow.grid$y.int
vv[vv>0]<-0
y.Adv.flow <- y.Adv.flow + as.matrix(vv * (
(1-AFDW.grid$y.int) * cbind(C.y.up,C,deparse.level = 0)
+ AFDW.grid$y.int * cbind(C,C.y.down,deparse.level = 0)))
}
x.flow <- x.Dif.flow + x.Adv.flow
y.flow <- y.Dif.flow + y.Adv.flow
z.Adv.flow <- 0
ww <- 0
if (masscons) {
ww <- -(flow.grid$x.int[-1, ] - flow.grid$x.int[-nrow(flow.grid$x.int), ] +
flow.grid$y.int[ ,-1] - flow.grid$y.int[ ,-ncol(flow.grid$y.int)] )
AFDW.z <- AFDW.grid$y.int[1]
if (any(ww >0) ) {
vv <- ww
vv[vv<0]<-0
z.Adv.flow <- z.Adv.flow + as.matrix(vv * (
AFDW.z * C + (1-AFDW.z) * C.z))
}
if (any (ww < 0)) {
vv <- ww
vv[vv>0]<-0
z.Adv.flow <- z.Adv.flow + as.matrix(vv * (
(1-AFDW.z) * C + AFDW.z * C.z))
}
}
## Impose boundary fluxes when needed
## Default boundary condition is no gradient
if (! is.null (F.x.up[1]))
x.flow[1,] <- F.x.up
if (! is.null (F.x.down[1]))
x.flow[nrow(x.flow),] <- F.x.down
if (! is.null (F.y.up[1]))
y.flow[,1] <- F.y.up
if (! is.null (F.y.down[1]))
y.flow[,ncol(y.flow)] <- F.y.down
## Calculate rate of change = flux gradient
dFdx <- - diff(x.flow)/ V
dFdy <- -t(diff(t(y.flow))/t(V))
dFdz <- - z.Adv.flow/V
if (!full.output) {
return (list (dC = dFdx + dFdy + dFdz, # Rate of change due to advective-diffuisve transport in each grid cell
F.x.up = x.flow[1,], # flux across upper boundary interface; positive = IN
F.x.down = x.flow[nrow(x.flow),], # flux across lower boundary interface; positive = OUT
F.y.up = y.flow[,1],
F.y.down = y.flow[,ncol(y.flow)],
F.z = z.Adv.flow ))
} else {
return (list (dC = dFdx + dFdy, # Rate of change in the centre of each grid cells
x.flow = x.flow, # flow across at the interface of each grid cell
y.flow = y.flow, # flow across at the interface of each grid cell
F.x.up = x.flow[1,],
F.x.down = x.flow[nrow(x.flow),],
F.y.up = y.flow[,1],
F.y.down = y.flow[,ncol(y.flow)],
F.z = z.Adv.flow ))
}
} # end tran.2D
