### R code from vignette source 'PDE.Rnw'

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


###################################################
### code chunk number 2: PDE.Rnw:175-176
###################################################
Grid <- setup.grid.1D(N = 1000, L = 10)


###################################################
### code chunk number 3: PDE.Rnw:182-184
###################################################
r  <- setup.prop.1D(grid = Grid, func = function(r) r)
r2 <- setup.prop.1D(grid = Grid, func = function(r) r^2)


###################################################
### code chunk number 4: PDE.Rnw:194-199
###################################################
library(ReacTran)
pde1D <-function(t, C, parms, A = 1)  {
  tran   <- tran.1D (C = C, A = A, D = D, C.down = Cext, dx = Grid)$dC
  list(tran - Q)  # the return value: rate of change
}


###################################################
### code chunk number 5: PDE.Rnw:204-207
###################################################
D    <- 1    # diffusion constant
Q    <- 1    # uptake rate
Cext <- 20


###################################################
### code chunk number 6: PDE.Rnw:232-237
###################################################
library(rootSolve)
Cartesian   <- steady.1D(y = runif(Grid$N),
  func = pde1D, parms = NULL, nspec = 1, A = 1)
Cylindrical <- steady.1D(y = runif(Grid$N),
  func = pde1D, parms = NULL, nspec = 1, A = r)


###################################################
### code chunk number 7: PDE.Rnw:239-243
###################################################
print(system.time(
  Spherical   <- steady.1D(y = runif(Grid$N),
    func = pde1D, parms = NULL, nspec = 1, A = r2)
))


###################################################
### code chunk number 8: pde
###################################################
plot(Cartesian, Cylindrical, Spherical, grid = Grid$x.mid, 
   main = "steady-state PDE", xlab = "x", ylab = "C",
   col = c("darkgreen", "blue", "red"), lwd = 3, lty = 1:3)

legend("bottomright", c("cartesian", "cylindrical", "spherical"),
  col = c("darkgreen", "blue", "red"), lwd = 3, lty = 1:3)


###################################################
### code chunk number 9: figpde
###################################################
plot(Cartesian, Cylindrical, Spherical, grid = Grid$x.mid, 
   main = "steady-state PDE", xlab = "x", ylab = "C",
   col = c("darkgreen", "blue", "red"), lwd = 3, lty = 1:3)

legend("bottomright", c("cartesian", "cylindrical", "spherical"),
  col = c("darkgreen", "blue", "red"), lwd = 3, lty = 1:3)


###################################################
### code chunk number 10: PDE.Rnw:275-278
###################################################
max(abs(Q/6/D*(r2$mid - 10^2) + Cext - Spherical$y))
max(abs(Q/4/D*(r2$mid - 10^2) + Cext - Cylindrical$y))
max(abs(Q/2/D*(r2$mid - 10^2) + Cext - Cartesian$y))


###################################################
### code chunk number 11: PDE.Rnw:285-291
###################################################
require(deSolve)
times <- seq(0, 100, by = 1)
system.time(
  out <- ode.1D(y = rep(1, Grid$N), times = times, func = pde1D,
    parms = NULL, nspec = 1, A = r2)
)


###################################################
### code chunk number 12: PDE.Rnw:298-299
###################################################
tail(out[, 1:4], n = 3)


###################################################
### code chunk number 13: yy
###################################################
image(out, grid = Grid$x.mid, xlab = "time, days", 
     ylab = "Distance, cm", main = "PDE", add.contour = TRUE)


###################################################
### code chunk number 14: PDE.Rnw:356-360
###################################################
dx    <- 0.2
xgrid <- setup.grid.1D(-100, 100, dx.1 = dx)
x     <- xgrid$x.mid
N     <- xgrid$N


###################################################
### code chunk number 15: PDE.Rnw:363-367
###################################################
uini  <- exp(-0.05 * x^2)
vini  <- rep(0, N)
yini  <- c(uini, vini)
times <- seq (from = 0, to = 50, by = 1)


###################################################
### code chunk number 16: PDE.Rnw:374-382
###################################################
wave <- function (t, y, parms) {
  u1 <- y[1:N]
  u2 <- y[-(1:N)]

  du1 <- u2
  du2 <- tran.1D(C = u1, C.up = 0, C.down = 0, D = 1, dx = xgrid)$dC
  return(list(c(du1, du2)))
}


###################################################
### code chunk number 17: PDE.Rnw:387-389
###################################################
out <- ode.1D(func = wave, y = yini, times = times, parms = NULL,
         nspec = 2, method = "ode45", dimens = N, names = c("u", "v"))


###################################################
### code chunk number 18: wave
###################################################
matplot.1D(out, which = "u", subset = time %in% seq(0, 50, by = 10), 
   type = "l", col = c("black", rep("darkgrey", 5)), lwd = 2,
   grid = x, xlim = c(-50,50))
legend("topright", lty = 1:6, lwd = 2, col = c("black", rep("darkgrey", 5)), 
       legend = paste("t = ",seq(0, 50, by = 10)))


###################################################
### code chunk number 19: wave
###################################################
matplot.1D(out, which = "u", subset = time %in% seq(0, 50, by = 10), 
   type = "l", col = c("black", rep("darkgrey", 5)), lwd = 2,
   grid = x, xlim = c(-50,50))
legend("topright", lty = 1:6, lwd = 2, col = c("black", rep("darkgrey", 5)), 
       legend = paste("t = ",seq(0, 50, by = 10)))


###################################################
### code chunk number 20: waveimg
###################################################
par(mar=c(0,0,0,0))
image(out, which = "u", method = "persp", main = "",
        border = NA, col = "lightblue", box = FALSE, 
        shade = 0.5, theta = 0, phi = 60)                                                     


###################################################
### code chunk number 21: PDE.Rnw:465-473
###################################################
require(ReacTran)
pde2D <- function (t, y, parms) {
  CONC <- matrix(nr = n, nc = n, y)
  Tran <- tran.2D(CONC, D.x = Dx, D.y = Dy, dx = dx,  dy = dy)
  dCONC <- Tran$dC - r * CONC
  dCONC[ii]<- dCONC[ii] + p
  return(list(as.vector(dCONC)))
}


###################################################
### code chunk number 22: PDE.Rnw:482-488
###################################################
n  <- 100
dy <- dx <- 1
Dy <- Dx <- 2
r  <- 0.001
p  <- runif(50)
ii <- trunc(cbind(runif(50)*n, runif(50)*n) + 1)


###################################################
### code chunk number 23: PDE.Rnw:496-502
###################################################
require(rootSolve)
Conc0 <- matrix(nr = n, nc = n, 10.)
print(system.time(
  ST <- steady.2D(y = Conc0, func = pde2D, parms = NULL, dimens = c(n, n),
   lrw = 600000)
))


###################################################
### code chunk number 24: D2
###################################################
image(ST, main = "steady-state 2-D PDE")


###################################################
### code chunk number 25: D2fig
###################################################
image(ST, main = "steady-state 2-D PDE")