https://github.com/cran/ReacTran
Raw File
Tip revision: 4945515d3d9eb03bb3f931bb41d2f124904f90bc authored by Karline Soetaert on 15 August 2017, 21:13:03 UTC
version 1.4.3.1
Tip revision: 4945515
ReacTran.R
### R code from vignette source 'ReacTran.rnw'

###################################################
### code chunk number 1: preliminaries
###################################################
library("ReacTran")
options(prompt = "> ")
options(width=75)


###################################################
### code chunk number 2: ReacTran.rnw:239-240
###################################################
 (grid <- setup.grid.1D(L = 100, dx.1 = 1, N = 50))


###################################################
### code chunk number 3: Grid
###################################################
 plot(grid)


###################################################
### code chunk number 4: grid
###################################################
 plot(grid)


###################################################
### code chunk number 5: ReacTran.rnw:293-296
###################################################
grid <- setup.grid.1D(L = 10, N = 100)
Db   <- setup.prop.1D(func = p.exp, grid = grid, 
                      y.0 = 5, y.inf = 0, x.L = 2)


###################################################
### code chunk number 6: ReacTran.rnw:305-309
###################################################
exp.inc <- function(x, y.0 = 1, y.inf = 0.5, x.L = 0, x.att = 1)
       return(1 - p.exp(x, y.0, y.inf, x.L, x.att))

VFsolid <- setup.prop.1D(func = exp.inc, grid = grid, y.0 = 0.9, y.inf = 0.7)


###################################################
### code chunk number 7: prop
###################################################
par(mfrow = c(1, 2))
plot(VFsolid, grid, xyswap = TRUE, type = "l", main = "1-porosity")
plot(Db, grid, xyswap = TRUE, type = "l", main = "Db")
par(mfrow = c(1, 1))


###################################################
### code chunk number 8: prop
###################################################
par(mfrow = c(1, 2))
plot(VFsolid, grid, xyswap = TRUE, type = "l", main = "1-porosity")
plot(Db, grid, xyswap = TRUE, type = "l", main = "Db")
par(mfrow = c(1, 1))


###################################################
### code chunk number 9: ReacTran.rnw:401-403
###################################################
tran.1D(C = 1:20, D = 0, flux.up = 1, v = 1, dx = 1)
tran.1D(C = 1:20, D = 0, flux.up = 1, v = 1, dx = 1, full.output = TRUE)


###################################################
### code chunk number 10: ReacTran.rnw:436-437
###################################################
parms <- c(F0 = 1, v = 1, k = 0.1, dx = 1)


###################################################
### code chunk number 11: ReacTran.rnw:439-452
###################################################
advModel <- function(t, C, parms) {
    
   with (as.list(parms), {

     Tran <- tran.1D(C = C, D = 0, flux.up = F0, v = v, dx = dx)
     Consumption <-  k*C
     dC   <- Tran$dC - Consumption
        
     return (list(dC = dC, Consumption= Consumption,
                  flux.up = Tran$flux.up, flux.down = Tran$flux.down))
      })

    }


###################################################
### code chunk number 12: ReacTran.rnw:469-471
###################################################
out <- steady.1D(func = advModel, y = runif(25), parms = parms,
                 nspec = 1, positive = TRUE)


###################################################
### code chunk number 13: ReacTran.rnw:492-493
###################################################
out


###################################################
### code chunk number 14: st1
###################################################
plot (out, xlab = "x", ylab = "Conc", main = "advection")


###################################################
### code chunk number 15: figst1
###################################################
plot (out, xlab = "x", ylab = "Conc", main = "advection")


###################################################
### code chunk number 16: ReacTran.rnw:525-526
###################################################
with (out, print(sum(Consumption)-(flux.up-flux.down)))


###################################################
### code chunk number 17: ReacTran.rnw:539-550
###################################################

Aggregate.Model <- function(time, O2, pars) {

  tran <- tran.1D(C = O2, C.down = C.ow.O2,
                  D = D.grid, A = A.grid,
                  VF = por.grid, dx = grid )

  reac <- - R.O2*(O2/(Ks+O2))
  return(list(dCdt = tran$dC + reac, reac = reac,
              flux.up = tran$flux.up, flux.down = tran$flux.down))
}


###################################################
### code chunk number 18: ReacTran.rnw:554-560
###################################################
C.ow.O2 <- 0.25     # concentration O2 water [micromol cm-3]
por     <- 0.8      # porosity
D       <- 400      # diffusion coefficient O2 [cm2 yr-1]
v       <- 0        # advective velocity [cm yr-1]
R.O2    <- 1000000  # O2 consumption rate [micromol cm-3 yr-1]
Ks      <- 0.005    # O2 saturation constant [micromol cm-3]


###################################################
### code chunk number 19: ReacTran.rnw:563-567
###################################################
R <- 0.025           # radius of the agggregate [cm]
N <- 100             # number of grid layers

grid <- setup.grid.1D(x.up = 0, L = R, N = N)


###################################################
### code chunk number 20: ReacTran.rnw:571-573
###################################################
por.grid <- setup.prop.1D(value = por, grid = grid)
D.grid <- setup.prop.1D(value = D, grid = grid)


###################################################
### code chunk number 21: ReacTran.rnw:583-586
###################################################
sphere.surf <- function (x)   4*pi*x^2

A.grid  <- setup.prop.1D(func = sphere.surf, grid = grid)


###################################################
### code chunk number 22: ReacTran.rnw:590-592
###################################################
O2.agg <- steady.1D (y = runif(N), func = Aggregate.Model, nspec = 1,
                     positive = TRUE, atol = 1e-10)


###################################################
### code chunk number 23: agg
###################################################
plot(O2.agg, grid = grid$x.mid, xlab = "distance from centre, cm", 
     ylab = "mmol/m3", 
     main = "Diffusion-reaction of O2 in a spherical aggregate")


###################################################
### code chunk number 24: figagg
###################################################
plot(O2.agg, grid = grid$x.mid, xlab = "distance from centre, cm", 
     ylab = "mmol/m3", 
     main = "Diffusion-reaction of O2 in a spherical aggregate")


###################################################
### code chunk number 25: ReacTran.rnw:620-622
###################################################
  O2.agg$flux.up
  O2.agg$flux.down


###################################################
### code chunk number 26: ReacTran.rnw:632-634
###################################################
 Volume <- A.grid$mid * grid$dx
(Consump <- - sum(O2.agg$reac * Volume * por.grid$mid))


###################################################
### code chunk number 27: ReacTran.rnw:638-639
###################################################
 (Fluxin <- - O2.agg$flux.down * A.grid$int[N+1])


###################################################
### code chunk number 28: ReacTran.rnw:644-645
###################################################
  Consump - Fluxin


###################################################
### code chunk number 29: ReacTran.rnw:727-738
###################################################
river.model <- function (t = 0, OC, pars = NULL)  {

  tran <- tran.volume.1D(C = OC, F.up = F.OC, F.lat = F.lat, 
                         Disp = Disp, flow = flow, V = Volume,
                         full.output = TRUE)
  reac <- - k*OC

  return(list(dCdt = tran$dC + reac,
            F.up = tran$F.up, F.down = tran$F.down,
            F.lat = tran$F.lat))
}


###################################################
### code chunk number 30: ReacTran.rnw:743-746
###################################################
nbox          <- 500                # number of grid cells
lengthEstuary <- 100000             # length of estuary [m]
BoxLength     <- lengthEstuary/nbox # [m]


###################################################
### code chunk number 31: ReacTran.rnw:751-756
###################################################
Distance      <- seq(BoxLength/2, by = BoxLength, len = nbox) # [m]

CrossArea <- 4000 + 72000 * Distance^5 /(Distance^5+50000^5)

Volume  <- CrossArea*BoxLength


###################################################
### code chunk number 32: ReacTran.rnw:761-763
###################################################
Disp    <- 1000   # m3/s, bulk dispersion coefficient
flow    <- 180    # m3/s, mean river flow


###################################################
### code chunk number 33: ReacTran.rnw:768-772
###################################################
F.OC    <- 180               # input organic carbon [mol s-1]
F.lat.0 <- F.OC              # lateral input organic carbon [mol s-1]

k       <- 10/(365*24*3600)  # decay constant organic carbon [s-1]


###################################################
### code chunk number 34: ReacTran.rnw:775-776
###################################################
F.lat <- rep(0, length.out = nbox)


###################################################
### code chunk number 35: ReacTran.rnw:781-783
###################################################
sol  <- steady.1D(y = runif(nbox), fun = river.model, nspec = 1,
                  atol = 1e-15, rtol = 1e-15, positive = TRUE)


###################################################
### code chunk number 36: ReacTran.rnw:786-791
###################################################
F.lat <- F.lat.0*dnorm(x = Distance/lengthEstuary,
                       mean = Distance[nbox/2]/lengthEstuary,
                       sd = 1/20, log = FALSE) /nbox
sol2 <- steady.1D(y = runif(nbox), fun = river.model, nspec = 1, 
                  atol = 1e-15, rtol = 1e-15, positive = TRUE)


###################################################
### code chunk number 37: ReacTran.rnw:799-801
###################################################
summary(sol)
summary(sol2)


###################################################
### code chunk number 38: est
###################################################
plot(sol, sol2, grid = Distance/1000, lwd = 2,
        main = "Organic carbon decay in an estuary",xlab = "distance [km]",
        ylab = "OC Concentration [mM]")
legend ("topright", col = 1:2, lwd = 2, c("baseline", "with lateral input"))


###################################################
### code chunk number 39: est
###################################################
plot(sol, sol2, grid = Distance/1000, lwd = 2,
        main = "Organic carbon decay in an estuary",xlab = "distance [km]",
        ylab = "OC Concentration [mM]")
legend ("topright", col = 1:2, lwd = 2, c("baseline", "with lateral input"))


###################################################
### code chunk number 40: ReacTran.rnw:824-825
###################################################
sum(sol$F.up) + sum(sol$F.lat) - sum(sol$F.down) - sum(sol$y*k*Volume)


###################################################
### code chunk number 41: ReacTran.rnw:886-894
###################################################
n     <- 100           # number of grid cells
dy    <- dx <- 100/n   # grid size

Dy    <- Dx <- 5   # diffusion coeff, X- and Y-direction
r     <- -0.02     # production/consumption rate
Bc    <- 300       # boundary concentration
irr   <- 20        # irrigation rate
vx    <- 1         # advection


###################################################
### code chunk number 42: ReacTran.rnw:898-899
###################################################
y  <- matrix(nrow = n, ncol = n, 0)


###################################################
### code chunk number 43: ReacTran.rnw:903-919
###################################################
Diff2D <- function (t, y, parms, N) {

  CONC <- matrix(nrow = N, ncol = N, y)
# Transport
  Tran    <-tran.2D(CONC, D.x = Dx, D.y = Dy, C.y.down = Bc,
                    dx = dx, dy = dy, v.x = vx)

# transport + reaction
  dCONC   <- Tran$dC + r*CONC

# Bioirrigation in a central spot
  mid <- N/2
  dCONC[mid, mid] <- dCONC[mid, mid]  + irr*(Bc - CONC[mid, mid])

  return (list(dCONC))
}


###################################################
### code chunk number 44: ReacTran.rnw:933-938
###################################################
print(system.time(
 std  <- steady.2D(func = Diff2D, y = as.vector(y), time = 0, N = n,
                   parms = NULL, lrw = 1000000, dimens = c(n, n),
                   nout = 0, positive = TRUE)
))


###################################################
### code chunk number 45: ReacTran.rnw:941-946
###################################################
times <- seq(0, 100, 5)
print(system.time(
  out2 <- ode.2D(func = Diff2D, y = as.vector(y), times = times, N = n,
                 parms = NULL, lrw = 10000000, dimens = c(n, n))
))


###################################################
### code chunk number 46: twod5
###################################################
mat <- matrix(nrow = n, ncol = n, subset(out2, time ==  20))
filled.contour(mat, zlim = c(0, Bc), color = femmecol,
               main = "after 20 time units")


###################################################
### code chunk number 47: twod
###################################################
mat <- matrix(nrow = n, ncol = n, subset(out2, time ==  20))
filled.contour(mat, zlim = c(0, Bc), color = femmecol,
               main = "after 20 time units")


###################################################
### code chunk number 48: twodst
###################################################
image (std, main = "steady-state", legend = TRUE)


###################################################
### code chunk number 49: twodst
###################################################
image (std, main = "steady-state", legend = TRUE)


back to top