From e853a9683900dba08b52b6c36dd0f59483bcb146 Mon Sep 17 00:00:00 2001 From: Marco De Lucia Date: Tue, 20 Dec 2022 22:53:21 +0100 Subject: [PATCH 1/5] Working out heter.diff scheme (not working yet) --- doc/ADI_scheme.org | 82 ++++++++++++++++++++++- scripts/Adi2D_Reference.R | 132 ++++++++++++++++++++++++++++++++++++-- 2 files changed, 208 insertions(+), 6 deletions(-) diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index 6529522..fbcbefb 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -5,11 +5,12 @@ #+OPTIONS: toc:nil -* Diffusion in 1D +* Homogeneous diffusion in 1D ** Finite differences with nodes as cells' centres -The 1D diffusion equation is: +The 1D diffusion equation for spatially constant diffusion +coefficients $\alpha$ is: #+NAME: eqn:1 \begin{align} @@ -304,3 +305,80 @@ form: #+LATEX: \clearpage + +* Heterogeneous diffusion + +If the diffusion coefficient $\alpha$ is spatially variable, [[eqn:1]] can +be rewritten: + +#+NAME: eqn:hetdiff +\begin{align} +\frac{\partial C }{\partial t} & = \frac{\partial}{\partial x} \left(\alpha(x) \frac{\partial C }{\partial x} \right) +\end{align} + +\noindent From the product rule for derivatives we obtain: + +#+NAME: eqn:product +\begin{equation} +\frac{\partial}{\partial x} \left(\alpha(x) \frac{\partial C }{\partial x}\right) = \frac{\partial \alpha}{\partial x} \cdot \frac{\partial C}{\partial x} + \alpha \frac{\partial^2 C }{\partial x^2} +\end{equation} + +\noindent Using a spatially centred second order finite difference +approximation at $x=x_i$ for both $\alpha$ and $C$, we have +#+NAME: eqn:hetdiff_fd +\begin{align} +\frac{\partial \alpha}{\partial x} \cdot \frac{\partial C}{\partial x} + + \alpha \frac{\partial^2 C }{\partial x^2} & \simeq \frac{\alpha_{i+1} - \alpha_{i-1}}{2\Delta x}\cdot\frac{C_{i+1} - C_{i-1}}{2\Delta x} + \alpha_i \frac{C_{i+1} - 2 C_{i} + C_{i-1}}{\Delta x^2} \\ \nonumber + & = \frac{1}{\Delta x^2} \frac{\alpha_{i+1} - \alpha_{i-1}}{4}(C_{i+1} - C_{i-1}) + \frac{\alpha_{i}}{\Delta x^2}(C_{i+1}-2C_i+C_{i-1})\\ \nonumber + & = \frac{1}{\Delta x^2} \left\{A C_{i+1} -2\alpha_i C_i + AC_{i-1})\right\} +\end{align} +\noindent having set +\[ A = \frac{\alpha_{i+1}-\alpha_{i-1}}{4} + \alpha_i \] + +\noindent In 2D the ADI scheme [[eqn:genADI]] with heterogeneous diffusion +coefficients can thus be written: + +#+NAME: eqn:genADI_het +\begin{equation} +\systeme{ + \displaystyle \frac{C^{t+1/2}_{i,j}-C^{t }_{i,j}}{\Delta t/2} = \displaystyle \frac{\partial}{\partial x} \left( \alpha^x_{i,j} \frac{\partial C^{t+1/2}_{i,j}}{\partial x}\right) + \frac{\partial}{\partial y} \left( \alpha^y_{i,j} \frac{\partial C^{t }_{i,j}}{\partial y}\right), + \displaystyle \frac{C^{t+1 }_{i,j}-C^{t+1/2}_{i,j}}{\Delta t/2} = \displaystyle \frac{\partial}{\partial x} \left( \alpha^x_{i,j} \frac{\partial C^{t+1/2}_{i,j}}{\partial x}\right) + \frac{\partial}{\partial y} \left( \alpha^y_{i,j} \frac{\partial C^{t+1}_{i,j}}{\partial y}\right) +} +\end{equation} + +\noindent We define for compactness $S_x=\frac{\Delta t}{2\Delta x^2}$ +and $S_y=\frac{\Delta t}{2\Delta y^2}$ and + +#+NAME: eqn:het_AB +\begin{equation} +\systeme{ + \displaystyle A_{i,j} = \displaystyle \frac{\alpha^x_{i+1,j} -\alpha^x_{i-1,j }}{4} + \alpha^x_{i,j}, + \displaystyle B_{i,j} = \displaystyle \frac{\alpha^y_{i, j+1}-\alpha^y_{i ,j-1}}{4} + \alpha^y_{i,j} +} +\end{equation} + +\noindent Plugging eq. ([[eqn:hetdiff_fd]]) into the first of equations +([[eqn:genADI_het]]) - so called "sweep by x" - and putting all implicit +terms (at time level $t+1/2$) on the left hand side we obtain: + +#+NAME: eqn:sweepX_het +\begin{equation}\displaystyle +\begin{split} +-S_x A_{i,j} C^{t+1/2}_{i+1,j} + (1 + 2S_x\alpha^x_{i,j})C^{t+1/2}_{i,j} - S_x A_{i,j}C^{t+1/2}_{i-1,j} = \\ + S_y B_{i,j} C^{t }_{i,j+1} + (1 - 2S_y\alpha^y_{i,j})C^{t }_{i,j} + S_y B_{i,j}C^{t }_{i,j-1} +\end{split} +\end{equation} + +\noindent In the same way for the second of eq. [[eqn:genADI_het]] we have: + +#+NAME: eqn:sweepY_het +\begin{equation}\displaystyle +\begin{split} +-S_y B_{i,j} C^{t+1 }_{i,j+1} + (1 + 2S_y\alpha^y_{i,j})C^{t+1 }_{i,j} - S_y B_{i,j}C^{t+1 }_{i,j-1} = \\ + S_x A_{i,j} C^{t+1/2}_{i+1,j} + (1 - 2S_x\alpha^x_{i,j})C^{t+1/2}_{i,j} + S_x A_{i,j}C^{t+1/2}_{i-1,j} +\end{split} +\end{equation} + +\noindent If the diffusion coefficients are constant, +$A_{i,j}=B_{i,j}=\alpha$ and the scheme reverts to the homogeneous +case. diff --git a/scripts/Adi2D_Reference.R b/scripts/Adi2D_Reference.R index e3c3fd6..f1d942c 100644 --- a/scripts/Adi2D_Reference.R +++ b/scripts/Adi2D_Reference.R @@ -1,3 +1,4 @@ +## Time-stamp: "Last modified 2022-12-20 21:17:32 delucia" ## Brutal implementation of 2D ADI scheme ## Square NxN grid with dx=dy=1 @@ -23,12 +24,12 @@ ADI <- function(n, dt, iter, alpha) { tmpY[i,] <- SweepByRow(i, resY, dt=dt, alpha=alpha) res <- t(tmpY) - out[[it]] <- res } + out[[it]] <- res + } return(out) } - ## Workhorse function to fill A, B and solve for a given *row* of the ## grid matrix SweepByRow <- function(i, field, dt, alpha) { @@ -48,12 +49,12 @@ SweepByRow <- function(i, field, dt, alpha) { B <- numeric(ncol(field)) ## We now distinguish the top and bottom rows - if (i == 1){ + if (i == 1) { ## top boundary, "i-1" doesn't exist or is at a ghost ## node/cell boundary (TODO) for (ii in seq_along(B)) B[ii] <- (-1 +2*Sy)*field[i,ii] - Sy*field[i+1,ii] - } else if (i == nrow(field)){ + } else if (i == nrow(field)) { ## bottom boundary, "i+1" doesn't exist or is at a ghost ## node/cell boundary (TODO) for (ii in seq_along(B)) @@ -122,3 +123,126 @@ plot(adi2[[length(adi2)]], ref2, log="xy", xlab="ADI", ylab="ode.2D (reference)" las=1, xlim=c(1E-15, 1), ylim=c(1E-15, 1)) abline(0,1) + +## Test heterogeneous scheme +ADIHet <- function(field, dt, iter, alpha) { + + if (!all.equal(dim(field), dim(alpha))) + stop("field and alpha are not matrix") + + ## now both field and alpha must be nx*ny matrices + nx <- ncol(field) + ny <- nrow(field) + dx <- dy <- 1 + + ## find out the center of the grid to apply conc=1 + cenx <- ceiling(nx/2) + ceny <- ceiling(ny/2) + field[cenx, ceny] <- 1 + + Aij <- Bij <- alpha + + for (i in seq(2,ncol(field)-1)) { + for (j in seq(2,nrow(field)-1)) { + Aij[i,j] <- (alpha[i+1,j]-alpha[i-1,j])/4 + alpha[i,j] + Bij[i,j] <- (alpha[i,j+1]-alpha[i,j-1])/4 + alpha[i,j] + } + } + + if (any(Aij<0) || any(Bij<0)) + stop("Aij or Bij are negative!") + + ## prepare containers for computations and outputs + tmpX <- tmpY <- res <- field + out <- vector(mode="list", length=iter) + + for (it in seq(1, iter)) { + for (i in seq(1, ny)) + tmpX[i,] <- SweepByRowHet(i, res, dt=dt, alpha=alpha, Aij, Bij) + + resY <- t(tmpX) + for (i in seq(1, nx)) + tmpY[i,] <- SweepByRowHet(i, resY, dt=dt, alpha=alpha, Bij, Aij) + + res <- t(tmpY) + out[[it]] <- res + } + + return(out) +} + + +## Workhorse function to fill A, B and solve for a given *row* of the +## grid matrix +SweepByRowHet <- function(i, field, dt, alpha, Aij, Bij) { + dx <- 1 ## fixed in our test + Sx <- Sy <- dt/2/dx/dx + + ## diagonal of A at once + A <- matrix(0, nrow(field), ncol(field)) + diag(A) <- 1+2*Sx*diag(alpha) + + ## adjacent diagonals "Sx" + for (ii in seq(1, nrow(field)-1)) { + A[ii+1, ii] <- -Sx*Aij[ii+1,ii] + A[ii, ii+1] <- -Sx*Aij[ii,ii+1] + } + + B <- numeric(ncol(field)) + + ## We now distinguish the top and bottom rows + if (i == 1) { + ## top boundary, "i-1" doesn't exist or is at a ghost + ## node/cell boundary (TODO) + for (ii in seq_along(B)) + B[ii] <- Sy*Bij[i+1,ii]*field[i+1,ii] + (1-2*Sy*Bij[i,ii])*field[i, ii] + } else if (i == nrow(field)) { + ## bottom boundary, "i+1" doesn't exist or is at a ghost + ## node/cell boundary (TODO) + for (ii in seq_along(B)) + B[ii] <- (1-2*Sy*Bij[i,ii])*field[i, ii] + Sy*Bij[i-1,ii]*field[i-1,ii] + + } else { + ## inner grid row, full expression + for (ii in seq_along(B)) + B[ii] <- Sy*Bij[i+1,ii]*field[i+1,ii] + (1-2*Sy*Bij[i,ii])*field[i, ii] + Sy*Bij[i-1,ii]*field[i-1,ii] + } + + x <- solve(A, B) + x +} + +## adi2 <- ADI(n=51, dt=10, iter=200, alpha=1E-3) +## ref2 <- DoRef(n=51, alpha=1E-3, dt=10, iter=200) + +n <- 51 +field <- matrix(0, n, n) +alphas <- matrix(1E-3*runif(n*n, 1,1.2), n, n) + +## for (i in seq(1,nrow(alphas))) +## alphas[i,] <- seq(1E-7,1E-3, length=n) + +#diag(alphas) <- rep(1E-2, n) + +adih1 <- ADIHet(field=field, dt=10, iter=100, alpha=alphas) +adi2 <- ADI(n=n, dt=10, iter=100, alpha=1E-3) + + +par(mfrow=c(1,3)) +image(adi2[[length(adi2)]]) +image(adih1[[length(adih1)]]) +points(0.5,0.5, col="red",pch=4) +plot(adih1[[length(adih1)]], adi2[[length(adi2)]], pch=4, log="xy") +abline(0,1) + + +sapply(adih1, sum) +sapply(adi2, sum) + +adi2 + + +par(mfrow=c(1,2)) +image(alphas) +image(adih1[[length(adih1)]]) +points(0.5,0.5, col="red",pch=4) From 7bd40639d56871c178fdaa3af76d0301b7b41e97 Mon Sep 17 00:00:00 2001 From: Marco De Lucia Date: Wed, 21 Dec 2022 13:55:18 +0100 Subject: [PATCH 2/5] MDL: direct discretization scheme, theory in doc and dirty implementation in scripts --- doc/ADI_scheme.org | 59 ++++++++++++++++++++- scripts/Adi2D_Reference.R | 105 +++++++++++++++++++++++++++++++++++++- 2 files changed, 161 insertions(+), 3 deletions(-) diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index fbcbefb..096a02a 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -316,6 +316,8 @@ be rewritten: \frac{\partial C }{\partial t} & = \frac{\partial}{\partial x} \left(\alpha(x) \frac{\partial C }{\partial x} \right) \end{align} +** Discretization of the equation using chain rule + \noindent From the product rule for derivatives we obtain: #+NAME: eqn:product @@ -381,4 +383,59 @@ terms (at time level $t+1/2$) on the left hand side we obtain: \noindent If the diffusion coefficients are constant, $A_{i,j}=B_{i,j}=\alpha$ and the scheme reverts to the homogeneous -case. +case. Problem with this discretization is that the terms in $A_{ij}$ +and $B_{ij}$ can be negative depending on the derivative of the +diffusion coefficient, resulting in unphysical values for the +concentrations. + +** Direct discretization + +As noted in literature (LeVeque and Numerical Recipes) a better way is +to discretize directly the physical problem ([[eqn:hetdiff]]) at points +halfway between grid points: + +\begin{align*} +\begin{cases} +\displaystyle \alpha(x_{i+1/2}) \frac{\partial C }{\partial x}(x_{i+1/2}) & \displaystyle = \alpha_{i+1/2} \left( \frac{C_{i+1} -C_{i}}{\Delta x} \right) \\ +\displaystyle \alpha(x_{i-1/2}) \frac{\partial C }{\partial x}(x_{i-1/2}) & \displaystyle = \alpha_{i-1/2} \left( \frac{C_{i} -C_{i-1}}{\Delta x} \right) +\end{cases} +\end{align*} + +\noindent A further differentiation gives us the centered +approximation of $\frac{\partial}{\partial x} \left(\alpha(x) +\frac{\partial C }{\partial x}\right)$: + +\begin{align*} +\frac{\partial}{\partial x} \left(\alpha(x) +\frac{\partial C }{\partial x}\right)(x_i) & \simeq \frac{1}{\Delta x}\left[\alpha_{i+1/2} \left( \frac{C_{i+1} -C_{i}}{\Delta x} \right) - \alpha_{i-1/2} \left( \frac{C_{i} -C_{i-1}}{\Delta x} \right) \right]\\ +&\displaystyle =\frac{1}{\Delta x^2} \left[ \alpha_{i+1/2}C_{i+1} - (\alpha_{i+1/2}+\alpha_{i-1/2}) C_{i} + \alpha_{i-1/2}C_{i-1}\right] +\end{align*} + +\noindent The ADI scheme with this approach becomes: + +#+NAME: eqn:genADI_hetdir +\begin{equation} +\left\{ +\begin{aligned} + \frac{C^{t+1/2}_{i,j}-C^{t }_{i,j}}{\Delta t/2} = & \frac{1}{\Delta x^2} \left[ \alpha_{i+1/2,j} C^{t+1/2}_{i+1,j} - (\alpha_{i+1/2,j}+\alpha_{i-1/2,j}) C^{t+1/2}_{i,j} + \alpha_{i-1/2,j} C^{t+1/2}_{i-1,j}\right] + \\ + & \frac{1}{\Delta y^2} \left[ \alpha_{i,j+1/2}C^{t}_{i,j+1} - (\alpha_{i,j+1/2}+\alpha_{i,j-1/2}) C^t_{i} + \alpha_{i,j-1/2}C^{t}_{i,j-1}\right]\\ +\frac{C^{t+1 }_{i,j}-C^{t+1/2}_{i,j}}{\Delta t/2} = & \frac{1}{\Delta x^2} \left[ \alpha_{i+1/2,j}C^{t+1/2}_{i+1,j} - (\alpha_{i+1/2,j}+\alpha_{i-1/2,j}) C_{i} + \alpha_{i-1/2,j}C^{t+1/2}_{i-1,j}\right] + \\ + & \frac{1}{\Delta y^2} \left[ \alpha_{i,j+1/2}C^{t+1}_{i,j+1} - (\alpha_{i,j+1/2}+\alpha_{i,j-1/2}) C^{t+1}_{i,j} + \alpha_{i,j-1/2}C^{t+1}_{i,j-1}\right] +\end{aligned} +\right. +\end{equation} + +\noindent Doing the usual algebra and separating implicit from +explicit terms, the two sweeps become: + +#+NAME: eqn:sweepX_hetdir +\begin{equation} +\left\{ +\begin{aligned} +-S_x \alpha^x_{i+1/2,j} C^{t+1/2}_{i+1,j} + (1 + S_x(\alpha^x_{i+1/2,j}+ \alpha^x_{i-1/2,j}))C^{t+1/2}_{i,j} - S_x \alpha^x_{i-1/2,j} C^{t+1/2}_{i-1,j} = & \\ + S_y \alpha^y_{i,j+1/2} C^{t }_{i,j+1} + (1 - S_y(\alpha^y_{i,j+1/2}+ \alpha^y_{i,j-1/2}))C^{t }_{i,j} + S_y \alpha^y_{i,j-1/2} C^{t }_{i,j-1} & \\[1em] +-S_y \alpha^y_{i,j+1/2} C^{t+1 }_{i,j+1} + (1 + S_y(\alpha^y_{i,j+1/2}+ \alpha^y_{i,j-1/2}))C^{t+1 }_{i,j} - S_y \alpha^y_{i,j-1/2} C^{t+1 }_{i,j-1} = & \\ + S_x \alpha^x_{i+1/2,j} C^{t+1/2}_{i+1,j} + (1 - S_x(\alpha^x_{i+1/2,j}+ \alpha^x_{i-1/2,j}))C^{t+1/2}_{i,j} + S_x \alpha^x_{i-1/2,j} C^{t+1/2}_{i-1,j} +\end{aligned} +\right. +\end{equation} diff --git a/scripts/Adi2D_Reference.R b/scripts/Adi2D_Reference.R index f1d942c..ffc8212 100644 --- a/scripts/Adi2D_Reference.R +++ b/scripts/Adi2D_Reference.R @@ -1,4 +1,4 @@ -## Time-stamp: "Last modified 2022-12-20 21:17:32 delucia" +## Time-stamp: "Last modified 2022-12-21 13:47:00 delucia" ## Brutal implementation of 2D ADI scheme ## Square NxN grid with dx=dy=1 @@ -124,7 +124,7 @@ plot(adi2[[length(adi2)]], ref2, log="xy", xlab="ADI", ylab="ode.2D (reference)" abline(0,1) -## Test heterogeneous scheme +## Test heterogeneous scheme, chain rule ADIHet <- function(field, dt, iter, alpha) { if (!all.equal(dim(field), dim(alpha))) @@ -246,3 +246,104 @@ par(mfrow=c(1,2)) image(alphas) image(adih1[[length(adih1)]]) points(0.5,0.5, col="red",pch=4) + + + + +## Test heterogeneous scheme, direct discretization +ADIHetDir <- function(field, dt, iter, alpha) { + + if (!all.equal(dim(field), dim(alpha))) + stop("field and alpha are not matrix") + + ## now both field and alpha must be nx*ny matrices + nx <- ncol(field) + ny <- nrow(field) + dx <- dy <- 1 + + ## find out the center of the grid to apply conc=1 + cenx <- ceiling(nx/2) + ceny <- ceiling(ny/2) + field[cenx, ceny] <- 1 + + ## prepare containers for computations and outputs + tmpX <- tmpY <- res <- field + out <- vector(mode="list", length=iter) + + for (it in seq(1, iter)) { + for (i in seq(2, ny-1)) { + Aij <- cbind((alpha[i,]+alpha[i-1,])/2, (alpha[i,]+alpha[i+1,])/2) + Bij <- cbind((alpha[,i]+alpha[,i-1])/2, (alpha[,i]+alpha[,i+1])/2) + tmpX[i,] <- SweepByRowHetDir(i, res, dt=dt, Aij, Bij) + } + resY <- t(tmpX) + for (i in seq(2, nx-1)) + tmpY[i,] <- SweepByRowHetDir(i, resY, dt=dt, Bij, Aij) + res <- t(tmpY) + out[[it]] <- res + } + + return(out) +} + + +## Direct discretization, Workhorse function to fill A, B and solve +## for a given *row* of the grid matrix +SweepByRowHetDir <- function(i, field, dt, Aij, Bij) { + dx <- 1 ## fixed in our test + Sx <- Sy <- dt/2/dx/dx + + ## diagonal of A at once + A <- matrix(0, nrow(field), ncol(field)) + diag(A) <- 1 + Sx*(Aij[,1]+Aij[,2]) + + ## adjacent diagonals "Sx" + for (ii in seq(1, nrow(field)-1)) { + A[ii+1, ii] <- -Sx*Aij[ii,1] + A[ii, ii+1] <- -Sx*Aij[ii,2] + } + + B <- numeric(ncol(field)) + + for (ii in seq_along(B)) + B[ii] <- Sy*Bij[ii,2]*field[i+1,ii] + (1 - Sy*(Bij[ii,1]+Bij[ii,2]))*field[i, ii] + Sy*Bij[ii,1]*field[i-1,ii] + + x <- solve(A, B) + x +} + +## adi2 <- ADI(n=51, dt=10, iter=200, alpha=1E-3) +## ref2 <- DoRef(n=51, alpha=1E-3, dt=10, iter=200) + +n <- 51 +field <- matrix(0, n, n) +alphas <- matrix(1E-3*runif(n*n, 1,1.2), n, n) + + +## for (i in seq(1,nrow(alphas))) +## alphas[i,] <- seq(1E-7,1E-3, length=n) + +#diag(alphas) <- rep(1E-2, n) + +adih <- ADIHetDir(field=field, dt=10, iter=100, alpha=alphas) +adi2 <- ADI(n=n, dt=10, iter=100, alpha=1E-3) + + +par(mfrow=c(1,3)) +image(adi2[[length(adi2)]]) +image(adih[[length(adih)]]) +points(0.5,0.5, col="red",pch=4) +plot(adih1[[length(adih1)]], adi2[[length(adi2)]], pch=4, log="xy") +abline(0,1) + + +sapply(adih, sum) +sapply(adi2, sum) + +adi2 + + +par(mfrow=c(1,2)) +image(alphas) +image(adih1[[length(adih1)]]) +points(0.5,0.5, col="red",pch=4) From 9815ebce9c7ed97d3c192a1f638aa54d5b5a5428 Mon Sep 17 00:00:00 2001 From: Marco De Lucia Date: Wed, 21 Dec 2022 18:54:15 +0100 Subject: [PATCH 3/5] Small fixes --- .gitignore | 1 + doc/ADI_scheme.org | 8 ++++---- scripts/Adi2D_Reference.R | 15 ++++++++++----- 3 files changed, 15 insertions(+), 9 deletions(-) diff --git a/.gitignore b/.gitignore index 1e841cf..cf9dc56 100644 --- a/.gitignore +++ b/.gitignore @@ -8,3 +8,4 @@ compile_commands.json .cache/ .ccls-cache/ /iwyu/ +.Rhistory diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index 096a02a..5cfa8b0 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -419,8 +419,8 @@ approximation of $\frac{\partial}{\partial x} \left(\alpha(x) \begin{aligned} \frac{C^{t+1/2}_{i,j}-C^{t }_{i,j}}{\Delta t/2} = & \frac{1}{\Delta x^2} \left[ \alpha_{i+1/2,j} C^{t+1/2}_{i+1,j} - (\alpha_{i+1/2,j}+\alpha_{i-1/2,j}) C^{t+1/2}_{i,j} + \alpha_{i-1/2,j} C^{t+1/2}_{i-1,j}\right] + \\ & \frac{1}{\Delta y^2} \left[ \alpha_{i,j+1/2}C^{t}_{i,j+1} - (\alpha_{i,j+1/2}+\alpha_{i,j-1/2}) C^t_{i} + \alpha_{i,j-1/2}C^{t}_{i,j-1}\right]\\ -\frac{C^{t+1 }_{i,j}-C^{t+1/2}_{i,j}}{\Delta t/2} = & \frac{1}{\Delta x^2} \left[ \alpha_{i+1/2,j}C^{t+1/2}_{i+1,j} - (\alpha_{i+1/2,j}+\alpha_{i-1/2,j}) C_{i} + \alpha_{i-1/2,j}C^{t+1/2}_{i-1,j}\right] + \\ - & \frac{1}{\Delta y^2} \left[ \alpha_{i,j+1/2}C^{t+1}_{i,j+1} - (\alpha_{i,j+1/2}+\alpha_{i,j-1/2}) C^{t+1}_{i,j} + \alpha_{i,j-1/2}C^{t+1}_{i,j-1}\right] +\frac{C^{t+1 }_{i,j}-C^{t+1/2}_{i,j}}{\Delta t/2} = & \frac{1}{\Delta y^2} \left[ \alpha_{i+1/2,j}C^{t+1/2}_{i+1,j} - (\alpha_{i+1/2,j}+\alpha_{i-1/2,j}) C_{i} + \alpha_{i-1/2,j}C^{t+1/2}_{i-1,j}\right] + \\ + & \frac{1}{\Delta x^2} \left[ \alpha_{i,j+1/2}C^{t+1}_{i,j+1} - (\alpha_{i,j+1/2}+\alpha_{i,j-1/2}) C^{t+1}_{i,j} + \alpha_{i,j-1/2}C^{t+1}_{i,j-1}\right] \end{aligned} \right. \end{equation} @@ -432,9 +432,9 @@ explicit terms, the two sweeps become: \begin{equation} \left\{ \begin{aligned} --S_x \alpha^x_{i+1/2,j} C^{t+1/2}_{i+1,j} + (1 + S_x(\alpha^x_{i+1/2,j}+ \alpha^x_{i-1/2,j}))C^{t+1/2}_{i,j} - S_x \alpha^x_{i-1/2,j} C^{t+1/2}_{i-1,j} = & \\ +-S_x \alpha^x_{i+1/2,j} C^{t+1/2}_{i+1,j} + (1 + S_x(\alpha^x_{i+1/2,j}+ \alpha^x_{i-1/2,j}))C^{t+1/2}_{i,j} - S_x \alpha^x_{i-1/2,j} C^{t+1/2}_{i-1,j} = \quad & \\ S_y \alpha^y_{i,j+1/2} C^{t }_{i,j+1} + (1 - S_y(\alpha^y_{i,j+1/2}+ \alpha^y_{i,j-1/2}))C^{t }_{i,j} + S_y \alpha^y_{i,j-1/2} C^{t }_{i,j-1} & \\[1em] --S_y \alpha^y_{i,j+1/2} C^{t+1 }_{i,j+1} + (1 + S_y(\alpha^y_{i,j+1/2}+ \alpha^y_{i,j-1/2}))C^{t+1 }_{i,j} - S_y \alpha^y_{i,j-1/2} C^{t+1 }_{i,j-1} = & \\ +-S_y \alpha^y_{i,j+1/2} C^{t+1 }_{i,j+1} + (1 + S_y(\alpha^y_{i,j+1/2}+ \alpha^y_{i,j-1/2}))C^{t+1 }_{i,j} - S_y \alpha^y_{i,j-1/2} C^{t+1 }_{i,j-1} = \qquad & \\ S_x \alpha^x_{i+1/2,j} C^{t+1/2}_{i+1,j} + (1 - S_x(\alpha^x_{i+1/2,j}+ \alpha^x_{i-1/2,j}))C^{t+1/2}_{i,j} + S_x \alpha^x_{i-1/2,j} C^{t+1/2}_{i-1,j} \end{aligned} \right. diff --git a/scripts/Adi2D_Reference.R b/scripts/Adi2D_Reference.R index ffc8212..97a476c 100644 --- a/scripts/Adi2D_Reference.R +++ b/scripts/Adi2D_Reference.R @@ -1,4 +1,4 @@ -## Time-stamp: "Last modified 2022-12-21 13:47:00 delucia" +## Time-stamp: "Last modified 2022-12-21 18:47:13 delucia" ## Brutal implementation of 2D ADI scheme ## Square NxN grid with dx=dy=1 @@ -319,21 +319,25 @@ n <- 51 field <- matrix(0, n, n) alphas <- matrix(1E-3*runif(n*n, 1,1.2), n, n) +alphas1 <- matrix(1E-5, n, 25) +alphas2 <- matrix(1E-5, n, 26) + +alphas <- cbind(alphas1, alphas2) ## for (i in seq(1,nrow(alphas))) ## alphas[i,] <- seq(1E-7,1E-3, length=n) #diag(alphas) <- rep(1E-2, n) -adih <- ADIHetDir(field=field, dt=10, iter=100, alpha=alphas) -adi2 <- ADI(n=n, dt=10, iter=100, alpha=1E-3) +adih <- ADIHetDir(field=field, dt=10, iter=200, alpha=alphas) +adi2 <- ADI(n=n, dt=10, iter=200, alpha=1E-5) par(mfrow=c(1,3)) image(adi2[[length(adi2)]]) image(adih[[length(adih)]]) points(0.5,0.5, col="red",pch=4) -plot(adih1[[length(adih1)]], adi2[[length(adi2)]], pch=4, log="xy") +plot(adih[[length(adih)]], adi2[[length(adi2)]], pch=4, log="xy") abline(0,1) @@ -345,5 +349,6 @@ adi2 par(mfrow=c(1,2)) image(alphas) -image(adih1[[length(adih1)]]) +points(0.5,0.5, col="red",pch=4) +image(adih[[length(adih)]]) points(0.5,0.5, col="red",pch=4) From e4ec0a11dad15249909404de38b5de7874fe4254 Mon Sep 17 00:00:00 2001 From: Marco De Lucia Date: Tue, 27 Dec 2022 15:37:41 +0100 Subject: [PATCH 4/5] MDL fixing docs --- doc/ADI_scheme.org | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index 5cfa8b0..e0d2abd 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -418,8 +418,8 @@ approximation of $\frac{\partial}{\partial x} \left(\alpha(x) \left\{ \begin{aligned} \frac{C^{t+1/2}_{i,j}-C^{t }_{i,j}}{\Delta t/2} = & \frac{1}{\Delta x^2} \left[ \alpha_{i+1/2,j} C^{t+1/2}_{i+1,j} - (\alpha_{i+1/2,j}+\alpha_{i-1/2,j}) C^{t+1/2}_{i,j} + \alpha_{i-1/2,j} C^{t+1/2}_{i-1,j}\right] + \\ - & \frac{1}{\Delta y^2} \left[ \alpha_{i,j+1/2}C^{t}_{i,j+1} - (\alpha_{i,j+1/2}+\alpha_{i,j-1/2}) C^t_{i} + \alpha_{i,j-1/2}C^{t}_{i,j-1}\right]\\ -\frac{C^{t+1 }_{i,j}-C^{t+1/2}_{i,j}}{\Delta t/2} = & \frac{1}{\Delta y^2} \left[ \alpha_{i+1/2,j}C^{t+1/2}_{i+1,j} - (\alpha_{i+1/2,j}+\alpha_{i-1/2,j}) C_{i} + \alpha_{i-1/2,j}C^{t+1/2}_{i-1,j}\right] + \\ + & \frac{1}{\Delta y^2} \left[ \alpha_{i,j+1/2}C^{t}_{i,j+1} - (\alpha_{i,j+1/2}+\alpha_{i,j-1/2}) C^t_{i,j} + \alpha_{i,j-1/2}C^{t}_{i,j-1}\right]\\ +\frac{C^{t+1 }_{i,j}-C^{t+1/2}_{i,j}}{\Delta t/2} = & \frac{1}{\Delta y^2} \left[ \alpha_{i+1/2,j}C^{t+1/2}_{i+1,j} - (\alpha_{i+1/2,j}+\alpha_{i-1/2,j}) C^t_{i,j} + \alpha_{i-1/2,j}C^{t+1/2}_{i-1,j}\right] + \\ & \frac{1}{\Delta x^2} \left[ \alpha_{i,j+1/2}C^{t+1}_{i,j+1} - (\alpha_{i,j+1/2}+\alpha_{i,j-1/2}) C^{t+1}_{i,j} + \alpha_{i,j-1/2}C^{t+1}_{i,j-1}\right] \end{aligned} \right. From cea941853bff2ae29feef928515c3db561cd7bc8 Mon Sep 17 00:00:00 2001 From: Marco De Lucia Date: Fri, 17 Mar 2023 10:57:07 +0100 Subject: [PATCH 5/5] Adi2D_Ref --- scripts/Adi2D_Reference.R | 26 +++++++++++++++++--------- 1 file changed, 17 insertions(+), 9 deletions(-) diff --git a/scripts/Adi2D_Reference.R b/scripts/Adi2D_Reference.R index 97a476c..9b0bdd1 100644 --- a/scripts/Adi2D_Reference.R +++ b/scripts/Adi2D_Reference.R @@ -1,4 +1,4 @@ -## Time-stamp: "Last modified 2022-12-21 18:47:13 delucia" +## Time-stamp: "Last modified 2023-01-05 17:52:55 delucia" ## Brutal implementation of 2D ADI scheme ## Square NxN grid with dx=dy=1 @@ -272,8 +272,8 @@ ADIHetDir <- function(field, dt, iter, alpha) { for (it in seq(1, iter)) { for (i in seq(2, ny-1)) { - Aij <- cbind((alpha[i,]+alpha[i-1,])/2, (alpha[i,]+alpha[i+1,])/2) - Bij <- cbind((alpha[,i]+alpha[,i-1])/2, (alpha[,i]+alpha[,i+1])/2) + Aij <- cbind(harm(alpha[i,], alpha[i-1,]), harm(alpha[i,], alpha[i+1,])) + Bij <- cbind(harm(alpha[,i], alpha[,i-1]), harm(alpha[,i], alpha[,i+1])) tmpX[i,] <- SweepByRowHetDir(i, res, dt=dt, Aij, Bij) } resY <- t(tmpX) @@ -286,6 +286,10 @@ ADIHetDir <- function(field, dt, iter, alpha) { return(out) } +harm <- function(x,y) 1/(1/x+1/y) + +harm(1,4) + ## Direct discretization, Workhorse function to fill A, B and solve ## for a given *row* of the grid matrix @@ -299,15 +303,17 @@ SweepByRowHetDir <- function(i, field, dt, Aij, Bij) { ## adjacent diagonals "Sx" for (ii in seq(1, nrow(field)-1)) { - A[ii+1, ii] <- -Sx*Aij[ii,1] - A[ii, ii+1] <- -Sx*Aij[ii,2] + A[ii+1, ii] <- -Sx*Aij[ii,2] # i-1/2 + A[ii, ii+1] <- -Sx*Aij[ii,1] # i+1/2 } B <- numeric(ncol(field)) for (ii in seq_along(B)) B[ii] <- Sy*Bij[ii,2]*field[i+1,ii] + (1 - Sy*(Bij[ii,1]+Bij[ii,2]))*field[i, ii] + Sy*Bij[ii,1]*field[i-1,ii] - + + lastA <<- A + lastB <<- B x <- solve(A, B) x } @@ -315,11 +321,11 @@ SweepByRowHetDir <- function(i, field, dt, Aij, Bij) { ## adi2 <- ADI(n=51, dt=10, iter=200, alpha=1E-3) ## ref2 <- DoRef(n=51, alpha=1E-3, dt=10, iter=200) -n <- 51 +n <- 5 field <- matrix(0, n, n) -alphas <- matrix(1E-3*runif(n*n, 1,1.2), n, n) +alphas <- matrix(1E-5*runif(n*n, 1,2), n, n) -alphas1 <- matrix(1E-5, n, 25) +alphas1 <- matrix(3E-5, n, 25) alphas2 <- matrix(1E-5, n, 26) alphas <- cbind(alphas1, alphas2) @@ -352,3 +358,5 @@ image(alphas) points(0.5,0.5, col="red",pch=4) image(adih[[length(adih)]]) points(0.5,0.5, col="red",pch=4) + +options(width=110)