15-859E Programming Assignment 2: Multigrid

Due Date: 11:59pm Wed. 21 Oct. 1998

Revision 2: 14 Oct. 98.

In this assignment, you'll solve a 2-D simulation problem using the multigrid method. The recommended problem is Poisson's Equation on a unit square, and working from Briggs' book. I'm providing a bit of starter code, but not nearly as much as in the previous assignment.

Tips

• If you choose to simulate Poisson's equation, see Briggs' equation 1.4. For variety, you might want to try various right hand sides (f).
• If you're more ambitious, you might try multigrid on another type of problem, say from computer vision, elasticity, heat diffusion, or fluid mechanics. In general, multigrid usually works well on elliptical partial differential equations with symmetric positive definite matrices. But don't get too carried away - I urge you to try multigrid on a square, power-of-two grid hierarchy before attempting new applications, irregular domains, or unstructured meshes. If you want to do new research using multigrid, you might want to save that and make it your project at the end of the semester.
• One way to get ambitious without as much risk is to do the basic assignment first, get it running fast enough (using optimization tips in Makefile plus code tuning and run it on a smallsh grid) that you can get a full solution in less than a second, then implement interactive control of the right hand side (f).
• I found that restriction with full weighting works converges far faster than injection, but you may want to do some experimentation. For a relaxation method, I recommend weighted Gauss-Seidel, but you may want to experiment.
• Note that it's never necessary to allocate space for the A matrix! In my implementation I had 2-D arrays for v, f, and r at power-of-two sizes, and those were the only large memory allocations. (Briggs doesn't talk about keeping an r grid, but I found that this simplified the restriction code).
• I'll assume in most of the following tips that you're programming in C++.
• Use double precision (64 bit) floating point, not single precision.
• The recommended starter code is in pubdir/src/multigrid, where pubdir = /afs/cs/project/classes-ph/859E/pub . Take a look at it here. simp.cc is a very short program showing how to create a grid (2-D array) of arbitrary size with SVL's "Mat" class, and how to time a program in software.
• If you choose Poisson's equation, it helps to make your grids (2^k+1)x(2^k+1) on a side, with the perimeter values of the v grid permanently set to zero, and only the center (2^k-1)x(2^k-1) sub-grid unknown. Again, see Briggs equation 1.4. The principal advantage of this is that the 3x3 neighborhood of any of the unknowns is within bounds, so for relaxation and restriction you won't need special case code to handle neighbors of an edge or corner grid node. The setup for other simulation problems may differ a bit.
• The small program levi.cc shows how to use a nifty argument parser I wrote. I find this very handy for adding command line options to a program. You could use it, for example, to set the grid size, number of iterations here and there (Briggs' parameters nu0, nu1, nu2), select which relaxation algorithm to use, set omega parameter for weighted Gauss-Seidel (if you use that), or turn on debugging. Source and documentation for the argument parser are here.
• I recommend you debug by working with small grids (9x9 or less, say) and printing them after each operation. To facilitate this, use the formatting features of printf (or ostream) to get things to line up, e.g.
0.00   1.00   2.00              0 1 2
90.33 140.67 100.05      not     90.33 140.67 100.05
1.00   2.00   5.00              1 2 5

The former was generated with printf(" %5.2f", m[i][j]), the latter with printf(" %g", m[i][j]).
• On an SGI, if you use SVL, debug with compiler options "-g -DVL_CHECKING" (the latter causes subscriptrange checking in SVL) and linker option "-lsvl.dbg" but compile for runs on big grids with "-O2" and without "-DVL_CHECKING" and link with "-lsvl". For me, this change made my program 5x faster. On Suns, the options are the same except for "-O" instead of "-O2". These issues are mentioned in Makefile.sgi and Makefile.sun.
• I found it helpful to have an option that runs the multigrid algorithm silently except for printing the relative residual (see above) of the final solution at the end. That way I could run it with different parameters (varying the number of relaxation steps, or switching relaxation methods, say) to find the combination that most quickly found an accurate solution. You could also make the whole thing interactive, if you wish.

15-859, Hierarchical Methods for Simulation
Paul Heckbert