Estimates the transport term (i.e. the rate of change of a concentration due to diffusion and advection) in a threedimensional rectangular model domain.
Do not use with too many boxes!
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16  tran.3D (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)

C 
concentration, expressed per unit volume, defined at the centre of each grid cell; Nx*Ny*Nz array [M/L3]. 
C.x.up 
concentration at upstream boundary in xdirection; matrix of dimensions Ny*Nz [M/L3]. 
C.x.down 
concentration at downstream boundary in xdirection; matrix of dimensions Ny*Nz [M/L3]. 
C.y.up 
concentration at upstream boundary in ydirection; matrix of dimensions Nx*Nz [M/L3]. 
C.y.down 
concentration at downstream boundary in ydirection; matrix of dimensions Nx*Nz [M/L3]. 
C.z.up 
concentration at upstream boundary in zdirection; matrix of dimensions Nx*Ny [M/L3]. 
C.z.down 
concentration at downstream boundary in zdirection; matrix of dimensions Nx*Ny [M/L3]. 
flux.x.up 
flux across the upstream boundary in xdirection, positive = INTO model domain; matrix of dimensions Ny*Nz [M/L2/T]. 
flux.x.down 
flux across the downstream boundary in xdirection, positive = OUT of model domain; matrix of dimensions Ny*Nz [M/L2/T]. 
flux.y.up 
flux across the upstream boundary in ydirection, positive = INTO model domain; matrix of dimensions Nx*Nz [M/L2/T]. 
flux.y.down 
flux across the downstream boundary in ydirection, positive = OUT of model domain; matrix of dimensions Nx*Nz [M/L2/T]. 
flux.z.up 
flux across the upstream boundary in zdirection, positive = INTO model domain; matrix of dimensions Nx*Ny [M/L2/T]. 
flux.z.down 
flux across the downstream boundary in zdirection, positive = OUT of model domain; matrix of dimensions Nx*Ny [M/L2/T]. 
a.bl.x.up 
transfer coefficient across the upstream boundary layer. in xdirection

a.bl.x.down 
transfer coefficient across the downstream boundary layer in xdirection;

a.bl.y.up 
transfer coefficient across the upstream boundary layer. in ydirection

a.bl.y.down 
transfer coefficient across the downstream boundary layer in ydirection;

a.bl.z.up 
transfer coefficient across the upstream boundary layer. in ydirection

a.bl.z.down 
transfer coefficient across the downstream boundary layer in zdirection;

D.grid 
diffusion coefficient defined on all grid cell interfaces. Should contain elements x.int, y.int, z.int, arrays with the values on the interfaces in x, y and zdirection, and with dimensions (Nx+1)*Ny*Nz, Nx*(Ny+1)*Nz and Nx*Ny*(Nz+1) respectively. [L2/T]. 
D.x 
diffusion coefficient in xdirection, defined on grid cell interfaces. One value, a vector of length (Nx+1), or a (Nx+1)* Ny *Nz array [L2/T]. 
D.y 
diffusion coefficient in ydirection, defined on grid cell interfaces. One value, a vector of length (Ny+1), or a Nx*(Ny+1)*Nz array [L2/T]. 
D.z 
diffusion coefficient in zdirection, defined on grid cell interfaces. One value, a vector of length (Nz+1), or a Nx*Ny*(Nz+1) array [L2/T]. 
v.grid 
advective velocity defined on all grid cell interfaces. Can be positive (downstream flow) or negative (upstream flow). Should contain elements x.int, y.int, z.int, arrays with the values on the interfaces in x, y and zdirection, and with dimensions (Nx+1)*Ny*Nz, Nx*(Ny+1)*Nz and Nx*Ny*(Nz+1) respectively. [L/T]. 
v.x 
advective velocity in the xdirection, defined on grid cell interfaces. Can be positive (downstream flow) or negative (upstream flow). One value, a vector of length (Nx+1), or a (Nx+1)*Ny*Nz array [L/T]. 
v.y 
advective velocity in the ydirection, defined on grid cell interfaces. Can be positive (downstream flow) or negative (upstream flow). One value, a vector of length (Ny+1), or a Nx*(Ny+1)*Nz array [L/T]. 
v.z 
advective velocity in the zdirection, defined on grid cell interfaces. Can be positive (downstream flow) or negative (upstream flow). One value, a vector of length (Nz+1), or a Nx*Ny*(Nz+1) array [L/T]. 
AFDW.grid 
weight used in the finite difference scheme for advection in the xdirection, defined on grid cell interfaces; backward = 1, centred = 0.5, forward = 0; default is backward. Should contain elements x.int, y.int, z.int, arrays with the values on the interfaces in x, y and zdirection, and with dimensions (Nx+1)*Ny*Nz, Nx*(Ny+1)*Nz and Nx*Ny*(Nz+1) respectively. []. 
AFDW.x 
weight used in the finite difference scheme for advection
in the xdirection, defined on grid cell interfaces; backward = 1,
centred = 0.5, forward = 0; default is backward.
One value, a vector of length (Nx+1),
a 
AFDW.y 
weight used in the finite difference scheme for advection
in the ydirection, defined on grid cell interfaces; backward = 1,
centred = 0.5, forward = 0; default is backward.
One value, a vector of length (Ny+1),
a 
AFDW.z 
weight used in the finite difference scheme for advection
in the zdirection, defined on grid cell interfaces; backward = 1,
centred = 0.5, forward = 0; default is backward.
One value, a vector of length (Nz+1),
a 
VF.grid 
Volume fraction. A list. Should contain elements x.int, y.int, z.int, arrays with the values on the interfaces in x, y and zdirection, and with dimensions (Nx+1)*Ny*Nz, Nx*(Ny+1)*Nz and Nx*Ny*(Nz+1) respectively. []. 
VF.x 
Volume fraction at the grid cell interfaces in the xdirection.
One value, a vector of length (Nx+1),
a 
VF.y 
Volume fraction at the grid cell interfaces in the ydirection.
One value, a vector of length (Ny+1),
a 
VF.z 
Volume fraction at the grid cell interfaces in the zdirection.
One value, a vector of length (Nz+1),
a 
A.grid 
Interface area, a list. Should contain elements x.int, y.int, z.int, arrays with the values on the interfaces in x, y and zdirection, and with dimensions (Nx+1)*Ny*Nz, Nx*(Ny+1)*Nz and Nx*Ny*(Nz+1) respectively. [L2]. 
A.x 
Interface area defined at the grid cell interfaces in
the xdirection. One value, a vector of length (Nx+1),
a 
A.y 
Interface area defined at the grid cell interfaces in
the ydirection. One value, a vector of length (Ny+1),
a 
A.z 
Interface area defined at the grid cell interfaces in
the zdirection. One value, a vector of length (Nz+1),
a 
dx 
distance between adjacent cell interfaces in the xdirection (thickness of grid cells). One value or vector of length Nx [L]. 
dy 
distance between adjacent cell interfaces in the ydirection (thickness of grid cells). One value or vector of length Ny [L]. 
dz 
distance between adjacent cell interfaces in the zdirection (thickness of grid cells). One value or vector of length Nz [L]. 
grid 
discretization grid, a list containing at least elements

full.check 
logical flag enabling a full check of the consistency
of the arguments (default = 
full.output 
logical flag enabling a full return of the output
(default = 
Do not use this with too large grid.
The boundary conditions are either
(1) zerogradient
(2) fixed concentration
(3) convective boundary layer
(4) fixed flux
This is also the order of priority. The zero gradient is the default, the fixed flux overrules all other.
a list containing:
dC 
the rate of change of the concentration C due to transport, defined in the centre of each grid cell, an array with dimension Nx*Ny*Nz [M/L3/T]. 
C.x.up 
concentration at the upstream interface in xdirection.
A matrix of dimension Ny*Nz [M/L3]. Only when 
C.x.down 
concentration at the downstream interface in xdirection.
A matrix of dimension Ny*Nz [M/L3]. Only when 
C.y.up 
concentration at the upstream interface in ydirection.
A matrix of dimension Nx*Nz [M/L3]. Only when 
C.y.down 
concentration at the downstream interface in ydirection.
A matrix of dimension Nx*Nz [M/L3]. Only when 
C.z.up 
concentration at the upstream interface in zdirection.
A matrix of dimension Nx*Ny [M/L3]. Only when 
C.z.down 
concentration at the downstream interface in zdirection.
A matrix of dimension Nx*Ny [M/L3]. Only when 
x.flux 
flux across the interfaces in xdirection of the grid cells.
A (Nx+1)*Ny*Nz array [M/L2/T]. Only when 
y.flux 
flux across the interfaces in ydirection of the grid cells.
A Nx*(Ny+1)*Nz array [M/L2/T]. Only when 
z.flux 
flux across the interfaces in zdirection of the grid cells.
A Nx*Ny*(Nz+1) array [M/L2/T]. Only when 
flux.x.up 
flux across the upstream boundary in xdirection, positive = INTO model domain. A matrix of dimension Ny*Nz [M/L2/T]. 
flux.x.down 
flux across the downstream boundary in xdirection, positive = OUT of model domain. A matrix of dimension Ny*Nz [M/L2/T]. 
flux.y.up 
flux across the upstream boundary in ydirection, positive = INTO model domain. A matrix of dimension Nx*Nz [M/L2/T]. 
flux.y.down 
flux across the downstream boundary in ydirection, positive = OUT of model domain. A matrix of dimension Nx*Nz [M/L2/T]. 
flux.z.up 
flux across the upstream boundary in zdirection, positive = INTO model domain. A matrix of dimension Nx*Ny [M/L2/T]. 
flux.z.down 
flux across the downstream boundary in zdirection, positive = OUT of model domain. A matrix of dimension Nx*Ny [M/L2/T]. 
Filip Meysman <filip.meysman@nioz.nl>, Karline Soetaert <karline.soetaert@nioz.nl>
Soetaert and Herman, a practical guide to ecological modelling  using R as a simulation platform, 2009. Springer
tran.cylindrical
, tran.spherical
for a discretisation of 3D transport equations in cylindrical and
spherical coordinates
tran.1D
, tran.2D
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40  ## =============================================================================
## Diffusion in 3D; imposed boundary conditions
## =============================================================================
diffusion3D < function(t, Y, par) {
yy < array(dim = c(n, n, n), data = Y) # vector to 3D array
dY < r * yy # consumption
BND < matrix(nrow = n, ncol = n, 1) # boundary concentration
dY < dY + tran.3D(C = yy,
C.x.up = BND, C.y.up = BND, C.z.up = BND,
C.x.down = BND, C.y.down = BND, C.z.down = BND,
D.x = Dx, D.y = Dy, D.z = Dz,
dx = dx, dy = dy, dz = dz, full.check = TRUE)$dC
return(list(dY))
}
# parameters
dy < dx < dz < 1 # grid size
Dy < Dx < Dz < 1 # diffusion coeff, X and Ydirection
r < 0.025 # consumption rate
n < 10
y < array(dim = c(n, n, n), data = 10.)
print(system.time(
ST3 < steady.3D(y, func = diffusion3D, parms = NULL,
pos = TRUE, dimens = c(n, n, n),
lrw = 2000000, verbose = TRUE)
))
pm < par(mfrow = c(1,1))
y < array(dim = c(n, n, n), data = ST3$y)
filled.contour(y[ , ,n/2], color.palette = terrain.colors)
# a selection in the xdirection
image(ST3, mfrow = c(2, 2), add.contour = TRUE, legend = TRUE,
dimselect = list(x = c(1, 4, 8, 10)))
par(mfrow = pm)

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.
All documentation is copyright its authors; we didn't write any of that.