You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
419 lines
14 KiB
419 lines
14 KiB
#============================================================================= |
|
# Compare PDEs to discrete Langrangian particle tracking simulation |
|
# |
|
# Script corresponds to the first application in the manuscript. |
|
# It can produce figure 1 from the manuscript. |
|
# |
|
# Coded by J. Rooze |
|
#============================================================================= |
|
|
|
library(deSolve) |
|
library(english) |
|
|
|
#create a grid in which space between nodes equals the jumping distance of molecules |
|
n_cells <- 50 #number of cells including boundaries |
|
l_domain <- 1 #length domain |
|
jump_distance <- l_domain / (n_cells-1) #jumping distance of particles |
|
jump_freq <- 1 #[1/time] #the frequency of successful jumps |
|
|
|
#Reactions |
|
R1 <- 0#10 # Production [conc / (time * location)] |
|
k <- 0# 1e-3 # Consumption rate constant [1/time] |
|
|
|
|
|
#functions below are slightly different than standard functions in moments package |
|
#b.c. sums are divided by n rather than (n-1), e.g. var = [sum(Xi-meanX)^2] / n |
|
|
|
#function to calculate variance |
|
moment_func <- function(X, m) { |
|
xi <- mean(X) |
|
var_X <- sum((X - xi)^2) / length(X) |
|
|
|
#treat special cases |
|
if (m == 1 ) { |
|
return(xi) |
|
} else if ((m < 2) || (round(m) != m)) { |
|
stop("Wrong input to moment calculation function") |
|
} else if (m == 2) { |
|
return(var_X) |
|
} else { |
|
return(sum((X - xi)^m) / (length(X) * sqrt(var_X)^m )) |
|
} |
|
} |
|
|
|
#Create populations of particles |
|
#make first 2 populations for the left and right boundaries |
|
particles_left_bnd <- c(rnorm(n = 1000, mean = 0, sd = 2 ), |
|
rnorm( 500, -0.5, 1 ), rnorm( 500, 0.5, 1 )) |
|
particles_right_bnd <- rnorm(n=4000, mean=1, sd=0.2 ) #standard Gaussian (skewness = 0, kurtosis = 3) |
|
#fill the grid interior with left boundary condition, and most right with with right population |
|
age_particles_init_grid <- lapply(1:n_cells, function(i_cell) { #as a list |
|
if (i_cell == n_cells) { |
|
particles_right_bnd |
|
} else |
|
particles_left_bnd |
|
}) |
|
age_particles_init <- c(rep(particles_left_bnd, n_cells -1), particles_right_bnd) #as vector |
|
|
|
|
|
################################################################################ |
|
## Setup and run particle tracking simulation |
|
################################################################################ |
|
|
|
#particle tracking |
|
n_init <- vapply(age_particles_init_grid, length, FUN.VALUE = n_cells) |
|
total_particles <- sum(n_init) |
|
xloc_particles <- rep(1:n_cells, n_init) |
|
age_particles <- age_particles_init |
|
|
|
#time-stepping |
|
dt <- 1 / jump_freq |
|
sim_time <- 5000 |
|
time <- 0 |
|
n_timesteps <- ceiling(sim_time/dt) |
|
output_times <- round( exp(seq(0, log(sim_time), length.out = 6 ))) |
|
|
|
|
|
#list to store results |
|
results <- vector(mode="list", length = length(output_times)+1) |
|
results[[1]] <- matrix(ncol = 2, c(xloc_particles, age_particles)) |
|
|
|
#simulate with particle tracking |
|
for(i_timestep in 1:n_timesteps) { |
|
#jump direction: -1 = left, 1 = right |
|
jumps <- sample(c(-1,1), replace=TRUE, size=total_particles) |
|
|
|
#accounting for jumps (all particles shift left or right) |
|
xloc_particles <- xloc_particles + jumps |
|
|
|
#enforce boundary conditions |
|
particles_at_left_bnd <- (xloc_particles == 1) |
|
particles_at_right_bnd <- (xloc_particles == n_cells) |
|
#remove particles at boundaries |
|
xloc_particles <- xloc_particles[!(particles_at_left_bnd | particles_at_right_bnd)] |
|
age_particles <- age_particles[!(particles_at_left_bnd | particles_at_right_bnd)] |
|
#place new particles at boundaries |
|
xloc_particles <- c(rep(1, length(particles_left_bnd)), #new particles left |
|
rep(n_cells, length(particles_right_bnd)), #new particles right |
|
xloc_particles) #particles within domain |
|
age_particles <- c(particles_left_bnd, particles_right_bnd, |
|
age_particles) |
|
|
|
#consumption reaction (concentration-dependent) |
|
if (!is.na(k) && k > 0) { |
|
for (i_cell in 2:(n_cells-1)) { |
|
#filter particles in specific cell |
|
particles_in_cell <- (xloc_particles == i_cell) |
|
#determine rate |
|
R2_cell <- round(k * sum(particles_in_cell)) |
|
|
|
if (R2_cell > 1) { |
|
#randomize sequence of particles to be removed |
|
toberemoved <- sample((1:length(xloc_particles))[particles_in_cell])[1:R2_cell] |
|
#remove concentration |
|
xloc_particles <- xloc_particles[-toberemoved] |
|
age_particles <- age_particles[-toberemoved] |
|
} |
|
} |
|
} |
|
|
|
#production reaction (not concentration-dependent) |
|
new_xloc_particles <- rep(2:(n_cells-1), each = R1) |
|
new_age_particles <- rep(0, length(new_xloc_particles)) |
|
#add the new particles to vectors |
|
xloc_particles <- c(xloc_particles, new_xloc_particles) |
|
age_particles <- c(age_particles, new_age_particles) |
|
|
|
#update number of total particles |
|
total_particles <- length(xloc_particles) |
|
|
|
#save result |
|
if (i_timestep %in% output_times) { |
|
results[[match(i_timestep, output_times) ]] <- matrix(ncol = 2, c(xloc_particles, age_particles)) |
|
print(i_timestep) |
|
} |
|
} |
|
|
|
#binning based on x location |
|
bin_particles <- function(result) { |
|
xloc <- result[,1] |
|
age <- result[,2] |
|
|
|
summary_grid <- matrix(nrow = 1+num_of_moments, ncol = n_cells) |
|
for (i_cell in 1:n_cells) { |
|
indices <- (xloc == i_cell) |
|
|
|
#check to avoid NA values when no particles are present |
|
if (sum(indices) == 0) { |
|
ages_icell <- c(0,0) #yields mu = 0 and sd = 0 |
|
} else |
|
ages_icell <- age[indices] |
|
|
|
|
|
summary_grid[1, i_cell] <- sum(indices) #number of particles |
|
for (i_moment in (1:num_of_moments) ) { #statistics |
|
summary_grid[1 + i_moment, i_cell] <- moment_func(ages_icell, m = i_moment) |
|
} |
|
} |
|
|
|
return(summary_grid) |
|
} |
|
|
|
|
|
################################################################################ |
|
## Setup and run simulation for partial differential equations |
|
################################################################################ |
|
|
|
#translate particle distribution to settings for pde |
|
#for pde grid |
|
nx <- n_cells |
|
particles_per_conc <- 1e3 |
|
num_of_moments <- 8 |
|
|
|
#boundary conditions (actually not used in pde, because boundaries included in grid) |
|
bc_mat <- matrix(ncol = 2, nrow = num_of_moments+1, 0) |
|
colnames(bc_mat) <- c("bc_left", "bc_right") |
|
bc_mat[1, ] <- c(length(particles_left_bnd), length(particles_right_bnd)) / particles_per_conc |
|
for (i_moment in 1:num_of_moments) { |
|
bc_mat[1+i_moment, ] <- vapply(list(particles_left_bnd, particles_right_bnd), |
|
moment_func, m = i_moment, FUN.VALUE = 2) |
|
} |
|
|
|
#initial conditions |
|
init_mat <- matrix(nrow = nx, ncol = num_of_moments+1, 0) |
|
init_mat[,1] <- approx(x = seq(0, 1, length=n_cells), |
|
y = vapply(age_particles_init_grid, length, FUN.VALUE = n_cells ), |
|
xout = seq(0, 1, length = nx ))$y / particles_per_conc |
|
for(i_moment in (1:num_of_moments) ) { |
|
init_mat[, 1+i_moment] <- approx(x = seq(0, 1, length=n_cells), |
|
y = vapply(age_particles_init_grid, moment_func, |
|
m = i_moment, FUN.VALUE = n_cells ), |
|
xout = seq(0, 1, length = nx ))$y * |
|
init_mat[,1] #multiply by concentration |
|
|
|
#account for difference between definitions moments and phi = (x - mu)^m / n |
|
if (i_moment > 2) { |
|
init_sd <- sqrt(init_mat[, 3] / init_mat[,1]) |
|
init_mat[, 1+i_moment] <- init_mat[, 1+i_moment] * init_sd^i_moment |
|
} |
|
} |
|
|
|
|
|
|
|
#simulate derived pde's, finite differences |
|
heat_pde <- function(t, y, parms) { |
|
with(as.list(parms), { |
|
|
|
#for uniform grid using arithmetic mean works |
|
if (length(D) == nx) { |
|
D_bnd <- 0.5 * (D[1:(nx-1)] + D[2:nx]) #length nx-1 |
|
} else if (length(D) == 1) { |
|
D_bnd <- rep(D, nx-1) |
|
} else |
|
stop("D vector in parms has the wrong length") |
|
|
|
|
|
#F = -D dC/dx |
|
flux <- - D_bnd * (y[2:nx] - y[1:(nx-1)]) / dx #length nx-1 |
|
|
|
#dC/dt = - dF/dx |
|
dCdt <- - (flux[2:(nx-1)] - flux[1:(nx-2)]) / dx #length nx-2 |
|
|
|
#first and last cell are fixed as boundary conditions, so there dC/dt=0 |
|
list(c(0,dCdt,0)) |
|
}) |
|
} |
|
|
|
coupled_pde <- function(t, y, parms) { |
|
with(as.list(parms), { |
|
#convert vector with state variables to matrix |
|
#order of columns: concentration, 1st moment, 2nd moment, etc. |
|
sv_mat <- matrix(nrow = nx, y) |
|
|
|
conc <- sv_mat[,1] |
|
mu <- sv_mat[,2] / conc |
|
var <- sv_mat[,3] / conc |
|
|
|
R1 <- c(0, rep(R1[1], nx-2), 0) |
|
R2 <- c(0, k * conc[-c(1, nx)], 0) #exclude boundaries |
|
aging <- 0 |
|
|
|
#matrix to store differentials |
|
dsvdt_mat <- matrix(nrow = nx, ncol = num_of_moments+1, 0) |
|
|
|
#change in concentration over time |
|
dsvdt_mat[,1] <- unlist( heat_pde(t, sv_mat[,1], parms) ) + R1 - R2 |
|
|
|
#d(mu * C)/dt |
|
dsvdt_mat[,2] <- unlist( heat_pde(t, sv_mat[,2], parms) ) + aging - mu * R2 |
|
|
|
#d(var * C)/dt |
|
varheat_term <- unlist( heat_pde(t, sv_mat[,3], parms) ) |
|
|
|
#calculate dmu/dx for left and right cell boundaries |
|
dmudx <- (mu[2:(nx)] - mu[1:(nx-1)])/dx |
|
#the gradient for interior cells, note that 2 multiplier from pde is in addition |
|
muterm <- c(0, D * conc[2:(nx-1)] * (dmudx[1:(nx-2)]^2 + dmudx[2:(nx-1)]^2), 0) |
|
|
|
#reaction terms |
|
|
|
#Production of material with age 0 |
|
dCmomentsdt_R1 <- matrix(nrow = nx, ncol = num_of_moments-1, 0) |
|
if (R1[2] != 0) { #test if production is turned on |
|
#for centralized moments mu0 = 1 and mu1=0 |
|
centralmoms <- sv_mat / conc |
|
centralmoms[,2] <- 0 |
|
|
|
#calculate the noncentral moments from central moments |
|
noncentralmoms <- sv_mat |
|
for (i_mom in 1:(num_of_moments+1) ) { |
|
n <- i_mom-1 |
|
k <- (1:i_mom)-1 |
|
noncentralmoms[,i_mom] <- vapply(1:nx, function(j_cell) { |
|
sum( |
|
choose(n, k) * centralmoms[j_cell, 1:i_mom] * mu[j_cell]^(n-k) |
|
) * conc[j_cell] |
|
}, FUN.VALUE = 1 ) |
|
} |
|
|
|
#create matrix to store results for var and higher moms |
|
dCmomentsdt <- matrix(nrow = nx, ncol = num_of_moments-1, 0) |
|
|
|
dmudt <- -mu * R1 / conc |
|
for (i_cell in 1:nx) { |
|
for (i_mom in 2:num_of_moments) { |
|
sumup <- 0 |
|
for (j in 0:i_mom) { |
|
mult <- choose(i_mom, j) |
|
|
|
mupowderiv <- j * mu[i_cell]^(j-1) * dmudt[i_cell] |
|
|
|
term <- (-1)^j * mult * mupowderiv * noncentralmoms[i_cell, i_mom-j+1] |
|
|
|
sumup <- sumup + term |
|
} |
|
dCmomentsdt_R1[i_cell, i_mom -1] <- sumup |
|
} |
|
} |
|
#respect boundary conditions |
|
dCmomentsdt_R1[1,] <- 0 |
|
dCmomentsdt_R1[nx,] <- 0 |
|
} |
|
|
|
#balance for variance |
|
varreacterm <- dCmomentsdt_R1[,1] - var * R2 |
|
dsvdt_mat[,3] <- varheat_term + muterm + varreacterm |
|
|
|
#calculate the derivatives for higher moments |
|
if (num_of_moments > 2) { |
|
for (i_moment in 3:num_of_moments) { |
|
#d(mean moment value * C)/dt |
|
#note that first column is for conc, 2nd for first moment, etc |
|
heat_term <- unlist( heat_pde(t, sv_mat[, 1+i_moment], parms) ) |
|
|
|
prevmoment <- sv_mat[,i_moment]/conc |
|
dprevmomentdx <- (prevmoment[2:(nx)] - prevmoment[1:(nx-1)])/dx #at cell boundaries |
|
dprevdmu <- 0.5 * (dprevmomentdx * dmudx)[1:(nx-2)] + |
|
0.5 * (dprevmomentdx * dmudx)[2:(nx-1)] #shifting to mid cells |
|
|
|
dmudx_term <- 2 * i_moment * D * conc * c(0, dprevdmu,0) |
|
|
|
#reaction terms |
|
currentmoment <- sv_mat[,i_moment+1] / conc |
|
reacterm <- dCmomentsdt_R1[,i_moment-1] - |
|
currentmoment * R2 |
|
dsvdt_mat[, 1+i_moment] <- heat_term + dmudx_term + reacterm |
|
} |
|
} |
|
|
|
list(c(dsvdt_mat)) |
|
}) |
|
} |
|
|
|
#simulate concentration |
|
parms <- c( |
|
D = jump_freq * jump_distance^2 * 0.5, |
|
nx = nx, |
|
dx = (jump_distance * (n_cells-1)) / (nx - 1), |
|
R1 = R1/particles_per_conc, |
|
k = k |
|
) |
|
|
|
pde_numerical_coupled <- ode(y = c(init_mat), func = coupled_pde, |
|
parms = parms, times = output_times, |
|
method = "vode") |
|
|
|
|
|
#plotting |
|
{ |
|
plot_rows <- floor(sqrt(num_of_moments+1)) |
|
plot_cols <- ceiling((num_of_moments+1)/plot_rows) |
|
par(mfrow = c(plot_rows, plot_cols)) |
|
|
|
#determine ylim for plotting |
|
y_lim_mat <- matrix(ncol = 2, nrow = num_of_moments+1, 0) |
|
for (i_output in 1:length(output_times)) { |
|
binned_result <- bin_particles( results[[i_output]] ) |
|
y_lim_mat[1,] <- c(min(binned_result[1, ]), max(binned_result[1, ])) |
|
for (j_prop in 2:(num_of_moments+1)) { |
|
y_lim_mat[j_prop,] <- c(min(binned_result[j_prop, ], y_lim_mat[j_prop,1]), |
|
max(binned_result[j_prop, ], y_lim_mat[j_prop,2])) |
|
} |
|
} |
|
|
|
#first plot: Number of Particles |
|
for (j_line in 1:length(output_times)) { |
|
plot_func <- plot |
|
if (j_line > 1) |
|
plot_func <- points |
|
|
|
binned_result <- bin_particles( results[[j_line]] ) |
|
|
|
plot_func(x = seq(0, l_domain, length.out = n_cells), y = binned_result[1,], |
|
ylim = y_lim_mat[1,]*c(0.9, 1.1), pch = 1, |
|
xlab = "x", ylab = "num. of particles", |
|
col = j_line) |
|
|
|
pde_conc <- pde_numerical_coupled[j_line, 2:(1+nx)] |
|
|
|
lines(x = seq(0, l_domain, length.out = nx), |
|
y = pde_conc * particles_per_conc, |
|
col = j_line, lwd = 1.5) |
|
} |
|
|
|
#plot the moments |
|
for (i_moment in 1:num_of_moments) { |
|
|
|
for (j_line in 1:length(output_times)) { |
|
plot_func <- plot |
|
if (j_line > 1) |
|
plot_func <- points |
|
|
|
binned_result <- bin_particles( results[[j_line]] ) |
|
|
|
plot_func(x = seq(0, l_domain, length.out = n_cells), y = binned_result[1+i_moment,], |
|
xlab = "x", ylab = paste(ordinal(i_moment), "moment"), |
|
col = j_line, ylim = y_lim_mat[1+i_moment,], pch = 1, lwd = 1.5) |
|
|
|
pde_conc <- unname(pde_numerical_coupled[j_line, 2:(1+nx)]) |
|
pde_moment <- unname(pde_numerical_coupled[j_line, (2+i_moment*nx):(1+nx*(i_moment+1))]/pde_conc) |
|
|
|
if (i_moment > 2) { |
|
pde_sd <- sqrt(unname(pde_numerical_coupled[j_line, (2+2*nx):(1+nx*(2+1))]) / pde_conc) |
|
pde_moment <- pde_moment / pde_sd^i_moment |
|
} |
|
|
|
|
|
lines(x = seq(0, l_domain, length.out = nx), |
|
y = pde_moment, |
|
col = j_line, lwd = 1.5) |
|
} |
|
} |
|
|
|
legend("topleft", legend = output_times, col = 1:j_line, pch = 1, title ='time', |
|
lty=1, lwd = 1.5) |
|
legend("topright", legend = c("particle track.", "pde soln"), |
|
pch = c(1,NA), lty = c(NA,1), lwd = c(NA, 1.5)) |
|
} |
|
|
|
|
|
|