MUDPACK: Multigrid Software for Elliptic Partial Differential Equations

(Fortran Code with OpenMP Directives for Shared Memory Parallelism)

by

John C. Adams
September 1999

John C. Adams is Available for Consulting

Version 5.0.1 of MUDPACK is identical to the previous version 5.0 except it fixes a problem in subroutine mud2sp: subroutine rmd2spp is now called instead of subroutine rcd2spp (there was one such instance of this incorrect call in mud2sp). Version 5.0.1 is incompatible with MUDPACK versions predating 4.0. For more details and parallel performance measurements, click on: Important information for users of earlier versions of MUDPACK software

This description consists of seven parts:

1. MUDPACK Introduction

2. MUDPACK Special Features

3. MUDPACK Solver Description

4. MUDPACK Solver Selection Flowchart

5. MUDPACK and Multigrid References

6. Obtaining MUDPACK Software

7. Acknowledgements

1. MUDPACK Introduction

MUDPACK is a collection of vectorized, portable, Fortran 77 subprograms, with a few Fortran90 extensions, for efficiently solving linear elliptic Partial Differential Equations (PDEs) using multigrid iteration. The most frequent Fortran90 extension used is the DO-END DO loop. Users may compile with any Fortran90 compiler, but Fortran77 compilers usually accept the code without complaint. In Fortran90 terminology, the source code is fixed-form, not free-form. OpenMP directives are included to enable shared memory parallelism. The package was created to make multigrid iteration available in user friendly form. The software is written in much the same format as the separable elliptic PDE package FISHPACK [5] (see FISHPACK). It extends the domain of solvable problems to include both separable and nonseparable PDEs. Detailed descriptions of earlier versions of MUDPACK are given [2,9].

Multigrid iteration (see [12,13,15,17,18,20]) combines classical iterative techniques, such as Gauss-Seidel line or point relaxation, with subgrid refinement procedures to yield a method superior to the iterative techniques alone. By iterating and transferring approximations and corrections at subgrid levels, a good initial guess and rapid convergence at the fine grid level can be achieved. Multigrid iteration requires less storage and computation than direct methods for nonseparable elliptic PDEs (e.g., see [7]) and is competitive with direct methods such as cyclic reduction [5,14,24,25] for separable equations. In particular, three-dimensional problems can often be handled at reasonable computational cost. Achieving optimal multigrid performance requires hand-tailored coding for certain problems. The generality of the equations solved by MUDPACK software may sometimes result in loss of efficiency. It is hoped that this is compensated for by the package's ease of use, applicability to a wide range of real problems (including those typically encountered in the atmospheric sciences at NCAR [8]), and avoidance of repeated "re-inventions of the wheel." Savings in human code development time can be at least as important as economic use of machine cycles. With careful selection of relaxation and multigrid parameters, optimal performance can often be realized using MUDPACK software. See [1,3,9,10] for a variety of problems where discretization level error (i.e., the same error that a direct method will produce) is reached in only one full multigrid cycle using MUDPACK solvers. Supercomputer performance from a decade ago [22] was measured with the examples in [1,7,9,10].

2. MUDPACK Special Features

* Solving Linear Elliptic PDEs in a Variety of Forms

These forms include real and complex, two- and three-dimensional, self-adjoint, separable and nonseparable, and PDEs with cross derivative terms (see part 3 of this file).

* Solving PDEs in Curvilinear Coordinate Systems

The solution regions are rectangular regions in the sense that the domain of each independent variable must be a bounded interval on the real line. This means that curvilinear coordinate systems such as spherical or cylindrical coordinates are acceptable. The codes are not restricted to Cartesian coordinates.

* Generating Second- and Fourth-order Approximations

Standard second-order finite difference approximations are generated on uniform grids superimposed on the solution region. These can be improved to fourth-order estimates using "deferred corrections" ([21,23]).

* Handling of General Boundary Conditions

Any combination of periodic, specified (Dirichlet), and mixed derivative boundary conditions is allowed. Some of the solvers allow oblique (non-normal) derivative boundary conditions.

* Ease of Input of the Continuous Problem

User defined input subroutines are the mechanisms for passing PDE coefficients and boundary conditions.

* Automatic Discretization of the Continuous Problem

The discretization is transparent to a user who only needs to supply the PDE, boundary conditions, and grid size information. Standard second-order finite difference formula are used to approximate first and second partial derivatives. The result is a linear block tridiagonal system of equations. More complex difference formula (asymmetric near boundaries [4]) are used with the fourth-order solvers. The coefficients multiplying the second partial derivatives in the PDE are adjusted during discretization at coarser grid levels if there are nonzero first-order coefficients which would destroy diagonal dominance. This is necessary to preserve convergence of iterative schemes.

* Use of Multigrid Iteration to Approximate the Discretization Equations

This is the essential feature of the MUDPACK software. It makes this complex collection of integrated numerical procedures available in friendly form.

* Flexibility in Choosing Grid Size

Second and fourth order approximations are generated on uniform l by m by n grids superimposed on boxes in three dimensions or l by m grids superimposed on rectangles in two dimensions. The grid sizes have the form:

l = p * 2(i-1) + 1

m = q * 2(j-1) + 1

n = r * 2(k-1) + 1

where p,q,r,i,j,k are positive integers. i,j,k > 0 determine the number and size of the subgrid levels employed by multigrid cycling. Values for p,q,r should be chosen as small as possible (typically 2,3 or 5) and values for i,j,k as large as possible within grid size requirements for efficient cycling. In particular, larger values for p,q or r can cause cause algorithm deterioration. For 2-d and 3-d nonseparable PDEs this can be bypassed by using one of the "hybrid" solvers described below.

Let G denote the finest l by m by n grid. In MUDPACK, multigrid cycling is implemented on the ascending chain of grids

G(1) < ... < G(s) < ... < G(t) = G

where t = max0(i,j,k) and each G(s) (s=1,...,t) has l(s) by m(s) by n(s) grid points given by:

l(s) = p * 2max0(i+s-t,1) + 1

m(s) = q * 2max0(j+s-t,1) + 1

n(s) = r * 2max0(k+s-t,1) + 1

When grid size requirements cannot be met with MUDPACK software (even with one of the hybrid solvers described below) then one option is to choose a grid which does satisfy the constraints which is as close as possible to the required grid and solve the problem there. The approximation can then be transferred to the required grid using multidimensional cubic interpolation.

* Selection of Multigrid Options

MUDPACK has options for implementing variants of multigrid iteration and default options for those preferring black box solvers. The default options (chosen for robustness) set cubic prolongation, fully weighted residual restriction, and W(2,1) cycling. The earlier version of MUDPACK described in [2,3] only allowed V(2,1) cycling with linear prolongation. This is still available as a possibly more efficient choice for certain problems.

* Selection of the Relaxation Method used within Multigrid Iteration

A relaxation menu is provided. It includes vectorized Gauss-Seidel schemes [16] on alternating points (red/black), lines (in any combination of directions) and planes (for three-dimensional anisotropic elliptic PDEs [26]). Choice of the correct relaxation method for a particular problem can be crucial. It depends on the relative grid and PDE coefficient size. Usually this can be pre-determined. Sometimes experimentation is required. Advice on method selection is given in the documentation.

* Availability of "hybrid" Multigrid/Direct Method Solvers

The certainty of direct methods is combined with the efficiency of multigrid iteration by providing "hybrid" solvers for 2-d and 3-d nonseparable PDEs. Gaussian elimination is used whenever the coarsest grid is encountered within multigrid cycling. This eliminates the usual constraint that the coarsest grid must have "few" points thus giving additional flexibility in choosing grid size. It also provides a way to compare approximations from multigrid and direct method solutions. The hybrid codes become full direct method solvers replacing the codes described in [6] if grid size arguments are chosen so that the coarse and fine grids coincide. Large storage and computational requirements make the use of the 3-d hybrid codes muh3,cuh3 as direct methods possible only for very coarse grids.

* Availability of Subroutines to Compute Residuals

Subroutines to compute fine grid residual after calling any of the second-order solvers are provided. The residual measures how well the current approximation satisfies the linear system of equations coming from the discretization. Residual ratios can be used to estimate the convergence rate of multigrid iteration.

* No Initial Guess Requirement

Unlike the case with classical iterative schemes, initial guesses are not necessary and should not be supplied unless they are very good (as, for example, when restarting multigrid iteration using an approximation generated earlier). Full multigrid cycling [12], beginning at the coarsest grid level, is used when there is no initial guess. Advice on how to use initial guesses within a time marching problem is given in the documentation.

* Non-initialization Calls

Redundant discretization and matrix factorization processes can (and should) be bypassed on recalls to the software. For example, this happens when only the right-hand side array has changed from a previous call or when more multigrid cycles are needed for additional accuracy.

* Error Control

Maximum relative error can be used to monitor convergence. Use of error control is optional and requires additional storage and computation.

* Flagging of Errors involving Input Parameters

This includes detection of singular and/or nonelliptic PDEs. Fatal and nonfatal errors are flagged.

* Output of Exact Minimal Work Space Requirements

This is especially important with three-dimensional problems where central memory is easily exhausted.

* Extensive Documentation and Test Programs

Users are encouraged to carefully read the documentation and execute the test program for the solver to be used. The next section provides links to documentation and fortran test program files.

3. MUDPACK Solver Description

Table 1 below lists all mudpack two- and three-dimensional, second and fourth order solvers for real and complex elliptic partial differential equations with and without cross derivative terms. Clicking on a solver will bring up its documentation file. Table 2 provides a list of the test and residual codes for each solver. The complete Fortran for these codes can be accessed by clicking on their name. The documentation for solvers and fortran code for the mudpack test programs and residual codes can also be downloaded as text files from this site by clicking on their names in the tables below. Section 6 of this document explains how the entire MUDPACK package (Fortran solvers and support files, test programs, and documentation) can be downloaded.

Table 1
An overview of MUDPACK Solvers
computation subprograms
2nd order/real 2D self-adjoint nonseparable mud2sa
2nd order/real 2D separable mud2sp
2nd order/real 2D nonseparable muh2, mud2
2nd order/real 2D with cross term muh2cr, mud2cr
4th order/real 2D separable mud24sp
4th order/real 2D nonseparable muh24, mud24
4th order/real 2D with cross term muh24cr, mud24cr
2nd order/real 3D self-adjoint nonseparable mud3sa
2nd order/real 3D separable mud3sp
2nd order/real 3D nonseparable muh3, mud3
2nd order/real 3D with cross term mud3cr
4th order/real 3D separable mud34sp
4th order/real 3D nonseparable mud34, muh34
2nd order/complex 2D separable cud2sp
2nd order/complex 2D nonseparable cuh2, cud2
2nd order/complex 2D with cross term cuh2cr, cud2cr
4th order/complex 2D separable cud24sp
4th order/complex 2D nonseparable cud24, cuh24
4th order/complex 2D with cross term cud24cr, cuh24cr
2nd order/complex 3D separable cud3sp
2nd order/complex 3D nonseparable cuh3, cud3
2nd order/complex 3D with cross term cud3cr
4th order/complex 3D separable cud34sp
4th order/complex 3D nonseparable cud34

Table 2
Solver test & residual codes
mud2sa tmud2sa
mud2sp tmud2sp, resm2sp
mud2,muh2 tmud2, muh2, resm2
mud2cr,muh2cr tmud2cr, tmuh2cr, resm2cr
mud24sp tmud24sp
mud24,muh24 tmud24, tmuh24
mud24cr,muh24cr tmud24cr, tmuh24cr
mud3sa tmud3sa
mud3sp tmud3sp, resm3sp
mud3,muh3 tmud3, tmuh3, resm3
mud3cr tmud3cr
mud34sp tmud34sp
mud34,muh34 tmud34, tmuh34
cud2sp tcud2sp, resc2sp
cud2,cuh2 tcud2, tcuh2, resc2
cud2cr,cuh2cr tcud2cr, tcuh2cr, resc2cr
cud24sp tcud24sp
cud24,cuh24 tcud24, tcuh24
cud24cr,cuh24cr tcud24cr, tcuh24cr
cud3sp tcud3sp, resc3sp
cud3,cuh3 tcud3, tcuh3, resc3
cud34sp tcud34sp
cud34 tcud34

4. MUDPACK Solver Selection

The following "flow chart" can be used in selecting the appropriate second-order software for the elliptic PDE to be solved:

(1) If the PDE is complex go to (9) else go to (2)

(2) If the PDE is three-dimensional go to (6) else go to (3)

(3) If the PDE is separable use mud2sp else go to (4)

(4) If the PDE has a cross derivative use muh2cr or mud2cr else go to (5)

(5) If the PDE is self-adjoint use mud2sa else use muh2 or mud2.

(6) If the PDE is separable use mud3sp else go to (7)

(7) If the PDE is self-adjoint use mud3sa else go to (8)

(8) If the PDE has cross derivatives use mud3cr else use muh3 or mud3.

(9) If the PDE is three dimensional go to (13) else go to (10)

(10) If the PDE is separable use cud2sp else go to (11)

(11) If the PDE has a cross derivative use cuh2cr or cud2cr else go to (12)

(12) Use cuh2 or cud2

(13) If the PDE is separable use cud3sp else use cuh3 or cud3.

Fourth-order solvers can improve the approximation if the corresponding second-order solver has reached discretization level error (i.e., the same error level that a direct method will reach) [1,3,10].

5. MUDPACK and Multigrid References

[1] J. Adams, "Multigrid Software for Elliptic Partial Differential Equations: MUDPACK," NCAR Technical Note-357+STR, Feb. 1991, 51 pages.

[2] J. Adams, "MUDPACK: Multigrid Fortran Software for the Efficient Solution of Linear Elliptic Partial Differential Equations," Applied Math. and Comput. Vol.34, Nov 1989, pp.113-146.

[3] J. Adams, "FMG Results with the Multigrid Software Package MUDPACK," Proceedings of the Fourth Copper Mountain Conference on Multigrid, SIAM, 1989, pp.1-12.

[4] J. Adams, "Fortran Subprograms for Finite Difference Formula," J. Comp. Phys.,Vol 26, Jan 1978, pp. 113-116.

[5] J. Adams, P. Swarztrauber, R. Sweet, "Efficient Fortran Subprograms for the Solution of Elliptic Partial Differential Equations," Elliptic Problem Solvers, Academic Press, 1982, pp.333-390.

[6] J. Adams, "New Software for Elliptic Partial Differential Equations," Computing Facility Notes 55, November 1978

[7] J. Adams, "Comparison of direct and iterative methods for approximating nonseparable elliptic PDEs at NCAR," SCD Computing News, Vol 10, Nov. 1989, pp.12-14.

[8] J. Adams, R. Garcia, B. Gross, J. Hack, D. Haidvogel, and V. Pizzo, "Applications of Multigrid Software in the Atmospheric Sciences," Monthly Weather Review,Vol. 120 # 7, July 1992, pp. 1447-1458.

[9] J. Adams, "MUDPACK: Multigrid Software for Linear Elliptic Partial Differential Equations," SCD UserDoc, Version 2.0, NCAR,February 1990.

[10] J. Adams, "Recent Enhancements in MUDPACK, A Multigrid Software Package for Elliptic Partial Differential Equations," Applied Math. and Comp., Vol. 43, May 1991, pp.79-94.

[11] J. Adams, "MUDPACK-2: Multigrid Software for Elliptic Partial Differential Equations on Uniform Grids with any Resolution," Applied. Math. and Comp., Vol. 53, February 1993, pp. 235-249.

[12] A. Brandt, "Multi-level Adaptive Solutions to Boundary Value Problems," Math. Comp., 31, 1977, pp.333-390.

[13] W. Briggs, "A Multigrid Tutorial," SIAM, Philadelphia,1987.

[14] B. Buzbee, G. Golub, and C. Nielson, "On direct methods for solving Poisson's equations," SIAM J. Numer. Anal., 7, 1970, pp.627-656.

[15] S. Fulton, R. Ciesielski, and W. Schubert, Multigrid methods for elliptic problems: a review. Monthly Weather Review, 114:943-959 (1986).

[16] W. Gentzsch, "Vectorization of Computer Programs with Applications to Computational Fluid Dynamics," Vieweg & Sohn, 1984 (246 pages).

[17] W. Hackbush and U. Trottenberg, "Multigrid Methods," Springer-Verlag, Berlin,1982.

[18] D. Jespersen, "Multigrid Methods for Partial Differential Equations." Studies in Numerical Analysis, Vol.24, MAA, 1984.

[19] J. Mandel and S, Parter, "On the Multigrid F-Cycle," Applied Math. and Comput., Vol 37, 1990, pp.19-36.

[20] S. McCormick, "Multigrid Methods," Vol 3 of SIAM Frontiers Series, SIAM, Philadelphia, 1987.

[21] V. Pereyra, "Highly Accurate Numerical Solution of Casilinear Elliptic Boundary-Value Problems in n Dimensions," Math. Comp., 24, 1970, pp.771-783.

[22] D. Sato, "PERFMON: The Cray Performance Monitor Utility," SCD UserDoc, Version 2.0, NCAR,March 1989.

[23] S. Schaffer, "Higher Order Multigrid Methods," Math. Comp., Vol 43, July 1984, pp. 89-115.

[24] P. Swarztrauber, "Fast Poisson Solvers," Studies in Numerical Analysis, Math. Assoc, America, 1985, pp. 319-369.

[25] R. Sweet, "A Parallel and Vector Variant of the Cyclic Reduction Algorithm," SIAM J. Sci. and Stat. Comp., Vol. 9, July 1988, pp. 761-766.

[26] C. Thole and U. Trottenberg, "Basic Smoothing Procedures for the Multigrid treatment of Elliptic 3D Operators," Applied Math. and Comp., 19, 1986, pp. 333-345.

6. Obtaining MUDPACK Software

MUDPACK solver documentation, test programs, and residual subroutines can be downloaded as text files from this web site (see Section 3 of this document).

The entire MUDPACK software package can now be downloaded after signing a UCAR licensing agreement which includes, among other provisions, that the software is not to be used for commercial purposes, modified, or distributed further.

DOWNLOAD MUDPACK SOFTWARE

When you order MUDPACK from NCAR:

You are assured of receiving original source code from its creator, so you avoid security concerns associated with pirated or "shared" software.

To order MUDPACK, complete the MUDPACK order form

The NCAR Technical Note, "MUltigriD Software for Elliptic Partial Differential Equations: MUDPACK" describing the earlier 1991 software version and reprints of the publications listed in Section 5 can be obtained from John C. Adams.

ACCESSING SOLVERS AT NCAR:

Binary libraries have been established at NCAR. If you compute at NCAR, please see this web page to learn which computers have the binary library, and the following web page to see how to load it:

http://www.scd.ucar.edu/docs/products

http://www.scd.ucar.edu/softlib/link.html

7. Acknowledgements

Steve McCormick introduced the author to the multigrid community and provided numerous helpful suggestions including the use of planar relaxation with the three-dimensional solvers. The importance of adjusting discretization coefficients at coarser grid levels for PDEs with nonzero first-order terms was pointed out by Klauss Steuben. Achi Brandt provided a complimentary foreword for the MUDPACK technical note [1]. A conversation with Achi Brandt affirmed that the default multigrid options in MUDPACK are a good choice and that the use of deferred corrections in obtaining fourth-order approximations with multigrid is a reasonable strategy. Dave Kennison of the NCAR graphics group provided the grid coarsening figure at the start of this document.


Return to beginning of this document