tran.2D {ReacTran} | R Documentation |
Estimates the transport term (i.e. the rate of change of a concentration due to diffusion and advection) in a two-dimensional model domain.
tran.2D ( C, C.x.up=C[1,], C.x.down=C[nrow(C),], C.y.up=C[,1], C.y.down=C[,ncol(C)], flux.x.up=NULL, flux.x.down=NULL, flux.y.up=NULL, flux.y.down=NULL, a.bl.x.up=NULL, C.bl.x.up=NULL, a.bl.x.down=NULL, C.bl.x.down=NULL, a.bl.y.up=NULL, C.bl.y.up=NULL, a.bl.y.down=NULL, C.bl.y.down=NULL, D.grid=NULL, D.x=NULL, D.y=D.x, v.grid=NULL, v.x=0, v.y=0, AFDW.grid=NULL, AFDW.x=1, AFDW.y=AFDW.x, VF.grid=NULL,VF.x=1, VF.y=VF.x, A.grid=NULL, A.x=1, A.y=1, grid=NULL, dx=NULL, dy=NULL, full.check = FALSE, full.output = FALSE)
C |
concentration, expressed per unit volume, defined at the centre of each grid cell; Nx*Ny matrix [M/L3]. |
C.x.up |
concentration at upstream boundary in x-direction; vector of length Ny [M/L3]. |
C.x.down |
concentration at downstream boundary in x-direction; vector of length Ny [M/L3]. |
C.y.up |
concentration at upstream boundary in y-direction; vector of length Nx [M/L3]. |
C.y.down |
concentration at downstream boundary in y-direction; vector of length Nx [M/L3]. |
flux.x.up |
flux across the upstream boundary in x-direction, positive = INTO model domain; vector of length Ny [M/L2/T]. |
flux.x.down |
flux across the downstream boundary in x-direction, positive = OUT of model domain; vector of length Ny [M/L2/T]. |
flux.y.up |
flux across the upstream boundary in y-direction, positive = INTO model domain; vector of length Nx [M/L2/T]. |
flux.y.down |
flux across the downstream boundary in y-direction, positive = OUT of model domain; vector of length Nx [M/L2/T]. |
a.bl.x.up |
transfer coefficient across the upstream boundary layer.
in x-direction;
Flux=a.bl.x.up*(C.bl.x.up-C[1,]) . One value [L/T].
|
C.bl.x.up |
concentration at the upstream boundary layer in x-direction; vector of length Ny [M/L3]. |
a.bl.x.down |
transfer coefficient across the downstream boundary
layer in x-direction;
Flux=a.bl.x.down*(C[Nx,]-C.bl.x.down) .
One value [L/T].
|
C.bl.x.down |
concentration at the downstream boundary layer in x-direction ; vector of length Ny [M/L3]. |
a.bl.y.up |
transfer coefficient across the upstream boundary layer.
in y-direction;
Flux=a.bl.y.up*(C.bl.y.up-C[,1]) . One value [L/T].
|
C.bl.y.up |
concentration at the upstream boundary layer in y-direction; vector of length Nx [M/L3]. |
a.bl.y.down |
transfer coefficient across the downstream boundary
layer in y-direction;
Flux=a.bl.y.down*(C[,Ny]-C.bl.y.down) .
One value [L/T].
|
C.bl.y.down |
concentration at the downstream boundary layer in y-direction ; vector of length Nx [M/L3]. |
D.grid |
diffusion coefficient defined on all grid cell
interfaces. A prop.2D list created by setup.prop.2D [L2/T].
|
D.x |
diffusion coefficient in x-direction, defined on grid cell
interfaces. One value, a vector of length (Nx+1),
a prop.1D list created by setup.prop.1D ,
or a (Nx+1)* Ny matrix [L2/T].
|
D.y |
diffusion coefficient in y-direction, defined on grid cell
interfaces. One value, a vector of length (Ny+1),
a prop.1D list created by setup.prop.1D ,
or a Nx*(Ny+1) matrix [L2/T].
|
v.grid |
advective velocity defined on all grid cell
interfaces. Can be positive (downstream flow) or negative (upstream flow).
A prop.2D list created by setup.prop.2D [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),
a prop.1D list created by setup.prop.1D ,
or a (Nx+1)*Ny matrix [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),
a prop.1D list created by setup.prop.1D ,
or a Nx*(Ny+1) matrix [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.
A prop.2D list created by setup.prop.2D [-].
|
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 matrix [-].
|
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) matrix [-].
|
VF.grid |
Volume fraction. A prop.2D list created by
setup.prop.2D [-].
|
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 matrix [-].
|
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) matrix [-].
|
A.grid |
Interface area. A prop.2D list created by
setup.prop.2D [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 matrix [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) matrix [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]. |
grid |
discretization grid, a list containing at least elements
dx , dx.aux , dy , dy.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).
|
The boundary conditions are either
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, a Nx*Ny matrix. [M/L3/T]. |
C.x.up |
concentration at the upstream interface in x-direction.
A vector of length Ny [M/L3]. Only when full.output = TRUE .
|
C.x.down |
concentration at the downstream interface in x-direction.
A vector of length Ny [M/L3]. Only when full.output = TRUE .
|
C.y.up |
concentration at the the upstream interface in y-direction.
A vector of length Nx [M/L3]. Only when full.output = TRUE .
|
C.y.down |
concentration at the downstream interface in y-direction.
A vector of length Nx [M/L3]. Only when full.output = TRUE .
|
x.flux |
flux across the interfaces in x-direction of the grid cells.
A (Nx+1)*Ny matrix [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) matrix [M/L2/T]. Only when full.output = TRUE .
|
flux.x.up |
flux across the upstream boundary in x-direction, positive = INTO model domain. A vector of length Ny [M/L2/T]. |
flux.x.down |
flux across the downstream boundary in x-direction, positive = OUT of model domain. A vector of length Ny [M/L2/T]. |
flux.y.up |
flux across the upstream boundary in y-direction, positive = INTO model domain. A vector of length Nx [M/L2/T]. |
flux.y.down |
flux across the downstream boundary in y-direction, positive = OUT of model domain. A vector of length Nx [M/L2/T]. |
Filip Meysman <f.meysman@nioo.knaw.nl>, Karline Soetaert <k.soetaert@nioo.knaw.nl>
Soetaert and Herman, 2009. a practical guide to ecological modelling - using R as a simulation platform. Springer
## ============================================================================= ## Testing the functions ## ============================================================================= # Parameters F <- 100 # input flux [micromol cm-2 yr-1] por <- 0.8 # constant porosity D <- 400 # mixing coefficient [cm2 yr-1] v <- 1 # advective velocity [cm yr-1] # Grid definition x.N <- 4 # number of cells in x-direction y.N <- 6 # number of cells in y-direction x.L <- 8 # domain size x-direction [cm] y.L <- 24 # domain size y-direction [cm] dx <- x.L/x.N # cell size x-direction [cm] dy <- y.L/y.N # cell size y-direction [cm] # Intial conditions C <- matrix(nrow=x.N, ncol=y.N, data=0, byrow=FALSE) # Boundary conditions: fixed concentration C.x.up <- rep(1, times=y.N) C.x.down <- rep(0, times=y.N) C.y.up <- rep(1, times=x.N) C.y.down <- rep(0, times=x.N) # Only diffusion tran.2D(full.output=TRUE, C=C, D.x=D, D.y=D, v.x=0, v.y=0, VF.x=por, VF.y=por, dx=dx, dy=dy, C.x.up=C.x.up, C.x.down=C.x.down, C.y.up=C.y.up,C.y.down=C.y.down) # Strong advection, backward (default), central and forward #finite difference schemes tran.2D(C=C, D.x=D, v.x=100*v, VF.x=por, dx=dx, dy=dy, C.x.up=C.x.up, C.x.down=C.x.down, C.y.up=C.y.up, C.y.down=C.y.down) tran.2D(AFDW.x=0.5, C=C, D.x=D, v.x=100*v, VF.x=por, dx=dx, dy=dy, C.x.up=C.x.up, C.x.down=C.x.down, C.y.up=C.y.up, C.y.down=C.y.down) tran.2D(AFDW.x=0, C=C, D.x=D, v.x=100*v, VF.x=por, dx=dx, dy=dy, C.x.up=C.x.up, C.x.down=C.x.down, C.y.up=C.y.up, C.y.down=C.y.down) # Boundary conditions: fixed fluxes flux.x.up <- rep(200, times=y.N) flux.x.down <- rep(-200, times=y.N) flux.y.up <- rep(200, times=x.N) flux.y.down <- rep(-200, times=x.N) tran.2D(C=C, D.x=D, v.x=0, VF.x=por, dx=dx, dy=dy, flux.x.up=flux.x.up, flux.x.down=flux.x.down, flux.y.up=flux.y.up, flux.y.down=flux.y.down) # Boundary conditions: convective boundary layer on all sides a.bl <- 800 # transfer coefficient C.bl.x.up <- rep(1, times=(y.N)) # fixed conc at boundary layer C.bl.y.up <- rep(1, times=(x.N)) # fixed conc at boundary layer tran.2D(full.output=TRUE, C=C, D.x=D, v.x=0, VF.x=por, dx=dx, dy=dy, C.bl.x.up=C.bl.x.up, a.bl.x.up=a.bl, C.bl.x.down=C.bl.x.up, a.bl.x.down=a.bl, C.bl.y.up=C.bl.y.up, a.bl.y.up=a.bl, C.bl.y.down=C.bl.y.up, a.bl.y.down=a.bl) # Runtime test with and without argument checking n.iterate <-1000 test1 <- function() { for (i in 1:n.iterate ) ST<-tran.2D(full.check=TRUE,C=C,D.x=D,v.x=0,VF.x=por, dx=dx,dy=dy,C.bl.x.up=C.bl.x.up,a.bl.x.up=a.bl,C.x.down=C.x.down) } system.time(test1()) test2 <- function() { for (i in 1:n.iterate ) ST<-tran.2D(full.output=TRUE,C=C,D.x=D,v.x=0,VF.x=por, dx=dx,dy=dy,C.bl.x.up=C.bl.x.up,a.bl.x.up=a.bl,C.x.down=C.x.down) } system.time(test2()) test3 <- function() { for (i in 1:n.iterate ) ST<-tran.2D(full.output=TRUE,full.check=TRUE,C=C,D.x=D,v.x=0, VF.x=por,dx=dx,dy=dy,C.bl.x.up=C.bl.x.up,a.bl.x.up=a.bl,C.x.down=C.x.down) } system.time(test3()) ## ============================================================================= ## A 2-D model with diffusion in x- and y direction and first-order ## consumption ## ============================================================================= N <- 51 # number of grid cells XX <- 10 # total size dy <- dx <- XX/N # grid size Dy <- Dx <- 0.1 # diffusion coeff, X- and Y-direction r <- 0.005 # consumption rate ini <- 1 # initial value at x=0 N2 <- ceiling(N/2) X <- seq (dx,by=dx,len=(N2-1)) X <- c(-rev(X),0,X) # The model equations Diff2D <- function (t,y,parms) { CONC <- matrix(nr=N,nc=N,y) dCONC <- tran.2D(CONC, D.x=Dx, D.y=Dy, dx=dx, dy=dy)$dC + r * CONC return (list(as.vector(dCONC))) } # initial condition: 0 everywhere, except in central point y <- matrix(nr=N,nc=N,data=0) y[N2,N2] <- ini # initial concentration in the central point... # solve for 10 time units times <- 0:10 out <- ode.2D (y=y, func=Diff2D, t=times, parms=NULL, dim = c(N,N), lrw = 160000) pm <- par (mfrow=c(2,2)) # Compare solution with analytical solution... for (i in seq(2,11,by=3)) { tt <- times[i] mat <- matrix(nr=N,nc=N,out[i,-1]) plot(X,mat[N2,],type="l",main=paste("time=",times[i]), ylab="Conc",col="red") ana <- ini*dx^2/(4*pi*Dx*tt)*exp(r*tt-X^2/(4*Dx*tt)) points(X,ana,pch="+") } legend ("bottom", col=c("red","black"), lty=c(1,NA), pch=c(NA,"+"), c("tran.2D","exact")) par("mfrow"=pm )