R: General Three-Dimensional Advective-Diffusive Transport
tran.3D
R Documentation
General Three-Dimensional Advective-Diffusive Transport
Description
Estimates the transport term (i.e. the rate of change of a concentration
due to diffusion and advection) in a three-dimensional rectangular
model domain.
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 x-direction;
matrix of dimensions Ny*Nz [M/L3].
C.x.down
concentration at downstream boundary in x-direction;
matrix of dimensions Ny*Nz [M/L3].
C.y.up
concentration at upstream boundary in y-direction;
matrix of dimensions Nx*Nz [M/L3].
C.y.down
concentration at downstream boundary in y-direction;
matrix of dimensions Nx*Nz [M/L3].
C.z.up
concentration at upstream boundary in z-direction;
matrix of dimensions Nx*Ny [M/L3].
C.z.down
concentration at downstream boundary in z-direction;
matrix of dimensions Nx*Ny [M/L3].
flux.x.up
flux across the upstream boundary in x-direction,
positive = INTO model domain; matrix of dimensions Ny*Nz [M/L2/T].
flux.x.down
flux across the downstream boundary in x-direction,
positive = OUT of model domain; matrix of dimensions Ny*Nz [M/L2/T].
flux.y.up
flux across the upstream boundary in y-direction,
positive = INTO model domain; matrix of dimensions Nx*Nz [M/L2/T].
flux.y.down
flux across the downstream boundary in y-direction,
positive = OUT of model domain; matrix of dimensions Nx*Nz [M/L2/T].
flux.z.up
flux across the upstream boundary in z-direction,
positive = INTO model domain; matrix of dimensions Nx*Ny [M/L2/T].
flux.z.down
flux across the downstream boundary in z-direction,
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 x-direction
Flux=a.bl.x.up*(C.x.up-C[1,,]). One value [L/T].
a.bl.x.down
transfer coefficient across the downstream boundary
layer in x-direction;
Flux=a.bl.x.down*(C[Nx,,]-C.x.down).
One value [L/T].
a.bl.y.up
transfer coefficient across the upstream boundary layer.
in y-direction
Flux=a.bl.y.up*(C.y.up-C[,1,]). One value [L/T].
a.bl.y.down
transfer coefficient across the downstream boundary
layer in y-direction;
Flux=a.bl.y.down*(C[,Ny,]-C.y.down).
One value [L/T].
a.bl.z.up
transfer coefficient across the upstream boundary layer.
in y-direction
Flux=a.bl.y.up*(C.y.up-C[,,1]). One value [L/T].
a.bl.z.down
transfer coefficient across the downstream boundary
layer in z-direction;
Flux=a.bl.z.down*(C[,,Nz]-C.z.down).
One value [L/T].
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 z-direction, and with dimensions
(Nx+1)*Ny*Nz, Nx*(Ny+1)*Nz and Nx*Ny*(Nz+1) respectively. [L2/T].
D.x
diffusion coefficient in x-direction, 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 y-direction, 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 z-direction, 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 z-direction, 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 x-direction, 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 y-direction, 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 z-direction, 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 x-direction, 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 z-direction, 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 x-direction, defined on grid cell interfaces; backward = 1,
centred = 0.5, forward = 0; default is backward.
One value, a vector of length (Nx+1),
a prop.1D list created by setup.prop.1D,
or a (Nx+1)*Ny*Nz array [-].
AFDW.y
weight used in the finite difference scheme for advection
in the y-direction, defined on grid cell interfaces; backward = 1,
centred = 0.5, forward = 0; default is backward.
One value, a vector of length (Ny+1),
a prop.1D list created by setup.prop.1D,
or a Nx*(Ny+1)*Nz array [-].
AFDW.z
weight used in the finite difference scheme for advection
in the z-direction, defined on grid cell interfaces; backward = 1,
centred = 0.5, forward = 0; default is backward.
One value, a vector of length (Nz+1),
a prop.1D list created by setup.prop.1D,
or a Nx*Ny*(Nz+1) array [-].
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 z-direction, 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 x-direction.
One value, a vector of length (Nx+1),
a prop.1D list created by setup.prop.1D,
or a (Nx+1)*Ny*Nz array [-].
VF.y
Volume fraction at the grid cell interfaces in the y-direction.
One value, a vector of length (Ny+1),
a prop.1D list created by setup.prop.1D,
or a Nx*(Ny+1)*Nz array [-].
VF.z
Volume fraction at the grid cell interfaces in the z-direction.
One value, a vector of length (Nz+1),
a prop.1D list created by setup.prop.1D,
or a Nx*Ny*(Nz+1) array [-].
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 z-direction, 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 x-direction. One value, a vector of length (Nx+1),
a prop.1D list created by setup.prop.1D,
or a (Nx+1)*Ny*Nz array [L2].
A.y
Interface area defined at the grid cell interfaces in
the y-direction. One value, a vector of length (Ny+1),
a prop.1D list created by setup.prop.1D,
or a Nx*(Ny+1)*Nz array [L2].
A.z
Interface area defined at the grid cell interfaces in
the z-direction. One value, a vector of length (Nz+1),
a prop.1D list created by setup.prop.1D,
or a Nx*Ny*(Nz+1) array [L2].
dx
distance between adjacent cell interfaces in the x-direction
(thickness of grid cells). One value or vector of length Nx [L].
dy
distance between adjacent cell interfaces in the y-direction
(thickness of grid cells). One value or vector of length Ny [L].
dz
distance between adjacent cell interfaces in the z-direction
(thickness of grid cells). One value or vector of length Nz [L].
grid
discretization grid, a list containing at least elements
dx, dx.aux, dy, dy.aux, dz, dz.aux
(see setup.grid.2D) [L].
full.check
logical flag enabling a full check of the consistency
of the arguments (default = FALSE; TRUE slows down
execution by 50 percent).
full.output
logical flag enabling a full return of the output
(default = FALSE; TRUE slows down execution by 20 percent).
Details
Do not use this with too large grid.
The boundary conditions are either
(1) zero-gradient
(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.
Value
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 x-direction.
A matrix of dimension Ny*Nz [M/L3]. Only when full.output = TRUE.
C.x.down
concentration at the downstream interface in x-direction.
A matrix of dimension Ny*Nz [M/L3]. Only when full.output = TRUE.
C.y.up
concentration at the upstream interface in y-direction.
A matrix of dimension Nx*Nz [M/L3]. Only when full.output = TRUE.
C.y.down
concentration at the downstream interface in y-direction.
A matrix of dimension Nx*Nz [M/L3]. Only when full.output = TRUE.
C.z.up
concentration at the upstream interface in z-direction.
A matrix of dimension Nx*Ny [M/L3]. Only when full.output = TRUE.
C.z.down
concentration at the downstream interface in z-direction.
A matrix of dimension Nx*Ny [M/L3]. Only when full.output = TRUE.
x.flux
flux across the interfaces in x-direction of the grid cells.
A (Nx+1)*Ny*Nz array [M/L2/T]. Only when full.output = TRUE.
y.flux
flux across the interfaces in y-direction of the grid cells.
A Nx*(Ny+1)*Nz array [M/L2/T]. Only when full.output = TRUE.
z.flux
flux across the interfaces in z-direction of the grid cells.
A Nx*Ny*(Nz+1) array [M/L2/T]. Only when full.output = TRUE.
flux.x.up
flux across the upstream boundary in x-direction,
positive = INTO model domain. A matrix of dimension Ny*Nz [M/L2/T].
flux.x.down
flux across the downstream boundary in x-direction,
positive = OUT of model domain. A matrix of dimension Ny*Nz [M/L2/T].
flux.y.up
flux across the upstream boundary in y-direction,
positive = INTO model domain. A matrix of dimension Nx*Nz [M/L2/T].
flux.y.down
flux across the downstream boundary in y-direction,
positive = OUT of model domain. A matrix of dimension Nx*Nz [M/L2/T].
flux.z.up
flux across the upstream boundary in z-direction,
positive = INTO model domain. A matrix of dimension Nx*Ny [M/L2/T].
flux.z.down
flux across the downstream boundary in z-direction,
positive = OUT of model domain. A matrix of dimension Nx*Ny [M/L2/T].
Author(s)
Filip Meysman <filip.meysman@nioz.nl>,
Karline Soetaert <karline.soetaert@nioz.nl>
References
Soetaert and Herman, a practical guide to ecological modelling - using R as
a simulation platform, 2009. Springer
See Also
tran.cylindrical, tran.spherical
for a discretisation of 3-D transport equations in cylindrical and
spherical coordinates
tran.1D, tran.2D
Examples
## =============================================================================
## Diffusion in 3-D; imposed boundary conditions
## =============================================================================
diffusion3D <- function(t, Y, par) {
yy <- array(dim = c(n, n, n), data = Y) # vector to 3-D 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 Y-direction
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 x-direction
image(ST3, mfrow = c(2, 2), add.contour = TRUE, legend = TRUE,
dimselect = list(x = c(1, 4, 8, 10)))
par(mfrow = pm)
Results
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(ReacTran)
Loading required package: rootSolve
Loading required package: deSolve
Attaching package: 'deSolve'
The following object is masked from 'package:graphics':
matplot
Loading required package: shape
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/ReacTran/tran.3D.Rd_%03d_medium.png", width=480, height=480)
> ### Name: tran.3D
> ### Title: General Three-Dimensional Advective-Diffusive Transport
> ### Aliases: tran.3D
> ### Keywords: utilities
>
> ### ** Examples
>
> ## =============================================================================
> ## Diffusion in 3-D; imposed boundary conditions
> ## =============================================================================
> diffusion3D <- function(t, Y, par) {
+
+ yy <- array(dim = c(n, n, n), data = Y) # vector to 3-D 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 Y-direction
> 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)
+ ))
[1] "Steady-state settings"
sparseType message
1 3D sparse 3-D jacobian, calculated internally
[1] "estimated number of nonzero elements: 6910"
[1] "estimated number of function calls: 1001"
[1] "number of species: 1"
[1] "dimensions: 10 10 10"
[1] "cyclic boundaries: 0 0 0"
mean residual derivative 4.73661e-07
[1] "precision at each steady state step"
[1] 1.105000e+01 4.736612e-07
[1] ""
[1] "--------------------"
[1] " Memory requirements"
[1] "--------------------"
par mess val
1 nnz the number of nonzero elements 6400
2 ngp the number of independent groups of state variables 12
3 nsp the length of the work array actually required. 82366
user system elapsed
0.080 0.292 0.371
>
> 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 x-direction
> image(ST3, mfrow = c(2, 2), add.contour = TRUE, legend = TRUE,
+ dimselect = list(x = c(1, 4, 8, 10)))
>
> par(mfrow = pm)
>
>
>
>
>
> dev.off()
null device
1
>