R/funcs_diffs_explicit_nolimF.R

In ahmad-alkadri/Rdiffsolver: A simple R package for solving diffusion equation with numerical methods.

#' 1D diffusion, boundary constant, FD explicit,
#' unbounded by F
#'
#' @description Function to solve diffusion equation
#' using finite difference method, forward euler or
#' explicit scheme, applied in the condition where:
#' a) the diffusion occurs inside a slab of material unidimensionally,
#' or only along one axis, b) the coefficient of diffusion is constant,
#' and c) the boundary conditions are constant.
#'
#' @usage mdfexdiffulibre(D,dt,dx,l,T,C_i=0,C_f=1)
#'
#' @param D Coefficient of diffusion, constant
#' @param dt different between each time steps
#' @param dx different between each space (x) steps
#' @param l half-length of the slab, usually in cm
#' @param T Total calculated diffusion time
#' @param C_i Initial concentration value inside the slab
#' @param C_f Dirichlet boundary condition,
#' final concentration coming from the outside of the slab,
#' with one or two elements. If there is only one element,
#' the two sides of the slab (x = -l and x = l) will have
#' the same C_f. If there are
#' two elements, the first element (C_f[1]) will be on
#' x = -l while the second one will be on x = l.
#'
#' @return A matrix with {round(T/dt,0)} number of row and
#' {round(L/dx,0)} number of column,
#' profiling the diffusion on slab.
#'
#' @examples
#' C_i = 0.00 # Initial concentration inside the slab
#' C_f = 1.00	# Final concentration coming from outside
#' D = 10^-7 # Coefficient of diffusion, cm^2/s
#' dt = 60 # difference between each time step
#' dx = 0.01 # space step in cm
#' l = 0.25 # half-thickness of the slab, in cm
#' T = 432000 # Total measured time in seconds (~5 days)
#' u <- mdfexdiffulibre(D,dt,dx,l,T,C_i,C_f)
#'
#' # Using plotly for plotting a contour plot
#' library(plotly)
#'
#' df.list <- list(x = seq(-l,l,length.out = ncol(u)),
#'                 y = seq(0,T,length.out = nrow(u)),
#'                 z = u)
#'
#' plot_ly() %>%
#'     add_contour(x = df.list$x, y = df.list$y, z = df.list$z) %>%
#'     layout(title = "Contour plot diffusion",
#'            xaxis = list(title = "x"),
#'            yaxis = list(title = "t (s)"))

mdfexdiffulibre <- function(D,dt,dx,l,T,C_i=0,C_f=1){

  if(length(C_f) > 2){
    stop("ERROR: C_f must consist of maximum 2 elements")
  }

  L = 2*l
  Nt = round(T/dt,0)
  # Mesh point in time
  t = seq(0, Nt*dt, length.out=Nt)
  F = D*dt/dx^2
  Nx = round(L/dx,0)
  # Mesh point in space
  x = seq(0,L,length.out=Nx)
  # Verify
  dx = x[2]-x[1]
  dt = t[2]-t[1]

  # Matrice vide pour u et u_n
  u = matrix(0,ncol=Nx,nrow=Nt)
  u_n	= seq(0,0,length.out=Nx)

  # Condition initiale
  for(i in 1:(Nx+1)){
    u_n[i] = C_i
  }

  # Compute
  for(n in 1:Nt){
    # Compute u at mesh points
    for(i in 2:Nx){
      u[n,i] = u_n[i]+F*(u_n[i-1]-2*u_n[i]+u_n[i+1])
    }
    # Insert boundary conditions
    u[n,1] = C_f[1] ; u[n,Nx] = C_f[length(C_f)]

    # Update u_n before next step
    u_n = u[n,]
  }

  return(u)

}