# R/tran.polar.R In ReacTran: Reactive Transport Modelling in 1d, 2d and 3d

#### Documented in tran.polar

```##==============================================================================
## Diffusive transport in polar coordinates (r, theta)
##==============================================================================

tran.polar <- function(C, C.r.up=NULL,    C.r.down=NULL,
C.theta.up=NULL,  C.theta.down=NULL,

flux.r.up=NULL,   flux.r.down=NULL,
flux.theta.up=NULL, flux.theta.down=NULL,

cyclicBnd = NULL,

D.r=1, D.theta=D.r,

r=NULL, theta=NULL,

full.output = FALSE )

{
# ------------------------------------------------------------------------------
# Initialisation
# ------------------------------------------------------------------------------

# a function to check the dimensionality of the system
checkDim <- function (dr, dtheta, nx="dr", ntheta="dtheta") {
if (length(dr) != 1 )
stop (nx,", should have length 1")
if (length(dtheta) != 1)
stop (ntheta,", should have length 1")
}

checkDim2 <- function (dr,dtheta,nx="dr",ntheta="dtheta") {
if (length(dr) != dimC[1]+1)
stop (nx,", should have length equal to dim(C)[1]+1 = ",dimC[1]+1," it is ",length(dr))

if (length(dtheta) != dimC[2]+1)
stop (ntheta,", should have length equal to dim(C)[2]+1 = ",dimC[2]+1," it is ",length(dtheta))
}

dimC <- dim (C)
Type <- 2

if (max(theta)> 2*pi)
stop ("theta should be < 2pi")
if (min(theta)< 0)
stop ("theta should be > 0 ")

checkDim2(r, theta,"r, the grid in x-direction",
"theta, the grid in theta-direction")
# central values
Nr       <- dimC[1]
r.c      <- 0.5*(r[-1]+r[-(Nr+1)])
dr       <- diff(r)
drint    <- c(r.c[1]-r[1],diff(r.c),r[Nr+1]-r.c[Nr])
divr     <- 1/r
divr[is.na(divr)] <- 0

Ntheta     <- dimC[2]
theta.c    <- 0.5*(theta[-1]+theta[-(Ntheta +1)])
dtheta     <- diff(theta)
dthetaint  <- c(theta.c[1]-theta[1],diff(theta.c),theta[Ntheta+1]-theta.c[Ntheta])

checkDim (D.r,D.theta,"D.r, the diffusion in x-direction ",
"D.theta, the diffusion in theta-direction ")

# check the dimensionality of the boundaries
if (! is.null(cyclicBnd)) {
# check if other boundaries not prescribed in this direction
}
if (! is.null(C.r.up)){
if( length(C.r.up) != 1 && length(C.r.up) != Ntheta)
stop ("'C.r.up' should have length 1 or equal to ", Ntheta)
} else if (! is.null(flux.r.up)) {
if( length(flux.r.up) != 1 && length(flux.r.up) != Ntheta)
stop ("'flux.r.up' should have length 1 or equal to ", Ntheta)
} else if (!1 %in% cyclicBnd)
C.r.up = C[1,]

if (! is.null(C.r.down)) {
if( length(C.r.down) != 1 && length(C.r.down) != Ntheta)
stop ("'C.r.down' should have length 1 or equal to ", Ntheta)
} else if (! is.null(flux.r.down)){
if(length(flux.r.down) != 1 && length(flux.r.down) != Ntheta)
stop ("'flux.r.down' should have length 1 or equal to ", Ntheta)
} else if (!1 %in% cyclicBnd)
C.r.down = C[Nr,]

if (! is.null(C.theta.up)) {
if ( length(C.theta.up) != 1 && length(C.theta.up) != Nr)
stop ("'C.theta.up' should have length 1 or equal to ", Nr)
} else if (! is.null(flux.theta.up)){
if (length(flux.theta.up) != 1 && length(flux.theta.up) != Nr)
stop ("'flux.theta.up' should have length 1 or equal to ", Nr)
} else if (!2 %in% cyclicBnd)
stop ("'flux.theta.up' OR 'C.theta.up' should be specified")

if (! is.null(C.theta.down)) {
if(length(C.theta.down) != 1 && length(C.theta.down) != Nr)
stop ("'C.theta.down' should have length 1 or equal to ", Nr)
} else if (! is.null(flux.theta.down)) {
if( length(flux.theta.down) != 1 && length(flux.theta.down) != Nr)
stop ("'flux.theta.down' should have length 1 or equal to ", Nr)
} else if (!2 %in% cyclicBnd)
stop ("'flux.theta.down' OR 'C.theta.down' should be specified")

# ------------------------------------------------------------------------------
# Function body: calculation
# ------------------------------------------------------------------------------

## Initialise 'fluxes' in all directions
Flux.r <- 0
Flux.theta <- 0

## First rewrite boundary values as "concentration"
## then perform diffusive transport
if (1 %in% cyclicBnd) {
C.r.up   <- (C[1,]*drint[Nr+1] + C[Nr,]*drint[1]) /(drint[1]+drint[Nr+1])
C.r.down <- C.r.up
}
if (is.null(C.r.up))
C.r.up <-  flux.r.up / D.r * drint[1]+ C[1,]
if (is.null(C.r.down))
C.r.down <- - flux.r.down / D.r * drint[Nr+1] + C[Nr,]

Flux.r <- D.r * (rbind(C.r.up,C, deparse.level = 0) -
rbind(C,C.r.down, deparse.level = 0))/drint

if (2 %in% cyclicBnd) {
C.theta.up   <- (C[,1]*dthetaint[Ntheta+1] + C[,Ntheta]*dthetaint[1]) /
(dthetaint[1]+dthetaint[Ntheta+1])
C.theta.down <- C.theta.up
}
if (is.null(C.theta.up))
C.theta.up <- flux.theta.up/D.theta * dthetaint[1]*r.c + C[,1]
if (is.null(C.theta.down))
C.theta.down <- - flux.theta.down /D.theta * dthetaint[Ntheta+1]*r.c + C[,Ntheta]

Flux.theta <- D.theta * (cbind(C.theta.up,C,deparse.level = 0) -
cbind(C,C.theta.down, deparse.level = 0))/
matrix(data= dthetaint,nrow=Nr,ncol=(Ntheta+1),byrow=TRUE) /r.c

## Calculate rate of change = flux gradient

dFdtheta <- 0.
dFdr   <- -diff(Flux.r * r)/dr/r.c
dFdtheta <- -t(diff(t(Flux.theta))/dtheta)/r.c

if (full.output)
return (list (dC = dFdr + dFdtheta ,                  # Rate of change
C.r.up     = C.r.up,
C.r.down   = C.r.down,
C.theta.up   = C.theta.up,
C.theta.down = C.theta.down,
r.flux   = Flux.r,
theta.flux = Flux.theta,
flux.r.up = Flux.r[1,],
flux.r.down = Flux.r[Nr+1,],
flux.theta.up = Flux.theta[,1],
flux.theta.down = Flux.theta[,Ntheta+1]

))
else
return (list (dC = dFdr + dFdtheta ,
flux.r.up = Flux.r[1,],
flux.r.down = Flux.r[Nr+1,],
flux.theta.up = Flux.theta[,1],
flux.theta.down = Flux.theta[,Ntheta+1]
))

} # end tran.polar
```

## Try the ReacTran package in your browser

Any scripts or data that you put into this service are public.

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