ME 540

Numerical Solutions Applied To Heat Transfer And Fluid Mechanics Problems


Domain Decomposition

To model fluid flow or heat transfer in complex spatial regions it is often useful to break the spatial domain into smaller regions. The primary focus is on solution within these subregions. Coupling between subregions becomes part of an iteration procedure. This approach also has the advantage that computations on subregions can be distributed to separate processors in a parallel computing environment. As an example suppose that a spatial domain has been split into six sub-domains. A grid has been established in each sub-domain and appropriate equations evaluated at all these points to produce a general linear system Ax=b. For this example I'll assume that sub-domain numbering and connectivity is such that a given sub-domain only interacts with sub-domains having the next higher and next lower indices. In this situation A can be written in block matrix form as:

If I denote xi as that portion of the solution vector evaluated within the ith sub-domain, then I can write an iteration scheme in the form:

S1x1(j+1) = b1 - C1,2 x2(j)

Si xi(j+1) = bi - Ci,i-1 xi-1(j) - Ci,i+1 xi+1(j) for i = 2,5

S6x6(j+1) = b6 - C6,5 x5(j)

You may recognize this as the being the same form as a point Jacobi iteration. It is generally called a Block Jacobi method. It has the strong advantage that each of the six block equations can be distributed to a separate processor. It has the disadvantage that information about a guess in a given sub-domain is only propagated one sub-domain per iteration. For example the initial guess at a solution in the first sub-domain does not impact the solution of the sixth sub-domain until the fifth iteration and effects of full feedback between the first and sixth sub-domains are not apparent until the 10th iteration. This significantly slows convergence, and eliminates error in a selective fashion. If we consider Fourier transforms of the error vector from iteration to iteration, components with the shortest wave length (cell to cell variation) are suppressed quickly. Those with the longest wavelengths are the most persistent. One solution is treat domain decomposition as part of a multigrid approach described below.


Multigrid Methods

The basic idea behind multigrid method is to deal with long wavelength errors by solving the problem on coarser grids. After defining the finest required grid resolutions, a sequence of successively coarser grids are defined for the same domain. Normally this is accomplished by deleting points in the even numbered planes (lines in 2-D) of each finer grid. This is continued until a grid is obtained where equations can be solved directly or through a tight convergence of a very inexpensive iteration. For example, start with the following grid with x and y spacing set to "d".

The next grid has mesh spacing equal to 2d

A final grid is defined with mesh spacing equal to 4d. A system of 9 equations is relatively inexpensive to solve with direct methods and no further refinement is required.

To proceed with the discussion, I'm going to introduce a subscript "h" to all matrix and vector objects to index the grid being used. When it is necessary to refer to operations between two grids "h" will indicate the finer and "H" the coarser.

The physical problem is modeled on each of the grids resulting in a series of linear systems:

Ahxh = bh

For my example h runs from 3 (finest) to 1 (coarsest). However, using this form of the equation won't get us too far. The object is to use the coarser grids to suppress error, so we need to transform the system to one dealing with error terms. At a given iteration and mesh define the solution error as:

eh(j) = xh - xh(j)

and the equation residual as:

rh(j) = bh - Ahxh(j)

Multiplying the error definition by Ah gives:

Ah eh(j) = Ah xh - Ah xh(j) = bh - Axh(j) = rh(j)

One form of the multigrid algorithm proceeds in the following series of steps.

1 At the finest grid level, generate some initial estimate of the solution to Ax=b and apply your favorite iteration method with a limited number of iterations or relatively loose convergence criterion to smooth the solution to Ax = b on that grid.
2 Now project rh(j) and perhaps eh(j) onto the next coarser grid. The simplest way to do this is a direct transfer of the values on the finer grid to the corresponding points on the coarser grid, but this does not produce the best results. A better solution is to also project values from adjacent fine grid points onto the course grid. If point l,m on the coarser grid is at the same spatial location as point i,j on the finer grid, then one such projection method is:

rH l,m = (4 rh i,j + 2(rh i+1,j +rh i-1,j +rh i,j+1 +rh i,j-1 ) + rh i+1,j+1 +rh i+1,j-1 +rh i-1,j-1 +rh i-1,j+1 )/16

3 Obtain an approximate solution to AH eH =rH on the coarser grid. If eh has been projected to eH, that is the initial estimate the coarser grid iteration. If not, the initial estimate is eH=0.
4 Interpolate the coarser grid results for error eH onto the finest grid.
5 Adjust your fine grid solution estimate based upon the new estimate of eh, and use your favorite iteration method again with a limited number of iterations. Check against a final convergence criterion. If no convergence go to step two to start another cycle.

Step three in the above algorithm follows one of two paths. If it uses the coarsest grid, then AH eH=rH is solved directly, and step 4 immediately follows. If not, a recursive pattern follows using a slight variation on the above five steps. The first and last steps are stated differently, because the coarser grids have no knowledge of the original right hand side "b", and a tight final convergence is not generally needed or desirable.

1 At intermediate grid levels, use your favorite iteration method with a limited number of iterations or relatively loose convergence criterion to smooth the solution to Ah eh =rh on that grid.
2 Now project rh(j) and perhaps eh(j) onto the next coarser grid (see the discussion of step 2 on the finest grid above).
3 Obtain an approximate solution to AH eH =rH on the coarser grid. If eh has been projected to eH, that is the initial estimate the coarser grid iteration. If not, the initial estimate is eH=0.
4 Interpolate the coarser grid results for error "eH" onto this grid to revise the estimate of eh on this grid.
5 Use your favorite iteration method again with a limited number of iterations or relatively loose convergence criterion to further improve the solution of Ah eh =rh on this grid. Move back up to the solution on the next finer grid.

Because the above algorithm proceeds directly from the finest grid down to the coarsest and back, it is often referred to as a V-cycle. For problems where tighter convergence on intermediate grids is desirable, the most efficient approach is not to tighten convergence in step 5 above. A better approach is to run through step 3 more than once. This results in a pattern of motion through the grids called a W-cycle.

One other variation on this method can be seen in the literature. An improved initial estimate can be generated by starting with a solution to Ax=b on the coarsest grid, and recursively using the result as an initial guess for a V-cycle solution on each of the successively finer grids.

Like other iteration procedures, multigrid method does not work well for all problems. In some cases it will fail where other iteration methods succeed. However, when it does work, it is normally more efficient than other options. The state of the art continues to improve with this method. When you need an iterative solver for large systems of equations, this should be your first attempt. You should also do a quick survey of the latest literature for improvements to details in the method. For more information on multigrid, and sample software, look at referenced sites on the class home page.


Other Lectures / Home


Created and Maintained by John Mahaffy : jhm@psu.edu