Fast Solvers For Complex Problems

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 finite-element shape. After a process called the FE-assembly, the PDE problem is transformed into a system of nonlinear or linear algebraic equations. For simplicity, assume we are solving a stationary (time-independent) 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 memory-intensive for larger 3D problems such as fluid-structure interactions or high-frequency 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 high-quality preconditioners, basically numerical methods that simplify A. However, a well-chosen 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 so-called V-cycle is the core of many multigrid methods. The GMG V-cycle 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, Gauss-Seidel, 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 coarse-mesh 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 V-cycle ends with postsmoothing at the finest mesh. The combined effect of smoothing and coarse-mesh 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 V-cycles are applied, or as a stand-alone linear solver where the GMG V-cycles 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 shape-function set. Users can produce the mesh hierarchy manually, or have the software generate it. The user selects standard pre and postsmoothers, specialized preconditioners, and coarse-mesh solvers in a general way from lists of solvers.

Many real-world 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 more-specialized SOR-vector preconditioner, intended for the vector Helmholtz equation as discretized with vector (or edge) elements, and the Vanka preconditioner, intended for the velocity-divergence 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 electromagnetic-wave 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 SOR-vector and SORU-vector 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 CPU-time comparison. An interesting number in this context is the break-even 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