Nothing
## -----------------------------------------------------------------------------
## GRADIENT functions for internal use of the ReacTran functions
## -----------------------------------------------------------------------------
## -----------------------------------------------------------------------------
## Binding one or two matrices to an array, to left and right
## -----------------------------------------------------------------------------
# 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))
}
# Differences of a 3-D array, no checks, internal use only
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")
}
##==============================================================================
## Transport in a three-dimensional finite difference grid
##==============================================================================
tran.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]],
flux.x.up = NULL, flux.x.down = NULL,
flux.y.up = NULL, flux.y.down = NULL,
flux.z.up = NULL, flux.z.down = NULL,
a.bl.x.up = NULL, a.bl.x.down = NULL,
a.bl.y.up = NULL, a.bl.y.down = NULL,
a.bl.z.up = NULL, a.bl.z.down = NULL,
D.grid = NULL, D.x = NULL, D.y = D.x, D.z = D.x,
v.grid = NULL, v.x = 0, v.y = 0, v.z = 0,
AFDW.grid = NULL, AFDW.x = 1, AFDW.y = AFDW.x, AFDW.z = AFDW.x,
VF.grid = NULL, VF.x = 1, VF.y = VF.x, VF.z = VF.x,
A.grid = NULL, A.x = 1, A.y = 1, A.z = 1,
grid = NULL, dx = NULL, dy = NULL, dz = NULL,
full.check = FALSE, full.output = FALSE)
{
if (is.null(grid))
if (is.null(dx) | is.null(dy) | is.null(dz))
stop("error: either 'grid' or ('dx' and 'dy' and 'dz') should be specified ")
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)
# DEFAULT INFILLING OF GRID PARAMETERS
# infilling of 3D grid
if (is.null(grid)) {
DX <- if (is.list(dx)) dx$dx else rep(dx,length.out=Nx)
DXaux <- if (is.list(dx)) dx$dx.aux else 0.5*(c(0,rep(dx,length.out=Nx))+
c(rep(dx,length.out=Nx),0))
DY <- if (is.list(dy)) dy$dx else rep(dy,length.out=Ny)
DYaux <- if (is.list(dy)) dy$dx.aux else 0.5*(c(0,rep(dy,length.out=Ny))+
c(rep(dy,length.out=Ny),0))
DZ <- if (is.list(dz)) dz$dx else rep(dz,length.out=Nz)
DZaux <- if (is.list(dz)) dz$dx.aux else 0.5*(c(0,rep(dz,length.out=Nz))+
c(rep(dz,length.out=Nz),0))
grid <- list(
dx = DX, dx.aux= DXaux,
dy = DY, dy.aux= DYaux,
dz = DZ, dz.aux= DZaux
)
}
#==============================================================================
# infilling of grids with x.int, y.int, z.int, x.mid, y.mid, z.mid needed
#==============================================================================
gridFill <- 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
G.grid$x.mid <- 0.5*(G.x[1:Nx,,]+G.x[2:(Nx+1),,])
} else if (length(G.x) == 1) {
G.grid$x.int <- array(data=G.x,dim=c(Nx+1,Ny,Nz))
G.grid$x.mid <- array(data=G.x,dim=c(Nx,Ny,Nz))
} else
stop (paste("error: ",Name,".x should be one element 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,".x array not of correct Nx, Ny+1 ,Nz dimensions", del=""))
G.grid$y.int <- G.y
G.grid$y.mid <- 0.5*(G.y[,1:Ny,]+G.y[,2:(Ny+1),])
} else if (length(G.y) == 1) {
G.grid$y.int <- array(data=G.y,dim=c(Nx,Ny+1,Nz))
G.grid$y.mid <- array(data=G.y,dim=c(Nx,Ny,Nz))
} else
stop (paste("error: ",Name,".y should be one element 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
G.grid$z.mid <- 0.5*(G.z[,1:Ny,]+G.z[,2:(Ny+1),])
} else if (length(G.z) == 1) {
G.grid$z.int <- array(data=G.z,dim=c(Nx,Ny,Nz+1))
G.grid$z.mid <- array(data=G.z,dim=c(Nx,Ny,Nz))
} else
stop (paste("error: ",Name,".z should be one element or an array", del=""))
G.grid
}
# Need this for VF and A (volume fraction and surface
if (is.null(VF.grid)) VF.grid <- gridFill(VF.x,VF.y,VF.z,"VF")
if (is.null(A.grid)) A.grid <- gridFill(A.x,A.y,A.z,"A")
#==============================================================================
# infilling of other 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(D.grid)) D.grid <- gridInt(D.x,D.y,D.z,"D")
if (is.null(v.grid)) v.grid <- gridInt(v.x,v.y,v.z,"v")
#==============================================================================
# INPUT CHECKS
#==============================================================================
if (full.check) {
# A function to test dimensions
TestDims <- function(var, N1, N2, varname, dimname) {
if (is.null(var)) return
if (length(var) == 1) return
if (any(dim(var) - c(N1,N2) !=0))
stop("error: ", varname, " should be of length 1 or a matrix with dim ", dimname)
}
## check dimensions of input concentrations
TestDims(C.x.up, Ny, Nz, "C.x.up", "Ny, Nz")
TestDims(C.x.down, Ny, Nz, "C.x.down", "Ny, Nz")
TestDims(C.y.up, Nx, Nz, "C.y.up", "Nx, Nz")
TestDims(C.y.down, Nx, Nz, "C.y.down", "Nx, Nz")
TestDims(C.z.up, Nx, Ny, "C.z.up", "Nx, Ny")
TestDims(C.z.down, Nx, Ny, "C.z.down", "Nx, Ny")
# check dimensions of input fluxes
TestDims(flux.x.up, Ny, Nz, "flux.x.up", "Ny, Nz")
TestDims(flux.x.down, Ny, Nz, "flux.x.down", "Ny, Nz")
TestDims(flux.y.up, Nx, Nz, "flux.y.up", "Nx, Nz")
TestDims(flux.y.down, Nx, Nz, "flux.y.down", "Nx, Nz")
TestDims(flux.z.up, Nx, Ny, "flux.z.up", "Nx, Ny")
TestDims(flux.z.down, Nx, Ny, "flux.z.down", "Nx, Ny")
## check input of grid
if (is.null(dx) && is.null(dy) && is.null(dz) && is.null(grid))
stop("error: dx, dy, dz and grid cannot be NULL at the same time")
gn <- names(grid)
if (! "dx" %in% gn)
stop("error: grid should be a list that contains 'dx' ")
if (! "dx.aux" %in% gn)
stop("error: grid should be a list that contains 'dx.aux' ")
if (! "dy" %in% gn)
stop("error: grid should be a list that contains 'dy' ")
if (! "dy.aux" %in% gn)
stop("error: grid should be a list that contains 'dy.aux' ")
if (! "dz" %in% gn)
stop("error: grid should be a list that contains 'dz' ")
if (! "dz.aux" %in% gn)
stop("error: grid should be a list that contains 'dz.aux' ")
if (is.null(grid$dx) || is.null(grid$dx.aux))
stop("error: the grid should be a list with (numeric) values for 'dx' and 'dx.aux' ")
if (is.null(grid$dy) || is.null(grid$dy.aux))
stop("error: the grid should be a list with (numeric) values for 'dy' and 'dy.aux' ")
if (is.null(grid$dz) || is.null(grid$dz.aux))
stop("error: the grid should be a list with (numeric) values for 'dz' and 'dz.aux' ")
if (any(grid$dx <= 0) || any(grid$dx.aux <= 0) )
stop("error: the grid distances dx and dx.aux should always be positive")
if (any(grid$dy <= 0) || any(grid$dy.aux <= 0) )
stop("error: the grid distances dy and dy.aux should always be positive")
if (any(grid$dz <= 0) || any(grid$dz.aux <= 0) )
stop("error: the grid distances dz and dz.aux should always be positive")
## 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 D.grid
if (is.null(D.x) && is.null(D.y) && is.null(D.grid))
stop("error: D.x, D.y, and D.grid cannot be NULL at the same time")
gn <- names(D.grid)
if (! "x.int" %in% gn)
stop("error: D.grid should be a list that contains 'x.int', the D values at the interfaces of the grid cells in x-direction")
if (! "y.int" %in% gn)
stop("error: D.grid should be a list that contains 'y.int', the D values at the interfaces of the grid cells in y-direction")
if (! "z.int" %in% gn)
stop("error: D.grid should be a list that contains 'z.int', the D values at the interfaces of the grid cells in z-direction")
if (is.null(D.grid$x.int))
stop("error: D.grid$x.int should be a list with (numeric) values")
if (is.null(D.grid$y.int))
stop("error: D.grid$y.int should be a list with (numeric) values")
if (is.null(D.grid$z.int))
stop("error: D.grid$z.int should be a list with (numeric) values")
if (any (D.grid$x.int < 0)||any (D.grid$y.int < 0)||any (D.grid$z.int < 0))
stop("error: the diffusion coefficient should always be positive")
## check input of v.grid
if (is.null(v.x) && is.null(v.y) && is.null(v.grid))
stop("error: v.x, v.y, and v.grid cannot be NULL at the same time")
gn <- names(v.grid)
if (! "x.int" %in% gn)
stop("error: v.grid should be a list that contains 'x.int', the velocity values at the interfaces of the grid cells in x-direction")
if (! "y.int" %in% gn)
stop("error: v.grid should be a list that contains 'y.int', the velocity values at the interfaces of the grid cells in y-direction")
if (! "z.int" %in% gn)
stop("error: v.grid should be a list that contains 'z.int', the velocity values at the interfaces of the grid cells in z-direction")
if (is.null(v.grid$x.int))
stop("error: the advective velocity v.grid$x.int should be a list with (numeric) values")
if (is.null(v.grid$y.int))
stop("error: the advective velocity v.grid$y.int should be a list with (numeric) values")
if (is.null(v.grid$z.int))
stop("error: the advective velocity v.grid$z.int should be a list with (numeric) values")
## check input of VF.grid
gn <- names(VF.grid)
if (! "x.int" %in% gn)
stop("error: VF.grid should be a list that contains 'x.int'")
if (! "y.int" %in% gn)
stop("error: VF.grid should be a list that contains 'y.int'")
if (! "z.int" %in% gn)
stop("error: VF.grid should be a list that contains 'z.int'")
if (! "x.mid" %in% gn)
stop("error: VF.grid should be a list that contains 'x.mid'")
if (! "y.mid" %in% gn)
stop("error: VF.grid should be a list that contains 'y.mid'")
if (! "z.mid" %in% gn)
stop("error: VF.grid should be a list that contains 'z.mid'")
if (is.null(VF.grid$x.int) || is.null(VF.grid$y.int)
|| is.null(VF.grid$x.mid) || is.null(VF.grid$y.mid)
|| is.null(VF.grid$z.mid) || is.null(VF.grid$z.int))
stop("error: VF should contain (numeric) values")
if (any (VF.grid$x.int < 0) || any (VF.grid$y.int < 0)
|| any (VF.grid$x.mid < 0) || any (VF.grid$y.mid < 0)
|| any (VF.grid$z.mid < 0) || any (VF.grid$z.int < 0))
stop("error: the VF values should always be positive")
## check input of A.grid
gn <- names(A.grid)
if (! "x.int" %in% gn)
stop("error: A.grid should be a list that contains 'x.int'")
if (! "y.int" %in% gn)
stop("error: A.grid should be a list that contains 'y.int'")
if (! "z.int" %in% gn)
stop("error: A.grid should be a list that contains 'z.int'")
if (! "x.mid" %in% gn)
stop("error: A.grid should be a list that contains 'x.mid'")
if (! "y.mid" %in% gn)
stop("error: A.grid should be a list that contains 'y.mid'")
if (! "z.mid" %in% gn)
stop("error: A.grid should be a list that contains 'z.mid'")
if (is.null(A.grid$x.int) || is.null(A.grid$y.int)
|| is.null(A.grid$x.mid) || is.null(A.grid$y.mid)
|| is.null(A.grid$z.mid) || is.null(A.grid$z.int))
stop("error: the VF should contain (numeric) values")
if (any (A.grid$x.int < 0) || any (A.grid$y.int < 0)
|| any (A.grid$x.mid < 0) || any (A.grid$y.mid < 0)
|| any (A.grid$z.mid < 0) || any (A.grid$z.int < 0))
stop("error: the A values should always be positive")
}
## FUNCTION BODY: CALCULATIONS
## Impose boundary flux at upstream x-boundary when needed
## Default boundary condition is no gradient
if (! is.null (flux.x.up[1])) {
nom <- flux.x.up + VF.grid$x.int[1,,]*(D.grid$x.int[1,,]/grid$dx.aux[1] +
(1-AFDW.grid$x.int[1,,])*v.grid$x.int[1,,])*C[1,,]
denom <- VF.grid$x.int[1,,]*(D.grid$x.int[1,,]/grid$dx.aux[1]+
AFDW.grid$x.int[1,,]*v.grid$x.int[1,,])
C.x.up <- nom/denom
}
## Impose boundary flux at downstream x-boundary when needed
## Default boundary condition is no gradient
if (! is.null (flux.x.down[1])) {
nom <- flux.x.down - VF.grid$x.int[(Nx+1),,]*(D.grid$x.int[(Nx+1),,]/
grid$dx.aux[Nx+1] + AFDW.grid$x.int[(Nx+1),,]*v.grid$x.int[(Nx+1),,])*C[Nx,,]
denom <- -VF.grid$x.int[(Nx+1),,]*(D.grid$x.int[(Nx+1),,]/grid$dx.aux[Nx+1]+
(1-AFDW.grid$x.int[(Nx+1),,])*v.grid$x.int[(Nx+1),,])
C.x.down <- nom/denom
}
# Impose boundary flux at upstream y-boundary when needed
# Default boundary condition is no gradient
if (! is.null (flux.y.up[1])) {
nom <- flux.y.up + VF.grid$y.int[,1,]*(D.grid$y.int[,1,]/grid$dy.aux[1] +
(1-AFDW.grid$y.int[,1,])*v.grid$y.int[,1,])*C[,1,]
denom <- VF.grid$y.int[,1,]*(D.grid$y.int[,1,]/grid$dy.aux[1]+
AFDW.grid$y.int[,1,]*v.grid$y.int[,1,])
C.y.up <- nom/denom
}
# Impose boundary flux at downstream y-boundary when needed
# Default boundary condition is no gradient
if (! is.null (flux.y.down[1])) {
nom <- flux.y.down - VF.grid$y.int[,(Ny+1),]*(D.grid$y.int[,(Ny+1),]/
grid$dy.aux[Ny+1] + AFDW.grid$y.int[,(Ny+1),]*v.grid$y.int[,(Ny+1),])*C[,Ny,]
denom <- -VF.grid$y.int[,(Ny+1),]*(D.grid$y.int[,(Ny+1),]/grid$dy.aux[Ny+1]+
(1-AFDW.grid$y.int[,(Ny+1),])*v.grid$y.int[,(Ny+1),])
C.y.down <- nom/denom
}
# Impose boundary flux at upstream z-boundary when needed
# Default boundary condition is no gradient
if (! is.null (flux.z.up[1])) {
nom <- flux.z.up + VF.grid$z.int[,,1]*(D.grid$z.int[,,1]/grid$dz.aux[1] +
(1-AFDW.grid$z.int[,,1])*v.grid$z.int[,,1])*C[,,1]
denom <- VF.grid$z.int[,,1]*(D.grid$z.int[,,1]/grid$dz.aux[1]+
AFDW.grid$z.int[,,1]*v.grid$z.int[,,1])
C.z.up <- nom/denom
}
# Impose boundary flux at downstream z-boundary when needed
# Default boundary condition is no gradient
if (! is.null (flux.z.down[1])) {
nom <- flux.z.down - VF.grid$z.int[,,(Nz+1)]*(D.grid$z.int[,,(Nz+1)]/
grid$dz.aux[Nz+1] + AFDW.grid$z.int[,,(Nz+1)]*v.grid$z.int[,,(Nz+1)])*C[,,Nz]
denom <- -VF.grid$z.int[,,(Nz+1)]*(D.grid$z.int[,,(Nz+1)]/grid$dz.aux[Nz+1]+
(1-AFDW.grid$z.int[,,(Nz+1)])*v.grid$z.int[,,(Nz+1)])
C.z.down <- nom/denom
}
## when upper boundary layer is present, calculate new C.x.up
if (!is.null(a.bl.x.up[1]) & !is.null(C.x.up[1])) {
nom <- a.bl.x.up*C.x.up + VF.grid$x.int[1,,]*(D.grid$x.int[1,,]/
grid$dx.aux[1] + (1-AFDW.grid$x.int[1,,])*v.grid$x.int[1,,])*C[1,,]
denom <- a.bl.x.up + VF.grid$x.int[1,,]*(D.grid$x.int[1,,]/grid$dx.aux[1]+
AFDW.grid$x.int[1,,]*v.grid$x.int[1,,])
C.x.up <- nom/denom
}
## when lower boundary layer is present, calculate new C.x.down
if (!is.null(a.bl.x.down[1]) & !is.null(C.x.down[1])) {
nom <- a.bl.x.down*C.x.down + VF.grid$x.int[(Nx+1),,]*(D.grid$x.int[(Nx+1),,]/
grid$dx.aux[(Nx+1)] + (1-AFDW.grid$x.int[(Nx+1),,])*
v.grid$x.int[(Nx+1),,])*C[Nx,,]
denom <- a.bl.x.down + VF.grid$x.int[(Nx+1),,]*(D.grid$x.int[(Nx+1),,]/
grid$dx.aux[(Nx+1)]+ AFDW.grid$x.int[(Nx+1),,]*v.grid$x.int[(Nx+1),,])
C.x.down <- nom/denom
}
## when upper y boundary layer is present, calculate new C.y.up
if (!is.null(a.bl.y.up[1]) & !is.null(C.y.up[1])) {
nom <- a.bl.y.up*C.y.up + VF.grid$y.int[,1,]*(D.grid$y.int[,1,]/
grid$dy.aux[1] + (1-AFDW.grid$y.int[,1,])*v.grid$y.int[,1,])*C[,1,]
denom <- a.bl.y.up + VF.grid$y.int[,1,]*(D.grid$y.int[,1,]/grid$dy.aux[1]+
AFDW.grid$y.int[,1,]*v.grid$y.int[,1,])
C.y.up <- nom/denom
}
## when lower y boundary layer is present, calculate new C.y.down
if (!is.null(a.bl.y.down[1]) & !is.null(C.y.down[1])) {
nom <- a.bl.y.down*C.y.down + VF.grid$y.int[,(Ny+1),]*
(D.grid$y.int[,(Ny+1),]/grid$dy.aux[(Ny+1)] +
(1-AFDW.grid$y.int[,(Ny+1),])*v.grid$y.int[,(Ny+1),])*C[,Ny,]
denom <- a.bl.y.down + VF.grid$y.int[,(Ny+1),]*(D.grid$y.int[,(Ny+1),]/
grid$dy.aux[(Ny+1)]+ AFDW.grid$y.int[,(Ny+1),]*v.grid$y.int[,(Ny+1),])
C.y.down <- nom/denom
}
## when upper z boundary layer is present, calculate new C.z.up
if (!is.null(a.bl.z.up[1]) & !is.null(C.z.up[1])) {
nom <- a.bl.z.up*C.z.up + VF.grid$z.int[,1,]*(D.grid$z.int[,1,]/
grid$dz.aux[1] + (1-AFDW.grid$z.int[,1,])*v.grid$z.int[,1,])*C[,1,]
denom <- a.bl.z.up + VF.grid$z.int[,1,]*(D.grid$z.int[,1,]/grid$dz.aux[1]+
AFDW.grid$z.int[,1,]*v.grid$z.int[,1,])
C.z.up <- nom/denom
}
## when lower z boundary layer is present, calculate new C.z.down
if (!is.null(a.bl.z.down[1]) & !is.null(C.z.down[1])) {
nom <- a.bl.z.down*C.z.down + VF.grid$z.int[,(Nz+1),]*
(D.grid$z.int[,,(Nz+1)]/grid$dz.aux[(Nz+1)] +
(1-AFDW.grid$z.int[,,(Nz+1)])*v.grid$z.int[,,(Nz+1)])*C[,,Nz]
denom <- a.bl.z.down + VF.grid$z.int[,,(Nz+1)]*(D.grid$z.int[,,(Nz+1)]/
grid$dz.aux[(Nz+1)]+ AFDW.grid$z.int[,,(Nz+1)]*v.grid$z.int[,,(Nz+1)])
C.z.down <- nom/denom
}
## Calculate diffusive part of the flux
# Nx = 10, Ny = 2, Nz = 3
# DX <- array(dim = c(Nx, Ny, Nz), dx.aux)
# DY <- aperm(array(dim = c(Ny,Nx,Nz), dy.aux),c(2,1,3))
# DZ <- aperm(array(dim=c(Nz,Ny,Nx),dz.aux),c(3,2,1))
x.Dif.flux <- -VF.grid$x.int * D.grid$x.int *
Diff3D(mbind(C.x.up, C, C.x.down, along=1), along=1)/
array(dim=c(Nx+1,Ny,Nz),data=grid$dx.aux)
y.Dif.flux <- -VF.grid$y.int * D.grid$y.int *
Diff3D(mbind(C.y.up, C, C.y.down, along=2), along=2)/
aperm(array(data=grid$dy.aux,dim=c(Ny+1,Nx,Nz)),c(2,1,3))
z.Dif.flux <- -VF.grid$z.int * D.grid$z.int *
Diff3D(mbind(C.z.up, C, C.z.down, along=3), along=3)/
aperm(array(data=grid$dz.aux,dim=c(Nz+1,Ny,Nx)),c(3,2,1))
## Calculate advective part of the flux
x.Adv.flux <- 0
if (any(v.grid$x.int >0) ) {
vv <- v.grid$x.int
vv[vv<0]<-0
x.Adv.flux <- x.Adv.flux + VF.grid$x.int * 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 (v.grid$x.int < 0)) {
vv <- v.grid$x.int
vv[vv>0]<-0
x.Adv.flux <- x.Adv.flux + VF.grid$x.int * 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.flux <- 0
if (any(v.grid$y.int >0) ) {
vv <- v.grid$y.int
vv[vv<0]<-0
y.Adv.flux <- y.Adv.flux + VF.grid$y.int * 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 (v.grid$y.int < 0)) {
vv <- v.grid$y.int
vv[vv>0]<-0
y.Adv.flux <- y.Adv.flux + VF.grid$y.int * 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.flux <- 0
if (any(v.grid$z.int >0) ) {
vv <- v.grid$z.int
vv[vv<0]<-0
z.Adv.flux <- z.Adv.flux + VF.grid$z.int * 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 (v.grid$z.int < 0)) {
vv <- v.grid$z.int
vv[vv>0]<-0
z.Adv.flux <- z.Adv.flux + VF.grid$z.int * 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.flux <- x.Dif.flux + x.Adv.flux
y.flux <- y.Dif.flux + y.Adv.flux
z.flux <- z.Dif.flux + z.Adv.flux
## Impose boundary fluxes when needed
## Default boundary condition is no gradient
if (! is.null (flux.x.up[1]))
x.flux[1,,] <- flux.x.up
if (! is.null (flux.x.down[1]))
x.flux[dim(x.flux)[1],,] <- flux.x.down
if (! is.null (flux.y.up[1]))
y.flux[,1,] <- flux.y.up
if (! is.null (flux.y.down[1]))
y.flux[,dim(y.flux)[2],] <- flux.y.down
if (! is.null (flux.z.up[1]))
z.flux[,,1] <- flux.z.up
if (! is.null (flux.z.down[1]))
z.flux[,,dim(z.flux)[3]] <- flux.z.down
## Calculate rate of change = flux gradient NOG DOEN
dFdx <- - (Diff3D(A.grid$x.int*x.flux,along=1)/ A.grid$x.mid/grid$dx) / VF.grid$x.mid
dFdy <- - (Diff3D(A.grid$y.int*y.flux,along=2)/ A.grid$y.mid/grid$dy) / VF.grid$y.mid
dFdz <- - (Diff3D(A.grid$z.int*z.flux,along=3)/ A.grid$z.mid/grid$dz) / VF.grid$z.mid
if (!full.output) {
return (list (dC = dFdx + dFdy + dFdz, # Rate of change due to advective-diffuisve transport in each grid cell
flux.x.up = x.flux[1,,], # flux across lower boundary interface; positive = IN
flux.x.down = x.flux[dim(x.flux)[1],,], # flux across lower boundary interface; positive = OUT
flux.y.up = y.flux[,1,], # flux across lower boundary interface; positive = IN
flux.y.down = y.flux[,dim(y.flux)[2],], # flux across lower boundary interface; positive = OUT
flux.z.up = z.flux[,,1], # flux across lower boundary interface; positive = IN
flux.z.down = z.flux[,,dim(z.flux)[3]])) # flux across lower boundary interface; positive = OUT
} else {
return (list (dC = dFdx + dFdy + dFdz, # Rate of change due to advective-diffuisve transport in each grid cell
C.x.up = C.x.up, # concentration at upper interface
C.x.down = C.x.down, # concentration at upper interface
C.y.up = C.y.up, # concentration at upper interface
C.y.down = C.y.down, # concentration at upper interface
C.z.up = C.z.up, # concentration at upper interface
C.z.down = C.z.down, # concentration at upper interface
x.flux = x.flux, # flux across at the interface of each grid cell
y.flux = y.flux, # flux across at the interface of each grid cell
z.flux = z.flux, # flux across at the interface of each grid cell
flux.x.up = x.flux[1,,], # flux across lower boundary interface; positive = IN
flux.x.down = x.flux[dim(x.flux)[1],,], # flux across lower boundary interface; positive = OUT
flux.y.up = y.flux[,1,], # flux across lower boundary interface; positive = IN
flux.y.down = y.flux[,dim(y.flux)[2],], # flux across lower boundary interface; positive = OUT
flux.z.up = z.flux[,,1], # flux across lower boundary interface; positive = IN
flux.z.down = z.flux[,,dim(z.flux)[3]])) # flux across lower boundary interface; positive = OUT
}
} # end tran.3D
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.