# COMSOL Blog

## On Solvers: Multigrid Methods

##### Valerio Marra | February 8, 2013

The differential equations that describe a real application admit an analytical solution only when several simplifying assumptions are made. The insights we gain from this approach are still valuable, but are not enough to confirm that our design is efficient or reduce the number of prototypes needed to reach a complete understanding of our application. This is why numerical solution methods are so important to us. They have been developed to overcome such limitations and allow us to represent our equations, on a grid obtained by partitioning the domain representing our application, in matrix form as , where is the vector solution. Once we have determined and , all we have to do is to “find” . It is important to note that the matrix form, both for stationary and transient applications, is obtained by finding a linear approximation to our equations. Nowadays, several solution methods are available. Here we will explore the ideas behind multigrid methods available in COMSOL Multiphysics.

### Solution Methods

Solution methods can be direct or iterative. Direct methods find an approximation of the solution by matrix factorization in a number of operations that depend on the number of unknowns. Factorization is expensive, but once it has been computed it is relatively inexpensive to solve for new right-hand sides . The approximate solution will be available only when all operations required by the factorization algorithm are executed.

Iterative methods begin with an approximated initial guess and then proceed to improve it by a succession of iterations. This means that, different from direct methods, we can stop the iterative algorithm at any iteration and have access to the approximate solution . Stopping the iterative process too early will result in an approximate solution with poor accuracy.

### Direct or Iterative?

It’s not easy to answer this question. To start, it depends on the application and the computer used. Since direct methods are both memory expensive and CPU time intensive, we prefer them for small to medium sized 2D and 3D applications. Conversely, iterative methods have a lower memory consumption and for large 3D applications outperform direct methods. It must be said that iterative methods are more difficult to tune and get working for matrices arising from multiphysics problems. As you can see a lot of different variables come into play when we have to make a decision about the best solver for the problem at hand. My suggestion is to use a simulation software that allows you to access the best-in-class solution methods. This means that you will be able to tackle your application with the right tools since your choice will be based on the physics involved and the computational resources available. It is indeed a good time to simulate your application since solution methods have improved at a rate even greater than the hardware, which exploded in capability in the past decades (see plot below). All the solution methods mentioned are available in COMSOL Multiphysics.

Rate of increase in performance of solution methods (solid) and hardware (dashed) vs. time. Legend: GE (Gaussian Elimination, direct method), GS (Gauss-Seidel, iterative method), SOR (Successive Over Relaxation, iterative method), CG (Conjugate Gradient, iterative method), MG (Multigrid, iterative method). Source: Report to the President of the United States titled “Computational Science: Ensuring Americaâ€™s Competitiveness” written by the Presidentâ€™s Information Technology Advisory Committee in 2005.

My colleagues developing the solvers available in COMSOL Multiphysics continually take advantage of these improvements and make sure that we offer you the highest-performance methods. As you can see from the plot above, multigrid methods are the ones that better stand the comparison with Moore’s law. Luckily, they are available in COMSOL Multiphysics. Since multigrid methods are the most powerful, I’m going to discuss the main idea behind these in the remainder of this blog post.

### Why Multigrid Methods are Needed

In order to introduce you to the basic ideas behind this solution method, I will present you with numerical experiments exposing the intrinsic limitations of iterative methods that multigrid methods are built on. For the sake of brevity, I will analyze only their qualitative behavior (I suggest you read A Multigrid Tutorial by William L. Briggs for more details on what is mentioned here).

Let’s start iterating with an approximated initial guess consisting of Fourier modes. They can be written as , where is the wavenumber (see plot below) and .

Fourier modes with . The mode consists of half sine waves on the domain. For small values of we have long smooth waves, while for larger values we have short oscillatory waves.

The qualitative plot below shows that convergence is slower for smooth waves and quicker for oscillatory waves. Smooth waves are detrimental for the performance of iterative methods; an intrinsic limitation of iterative methods. A typical approximated initial guess will contain several Fourier modes: convergence will start to slow down after the first iterations. That is because oscillatory components are efficiently eliminated from the error, while the smooth ones prevail and are left almost unchanged at every iteration, in other words convergence is stalling.

Log of the maximum norm of the algebraic error plotted vs. the iteration number for three different initial guesses
().

Moreover, and ironically, decreasing the grid spacing will not improve the accuracy of the solution. Having defined as the grid spacing, a Fourier mode with wavelength less than (oscillatory wave) will appear on the grid as having a wavelength greater than (smooth wave) because of spatial aliasing. In the regions of the domain representing our application where we want more accuracy, decreasing the grid spacing will actually worsen the performance of iterative methods. This is another intrinsic limitation.

Multigrid methods overcome these limitations by adapting iterative methods to make them efficient for all Fourier modes. This will be discussed in the second part of this post that will be published next week, so stay tuned!