**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 **A**x=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 x_{i }as that portion of the solution vector evaluated within the ith sub-domain, then I can write
an iteration scheme in the form:

**S**_{1}x_{1}^{(j+1) }= b_{1} - **C**_{1,2} x_{2}^{(j)}

**S**_{i} x_{i}^{(j+1) }= b_{i} - **C**_{i,i-1} x_{i-1}^{(j)} - **C**_{i,i+1} x_{i+1}^{(j)} for i = 2,5

**S**_{6}x_{6}^{(j+1) }= b_{6} - **C**_{6,5} x_{5}^{(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 10^{th} 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.

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 physical problem is modeled on each of the grids resulting in a series of linear systems:

**A**_{h}x_{h} = b_{h}

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:

e_{h}^{(j)} = x_{h} - x_{h}^{(j)}

and the equation residual as:

r_{h}^{(j)} = b_{h} - **A _{}**hx

Multiplying the error definition by A_{h} gives:

_{}^{} **A**_{h} e_{h}^{(j)} = **A**_{h} x_{h} - **A**_{h} x_{h}^{(j)} = b_{h} - **A**x_{h}^{(j)} = r_{h}^{(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 r_{h}^{(j)} and perhaps e_{h}^{(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:
r |

3 | Obtain an approximate solution to A_{H} e_{H} =r_{H} on the coarser grid. If e_{h}
has been projected to e_{H}, that is the initial estimate the coarser grid
iteration. If not, the initial estimate is e_{H}=0. |

4 | Interpolate the coarser grid results for error e_{H} onto the finest grid. |

5 | Adjust your fine grid solution estimate based upon the new estimate of
e_{h}, 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
**A**_{H} e_{H}=r_{H} 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 A_{h} e_{h} =r_{h} on that grid. |

2 | Now project r_{h}^{(j)} and perhaps e_{h}^{(j)} onto the next coarser grid (see the
discussion of step 2 on the finest grid above). |

3 | Obtain an approximate solution to A_{H} e_{H} =r_{H} on the coarser grid. If e_{h}
has been projected to e_{H}, that is the initial estimate the coarser grid
iteration. If not, the initial estimate is e_{H}=0. |

4 | Interpolate the coarser grid results for error "e_{H}" onto this grid to revise
the estimate of e_{h} 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 A_{h} e_{h} =r_{h} 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 **A**x=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.

_{}