Working out heter.diff scheme (not working yet)

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.

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)