EC526-Assignment 6 Solved

GOAL: MG (multigrid) and FFT (fast Fourier transform) are recursive divide and conquer methods that are the bedrock of fast algorithms in both computer science tasks and numerical methods for high performance computing. They work magically but to appreciate them you just need to code them. Their magic is easier to see in practice than to explain: Learn by doing!
Test Problem: Laplace Solver with Point Charge
Let’s start again with the simple iterative solvers Jacobi, Gauss Seidel and Red/Black (or even/odd) iterations in 1D. We want to compare them in 1D and 2D with multigrid and direct solve be Fourier Transforms! Remember all these methods must give the same solution when they converge, so you can use one method to debug the next. Also in 1D the exact solution is known even on a grid with spacing h, so there is no confusion of what the answer is. See the figure in HW6code. (Using a simple exact solution to debug a code is perhaps the most important debugging tool used by experts but too often not taught in class rooms.) The 1D Laplace equation for a function φ(x) of a real variable x ∈ [a,b] and put it on a grid with spacing h is:

                                                   ]                                (1)

Let’s take this to model temperature in the cold night with a heater in the middle of the 1D room. The walls are set to temperatures T1 and T2. For simplicity let’s take a = −1,b = 1 with a grid from i = −N,···N with h = 1/N. Let’s set the walls to absolute zero T1 = 0 and T2 = 0 with a single heater in the middle of the room! (I guess you are in outer space in a linear room. The matrix A T = b equation is

                                                                                                          (2)

with T[N] = T[−N] = 0. Actually we know the solution even on a grid with spacing h = 1/N

!

                                                                       T[i] = (N − |i|)(hα/2)     for     x = ih                                                (3)

Since this is exact, we can even look at for h → 0.

                                                                                           )                                                              (4)

whatever “δ(x)” means. The rule is the integral over δ(x) is 1. Let’s use this as the definition. (Physicists and engineers cheat this way all the time. Actually it is not cheating – it is smart. [1]) Now with h = 1/N our solution is T = h(N − |i|)(α/2) → (1 − |x|)(α/2) is independent of h so it is also the continuum solution when we take h → 0 or N → ∞. It has zero second derivative for all x 6= 0. What happens at x = 0? Lets check the solution by comparing the LHS and RHS of Eq. 4 near x = 0,

                                                     LHS      = 

                                                     RHS      =                                                                                   (5)

for 

0.1       Doubling Trick
The code on GitHUB has periodic boundary conditions, stabilized by a small mass term instead of fixed boundaries or a source. As a convenience we introduce a doubling trick to replace the boundary condition with T = 0 by an antisymmetric periodic problem. This avoids using boundary conditions with a doubled periodic lattice with images charges. Helpful to simplify our MG and FFT exercise. The exact solution we gave above strictly is the limit to zero mass. Try it and you will see.

The problem we described above has fixed boundary conditions so the value of T[i] at i = 0 and i = 2N are not variables. So there are really are 2N − 1 free variables. We would like that to have 2n variables to do the Multigrid or FFT so let’s be creative but with 2N − 1 = 2n this is impossible.

There are lots of methods to deal with boundary condition BUT to make life easy (always a good idea at first) let’s turn this into a periodic problem. This is easy. Just double the size of the problem and make it odd with reflection around the boundaries. To do this take double the 2N + 1 points to N0 = 4N + 2 − 2 = 4N points and make the system anti-periodic. ( We have - 2 because when you joint the two parts the T = 0 ends are identified.) Namely put a positive source at i = N again and a negative source at i = −N which is equivalent mod N0 to a negative source at i = 3N (i.e. −NmodN0 = 3N). Now the solution automatically will vanish at i = 0 and i = 2N by symmetry so we can solve this problem with multigrid and FFTs on a periodic grid of size N0 = 4N. The periodic problem has, for doubleT[N0],

                                                               (6)

defining wrap around conditions: T[i+1] ≡ T[(i+1)%N0] and T[i−1] ≡ T[(i−1+N0)%N0] for i = 0,··· ,N0 − 1. We assume powers of 2 (N0 = 2n). m2 is a small number that will not change the solution very much. In your code you can set h = 1. Sept we have changed the parameter N for convenience in the next part.

Part I Recursive Multigrid Solver
The exercise is to use a point source into the Laplace problem in both 1D and 2D and solve for the solution using a simple iteration and a Multigrid program. First solve this problem with a simple Jacobi, Gauss Seidel and Red/Black iteration with fixed boundary condition as a test program. Actually you already have this code form HW5. To give you an easy

The problem is to program this with multigrid and compare the iteration number for large N and scaling with respect to N relative to the single level algorithms. To implement Multigrid start on GitHub there are both 1D and 2D multigrid codes with periodic boundary condition but now point sources in HW6code/MG and the doubling trick to mimic the fixed boundary conditions. If the Multigrid is really working you can take it almost to zero. (Note in 2D you double twice. Once on each axes. e.g a “quadrupling trick” ?)

.The code is recursive and has 3 basic routines, described in class and the lecture notes. The top level has h = 1 and N0 = 2n grid points. The problem is to program this with multigrid and compare the iteration number for large N and scaling with respect to N relative to the single level algorithms aboveThis is called level 0. Below it are levels with spacing 2levelh and N0/2level grid points. Note that there 3 crucial function:

Proj(double *rHat, double r, int level);                                             // Project level to level +1

Inter(double *error,double *errorHat,int level);// Interpolate level to level -1

Iterate(double *phi_new, *phi, *b , level);                                     // Iterate Once on level
Coding Exercise #1: Point source with Images
Now for the problem you should slightly modify the MG code provided in both 1D and 2D, In 2D you may simplify the exercise still further with just one point source in the middle and periodic boundary conditions in both x and y directions. (For further reference possibly for you project see Sec 10.2 in class notes for fixed boundary conditions.) The program should stop iterating when each method has reached single precision convergence: