##============================================================================== ## 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