BVP reaction transport

## =============================================================================
## R-Code to solve example 11.9 from the book
## K. Soetaert, J. Cash and F. Mazzia, 2012.  
## Solving differential equations in R. UseR, Springer, 248 pp. 
## http://www.springer.com/statistics/computational+statistics/book/978-3-642-28069-6.
## implemented by Karline Soetaert
## =============================================================================

## Steady-state solution of a reaction-transport problem
library(ReacTran)

N    <- 1000
Grid <- setup.grid.1D(N = N, L = 100000)
v <- 1000; D <- 1e7; O2s <- 300
NH3in <- 500; O2in <- 100; NO3in <- 50
r <- 0.1; k <- 1.; p <- 0.1

Estuary <- function(t, y, parms)  {
 NH3 <- y[1:N]
 NO3 <- y[(N+1):(2*N)]
 O2  <- y[(2*N+1):(3*N)]
 tranNH3<- tran.1D (C = NH3, D = D, v = v,
             C.up = NH3in, C.down = 10,  dx = Grid)$dC
 tranNO3<- tran.1D (C = NO3, D = D, v = v,
             C.up = NO3in,  C.down = 30,  dx = Grid)$dC
 tranO2 <- tran.1D (C = O2 , D = D, v = v,
             C.up = O2in, C.down = 250, dx = Grid)$dC

 reaeration <- p * (O2s - O2)
 r_nit      <- r * O2 / (O2 + k) * NH3

 dNH3    <- tranNH3 - r_nit
 dNO3    <- tranNO3 + r_nit
 dO2     <- tranO2  - 2 * r_nit + reaeration

 list(c(  dNH3, dNO3, dO2 ))
}

print(system.time(
std <- steady.1D(y = runif(3 * N), parms = NULL,
             names=c("NH3", "NO3", "O2"),
             func = Estuary, dimens = N,
             positive = TRUE)
))

NH3in <- 100
std2 <- steady.1D(y = runif(3 * N), parms = NULL,
             names=c("NH3", "NO3", "O2"),
             func = Estuary, dimens = N,
             positive = TRUE)

plot(std, std2, grid = Grid$x.mid, ylab = "mmol/m3", lwd = 2,
     xlab = "m", mfrow = c(1,3), col = c("black", "red"))
legend("bottomright", lty = 1:2, title = "NH3in", lwd = 2,
      legend = c(500, 100), col = c("black", "red"))