Jacob Yström
Group leader for Solvers at the Development Dept.
Comsol Inc.
Stockholm, Sweden
Edited by Leslie Gordon

The models use ordinary differential equations (ODEs) and partial differential equations (PDEs). The software subdivides the particular physical system, represented by a geometric model, into a mesh, or a finite number of smaller, discreet elements.
Instead of directly solving PDEs over each region, the FE method first numerically approximates the unknowns of interest (for example, temperature, frequency, and velocity) using predetermined mathematical functions expressed interms of the geometrical locations (nodes) used to define the finiteelement shape. After a process called the FEassembly, the PDE problem is transformed into a system of nonlinear or linear algebraic equations. For simplicity, assume we are solving a stationary (timeindependent) PDE problem.
Algorithms known as solvers then come into play. They further divide the problem into manageable chunks. An important such chunk is to solve large systems of linear algebraic equations, generally in the form Ax = b, where A represents a sparse matrix (most of the entries are zero), b represents known quantities that, for example, come from source terms and boundary forces, and x represents what the user wants to find out, or the degrees of freedom.
Linear system solvers come in a few different flavors. For example, what are known as direct solvers solve problems with Gauss elimination. Direct solvers are generally fine for 2D problems such as predicting the propagation of cracks caused by a high concentration of stresses. But the solvers are too memoryintensive for larger 3D problems such as fluidstructure interactions or highfrequency wave problems. Such problems easily involve millions of degrees of freedoms.
These complex problems demand iterative solvers, which start from a suitable initial guess and then calculate successive updates that ultimately converge (approach a certain value) to a solution. Iterative solvers are more difficult to use because they require highquality preconditioners, basically numerical methods that simplify A. However, a wellchosen preconditioner reduces the number of iterations required and, therefore,saves computing time.
Among the many available preconditioners, the geometric multigrid (GMG) technique is widely used because it efficiently handles a large class of problems, particularly those in structural mechanics, fluid flow, and electromagnetics. Analysis software such as Comsol Multiphysics supports the GMG method.
BASIC GEOMETRIC MULTIGRID CONCEPT
The GMG algorithm jumps between different linear systems representing coarser and finer meshes of the underlying PDEs. The idea is to perform a fraction of the computations on a fine mesh while solving the full system under study on a coarser one. This provides an iterative algorithm that is fast and memory efficient.
The socalled Vcycle is the core of many multigrid methods. The GMG Vcycle starts with a fine mesh and uses a simple iterative solver to apply what's called presmoothing to the residual function associated with the error in the approximation.
This smoothing often only takes a few iterations using the standard Jacobi, GaussSeidel, or successive overrelaxation (SOR) method. Each of these would split the matrix a different way to solve the linear system. The software then maps the result to a coarser mesh using an operator called a restriction. This smoothing and restriction procedure repeats until the algorithm reaches the model's coarsest mesh.
Here, the GMG method uses another kind of solver — a coarsemesh solver — often a direct linear solver. The software then applies an operator called a prolongation to map the coarse mesh's corrected residual to a finer mesh. Then simple iterative solvers are applied, a process known as postsmoothing. The prolongation and smoothing procedure is repeated until the Vcycle ends with postsmoothing at the finest mesh. The combined effect of smoothing and coarsemesh correction substantially reduces the residuals, which, in turn, allows handling larger problems.
In Comsol Multiphysics, the GMG technique can be used as a preconditioner where a fixed number of GMG Vcycles are applied, or as a standalone linear solver where the GMG Vcycles are repeated until convergence. The user applies the GMG cycle to a hierarchy of mesh cases, where each consists of a mesh and a corresponding shapefunction set. Users can produce the mesh hierarchy manually, or have the software generate it. The user selects standard pre and postsmoothers, specialized preconditioners, and coarsemesh solvers in a general way from lists of solvers.
Many realworld problems are handled successfully with the GMG method and by using standard preconditioners as smoothers. Yet there are problems that don't lend themselves well to using standard preconditioners. Comsol Multiphysics also supplies the morespecialized SORvector preconditioner, intended for the vector Helmholtz equation as discretized with vector (or edge) elements, and the Vanka preconditioner, intended for the velocitydivergence formulation of the incompressible NavierStokes equations.
MULTIGRID AND DIRECT SOLVERS
To compare the performance of multigrid with direct solvers, take as an example the model of a 3D electromagneticwave problem, the Helmholtz equation for a corner cube reflector. The large system of linear equations that results from applying FE analysis with linear vector elements is complex and symmetric, but not Hermitian (in this case the matrix is real symmetric except for the main diagonal where complex numbers are encountered). A direct solver such as SPOOLES can take advantage of such symmetry and work on half the matrix, whereas UMFPACK, another direct solver, must work on the whole matrix.
To use the GMRES iterative method with a GMG preconditioner (GMRES/MG) in Comsol Multiphysics, a user would first select the appropriate smoothers. In this example, say, the user chooses the SORvector and SORUvector algorithms for the pre and postsmoothers. The user manually generates a hierarchy of three nested meshes by starting with a coarse mesh and then applying the Regular Refinement feature twice in sequence. The procedure does not require tuning of the restarted GMRES/MG method, which typically converges after 20 to 25 iterations.
Tests have shown that, for a given amount of time, the GMRES/GMG method solves far larger problems than SPOOLES or UMFPACK can handle. Plotting the number of calculated degrees of freedom against CPU memory usage for the GMRES/GMG method shows the approximate linear behavior of the resulting curve. Similar results hold for a CPUtime comparison. An interesting number in this context is the breakeven point, where solution time and memory usage are better for the GMRES/GMG method than with the direct solvers. The number is about 50,000 complex degrees of freedom. Thus, GMRES/GMG is faster and more memory eco nomical than direct solvers for problems with greater numbers of degrees of freedom.
MAKE CONTACT
Comsol Inc.,comsol.com