tran.2D: General Two-Dimensional Advective-Diffusive Transport

Description Usage Arguments Details Value Note Author(s) References See Also Examples

View source: R/tran.2D.R

Description

Estimates the transport term (i.e. the rate of change of a concentration due to diffusion and advection) in a two-dimensional model domain.

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
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, a.bl.x.down = NULL, 
         a.bl.y.up = NULL, a.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)

Arguments

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.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].

D.grid

diffusion coefficient defined on all grid cell interfaces. A prop.2D list created by setup.prop.2D [L2/T]. See last example for creating spatially-varying diffusion coefficients.

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- and y- 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).

Details

The boundary conditions are either

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, 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].

Note

It is much more efficient to use the grid input rather than vectors or single numbers.

Thus: to optimise the code, use setup.grid.2D to create the grid, and use setup.prop.2D to create D.grid, v.grid, AFDW.grid, VF.grid, and A.grid, even if the values are 1 or remain constant.

There is no provision (yet) to deal with cross-diffusion. Set D.x and D.y different only if cross-diffusion effects are unimportant.

Author(s)

Filip Meysman <[email protected]>, Karline Soetaert <[email protected]>

References

Soetaert and Herman, 2009. a practical guide to ecological modelling - using R as a simulation platform. Springer

See Also

tran.polar for a discretisation of 2-D transport equations in polar coordinates

tran.1D, tran.3D

Examples

  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
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
## =============================================================================
## 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(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, full.output = TRUE)

# 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.x.up <- rep(1, times = (y.N)) # fixed conc at boundary layer
C.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.x.up   = C.x.up, a.bl.x.up = a.bl, 
  C.x.down = C.x.up, a.bl.x.down = a.bl, 
  C.y.up   = C.y.up, a.bl.y.up = a.bl,
  C.y.down = C.y.up, a.bl.y.down = a.bl)

# Runtime test with and without argument checking

n.iterate <-500

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.x.up = C.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.x.up = C.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.x.up = C.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 - unefficient implementation
## =============================================================================

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(nrow = N, ncol = N, y)
 dCONC <- tran.2D(CONC, D.x = Dx, D.y = Dy, dx = dx, dy = dy)$dC + r * CONC

 return (list(dCONC))

}

# initial condition: 0 everywhere, except in central point
y <- matrix(nrow = N, ncol = 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(nrow = N, ncol = N, 
                  data = subset(out, time == tt))
  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 )



## =============================================================================
## A 2-D model with diffusion in x- and y direction and first-order
## consumption - more efficient implementation, specifying ALL 2-D grids
## =============================================================================

N     <- 51          # number of grid cells
Dy    <- Dx <- 0.1   # diffusion coeff, X- and Y-direction
r     <- 0.005       # consumption rate
ini   <- 1           # initial value at x=0

x.grid    <- setup.grid.1D(x.up = -5, x.down = 5, N = N)
y.grid    <- setup.grid.1D(x.up = -5, x.down = 5, N = N)
grid2D    <- setup.grid.2D(x.grid, y.grid)

D.grid    <- setup.prop.2D(value = Dx, y.value = Dy, grid = grid2D)
v.grid    <- setup.prop.2D(value = 0, grid = grid2D)
A.grid    <- setup.prop.2D(value = 1, grid = grid2D)
AFDW.grid <- setup.prop.2D(value = 1, grid = grid2D)
VF.grid   <- setup.prop.2D(value = 1, grid = grid2D)

# The model equations - using the grids

Diff2Db <- function (t, y, parms)  {

   CONC  <- matrix(nrow = N, ncol = N, data = y)

   dCONC <- tran.2D(CONC, grid = grid2D, D.grid = D.grid, 
      A.grid = A.grid, VF.grid = VF.grid, AFDW.grid = AFDW.grid, 
      v.grid = v.grid)$dC + r * CONC
  
  return (list(dCONC))
}

# initial condition: 0 everywhere, except in central point
y <- matrix(nrow = N, ncol = N, data = 0)
y[N2,N2] <- ini  # initial concentration in the central point...

# solve for 8 time units
times <- 0:8
outb <- ode.2D (y = y, func = Diff2Db, t = times, parms = NULL,
                dim = c(N, N), lrw = 160000)

image(outb, ask = FALSE, mfrow = c(3, 3), main = paste("time", times))

## =============================================================================
## Same 2-D model, but now with spatially-variable diffusion coefficients
## =============================================================================

N     <- 51          # number of grid cells
r     <- 0.005       # consumption rate
ini   <- 1           # initial value at x=0
N2    <- ceiling(N/2)

D.grid <- list()

# Diffusion on x-interfaces
D.grid$x.int <- matrix(nrow = N+1, ncol = N, data = runif(N*(N+1)))

# Diffusion on y-interfaces
D.grid$y.int <- matrix(nrow = N, ncol = N+1, data = runif(N*(N+1)))

dx <- 10/N
dy <- 10/N

# The model equations

Diff2Dc <- function (t, y, parms)  {

   CONC  <- matrix(nrow = N, ncol = N, data = y)

   dCONC <- tran.2D(CONC, dx = dx, dy = dy, D.grid = D.grid)$dC + r * CONC

  return (list(dCONC))
}

# initial condition: 0 everywhere, except in central point
y <- matrix(nrow = N, ncol = N, data = 0)
y[N2, N2] <- ini  # initial concentration in the central point...

# solve for 8 time units
times <- 0:8
outc <- ode.2D (y = y, func = Diff2Dc, t = times, parms = NULL,
                dim = c(N, N), lrw = 160000)

outtimes <- c(1, 3, 5, 7)
image(outc, ask = FALSE, mfrow = c(2, 2), main = paste("time", outtimes),
      legend = TRUE, add.contour = TRUE, subset = time %in% outtimes)

ReacTran documentation built on Aug. 15, 2017, 3:01 p.m.