tran.volume.3D.R
##==============================================================================
## Take differences of an array
##==============================================================================
diff3D <- function(Array, along) {
dimens <- dim(Array)
if (length(dimens) != 3)
stop("'Array' should be an array with dimension = 3")
if (along == 1)
return(Array[2:dimens[1], , ] - Array[1:(dimens[1]-1), , ])
else if (along == 2)
return(Array[, 2:dimens[2], ] - Array[, 1:(dimens[2]-1), ])
else if (along == 3)
return(Array[, , 2:dimens[3]] - Array[, , 1:(dimens[3]-1)])
else
stop("'along' should be 1, 2, or 3")
}
# function to bind two matrices to an array, to left and right
Mbind <- function (Matleft = NULL, Array, Matright = NULL, along = 1) {
dimens <- dim(Array)
if (length(dimens) != 3)
stop("'Array' should be an array with dimension = 3")
if (along < 1 | along > 3)
stop("'along' should be 1, 2, or 3")
if (! is.null(Matright)) {
if (sum(abs(dimens[-along] - dim(Matright))) != 0)
stop("'Matright' not compatible with Array")
dimens[along] <- dimens[along] + 1
Array <- mbindr(Array, Matright, along)
}
if (! is.null(Matleft)) {
if (sum(abs(dimens[-along] - dim(Matleft))) != 0)
stop("'Matleft' not compatible with Array")
dimens[along] <- dimens[along] + 1
Array <- mbindl(Matleft, Array, along)
}
Array
}
# function to bind two matrices to an array, to left and right
# NO error checking
mbind <- function (Mat1, Array, Mat2, along = 1) {
dimens <- dim(Array)
dimens[along] <- dimens[along] + 2
if (along == 3)
array(dim = dimens, data = c(Mat1, Array, Mat2))
else if (along == 1)
aperm(array(dim = dimens[c(3, 2, 1)],
data = c(t(Mat1), aperm(Array, c(3, 2, 1)), t(Mat2))),
c(3, 2, 1))
else if (along == 2)
aperm(array(dim = dimens[c(1, 3, 2)],
data = c(Mat1, aperm(Array, c(1, 3, 2)), Mat2)),
c(1, 3, 2))
}
# function to bind a matrix to an array on the left
mbindl <- function (Mat1, Array, along = 1) {
dimens <- dim(Array)
dimens[along] <- dimens[along]+1
if (along == 3)
array(dim = dimens, data = c(Mat1, Array))
else if (along == 1)
aperm(array(dim = dimens[c(3, 2, 1)],
data = c(t(Mat1), aperm(Array, c(3, 2, 1)))),
c(3, 2, 1))
else if (along == 2)
aperm(array(dim = dimens[c(1, 3, 2)],
data = c(Mat1, aperm(Array, c(1, 3, 2)))),
c(1, 3, 2))
}
# function to bind a matrix to an array on the right
mbindr <- function (Array, Mat2, along = 1) {
dimens <- dim(Array)
dimens[along] <- dimens[along]+1
if (along == 3)
array(dim = dimens, data = c(Array, Mat2))
else if (along == 1)
aperm(array(dim = dimens[c(3, 2, 1)],
data = c(aperm(Array, c(3, 2, 1)), t(Mat2))),
c(3, 2, 1))
else if (along == 2)
aperm(array(dim = dimens[c(1, 3, 2)],
data = c(aperm(Array, c(1, 3, 2)), Mat2)),
c(1, 3, 2))
}
##==============================================================================
## Transport in a three-dimensional finite difference grid
##==============================================================================
tran.volume.3D <- function(C,
C.x.up = C[1, , ], C.x.down = C[dim(C)[1], ,],
C.y.up = C[, 1, ], C.y.down = C[,dim(C)[2],],
C.z.up = C[, , 1], C.z.down = C[,,dim(C)[3]],
F.x.up = NULL, F.x.down = NULL,
F.y.up = NULL, F.y.down = NULL,
F.z.up = NULL, F.z.down = NULL,
Disp.grid = NULL, Disp.x = NULL, Disp.y = Disp.x, Disp.z = Disp.x,
flow.grid = NULL, flow.x = 0, flow.y = 0, flow.z = 0,
AFDW.grid = NULL, AFDW.x = 1, AFDW.y = AFDW.x, AFDW.z = AFDW.x,
V = NULL, full.check = FALSE, full.output = FALSE)
{
Nx <- dim(C)[1]
Ny <- dim(C)[2]
Nz <- dim(C)[3]
if(length (C.x.up) == 1) C.x.up <- matrix(nrow=Ny,ncol=Nz,C.x.up)
if(length (C.x.down) == 1) C.x.down <- matrix(nrow=Ny,ncol=Nz,C.x.down)
if(length (C.y.up) == 1) C.y.up <- matrix(nrow=Nx,ncol=Nz,C.y.up)
if(length (C.y.down) == 1) C.y.down <- matrix(nrow=Nx,ncol=Nz,C.y.down)
if(length (C.z.up) == 1) C.z.up <- matrix(nrow=Nx,ncol=Ny,C.z.up)
if(length (C.z.down) == 1) C.z.down <- matrix(nrow=Nx,ncol=Ny,C.z.down)
if (is.null(V))
stop("volume of each grid cell, 'V' should be specified ")
#==============================================================================
# infilling of grids with only x.int and y.int needed
#==============================================================================
gridInt <- function(G.x,G.y,G.z,Name) { # define a function first
# check if G.x and G.y and G.z is not NULL
if (is.null(G.x) | is.null(G.y) | is.null(G.z))
stop( (paste("error: ",Name,"and (",Name,".x and", Name,".y and",
Name,".z) cannot be NULL at the same time", del="")))
G.grid <- list()
# infilling of x-array
if (is.array(G.x)) {
if (sum(abs(dim(G.x) - c(Nx+1,Ny,Nz)))!=0)
stop (paste("error: ",Name,".x array not of correct (Nx+1) Ny ,Nz dimensions", del=""))
G.grid$x.int <- G.x
} else if (length(G.x) == 1) {
G.grid$x.int <- array(data=G.x,dim=c(Nx+1,Ny,Nz))
} else if (length(G.x) == Nx+1) {
G.grid$x.int <- array(data=G.x,dim=c(Nx+1,Ny,Nz))
} else
stop (paste("error: ",Name,".x should be one element, a vector, or an array", del=""))
# infilling of y-array
if (is.array(G.y)) {
if (sum(abs(dim(G.y) - c(Nx,Ny+1,Nz)))!=0)
stop (paste("error: ",Name,".y array not of correct Nx, Ny+1 ,Nz dimensions", del=""))
G.grid$y.int <- G.y
} else if (length(G.y) == 1) {
G.grid$y.int <- array(data=G.y,dim=c(Nx,Ny+1,Nz))
} else if (length(G.y) == Ny+1) {
G.grid$y.int <- aperm(array(data=G.y,dim=c(Ny+1,Nx,Nz)),c(2,1,3))
} else
stop (paste("error: ",Name,".y should be one element, a vector or an array", del=""))
# infilling of z-array
if (is.array(G.z)) {
if (sum(abs(dim(G.z) - c(Nx,Ny,Nz+1)))!=0)
stop (paste("error: ",Name,".z array not of correct Nx, Ny ,Nz+1 dimensions", del=""))
G.grid$z.int <- G.z
} else if (length(G.z) == 1) {
G.grid$z.int <- array(data=G.z,dim=c(Nx,Ny,Nz+1))
} else if (length(G.z) == Nz+1) {
G.grid$z.int <- aperm(array(data=G.z,dim=c(Nz+1,Ny,Nx)),c(3,2,1))
} else
stop (paste("error: ",Name,".z should be one element, a vector or an array", del=""))
G.grid
}
# Need this for AFDW , D and v
if (is.null(AFDW.grid)) AFDW.grid <- gridInt(AFDW.x,AFDW.y,AFDW.z,"AFDW")
if (is.null(Disp.grid)) Disp.grid <- gridInt(Disp.x,Disp.y,Disp.z,"Disp")
if (is.null(flow.grid)) flow.grid <- gridInt(flow.x,flow.y,flow.z,"flow")
#==============================================================================
# INPUT CHECKS
#==============================================================================
if (full.check) {
## check dimensions of input concentrations
if (!is.null(C.x.up)) {
if (!((length(C.x.up)==1) || (sum(abs(dim(C.x.up) - c(Ny,Nz)))==0)))
stop("error: C.x.up should be of length 1 or a matrix with dim Ny, Nz")
}
if (!is.null(C.x.down)) {
if (!((length(C.x.down)==1) || (sum(abs(dim(C.x.down) - c(Ny,Nz)))==0)))
stop("error: C.x.down should be of length 1 or a matrix with dim Ny, Nz")
}
if (!is.null(C.y.up)) {
if (!((length(C.y.up)==1) || (sum(abs(dim(C.y.up) - c(Nx,Nz)))==0)))
stop("error: C.y.up should be of length 1 or a matrix with dim Nx, Nz")
}
if (!is.null(C.y.down)) {
if (!((length(C.y.down)==1) || (sum(abs(dim(C.y.down) - c(Nx,Nz)))==0)))
stop("error: C.y.down should be of length 1 or a matrix with dim Nx, Nz")
}
if (!is.null(C.z.up)) {
if (!((length(C.z.up)==1) || (sum(abs(dim(C.z.up) - c(Nx,Ny)))==0)))
stop("error: C.z.up should be of length 1 or a matrix with dim Nx, Ny")
}
if (!is.null(C.z.down)) {
if (!((length(C.z.down)==1) || (sum(abs(dim(C.z.down) - c(Nx,Ny)))==0)))
stop("error: C.z.down should be of length 1 or a matrix with dim Nx, Ny")
}
# check dimensions of input fluxes
if (!is.null(F.x.up)) {
if (!((length(F.x.up)==1) || (sum(abs(dim(F.x.up) - c(Ny,Nz)))!=0)))
stop("error: F.x.up should be of length 1 or a matrix with dim Ny, Nz")
}
if (!is.null(F.x.down)) {
if (!((length(F.x.down)==1) || (sum(abs(dim(F.x.down) - c(Ny,Nz)))!=0)))
stop("error: F.x.down should be of length 1 or a matrix with dim Ny, Nz")
}
if (!is.null(F.y.up)) {
if (!((length(F.y.up)==1) || (sum(abs(dim(F.y.up) - c(Nx,Nz)))!=0)))
stop("error: F.y.up should be of length 1 or a matrix with dim Nx, Nz")
}
if (!is.null(F.y.down)) {
if (!((length(F.y.down)==1) || (sum(abs(dim(F.y.down) - c(Nx,Nz)))!=0)))
stop("error: F.y.down should be of length 1 or a matrix with dim Nx, Nz")
}
if (!is.null(F.z.up)) {
if (!((length(F.z.up)==1) || (sum(abs(dim(F.z.up) - c(Nx,Ny)))!=0)))
stop("error: F.z.up should be of length 1 or a matrix with dim Nx, Ny")
}
if (!is.null(F.z.down)) {
if (!((length(F.z.down)==1) || (sum(abs(dim(F.z.down) - c(Nx,Ny)))!=0)))
stop("error: F.z.down should be of length 1 or a matrix with dim Nx, Ny")
}
## check input of volumes
if (any(V <= 0))
stop("error: the volumes should always be positive")
if (dim(V)[1] != Nx || dim(V)[2] != Ny || dim(V)[3] != Nz)
stop("error: the dimension of 'V' should be = dimension of 'C'")
## check input of AFDW.grid
if (is.null(AFDW.x) && is.null(AFDW.y) && is.null(AFDW.z) && is.null(AFDW.grid))
stop("error: AFDW.x, AFDW.y, AFDW.z and AFDW.grid cannot be NULL at the same time")
gn <- names(AFDW.grid)
if (! "x.int" %in% gn)
stop("error: AFDW.grid should be a list that contains 'x.int', the 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 values at the interfaces of the grid cells in y-direction")
if (! "z.int" %in% gn)
stop("error: AFDW.grid should be a list that contains 'z.int', the values at the interfaces of the grid cells in z-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 (is.null(AFDW.grid$z.int))
stop("error: AFDW.grid$z.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")
if (any (AFDW.grid$z.int < 0)||any (AFDW.grid$z.int > 1))
stop("error: the AFDW should range between 0 and 1")
## check input of Disp.grid
if (is.null(Disp.x) && is.null(Disp.y) && is.null(Disp.grid))
stop("error: Disp.x, Disp.y, and Disp.grid cannot be NULL at the same time")
gn <- names(Disp.grid)
if (! "x.int" %in% gn)
stop("error: DDisp.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 (! "z.int" %in% gn)
stop("error: Disp.grid should be a list that contains 'z.int', the Disp values at the interfaces of the grid cells in z-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 (is.null(Disp.grid$z.int))
stop("error: Disp.grid$z.int should be a list with (numeric) values")
if (any (Disp.grid$x.int < 0)||any (Disp.grid$y.int < 0)||any (Disp.grid$z.int < 0))
stop("error: the diffusion coefficient should always be positive")
## check input of flow.grid
if (is.null(flow.x) && is.null(flow.y) && is.null(flow.grid))
stop("error: flow.x, flow.y, and flow.grid cannot be NULL at the same time")
gn <- names(flow.grid)
if (! "x.int" %in% gn)
stop("error: flow.grid should be a list that contains 'x.int', the flow values 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 flow values at the interfaces of the grid cells in y-direction")
if (! "z.int" %in% gn)
stop("error: flow.grid should be a list that contains 'z.int', the flow values at the interfaces of the grid cells in z-direction")
if (is.null(flow.grid$x.int))
stop("error: flow.grid$x.int should be a list with (numeric) values")
if (is.null(flow.grid$y.int))
stop("error: flow.grid$y.int should be a list with (numeric) values")
if (is.null(flow.grid$z.int))
stop("error: flow.grid$z.int should be a list with (numeric) values")
}
## FUNCTION BODY: CALCULATIONS
## Calculate diffusive part of the flow
x.Dif.flow <- -Disp.grid$x.int *
diff3D(mbind(C.x.up, C, C.x.down, along=1), along=1)
y.Dif.flow <- -Disp.grid$y.int *
diff3D(mbind(C.y.up, C, C.y.down, along=2), along=2)
z.Dif.flow <- -Disp.grid$z.int *
diff3D(mbind(C.z.up, C, C.z.down, along=3), along=3)
## 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 + vv * (
AFDW.grid$x.int * mbindl (C.x.up, C, along = 1)
+ (1-AFDW.grid$x.int) * mbindr (C, C.x.down, along = 1))
}
if (any (flow.grid$x.int < 0)) {
vv <- flow.grid$x.int
vv [vv > 0] <- 0
x.Adv.flow <- x.Adv.flow + vv * (
(1-AFDW.grid$x.int) * mbindl(C.x.up, C, along = 1)
+ AFDW.grid$x.int * mbindr(C, C.x.down, along = 1))
}
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 + vv * (
AFDW.grid$y.int * mbindl(C.y.up , C, along = 2)
+ (1-AFDW.grid$y.int) * mbindr(C, C.y.down, along = 2))
}
if (any (flow.grid$y.int < 0)) {
vv <- flow.grid$y.int
vv[vv>0]<-0
y.Adv.flow <- y.Adv.flow + vv * (
(1-AFDW.grid$y.int) * mbindl(C.y.up , C, along = 2)
+ AFDW.grid$y.int * mbindr(C, C.y.down, along = 2))
}
z.Adv.flow <- 0
if (any(flow.grid$z.int > 0 ) ) {
vv <- flow.grid$z.int
vv[vv<0]<-0
z.Adv.flow <- z.Adv.flow + vv * (
AFDW.grid$z.int * mbindl(C.z.up, C, along = 3)
+ (1-AFDW.grid$z.int) * mbindr(C, C.z.down, along = 3))
}
if (any (flow.grid$z.int < 0)) {
vv <- flow.grid$z.int
vv[vv>0]<-0
z.Adv.flow <- z.Adv.flow + vv * (
(1-AFDW.grid$z.int) * mbindl(C.z.up, C, along = 3)
+ AFDW.grid$z.int * mbindr(C, C.z.down, along = 3))
}
x.flow <- x.Dif.flow + x.Adv.flow
y.flow <- y.Dif.flow + y.Adv.flow
z.flow <- z.Dif.flow + z.Adv.flow
## 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[dim(x.flow)[1],,] <- 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[,dim(y.flow)[2],] <- F.y.down
if (! is.null (F.z.up[1]))
z.flow[,,1] <- F.z.up
if (! is.null (F.z.down[1]))
z.flow[,,dim(z.flow)[3]] <- F.z.down
## Calculate rate of change = flow gradient NOG DOEN
dFdx <- - (diff3D(x.flow,along=1))/ V
dFdy <- - (diff3D(y.flow,along=2))/ V
dFdz <- - (diff3D(z.flow,along=3))/ 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,,], # flow across lower boundary interface; positive = IN
F.x.down = x.flow[dim(x.flow)[1],,], # flow across lower boundary interface; positive = OUT
F.y.up = y.flow[,1,], # flow across lower boundary interface; positive = IN
F.y.down = y.flow[,dim(y.flow)[2],], # flow across lower boundary interface; positive = OUT
F.z.up = z.flow[,,1], # flow across lower boundary interface; positive = IN
F.z.down = z.flow[,,dim(z.flow)[3]])) # flow across lower boundary interface; positive = OUT
} else {
return (list (dC = dFdx + dFdy + dFdz, # Rate of change due to transport in each grid cell
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
z.flow = z.flow, # flow across at the interface of each grid cell
F.x.up = x.flow[1,,], # flow across lower boundary interface; positive = IN
F.x.down = x.flow[dim(x.flow)[1],,], # flow across lower boundary interface; positive = OUT
F.y.up = y.flow[,1,], # flow across lower boundary interface; positive = IN
F.y.down = y.flow[,dim(y.flow)[2],], # flow across lower boundary interface; positive = OUT
F.z.up = z.flow[,,1], # flow across lower boundary interface; positive = IN
F.z.down = z.flow[,,dim(z.flow)[3]])) # flow across lower boundary interface; positive = OUT
}
} # end tran.3D
