Solving the 2D Heat Equation Using Finite Differencing on Complex Geometry¶

Completed as a requirement for CS 555 with Professor Andreas Kloeckner, this project solves the heat equation

$$u_t = \Delta u$$

on a complex two-dimentional domain; a room with with walls (Neumann boundary condition), windows (Dirichlet boundary condition) and radiators (heat sources).

In Part 1, we build up our grid, develope a scheme to handle boundary conditions, preform von Neumann stability analysis to choose a time step for our forward Euler solver, and then use an explicit time integrator to solve our equation.

In Part 2, we use a second-order accurate diagonally implicit time integrator, solve for a steady state solution, and implement the Jacobi and damped Jacobi method to solve for our steady state solution.

In Part 3, we implemet a V-Cycle Geometric Multigrid up to a 4-factor coarsening of our mesh, using restriction, prolongation, and Laplace matrices.

Part 1: An Explicit Solver¶

Setting up the grid¶

Building the Laplacian Matrix using with second order centered finite differences¶

Interior Regions¶

The finite difference scheme for the inerior region of the domain can be written as $$\begin{align*} \Delta \tilde{u}_i^j &= \frac{u_{i+1}^j - 2u_i^{j} + u_{i-1}^j}{h_x^2} + \frac{u_{i}^{j+1} - 2u_i^j + u_{i}^{j-1}}{h_y^2} \\ &= \left( \frac{u_{i+1}^{j} - u_{i}^{j}}{h_x^2} \right) + \left( \frac{-u_{i}^{j} + u_{i-1}^{j}}{h_x^2} \right) + \left( \frac{u_{i}^{j+1} - u_{i}^{j}}{h_y^2} \right) + \left( \frac{-u_{i}^{j} + u_{i}^{j-1}}{h_y^2} \right) \end{align*}$$

Nuemann Boundary Conditions¶

¶

The Nuemann boundary conditions give us that

$$\vec{n} \cdot \nabla u(x,t) = 0$$

When $u_{k-1}$ is in the wall, we use the left stencil D[0] = [-1.5 2. -0.5]

$$\frac{-3u_{k-1} + 4u_{k} -u_{k+1}}{2h} = 0$$

and we solve for $u_{k-1}$ such that

$$u_{k-1} = \frac{1}{3}(4u_{k} - u_{k+1})$$

Similarly, when $u_{k+1}$ is in the wall, we use the right stencil D[-1] = [ 0.5 -2. 1.5]

$$\frac{u_{k-1} - 4u_{k} + 3u_{k+1}}{2h} = 0$$

and we solve for $u_{k+1}$ such that

$$u_{k+1} = \frac{1}{3}(4u_{k} - u_{k-1})$$

Dirichlet boundary condition¶

This is simple enought to implement, we just do nothing to the neighbor when it is at a window

Testing Convergence of the Laplacian¶

Testing boundary conditions¶

To find the error in the infinity norm, we first create a test equation using the method of manufactured solutions. We do this on a square domain such that the the sides are Dirchle and the top and bottom are both Neumann. The we get the following boundary conditions: $$\begin{align*} \vec{n} \cdot \nabla u(x,0,t) &= 0 \\ \vec{n} \cdot \nabla u(x,1,t) &= 0 \\ u(0,y,t) &= 0 \\ u(1,y,t) &= 0 \end{align*}$$ on the grid $x \in [0,1], y \in [0,1]$

We use $u(x,y,t) = \sin(\alpha x ) \cos(\beta y) \cos(t)$ as our manufatured solution, and solving for $\alpha$ and $\beta$ on with our boudary equations yeilds $\alpha = \beta = \pi$. From here, we compute $$\Delta u = -2 \pi^2 \sin(\pi x ) \cos(\pi y) \cos(t)$$

We take the $\ell^\infty$ error of $\Delta_h u_h - \Delta u$ (the max absolute value of our numerical solution minus our true solution) on coarser and coarser grids, keeping $t$ fixed and arbitrary, plotting against $\mathcal{O}(h_x^n)$.

The results of this plot show that our Laplacian error has $\mathcal{O}(h_x^2)$ tuncation error, which is what we were aiming for.

Forward Euler on the system¶

Foreward Euler is given by $$\frac{u_{t+1} - u_t}{h_t} = \left(\frac{u_{i+1}^j - 2u_i^{j} + u_{i-1}^j}{h_x^2} + \frac{u_{i}^{j+1} - 2u_i^j + u_{i}^{j-1}}{h_y^2} \right) u_t$$ We preform Von Neummann stability analysis where $$\begin{align*} u(\theta, \varphi) &= \frac{ht}{hx^2}e^{-i\theta} - 2 \frac{ht}{hx^2} + \frac{ht}{hx^2}e^{i\theta} + \frac{ht}{hy^2}e^{-i \varphi} - 2 \frac{ht}{hy^2} + \frac{ht}{hy^2} e ^{i \varphi} \\ &= \frac{ht}{hx^2} 2\cos(\theta) - 2 \frac{ht}{hx^2} + \frac{ht}{hy^2} 2\cos(\varphi)- 2 \frac{ht}{hy^2} \\ p(\theta, \varphi) &= 1 \end{align*}$$ and thus $$\begin{align*} |s(\theta, \varphi)|^2 &= \left|2\frac{ht}{hx^2} (1-\cos(\theta)) + 2\frac{ht}{hy^2} (1-\cos(\varphi)) \right|^2 \end{align*}$$ It is clear to see that this equation is at a maximum when $\theta = \varphi = \pi$ Thus, $$\begin{align*} |s(\pi + n, \pi + m)|^2 \leq 1 & \implies \left|4\frac{ht}{hx^2} + 4\frac{ht}{hy^2} \right|^2 \leq 1 \\ &\implies 4 \frac{h_t}{\min \{ h_x^2, h_y^2 \} } \leq 1 \\ &\implies h_t \leq \frac{\min\{h_x^2, h_y^2 \}}{4} \end{align*}$$

We now asses any benefit to using an explicit time integrator that has an order of accuracy higher than first. Take a time itterator, noting that $\mathcal{O}(h_t) = \mathcal{O}(h_x^2)$ with the timestep taken above. $$u(t)' = u(t) + \Delta u(t) + \mathcal{O}(h_t^2)$$ But $\Delta u(t)$ has $\mathcal{O}(h_x^2) = \mathcal{O}(h_t)$ trucation error, and thus any explicit time intergrator would suffer from $\mathcal{O}(h_t)$ truncation error, and an explicit scheme would never achive more that first order accuracy.

Part 2: Implicit Solver, Steady-State Solution¶

second-order accurate diagonally implicit time integrator: the diagonaly implicit midpoint method $$u_{n+1} = u_n + k \left( f \left( \frac{1}{2} \left(y_n + y_{n+1}\right) \right) \right)$$ $$u_{n+1} = u_n + \frac{k}{2} \left( \Delta u_n + r(t) + \Delta u_{n+1} + r(t) \right)$$ which can be written as the solveable system $$\left(I - \frac{k}{2} \Delta \right)u_{n+1} = u_n + \frac{k}{2} \Delta u_n + 2r(t)$$

Steady State Solution¶

Our equation enters a steady state when $u_t = 0$. Thus, our equation becomes $$\begin{align*} 0 &= \Delta u + r(x) \\ \Delta u &= -r(x) \end{align*}$$ From here, we can use a sparse solver such as scipy.linalg.factorized to solve for the steady state

Jacobi Method¶

The Jacobi Method is derived from

$$\begin{align*} A\textbf{x} &= b \\ A &= D + M \\ (D + M)x &= b \\ Dx &= -Mx + b \\ \textbf{x}_{k+1} &= D^{-1} \left( -M \textbf{x} + b \right) \end{align*}$$

We apply this to the steady state system

$$\begin{align*} \Delta u &= -r(x) \\ u_{k+1} &= D^{-1} \left( -Mu -r(x) \right) \end{align*}$$

Running Jacobi for 300 itterations gives a solution with bounded and rappid occilations between values. First, we notice that $-\Delta$ has negative eigenvalues, as $\Delta$ is positive definite. When we view Jacobi as an explicit time integrator, this begins to make sence. First, we justify that Jacobi CAN, in fact, be viewed as an explicit time integrator. Take a general explicit time inegrator for the heat equation $$u_{n+1} = (I + h \cdot \Delta)u_n$$ and take Jacobi $$u_{n+1} = D^{-1}(-Mu_n)$$ Where $M = \Delta - D$. Thus, our equation becomes $$\begin{align*} u_{n+1} &= D^{-1}(-Mu_n) \\ &= D^{-1}(-\Delta + D)u_n \\ &= (-D^{-1} \Delta + I)u_n \end{align*}$$ Which is an explicit time integrator. We know that $-D^{-1} \Delta$ has negative eigenvalues from above. We know from stability analysis that our system is stable within the range of $-1 \to 0$ on the real axis of the imaginary plane. But we also have that the closer our eigenvalues get to $-1$, the more our solution occilates (while remaining stable). This occilation seems to be exactly the occilation we see in our Jacobi solution.

Damped Jacobi Method¶

Let $\tilde{u}_{k+1}$ be the iterate computed by jacobi, then let $$u_{k+1} = \alpha u_k + (1-\alpha) \tilde{u}_{k+1}$$

This damped Jacobi seems to act as a smoother on the overall output od the jacobi function. This may result from the fact that we are only applying a small percentage of jacobi at every step, so the osscialtions we see in undamped jacobi do not come into play. From imperical tests, it seems to me that the more I decrease the alpha term, the less-smoth the immage gets, as more jacobi is being applied.

Part 3: Geometric Multigrid¶

Constructing the restriction operator: picking values off the fine mesh a certian points¶

Constructing the prolongation operator¶

Let $u_f$ be a point on the fine grid and $u_c$ be a point on the coasrse grid. We now define how to construct the prolongation opperator.

The easiest case to deal with is when $u_f$ corresponds to $u_c$, and no interpolation is nessasary.

Next, we consider a point $u_f$ that is on one of the coarse grid axes shuch that $u_{f+1}$ and $u_{f-1}$ is on the coarse grid. We interpolate to second order such that

$$u_f = \frac{u_{f+1} + u_{f-1}}{2}$$

We now consider when $u_f$ is not on any coarse axis. We again interpolate in the second order such that

$$u_f = \frac{u_f^{NW} + u_f^{SW} + u_f^{SE} + u_f^{NE}}{4}$$

where $u_f^{NW}$ is the north-west neighbor.

Prolongation boundary condition handeling¶

As with the leplacian, the Dirichlet boundary conditions are easy to handle. If $u_{f+1}$ is in a wall, set $u_{f+1} = 0$.

For the Neumann boundary conditions, we first consider a point that is on a coarse axis. WLOG, assume $u_{f+1}$ is in a wall. Then we use the second order first derivative edge stencil

$$u_{f+1} = \frac{1}{3}(4 u_f - u_{f-1})$$

and our interpolation becomes

$$\begin{align*} u_f = \frac{u_{f+1} + u_{f-1}}{2} &\implies u_f = \frac{\frac{1}{3}(4 u_f - u_{f-1}) + u_{f-1}}{2} \\ &\implies (2 -\frac{4}{3})u_f = (1-\frac{1}{3})u_{f-1} \\ &\implies u_f = \left(\frac{1}{2} + \frac{1}{2} \right)u_{f-1} \end{align*}$$

We now consider the case where we interpolate a point that is not on any coarse axis. WLOG, assume $u_f^{NW}$ is in the wall. Then

$$u_{NW} = \frac{1}{3}(4 u_f - u_{f}^{SE})$$

and our interpolation becomes

$$\begin{align*} 4u_f = \frac{1}{3}(4 u_f - u_{f}^{SE}) + u_f^{SW} + u_f^{SE} + u_f^{NE} &\implies 4 \left(1 - \frac{1}{3} \right)u_f = u_f^{SW} + \left(1-\frac{1}{3} \right)u_f^{SE} + u_f^{NE} \\ &\implies \frac{8}{3}u_f = u_f^{SW} + \frac{2}{3}u_f^{SE} + u_f^{NE} \\ &\implies u_f = \left( \frac{1}{4} + \frac{1}{8} \right) (u_f^{SW} + u_f^{NE}) + \frac{1}{4}u_f^{SE} \\ \end{align*}$$

Now we consider the diagnoal case where two neighbors are in the wall. WLOG, let $u_f^{NW}$ and $u_f^{NE}$ be in the wall. Then we have

$$u_f = \left( \frac{1}{4} + \frac{1}{8} \right) (u_f^{SW} + u_f^{NE}) + \frac{1}{4}u_f^{SE}$$

and

$$u_f^{NE} = \frac{1}{3}(4 u_f - u_{f}^{SW})$$

thus

$$u_f = \frac{1}{2} (u_f^{SW} + u_f^{SE})$$

Finaly, we consider the case where only one neighbor is in the domain. WLOG, assume $u_f^{NE}$ is the only neighbor in the domain. We use the diagonal line interpolation

$$u_f = \frac{1}{2}(u_f^{NE -1} + u_f^{NE})$$

Then we take the edge stencil

$$u_f^{NE - 1} = \frac{1}{3} \left( 4 u_f - u^{NE} \right)$$

and we get

$$u_f = \left( \frac{1}{4} + \frac{3}{4} \right)u^{NE}$$

Testing convergence of the prolongation opperator¶

We test convergence using the same test domain and solution as with the Leplacian convergence test, with Dirichlet boundary's on top and Neumann boundaries on the side

We will plot $||u - P_2 ~P_4~R_4 ~R_2 ~u||_\infty$ versus $h_x$ where $P$ is the prolongation matrix.

Multigrid¶

Error and Blow-Up¶

In the calculations bellow, the multigrid proccess blows-up. I believe that this is a resul of the prolongation opperator failing to interpolate along the wall and window that is not 90 degrees. The plot bellow shows the error in in the steady state solution after one cycle of prolongation and restriction. The error seems to be localized to the stair-step window boundary. I am not quite sure what could be the cause of this error, as I do not believe that wall moving could happen after one interpolation cycle, though I am not sure that I have exausted every scenario on my test grid. Additionaly, I do not believe that my corner stencil is wrong as my other corners do not have error. My method to begug this would be to more extensivly exaust all possible scenarios on a test domain and extract the results to my problem.

As noted above, my residual does not decrease in the expected manner. Again, I believe that this is do to error in my interpolation near the stair-step boundary.

Animation of Heat Dispersion on a Course Grid¶