From c7efae395fd0603ed9b8aaed69ad56db4713bb3c Mon Sep 17 00:00:00 2001 From: Marco De Lucia Date: Fri, 21 Jul 2023 13:20:37 +0200 Subject: [PATCH] MDL: added HetDiff.R script to generate and visualize reference simulations in R --- scripts/Adi2D_Reference.R | 16 ++- scripts/HetDiff.R | 227 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 237 insertions(+), 6 deletions(-) create mode 100644 scripts/HetDiff.R diff --git a/scripts/Adi2D_Reference.R b/scripts/Adi2D_Reference.R index 947ddc8..7636b01 100644 --- a/scripts/Adi2D_Reference.R +++ b/scripts/Adi2D_Reference.R @@ -1,4 +1,4 @@ -## Time-stamp: "Last modified 2023-05-11 17:31:41 delucia" +## Time-stamp: "Last modified 2023-07-20 15:37:25 delucia" ## Brutal implementation of 2D ADI scheme ## Square NxN grid with dx=dy=1 @@ -286,7 +286,11 @@ ADIHetDir <- function(field, dt, iter, alpha) { return(out) } -harm <- function(x,y) 1/(1/x+1/y) +harm <- function(x,y) { + if (length(x) != 1 || length(y) != 1) + stop("x & z have different lengths") + 2/(1/x+1/y) +} harm(1,4) @@ -428,10 +432,10 @@ FTCS_2D <- function(field, dt, iter, alpha) { ## res[i,j] *(mean(alpha[i,j+1],alpha[i,j])+mean(alpha[i,j-1],alpha[i,j])) + ## res[i,j-1]*mean(alpha[i,j-1],alpha[i,j])) tmp[i,j] <- res[i,j] + - + tsteps[innerit]/dx/dx * ((res[i+1,j]-res[i,j]) * mean(alpha[i+1,j],alpha[i,j]) - - (res[i,j]-res[i-1,j]) * mean(alpha[i-1,j],alpha[i,j])) + - + tsteps[innerit]/dx/dx * ((res[i,j+1]-res[i,j]) * mean(alpha[i,j+1],alpha[i,j]) - - (res[i,j]-res[i,j-1]) * mean(alpha[i,j-1],alpha[i,j])) + + tsteps[innerit]/dx/dx * ((res[i+1,j]-res[i,j]) * harm(alpha[i+1,j],alpha[i,j]) - + (res[i,j]-res[i-1,j]) * harm(alpha[i-1,j],alpha[i,j])) + + + tsteps[innerit]/dx/dx * ((res[i,j+1]-res[i,j]) * harm(alpha[i,j+1],alpha[i,j]) - + (res[i,j]-res[i,j-1]) * harm(alpha[i,j-1],alpha[i,j])) } } ## swap back tmp to res for the next inner iteration diff --git a/scripts/HetDiff.R b/scripts/HetDiff.R new file mode 100644 index 0000000..e821e09 --- /dev/null +++ b/scripts/HetDiff.R @@ -0,0 +1,227 @@ +## Time-stamp: "Last modified 2023-07-21 11:23:08 delucia" + +library(ReacTran) +library(deSolve) +options(width=114) + +## harmonic mean +harm <- function(x,y) { + if (length(x) != 1 || length(y) != 1) + stop("x & y have different lengths") + 2/(1/x+1/y) +} + +## harm(0, 1) ## 0 +## harm(1, 2) ## 0 + + +############# Providing coeffs on the interfaces +N <- 11 # number of grid cells +ini <- 1 # initial value at x=0 +N2 <- ceiling(N/2) +L <- 10 # domain side + +## Define diff.coeff per cell, in 4 quadrants +alphas <- matrix(0, N, N) +alphas[1:N2, 1:N2] <- 1 +alphas[1:N2, seq(N2+1,N)] <- 0.1 +alphas[seq(N2+1,N), 1:N2] <- 0.01 +alphas[seq(N2+1,N), seq(N2+1,N)] <- 0.001 + +image(log10(alphas), col=heat.colors(4)) + +r180 <- function(x) { + xx <- rev(x) + dim(xx) <- dim(x) + xx +} +mirror <- function (x) { + xx <- as.data.frame(x) + xx <- rev(xx) + xx <- as.matrix(xx) + xx +} + +array_to_LaTeX <- function(arr) { + rows <- apply(arr, MARGIN=1, paste, collapse = " & ") + matrix_string <- paste(rows, collapse = " \\\\ ") + return(paste("\\begin{bmatrix}", matrix_string, "\\end{bmatrix}")) +} + + +cat(array_to_LaTeX(mirror(r180(alphas)))) + + + +r180(alphas) + +filled.contour(log10(alphas), col=terrain.colors(4), nlevels=4) + +cmpharm <- function(x) { + y <- c(0, x, 0) + ret <- numeric(length(x)+1) + for (i in seq(2, length(y))) { + ret[i-1] <- harm(y[i], y[i-1]) + } + ret +} + +## Construction of the 2D grid +x.grid <- setup.grid.1D(x.up = 0, L = L, N = N) +y.grid <- setup.grid.1D(x.up = 0, L = L, N = N) +grid2D <- setup.grid.2D(x.grid, y.grid) + +D.grid <- list() + +# Diffusion on x-interfaces +D.grid$x.int <- apply(alphas, 1, cmpharm) + +# Diffusion on y-interfaces +## matrix(nrow = N, ncol = N+1, data = rep(c(rep(1E-1, 50),5.E-1,rep(1., 50)), 100) ) +D.grid$y.int <- t(apply(alphas, 2, cmpharm)) + +dx <- L/N +dy <- L/N + +# The model equations +Diff2Dc <- function(t, y, parms) { + CONC <- matrix(nrow = N, ncol = N, data = y) + dCONC <- tran.2D(CONC, dx = dx, dy = dy, D.grid = D.grid)$dC + return(list(dCONC)) +} + +## initial condition: 0 everywhere, except in central point +y <- matrix(nrow = N, ncol = N, data = 0) +y[N2, N2] <- ini # initial concentration in the central point... + +## solve for 10 time units +times <- 0:10 +outc <- ode.2D(y = y, func = Diff2Dc, t = times, parms = NULL, + dim = c(N, N), lrw = 1860000) + +outtimes <- c(1, 4, 7, 10) + +cairo_pdf("deSolve_AlphaHet1.pdf", family="serif", width=12, height=12) +image(outc, ask = FALSE, mfrow = c(2, 2), main = paste("time", outtimes), + legend = TRUE, add.contour = FALSE, subset = time %in% outtimes, xlab="",ylab="", axes=FALSE) +dev.off() + +outc + + + + + +#################### 2D visualization + +## Version of Rmufits::PlotCartCellData with the ability to fix the +## "breaks" for color coding of 2D simulations +Plot2DCellData <- function(data, grid, nx, ny, contour = TRUE, + nlevels = 12, breaks, palette = "heat.colors", + rev.palette = TRUE, scale = TRUE, plot.axes=TRUE, ...) { + if (!missing(grid)) { + xc <- unique(sort(grid$cell$XCOORD)) + yc <- unique(sort(grid$cell$YCOORD)) + nx <- length(xc) + ny <- length(yc) + if (!length(data) == nx * ny) + stop("Wrong nx, ny or grid") + } else { + xc <- seq(1, nx) + yc <- seq(1, ny) + } + z <- matrix(round(data, 6), ncol = nx, nrow = ny, byrow = TRUE) + pp <- t(z[rev(seq(1, nrow(z))), ]) + + if (missing(breaks)) { + breaks <- pretty(data, n = nlevels) + } + + breakslen <- length(breaks) + colors <- do.call(palette, list(n = breakslen - 1)) + if (rev.palette) + colors <- rev(colors) + if (scale) { + par(mfrow = c(1, 2)) + nf <- layout(matrix(c(1, 2), 1, 2, byrow = TRUE), widths = c(4, + 1)) + } + par(las = 1, mar = c(5, 5, 3, 1)) + image(xc, yc, pp, xlab = "X [m]", ylab = "Y[m]", las = 1, asp = 1, + breaks = breaks, col = colors, axes = FALSE, ann=plot.axes, + ...) + + if (plot.axes) { + axis(1) + axis(2) + } + if (contour) + contour(unique(sort(xc)), unique(sort(yc)), pp, breaks = breaks, + add = TRUE) + if (scale) { + par(las = 1, mar = c(5, 1, 5, 5)) + PlotImageScale(data, breaks = breaks, add.axis = FALSE, + axis.pos = 4, col = colors) + axis(4, at = breaks) + } + invisible(pp) +} + +PlotImageScale <- function(z, zlim, col = grDevices::heat.colors(12), breaks, + axis.pos = 1, add.axis = TRUE, ...) { + if (!missing(breaks)) { + if (length(breaks) != (length(col) + 1)) { + stop("must have one more break than colour") + } + } + if (missing(breaks) & !missing(zlim)) { + breaks <- seq(zlim[1], zlim[2], length.out = (length(col) + 1)) + } + if (missing(breaks) & missing(zlim)) { + zlim <- range(z, na.rm = TRUE) + zlim[2] <- zlim[2] + c(zlim[2] - zlim[1]) * (0.001) + zlim[1] <- zlim[1] - c(zlim[2] - zlim[1]) * (0.001) + breaks <- seq(zlim[1], zlim[2], length.out = (length(col) + 1)) + } + + poly <- vector(mode = "list", length(col)) + for (i in seq(poly)) { + poly[[i]] <- c(breaks[i], breaks[i + 1], breaks[i + 1], + breaks[i]) + } + if (axis.pos %in% c(1, 3)) { + ylim <- c(0, 1) + xlim <- range(breaks) + } + if (axis.pos %in% c(2, 4)) { + ylim <- range(breaks) + xlim <- c(0, 1) + } + plot(1, 1, t = "n", ylim = ylim, xlim = xlim, axes = FALSE, + xlab = "", ylab = "", xaxs = "i", yaxs = "i", ...) + for (i in seq(poly)) { + if (axis.pos %in% c(1, 3)) { + polygon(poly[[i]], c(0, 0, 1, 1), col = col[i], border = NA) + } + if (axis.pos %in% c(2, 4)) { + polygon(c(0, 0, 1, 1), poly[[i]], col = col[i], border = NA) + } + } + box() + if (add.axis) { + axis(axis.pos) + } +} + +cairo_pdf("AlphaHet1.pdf", family="serif", width=8) +par(mar = c(1,1,1,1)) +Plot2DCellData(log10(mirror(alphas)), nx=N, ny=N, nlevels=5, palette = terrain.colors, contour=FALSE, plot.axes=FALSE, + scale = F, + main=expression(log[10](alpha))) +text(3,8,"1") +text(8,8,"0.1") +text(3,3,"0.01") +text(8,3,"0.001") +# title("Diff. Coefficients (log10)") +dev.off() +