From 9815ebce9c7ed97d3c192a1f638aa54d5b5a5428 Mon Sep 17 00:00:00 2001 From: Marco De Lucia Date: Wed, 21 Dec 2022 18:54:15 +0100 Subject: [PATCH] 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)