# V-cycle Multigrid for 2D transient heat transfer on a square plate using finite difference

I'm currently developing a program to solve 2D transient state heat conduction on a square plate using the V-cycle multigrid. Althought my program is able to reach the steady state solution, it's computational time is longer then just running the problem using Gauss-seidel method.

Problem case: 0.1 m by 0.1 m square plate with fixed temperatures. Top: 20°C and the other three sides to be 40°C. Assuming constant material properties, no internal heat generation, and equal grid length $\Delta x=\Delta y$.

Parameters: thermal diffusion ($\alpha=23.1\times 10^{-6}$) (using steel for now as a guide), $\Delta t=0.01$

Methodology : implicit finite difference method.

$$T_2(i,j)=\frac{T_1 + F_O (T_2 (i-1, j) + T_2(i+1, j) + T_2(i, j+1) + T_2(i,j-1)}{1 + 4 F_O}$$

where $F_O$ is the fourier number $\alpha \Delta t/ \Delta x^2$.

• Step 1: Pre-smoothing- using the implicit finitie difference method shown above. The problem is smoothen by 2 cycles of red-balck gauss seidel.

• Step 2: compute residual- the residual is computed using

$$\operatorname{res}=\frac{T_1 + F_O [T_2(i-1,j) + T_2(i+1,j) + T_2(i,j+1) + T_2(i, j-1)]}{1 + 4 F_O} - T_2(i,j)$$

• Step 3: The residual is restricted using full weightage

• Step 4: The residual equation is solved with the intial guess of the error to be zero. I solve it by using $$\operatorname{error}(i,j)= \frac{\operatorname{res}(i,j) + F_O [\operatorname{error} (i-1,j) + \operatorname{error}(i+1,j) + \operatorname{error}(i,j-1)+ \operatorname{error}(i,j+1)]}{1 + 4 F_O}$$

I'm not sure if I used the correct equation to solve the residual equation and as the gird size decrease $\Delta x$ will increase therefore the Fourier number ($F_O$) will change . So I believe this maybe the problem to my program but i'm not sure as there are not alot of information regarding this. (I have been looking for a few weeks). I will really appreciate the help from you guys. If anyone needs any other information, please comment.

• What are your boundary conditions and are you treating them properly? Is your coarse grid operator equal to RAI where I is the interpolation operator, R is the restriction operator, and A is the original fine grid operator? Mar 17 '17 at 17:25
• hi nukeguy, My boundary conditions are fixed temperatures on all four side. The top is 20°C and the other three are 40°C. OH my... I have totally forgot about galerkin's coarse grid operator. I can use that to obtain my coarse grid operator which may solve my fourier number problem.. I will try it tomorrow and update you with the results Thanks so much for the reminder nukeguy!! Mar 17 '17 at 17:45
• How many levels is your multigrid? Can you plot convergence rate easily for us? Is the code available on github? I think you need to specifically and quantitatively describe the issues that you're having in order for people to help diagnose your problem. Mar 18 '17 at 2:01
• Hi Charlie, The coarsest grid for my multigrid is 3 X 3. I have plotted the error reduction rate for you (the G in the graph represents the grid size). Sorry, but the code is not on github. Currently, the main problem I'm facing is that the computational time for the Multigrid method is longer as compared to a single grid (Gauss-Seidel) which does not make sense to me. But the result obtained from both methods are the same. I assume that the coarse grid operator i am currently using to solve my residual equation is wrong.(Shown on the post). Mar 18 '17 at 9:58
• Good updates so far.. are you using cell center or cell corner data? (Or cell face/cell edge)? It looks like MG seems to outperform GS at the finest mesh. This seems correct since data on the coarser grids will have higher frequencies than the finer grids, so MG will be less likely to outperform GS. The plateau for the orange and grey curves look like a warning sign though. Mar 31 '17 at 5:25

Regarding your issue of multigrid taking longer than Gauss-Seidel:

1. It could be that you didn't code the interpolation/restriction operations very efficiently. Do you make sure to take full advantage of the sparsity of the matrices? We write the coarse grid operator as something like RAI, but you shouldn't multiply out the matrices directly. All three are sparse matrices and you can save a lot of multiplications by being smart about it.

2. In your figure, multigrid converges in fewer iterations than G-S, so it is doing something -- in most cases, you need less than half as many iterations with multigrid. And it converges to the right answer, so that leads me to believe you've probably coded it correctly. However, it doesn't seem like G-S is taking very many iterations at all. 20 iterations isn't a lot for G-S. The benefit you get from multigrid will be minimal at best since G-S isn't really being all that slow for this problem. The point of using multigrid is to accelerate smoothers like G-S in cases where they converge really slowly (i.e., in cases with low frequency error modes). Since there isn't much to accelerate, the overhead of setting up the multigrid stuff is probably comparable to, if not greater than, the speedup you get from multigrid...

One big benefit with multigrid is that the number of iterations is usually independent of (or at least really insensitive) the problem size. The same cannot be said for G-S. So you could try bigger problems to see if you get a speedup, and you could try seeding the problem with something that you know will give a low frequency error mode (which G-S should struggle with).

Also, a side comment: unless you're actually computing with 43 digits of precision, your results are meaningless below a tolerance of $10^{-15}$ or so. I suggest cutting that part out of your plot. In all practical applications, I don't think you'd need to go much beyond $10^{-12}$.

• Hi Nukeguy, Although it's being 9 months since I posted this (completely forgotten I posted here), I would like to thank you for taking the time to write this answer. Your second point was spot-on, my problem statement was too simple to require multigrid and the overhead setting up of the multigrid indeed slowed the computation time. Nov 25 '17 at 10:54

Your residual calculation is incorrect. If $[A]\{u\}=\{b\}$ is your system of linear equations, the residual is $res=\{b\}-[A]\{u_k\}$,where $u_k$ is the approximation.

In your case it should be calculated as: $$\operatorname{res}=T_1 + F_O [T_2(i-1,j) + T_2(i+1,j) + T_2(i,j+1) + T_2(i, j-1)] - (1 + 4 F_O) T_2(i,j)$$

As you can see, it just needed to be multiplied by $(1 + 4 F_O)$